diff --git a/build.py b/build.py index 5a748a10f4..89af468e38 100755 --- a/build.py +++ b/build.py @@ -69,10 +69,10 @@ ) parser.add_argument( - "--ci-bench", - action=argparse.BooleanOptionalAction, - default=False, - help="Run the benchmarking script that is run in CI (default is --no-ci-bench)", + "--ci-bench", + action=argparse.BooleanOptionalAction, + default=False, + help="Run the benchmarking script that is run in CI (default is --no-ci-bench)", ) args = parser.parse_args() @@ -303,25 +303,46 @@ def run_python_integration_tests(cwd, interpreter): def run_ci_historic_benchmark(): branch = "main" output = subprocess.check_output( - ["git", "rev-list", "--since=1 week ago", "--pretty=format:%ad__%h", "--date=short", branch] + [ + "git", + "rev-list", + "--since=1 week ago", + "--pretty=format:%ad__%h", + "--date=short", + branch, + ] ).decode("utf-8") - print('\n'.join([line for i, line in enumerate(output.split('\n')) if i % 2 == 1])) + print("\n".join([line for i, line in enumerate(output.split("\n")) if i % 2 == 1])) output = subprocess.check_output( - ["git", "rev-list", "--since=1 week ago", "--pretty=format:%ad__%h", "--date=short", branch] + [ + "git", + "rev-list", + "--since=1 week ago", + "--pretty=format:%ad__%h", + "--date=short", + branch, + ] ).decode("utf-8") - date_and_commits = [line for i, line in enumerate(output.split('\n')) if i % 2 == 1] + date_and_commits = [line for i, line in enumerate(output.split("\n")) if i % 2 == 1] for date_and_commit in date_and_commits: print("benching commit", date_and_commit) result = subprocess.run( - ["cargo", "criterion", "--message-format=json", "--history-id", date_and_commit], + [ + "cargo", + "criterion", + "--message-format=json", + "--history-id", + date_and_commit, + ], capture_output=True, - text=True + text=True, ) with open(f"{date_and_commit}.json", "w") as f: f.write(result.stdout) + if build_pip: step_start("Building the pip package") @@ -514,6 +535,7 @@ def run_ci_historic_benchmark(): "ipykernel", "nbconvert", "pandas", + "qutip", "qiskit>=1.3.0,<2.0.0", ] subprocess.run(pip_install_args, check=True, text=True, cwd=root_dir, env=pip_env) diff --git a/library/chemistry/qsharp.json b/library/chemistry/qsharp.json new file mode 100644 index 0000000000..394a533cb1 --- /dev/null +++ b/library/chemistry/qsharp.json @@ -0,0 +1,42 @@ +{ + "lints": [ + { + "lint": "deprecatedNewtype", + "level": "error" + }, + { + "lint": "deprecatedFunctionConstructor", + "level": "error" + }, + { + "lint": "deprecatedWithOperator", + "level": "error" + }, + { + "lint": "deprecatedDoubleColonOperator", + "level": "error" + }, + { + "lint": "deprecatedSet", + "level": "error" + }, + { + "lint": "needlessParens", + "level": "error" + } + ], + "files": [ + "src/Tests.qs", + "src/Utils.qs", + "src/JordanWigner/Convenience.qs", + "src/JordanWigner/ConvenienceVQE.qs", + "src/JordanWigner/JordanWignerBlockEncoding.qs", + "src/JordanWigner/JordanWignerClusterOperatorEvolutionSet.qs", + "src/JordanWigner/JordanWignerEvolutionSet.qs", + "src/JordanWigner/JordanWignerOptimizedBlockEncoding.qs", + "src/JordanWigner/JordanWignerVQE.qs", + "src/JordanWigner/OptimizedBEOperator.qs", + "src/JordanWigner/StatePreparation.qs", + "src/JordanWigner/Utils.qs" + ] +} diff --git a/library/chemistry/src/JordanWigner/Convenience.qs b/library/chemistry/src/JordanWigner/Convenience.qs new file mode 100644 index 0000000000..491be6289e --- /dev/null +++ b/library/chemistry/src/JordanWigner/Convenience.qs @@ -0,0 +1,107 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. + +export + TrotterStepOracle, + QubitizationOracle, + OptimizedQubitizationOracle; + +import JordanWigner.JordanWignerBlockEncoding.JordanWignerBlockEncodingGeneratorSystem; +import JordanWigner.JordanWignerEvolutionSet.JordanWignerFermionEvolutionSet; +import JordanWigner.JordanWignerEvolutionSet.JordanWignerGeneratorSystem; +import JordanWigner.JordanWignerOptimizedBlockEncoding.JordanWignerOptimizedBlockEncoding; +import JordanWigner.JordanWignerOptimizedBlockEncoding.PauliBlockEncoding; +import JordanWigner.JordanWignerOptimizedBlockEncoding.QuantumWalkByQubitization; +import JordanWigner.Utils.JordanWignerEncodingData; +import JordanWigner.Utils.MultiplexOperationsFromGenerator; +import JordanWigner.Utils.TrotterSimulationAlgorithm; +import Std.Convert.IntAsDouble; +import Std.Math.Ceiling; +import Std.Math.Lg; +import Utils.EvolutionGenerator; + +// Convenience functions for performing simulation. + +/// # Summary +/// Returns Trotter step operation and the parameters necessary to run it. +/// +/// # Input +/// ## jwHamiltonian +/// Hamiltonian described by `JordanWignerEncodingData` format. +/// ## trotterStepSize +/// Step size of Trotter integrator. +/// ## trotterOrder +/// Order of Trotter integrator. +/// +/// # Output +/// A tuple where: `Int` is the number of qubits allocated, +/// `Double` is `1.0/trotterStepSize`, and the operation +/// is the Trotter step. +function TrotterStepOracle(jwHamiltonian : JordanWignerEncodingData, trotterStepSize : Double, trotterOrder : Int) : (Int, (Double, (Qubit[] => Unit is Adj + Ctl))) { + let generatorSystem = JordanWignerGeneratorSystem(jwHamiltonian.Terms); + let evolutionGenerator = new EvolutionGenerator { EvolutionSet = JordanWignerFermionEvolutionSet(), System = generatorSystem }; + let simulationAlgorithm = TrotterSimulationAlgorithm(trotterStepSize, trotterOrder); + let oracle = simulationAlgorithm(trotterStepSize, evolutionGenerator, _); + let nTargetRegisterQubits = jwHamiltonian.NumQubits; + let rescaleFactor = 1.0 / trotterStepSize; + return (nTargetRegisterQubits, (rescaleFactor, oracle)); +} + + +function QubitizationOracleSeperatedRegisters(jwHamiltonian : JordanWignerEncodingData) : ((Int, Int), (Double, ((Qubit[], Qubit[]) => Unit is Adj + Ctl))) { + let generatorSystem = JordanWignerBlockEncodingGeneratorSystem(jwHamiltonian.Terms); + let (nTerms, genIdxFunction) = generatorSystem!; + let (oneNorm, blockEncodingReflection) = PauliBlockEncoding(generatorSystem); + let nTargetRegisterQubits = jwHamiltonian.NumQubits; + let nCtrlRegisterQubits = Ceiling(Lg(IntAsDouble(nTerms))); + return ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, QuantumWalkByQubitization(blockEncodingReflection))); +} + +/// # Summary +/// Returns Qubitization operation and the parameters necessary to run it. +/// +/// # Input +/// ## jwHamiltonian +/// Hamiltonian described by `JordanWignerEncodingData` format. +/// +/// # Output +/// A tuple where: `Int` is the number of qubits allocated, +/// `Double` is the one-norm of Hamiltonian coefficients, and the operation +/// is the Quantum walk created by Qubitization. +function QubitizationOracle(jwHamiltonian : JordanWignerEncodingData) : (Int, (Double, (Qubit[] => Unit is Adj + Ctl))) { + let ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, oracle)) = QubitizationOracleSeperatedRegisters(jwHamiltonian); + let nQubits = nCtrlRegisterQubits + nTargetRegisterQubits; + return (nQubits, (oneNorm, MergeTwoRegisters(oracle, nTargetRegisterQubits, _))); +} + + +operation MergeTwoRegisters(oracle : ((Qubit[], Qubit[]) => Unit is Adj + Ctl), nSystemQubits : Int, allQubits : Qubit[]) : Unit is Adj + Ctl { + oracle(allQubits[nSystemQubits..Length(allQubits) - 1], allQubits[0..nSystemQubits - 1]); +} + + +function OptimizedQubitizationOracleSeperatedRegisters(jwHamiltonian : JordanWignerEncodingData, targetError : Double) : ((Int, Int), (Double, ((Qubit[], Qubit[]) => Unit is Adj + Ctl))) { + let ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, blockEncodingReflection)) = JordanWignerOptimizedBlockEncoding(targetError, jwHamiltonian.Terms, jwHamiltonian.NumQubits); + return ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, QuantumWalkByQubitization(blockEncodingReflection))); +} + + +/// # Summary +/// Returns T-count optimized Qubitization operation +/// and the parameters necessary to run it. +/// +/// # Input +/// ## jwHamiltonian +/// Hamiltonian described by `JordanWignerEncodingData` format. +/// ## targetError +/// Error of auxillary state preparation step. +/// +/// # Output +/// A tuple where: `Int` is the number of qubits allocated, +/// `Double` is the one-norm of Hamiltonian coefficients, and the operation +/// is the Quantum walk created by Qubitization. +function OptimizedQubitizationOracle(jwHamiltonian : JordanWignerEncodingData, targetError : Double) : (Int, (Double, (Qubit[] => Unit is Adj + Ctl))) { + let ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, oracle)) = OptimizedQubitizationOracleSeperatedRegisters(jwHamiltonian, targetError); + let nQubits = nCtrlRegisterQubits + nTargetRegisterQubits; + return (nQubits, (oneNorm, MergeTwoRegisters(oracle, nTargetRegisterQubits, _))); +} diff --git a/library/chemistry/src/JordanWigner/ConvenienceVQE.qs b/library/chemistry/src/JordanWigner/ConvenienceVQE.qs new file mode 100644 index 0000000000..ffb0872490 --- /dev/null +++ b/library/chemistry/src/JordanWigner/ConvenienceVQE.qs @@ -0,0 +1,92 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. + +export + EstimateEnergyWrapper, + EstimateEnergy; + +import JordanWigner.JordanWignerEvolutionSet.JordanWignerGeneratorSystem; +import JordanWigner.JordanWignerVQE.EstimateTermExpectation; +import JordanWigner.JordanWignerVQE.ExpandedCoefficients; +import JordanWigner.JordanWignerVQE.MeasurementOperators; +import JordanWigner.StatePreparation.PrepareTrialState; +import JordanWigner.Utils.JordanWignerEncodingData; +import JordanWigner.Utils.JordanWignerInputState; +import JordanWigner.Utils.JWOptimizedHTerms; + +/// # Summary +/// Estimates the energy of the molecule by summing the energy contributed by the individual Jordan-Wigner terms. +/// This convenience wrapper takes input in raw data types and converts them to the Jordan-Wigner encoding data type. +/// +/// # Description +/// This operation implicitly relies on the spin up-down indexing convention. +/// +/// # Input +/// ## jwHamiltonian +/// The Jordan-Wigner Hamiltonian, represented in simple data types rather than a struct. +/// ## nSamples +/// The number of samples to use for the estimation of the term expectations. +/// +/// # Output +/// The estimated energy of the molecule +operation EstimateEnergyWrapper(jwHamiltonian : (Int, ((Int[], Double[])[], (Int[], Double[])[], (Int[], Double[])[], (Int[], Double[])[]), (Int, ((Double, Double), Int[])[]), Double), nSamples : Int) : Double { + let (nQubits, jwTerms, inputState, energyOffset) = jwHamiltonian; + let (hterm0, hterm1, hterm2, hterm3) = jwTerms; + let jwTerms = new JWOptimizedHTerms { HTerm0 = hterm0, HTerm1 = hterm1, HTerm2 = hterm2, HTerm3 = hterm3 }; + let (inputState1, inputState2) = inputState; + mutable jwInputState = []; + for entry in inputState2 { + let (amp, idicies) = entry; + jwInputState += [new JordanWignerInputState { Amplitude = amp, FermionIndices = idicies }]; + } + let inputState = (inputState1, jwInputState); + let jwHamiltonian = new JordanWignerEncodingData { + NumQubits = nQubits, + Terms = jwTerms, + InputState = inputState, + EnergyOffset = energyOffset + }; + return EstimateEnergy(jwHamiltonian, nSamples); +} + +/// # Summary +/// Estimates the energy of the molecule by summing the energy contributed by the individual Jordan-Wigner terms. +/// +/// # Description +/// This operation implicitly relies on the spin up-down indexing convention. +/// +/// # Input +/// ## jwHamiltonian +/// The Jordan-Wigner Hamiltonian. +/// ## nSamples +/// The number of samples to use for the estimation of the term expectations. +/// +/// # Output +/// The estimated energy of the molecule +operation EstimateEnergy(jwHamiltonian : JordanWignerEncodingData, nSamples : Int) : Double { + // Initialize return value + mutable energy = 0.; + + // Unpack information and qubit Hamiltonian terms + let (nQubits, jwTerms, inputState, energyOffset) = jwHamiltonian!; + + // Loop over all qubit Hamiltonian terms + let (nTerms, indexFunction) = (JordanWignerGeneratorSystem(jwTerms))!; + + for idxTerm in 0..nTerms - 1 { + let term = indexFunction(idxTerm); + let ((idxTermType, coeff), idxFermions) = term!; + let termType = idxTermType[0]; + + let ops = MeasurementOperators(nQubits, idxFermions, termType); + let coeffs = ExpandedCoefficients(coeff, termType); + + // The private wrapper enables fast emulation during expectation estimation + let inputStateUnitary = PrepareTrialState(inputState, _); + + let jwTermEnergy = EstimateTermExpectation(inputStateUnitary, ops, coeffs, nQubits, nSamples); + energy += jwTermEnergy; + } + + return energy + energyOffset; +} diff --git a/library/chemistry/src/JordanWigner/JordanWignerBlockEncoding.qs b/library/chemistry/src/JordanWigner/JordanWignerBlockEncoding.qs new file mode 100644 index 0000000000..2e6652e18d --- /dev/null +++ b/library/chemistry/src/JordanWigner/JordanWignerBlockEncoding.qs @@ -0,0 +1,211 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. + +export JordanWignerBlockEncodingGeneratorSystem; + +import JordanWigner.Utils.JWOptimizedHTerms; +import JordanWigner.Utils.RangeAsIntArray; +import Std.Arrays.IndexRange; +import Utils.GeneratorIndex; +import Utils.GeneratorSystem; +import Utils.HTermToGenIdx; +import Utils.IsNotZero; + +// This block encoding for qubitization runs off data optimized for a Jordan–Wigner encoding. +// This collects terms Z, ZZ, PQandPQQR, hpqrs separately. +// This only apples the needed hpqrs XXXX XXYY terms. + +// Convention for GeneratorIndex = ((Int[],Double[]), Int[]) +// We index single Paulis as 0 for I, 1 for X, 2 for Y, 3 for Z. +// We index Pauli strings with arrays of integers e.g. a = [3,1,1,2] for ZXXY. +// We assume the zeroth element of Double[] is the angle of rotation +// We index the qubits that Pauli strings act on with arrays of integers e.g. q = [2,4,5,8] for Z_2 X_4 X_5, Y_8 +// An example of a Pauli string GeneratorIndex is thus ((a,b), q) + +// Consider the Hamiltonian H = 0.1 XI + 0.2 IX + 0.3 ZY +// Its GeneratorTerms are (([1],b),[0]), 0.1), (([1],b),[1]), 0.2), (([3,2],b),[0,1]), 0.3). + +/// # Summary +/// Converts a `GeneratorIndex` describing a Z term to +/// an expression `GeneratorIndex[]` in terms of Paulis. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a Z term. +/// +/// # Output +/// `GeneratorIndex[]` expressing Z term as Pauli terms. +function ZTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] { + let (_, coeff) = term.Term; + return [new GeneratorIndex { Term = ([3], coeff), Subsystem = term.Subsystem }]; +} + +/// # Summary +/// Converts a `GeneratorIndex` describing a ZZ term to +/// an expression `GeneratorIndex[]` in terms of Paulis. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a ZZ term. +/// +/// # Output +/// `GeneratorIndex[]` expressing ZZ term as Pauli terms. +function ZZTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] { + let (_, coeff) = term.Term; + return [new GeneratorIndex { Term = ([3, 3], coeff), Subsystem = term.Subsystem }]; +} + +/// # Summary +/// Converts a `GeneratorIndex` describing a PQ term to +/// an expression `GeneratorIndex[]` in terms of Paulis +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a PQ term. +/// +/// # Output +/// `GeneratorIndex[]` expressing PQ term as Pauli terms. +function PQTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] { + let (_, coeff) = term.Term; + let newCoeff = [coeff[0]]; + let qubitPidx = term.Subsystem[0]; + let qubitQidx = term.Subsystem[1]; + let qubitIndices = RangeAsIntArray(qubitPidx..qubitQidx); + return [ + new GeneratorIndex { Term = (([1] + Repeated(3, Length(qubitIndices) - 2)) + [1], newCoeff), Subsystem = qubitIndices }, + new GeneratorIndex { Term = (([2] + Repeated(3, Length(qubitIndices) - 2)) + [2], newCoeff), Subsystem = qubitIndices } + ]; +} + +/// # Summary +/// Converts a `GeneratorIndex` describing a PQ or PQQR term to +/// an expression `GeneratorIndex[]` in terms of Paulis +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a PQ or PQQR term. +/// +/// # Output +/// `GeneratorIndex[]` expressing PQ or PQQR term as Pauli terms. +function PQandPQQRTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] { + let (_, coeff) = term.Term; + let newCoeff = [coeff[0]]; + + if Length(term.Subsystem) == 2 { + return PQTermToPauliGenIdx(term); + } else { + let qubitPidx = term.Subsystem[0]; + let qubitQidx = term.Subsystem[1]; + let qubitRidx = term.Subsystem[3]; + + if (qubitPidx < qubitQidx and qubitQidx < qubitRidx) { + + // Apply XZ..ZIZ..ZX + let qubitIndices = RangeAsIntArray(qubitPidx..qubitQidx - 1) + RangeAsIntArray(qubitQidx + 1..qubitRidx); + return [ + new GeneratorIndex { Term = (([1] + Repeated(3, Length(qubitIndices) - 2)) + [1], newCoeff), Subsystem = qubitIndices }, + new GeneratorIndex { Term = (([2] + Repeated(3, Length(qubitIndices) - 2)) + [2], newCoeff), Subsystem = qubitIndices } + ]; + } else { + + // Apply ZI..IXZ..ZX or XZ..ZXI..IZ + let qubitIndices = RangeAsIntArray(qubitPidx..qubitRidx) + [qubitQidx]; + return [ + new GeneratorIndex { Term = (([1] + Repeated(3, Length(qubitIndices) - 3)) + [1, 3], newCoeff), Subsystem = qubitIndices }, + new GeneratorIndex { Term = (([2] + Repeated(3, Length(qubitIndices) - 3)) + [2, 3], newCoeff), Subsystem = qubitIndices } + ]; + } + } +} + +/// # Summary +/// Converts a `GeneratorIndex` describing a PQRS term to +/// an expression `GeneratorIndex[]` in terms of Paulis +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a PQRS term. +/// +/// # Output +/// `GeneratorIndex[]` expressing PQRS term as Pauli terms. +function V0123TermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] { + let (_, v0123) = term.Term; + let qubitsPQ = term.Subsystem[0..1]; + let qubitsRS = term.Subsystem[2..3]; + let qubitsPQJW = RangeAsIntArray(qubitsPQ[0] + 1..qubitsPQ[1] - 1); + let qubitsRSJW = RangeAsIntArray(qubitsRS[0] + 1..qubitsRS[1] - 1); + let ops = [[1, 1, 1, 1], [1, 1, 2, 2], [1, 2, 1, 2], [2, 1, 1, 2], [2, 2, 2, 2], [2, 2, 1, 1], [2, 1, 2, 1], [1, 2, 2, 1]]; + mutable genIdxes = Repeated(new GeneratorIndex { Term = ([0], [0.0]), Subsystem = [0] }, 8); + mutable nonZero = 0; + + for idxOp in IndexRange(ops) { + if (IsNotZero(v0123[idxOp % 4])) { + let newCoeff = [v0123[idxOp % 4]]; + genIdxes w/= nonZero <- new GeneratorIndex { + Term = (ops[idxOp] + Repeated(3, Length(qubitsPQJW) + Length(qubitsRSJW)), newCoeff), + Subsystem = ((qubitsPQ + qubitsRS) + qubitsPQJW) + qubitsRSJW + }; + nonZero = nonZero + 1; + } + } + + return genIdxes[0..nonZero - 1]; +} + +/// # Summary +/// Converts a Hamiltonian described by `JWOptimizedHTerms` +/// to a `GeneratorSystem` expressed in terms of the Pauli +/// `GeneratorIndex`. +/// +/// # Input +/// ## data +/// Description of Hamiltonian in `JWOptimizedHTerms` format. +/// +/// # Output +/// Representation of Hamiltonian as `GeneratorSystem`. +function JordanWignerBlockEncodingGeneratorSystem(data : JWOptimizedHTerms) : GeneratorSystem { + let (ZData, ZZData, PQandPQQRData, h0123Data) = data!; + mutable genIdxes = Repeated( + new GeneratorIndex { Term = ([0], [0.0]), Subsystem = [0] }, + ((Length(ZData) + Length(ZZData)) + 2 * Length(PQandPQQRData)) + 8 * Length(h0123Data) + ); + mutable startIdx = 0; + + for idx in IndexRange(ZData) { + // Array of Arrays of Length 1 + genIdxes w/= idx <- (ZTermToPauliGenIdx(HTermToGenIdx(ZData[idx], [0])))[0]; + } + + startIdx = Length(ZData); + + for idx in IndexRange(ZZData) { + // Array of Arrays of Length 1 + genIdxes w/= startIdx + idx <- (ZZTermToPauliGenIdx(HTermToGenIdx(ZZData[idx], [1])))[0]; + } + + startIdx = startIdx + Length(ZZData); + + for idx in IndexRange(PQandPQQRData) { + + // Array of Arrays of Length 2 + let genArr = PQandPQQRTermToPauliGenIdx(HTermToGenIdx(PQandPQQRData[idx], [2])); + genIdxes w/= startIdx + 2 * idx <- genArr[0]; + genIdxes w/= (startIdx + 2 * idx) + 1 <- genArr[1]; + } + + startIdx = startIdx + 2 * Length(PQandPQQRData); + mutable finalIdx = startIdx; + + for idx in 0..Length(h0123Data) - 1 { + // Array of Arrays of Length up to 8 + let genArr = V0123TermToPauliGenIdx(HTermToGenIdx(h0123Data[idx], [3])); + + for idx0123 in IndexRange(genArr) { + genIdxes w/= finalIdx <- genArr[idx0123]; + finalIdx += 1; + } + } + + let genIdxes = genIdxes[0..finalIdx - 1]; + return new GeneratorSystem { NumEntries = finalIdx, EntryAt = idx -> genIdxes[idx] }; +} diff --git a/library/chemistry/src/JordanWigner/JordanWignerClusterOperatorEvolutionSet.qs b/library/chemistry/src/JordanWigner/JordanWignerClusterOperatorEvolutionSet.qs new file mode 100644 index 0000000000..b7ecdd9fe3 --- /dev/null +++ b/library/chemistry/src/JordanWigner/JordanWignerClusterOperatorEvolutionSet.qs @@ -0,0 +1,252 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. + +export + JordanWignerClusterOperatorEvolutionSet, + JordanWignerClusterOperatorGeneratorSystem; + +import JordanWigner.Utils.JordanWignerInputState; +import Std.Arrays.IndexRange; +import Std.Math.Max; +import Std.Math.Min; +import Utils.GeneratorIndex; +import Utils.GeneratorSystem; + +/// # Summary +/// Computes Z component of Jordan–Wigner string between +/// fermion indices in a fermionic operator with an even +/// number of creation / annihilation operators. +/// +/// # Input +/// ## nFermions +/// The Number of fermions in the system. +/// ## idxFermions +/// fermionic operator indices. +/// +/// # Output +/// Bitstring `Bool[]` that is `true` where a `PauliZ` should be applied. +/// +/// # Example +/// ```qsharp +/// let bitString = ComputeJordanWignerBitString(6, [0,1,2,6]) ; +/// // bitString is [false, false, false ,true, true, true, false]. +/// ``` +function ComputeJordanWignerBitString(nFermions : Int, idxFermions : Int[]) : Bool[] { + if Length(idxFermions) % 2 != 0 { + fail $"ComputeJordanWignerString failed. `idxFermions` must contain an even number of terms."; + } + + mutable zString = [false, size = nFermions]; + for fermionIdx in idxFermions { + if fermionIdx >= nFermions { + fail $"ComputeJordanWignerString failed. fermionIdx {fermionIdx} out of range."; + } + for idx in 0..fermionIdx { + zString w/= idx <- not zString[idx]; + } + } + + for fermionIdx in idxFermions { + zString w/= fermionIdx <- false; + } + return zString; +} + +// Identical to `ComputeJordanWignerBitString`, except with the map +// false -> PauliI and true -> PauliZ +function ComputeJordanWignerPauliZString(nFermions : Int, idxFermions : Int[]) : Pauli[] { + let bitString = ComputeJordanWignerBitString(nFermions, idxFermions); + mutable pauliString = Repeated(PauliI, Length(bitString)); + for idx in IndexRange(bitString) { + if bitString[idx] { + pauliString w/= idx <- PauliZ; + } + } + return pauliString; +} + +// Identical to `ComputeJordanWignerPauliZString`, except that some +// specified elements are substituted. +function ComputeJordanWignerPauliString(nFermions : Int, idxFermions : Int[], pauliReplacements : Pauli[]) : Pauli[] { + mutable pauliString = ComputeJordanWignerPauliZString(nFermions, idxFermions); + + for idx in IndexRange(idxFermions) { + let idxFermion = idxFermions[idx]; + let op = pauliReplacements[idx]; + pauliString w/= idxFermion <- op; + } + + return pauliString; +} + +/// # Summary +/// Applies time-evolution by a cluster operator PQ term described by a `GeneratorIndex`. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a cluster operator PQ term. +/// ## stepSize +/// Duration of time-evolution. +/// ## qubits +/// Qubits of Hamiltonian. +operation ApplyJordanWignerClusterOperatorPQTerm(term : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + let ((idxTermType, coeff), idxFermions) = term!; + let p = idxFermions[0]; + let q = idxFermions[1]; + if p == q { + fail $"Unitary coupled-cluster PQ failed: indices {p}, {q} must be distinct"; + } + let angle = 0.5 * coeff[0] * stepSize; + let ops = [[PauliX, PauliY], [PauliY, PauliX]]; + let signs = [+ 1.0, -1.0]; + + for i in IndexRange(ops) { + let pauliString = ComputeJordanWignerPauliString(Length(qubits), idxFermions, ops[i]); + Exp(pauliString, signs[i] * angle, qubits); + } +} + + +/// # Summary +/// Applies time-evolution by a cluster operator PQRS term described by a `GeneratorIndex`. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a cluster operator PQRS term. +/// ## stepSize +/// Duration of time-evolution. +/// ## qubits +/// Qubits of Hamiltonian. +operation ApplyJordanWignerClusterOperatorPQRSTerm(term : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + let ((idxTermType, coeff), idxFermions) = term!; + let p = idxFermions[0]; + let q = idxFermions[1]; + let r = idxFermions[2]; + let s = idxFermions[3]; + let angle = 0.125 * coeff[0] * stepSize; + + if p == q or p == r or p == s or q == r or q == s or r == s { + fail ($"Unitary coupled-cluster PQRS failed: indices {p}, {q}, {r}, {s} must be distinct"); + } + + let x = PauliX; + let y = PauliY; + + let ops = [[y, y, x, y], [x, x, x, y], [x, y, y, y], [y, x, y, y], [x, y, x, x], [y, x, x, x], [y, y, y, x], [x, x, y, x]]; + let (sortedIndices, signs, globalSign) = JordanWignerClusterOperatorPQRSTermSigns([p, q, r, s]); + + for i in IndexRange(ops) { + let pauliString = ComputeJordanWignerPauliString(Length(qubits), sortedIndices, ops[i]); + Exp(pauliString, globalSign * signs[i] * angle, qubits); + } +} + +function JordanWignerClusterOperatorPQRSTermSigns(indices : Int[]) : (Int[], Double[], Double) { + let p = indices[0]; + let q = indices[1]; + let r = indices[2]; + let s = indices[3]; + mutable sorted = [0, size = 4]; + mutable signs = [0.0, size = 8]; + mutable sign = 1.0; + + if (p > q) { + sign = sign * -1.0; + } + if (r > s) { + sign = sign * -1.0; + } + if (Min([p, q]) > Min([r, s])) { + sign = sign * -1.0; + sorted = [Min([r, s]), Max([r, s]), Min([p, q]), Max([p, q])]; + } else { + sorted = [Min([p, q]), Max([p, q]), Min([r, s]), Max([r, s])]; + } + // sorted is now in the order + // [p`,q`,r`,s`], where p` (r,s) + if (q1 < r1) { + // p1 < q1 < r1 < s1 + return ([p1, q1, r1, s1], [1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0], sign); + } + // Case interleaved + elif (q1 > r1 and q1 < s1) { + // p1 < r1 < q1 < s1 + return ([p1, r1, q1, s1], [-1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0], sign); + } + // Case contained + elif (q1 > r1 and q1 > s1) { + // p1 < r1 < s1 < q1 + return ([p1, r1, s1, q1], [1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0], sign); + } else { + fail ("Completely invalid cluster operator specified."); + } +} + +/// # Summary +/// Converts a Hamiltonian described by `JWOptimizedHTerms` +/// to a `GeneratorSystem` expressed in terms of the +/// `GeneratorIndex` convention defined in this file. +/// +/// # Input +/// ## data +/// Description of Hamiltonian in `JWOptimizedHTerms` format. +/// +/// # Output +/// Representation of Hamiltonian as `GeneratorSystem`. +function JordanWignerClusterOperatorGeneratorSystem(data : JordanWignerInputState[]) : GeneratorSystem { + new GeneratorSystem { NumEntries = Length(data), EntryAt = JordanWignerStateAsGeneratorIndex(data, _) } +} + +function JordanWignerStateAsGeneratorIndex(data : JordanWignerInputState[], idx : Int) : GeneratorIndex { + let ((real, imaginary), idxFermions) = data[idx]!; + + if Length(idxFermions) == 2 { + // PQ term + new GeneratorIndex { Term = ([0], [real]), Subsystem = idxFermions } + } elif Length(idxFermions) == 4 { + // PQRS term + new GeneratorIndex { Term = ([2], [real]), Subsystem = idxFermions } + } else { + // Any other term in invalid + new GeneratorIndex { Term = ([-1], [0.0]), Subsystem = [0] } + } +} + +/// # Summary +/// Represents a dynamical generator as a set of simulatable gates and an +/// expansion in the JordanWigner basis. +/// +/// # Input +/// ## generatorIndex +/// A generator index to be represented as unitary evolution in the JordanWigner. +/// ## stepSize +/// Dummy variable to match signature of simulation algorithms. +/// ## qubits +/// Register acted upon by time-evolution operator. +operation JordanWignerClusterOperatorImpl(generatorIndex : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + let ((idxTermType, idxDoubles), idxFermions) = generatorIndex!; + let termType = idxTermType[0]; + + if termType == 0 { + ApplyJordanWignerClusterOperatorPQTerm(generatorIndex, stepSize, qubits); + } elif termType == 2 { + ApplyJordanWignerClusterOperatorPQRSTerm(generatorIndex, stepSize, qubits); + } +} + +/// # Summary +/// Represents a dynamical generator as a set of simulatable gates and an +/// expansion in the JordanWigner basis. +/// +/// # Output +/// An evolution set function that maps a `GeneratorIndex` for the JordanWigner basis to +/// an evolution unitary operation. +function JordanWignerClusterOperatorEvolutionSet() : GeneratorIndex -> (Double, Qubit[]) => Unit is Adj + Ctl { + generatorIndex -> (stepSize, qubits) => JordanWignerClusterOperatorImpl(generatorIndex, stepSize, qubits) +} diff --git a/library/chemistry/src/JordanWigner/JordanWignerEvolutionSet.qs b/library/chemistry/src/JordanWigner/JordanWignerEvolutionSet.qs new file mode 100644 index 0000000000..cdcbabcec5 --- /dev/null +++ b/library/chemistry/src/JordanWigner/JordanWignerEvolutionSet.qs @@ -0,0 +1,240 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +export + JordanWignerGeneratorSystem, + JordanWignerFermionEvolutionSet; + +import JordanWigner.Utils.JWOptimizedHTerms; +import Std.Arrays.IndexRange; +import Std.Arrays.Subarray; +import Utils.GeneratorIndex; +import Utils.GeneratorSystem; +import Utils.HTermsToGenSys; +import Utils.IsNotZero; + +// This evolution set runs off data optimized for a Jordan–Wigner encoding. +// This collects terms Z, ZZ, PQandPQQR, hpqrs separately. +// This only apples the needed hpqrs XXXX XXYY terms. +// Operations here are expressed in terms of Exp([...]) + +// Convention for GeneratorIndex = ((Int[],Double[]), Int[]) +// We index single Paulis as 0 for I, 1 for X, 2 for Y, 3 for Z. +// We index Pauli strings with arrays of integers e.g. a = [3,1,1,2] for ZXXY. +// We assume the zeroth element of Double[] is the angle of rotation +// We index the qubits that Pauli strings act on with arrays of integers e.g. q = [2,4,5,8] for Z_2 X_4 X_5, Y_8 +// An example of a Pauli string GeneratorIndex is thus ((a,b), q) + +// Consider the Hamiltonian H = 0.1 XI + 0.2 IX + 0.3 ZY +// Its GeneratorTerms are (([1],b),[0]), 0.1), (([1],b),[1]), 0.2), (([3,2],b),[0,1]), 0.3). + +/// # Summary +/// Applies time-evolution by a Z term described by a `GeneratorIndex`. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a Z term. +/// ## stepSize +/// Duration of time-evolution. +/// ## qubits +/// Qubits of Hamiltonian. +operation ApplyJordanWignerZTerm(term : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + let (_, coeff) = term.Term; + let angle = (1.0 * coeff[0]) * stepSize; + let qubit = qubits[term.Subsystem[0]]; + Exp([PauliZ], angle, [qubit]); +} + + +/// # Summary +/// Applies time-evolution by a ZZ term described by a `GeneratorIndex`. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a ZZ term. +/// ## stepSize +/// Duration of time-evolution. +/// ## qubits +/// Qubits of Hamiltonian. +operation ApplyJordanWignerZZTerm(term : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + let (_, coeff) = term.Term; + let angle = (1.0 * coeff[0]) * stepSize; + let qubitsZZ = Subarray(term.Subsystem[0..1], qubits); + Exp([PauliZ, PauliZ], angle, qubitsZZ); +} + + +/// # Summary +/// Applies time-evolution by a PQ term described by a `GeneratorIndex`. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a PQ term. +/// ## stepSize +/// Duration of time-evolution. +/// ## extraParityQubits +/// Optional parity qubits that flip the sign of time-evolution. +/// ## qubits +/// Qubits of Hamiltonian. +operation ApplyJordanWignerPQTerm(term : GeneratorIndex, stepSize : Double, extraParityQubits : Qubit[], qubits : Qubit[]) : Unit is Adj + Ctl { + let ((idxTermType, coeff), idxFermions) = term!; + let angle = (1.0 * coeff[0]) * stepSize; + let qubitsPQ = Subarray(idxFermions[0..1], qubits); + let qubitsJW = qubits[idxFermions[0] + 1..idxFermions[1] - 1]; + let ops = [[PauliX, PauliX], [PauliY, PauliY]]; + let padding = Repeated(PauliZ, Length(qubitsJW) + Length(extraParityQubits)); + + for op in ops { + Exp(op + padding, angle, qubitsPQ + qubitsJW + extraParityQubits); + } +} + + +/// # Summary +/// Applies time-evolution by a PQ or PQQR term described by a `GeneratorIndex`. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a PQ or PQQR term. +/// ## stepSize +/// Duration of time-evolution. +/// ## qubits +/// Qubits of Hamiltonian. +operation ApplyJordanWignerPQandPQQRTerm(term : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + let ((idxTermType, coeff), idxFermions) = term!; + let angle = (1.0 * coeff[0]) * stepSize; + let qubitQidx = idxFermions[1]; + + // For all cases, do the same thing: + // p < r < q (1/4)(1-Z_q)(Z_{r-1,p+1})(X_p X_r + Y_p Y_r) (same as Hermitian conjugate of r < p < q) + // q < p < r (1/4)(1-Z_q)(Z_{r-1,p+1})(X_p X_r + Y_p Y_r) + // p < q < r (1/4)(1-Z_q)(Z_{r-1,p+1})(X_p X_r + Y_p Y_r) + + // This amounts to applying a PQ term, followed by same PQ term after a CNOT from q to the parity bit. + if Length(idxFermions) == 2 { + let termPR0 = new GeneratorIndex { Term = (idxTermType, [1.0]), Subsystem = idxFermions }; + ApplyJordanWignerPQTerm(termPR0, angle, [], qubits); + } else { + if idxFermions[0] < qubitQidx and qubitQidx < idxFermions[3] { + let termPR1 = new GeneratorIndex { Term = (idxTermType, [1.0]), Subsystem = [idxFermions[0], idxFermions[3] - 1] }; + let excludingQ = if qubitQidx > 0 { qubits[0..qubitQidx-1] + qubits[qubitQidx + 1...] } else { qubits[1...] }; + ApplyJordanWignerPQTerm(termPR1, angle, [], excludingQ); + } else { + let termPR1 = new GeneratorIndex { Term = (idxTermType, [1.0]), Subsystem = [0, idxFermions[3] - idxFermions[0]] }; + ApplyJordanWignerPQTerm(termPR1, angle, [qubits[qubitQidx]], qubits[idxFermions[0]..idxFermions[3]]); + } + } +} + + +/// # Summary +/// Applies time-evolution by a PQRS term described by a given index. +/// +/// # Input +/// ## term +/// The index representing a PQRS term to be applied. +/// ## stepSize +/// Duration of time-evolution. +/// ## qubits +/// Qubits to apply the given term to. +operation ApplyJordanWigner0123Term(term : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + let ((idxTermType, v0123), idxFermions) = term!; + let angle = stepSize; + let qubitsPQ = Subarray(idxFermions[0..1], qubits); + let qubitsRS = Subarray(idxFermions[2..3], qubits); + let qubitsPQJW = qubits[idxFermions[0] + 1..idxFermions[1] - 1]; + let qubitsRSJW = qubits[idxFermions[2] + 1..idxFermions[3] - 1]; + let ops = [[PauliX, PauliX, PauliX, PauliX], [PauliX, PauliX, PauliY, PauliY], [PauliX, PauliY, PauliX, PauliY], [PauliY, PauliX, PauliX, PauliY], [PauliY, PauliY, PauliY, PauliY], [PauliY, PauliY, PauliX, PauliX], [PauliY, PauliX, PauliY, PauliX], [PauliX, PauliY, PauliY, PauliX]]; + + for idxOp in IndexRange(ops) { + if (IsNotZero(v0123[idxOp % 4])) { + Exp(ops[idxOp] + Repeated(PauliZ, Length(qubitsPQJW) + Length(qubitsRSJW)), angle * v0123[idxOp % 4], ((qubitsPQ + qubitsRS) + qubitsPQJW) + qubitsRSJW); + } + } +} + + +/// # Summary +/// Converts a Hamiltonian described by `JWOptimizedHTerms` +/// to a `GeneratorSystem` expressed in terms of the +/// `GeneratorIndex` convention defined in this file. +/// +/// # Input +/// ## data +/// Description of Hamiltonian in `JWOptimizedHTerms` format. +/// +/// # Output +/// Representation of Hamiltonian as `GeneratorSystem`. +function JordanWignerGeneratorSystem(data : JWOptimizedHTerms) : GeneratorSystem { + let (ZData, ZZData, PQandPQQRData, h0123Data) = data!; + let ZGenSys = HTermsToGenSys(ZData, [0]); + let ZZGenSys = HTermsToGenSys(ZZData, [1]); + let PQandPQQRGenSys = HTermsToGenSys(PQandPQQRData, [2]); + let h0123GenSys = HTermsToGenSys(h0123Data, [3]); + let sum = AddGeneratorSystems(ZGenSys, ZZGenSys); + let sum = AddGeneratorSystems(sum, PQandPQQRGenSys); + return AddGeneratorSystems(sum, h0123GenSys); +} + +/// # Summary +/// Adds two `GeneratorSystem`s to create a new `GeneratorSystem`. +/// +/// # Input +/// ## generatorSystemA +/// The first `GeneratorSystem`. +/// ## generatorSystemB +/// The second `GeneratorSystem`. +/// +/// # Output +/// A `GeneratorSystem` representing a system that is the sum of the +/// input generator systems. +function AddGeneratorSystems(generatorSystemA : GeneratorSystem, generatorSystemB : GeneratorSystem) : GeneratorSystem { + let nTermsA = generatorSystemA.NumEntries; + let nTermsB = generatorSystemB.NumEntries; + let generatorIndexFunction = idx -> { + if idx < nTermsA { + return generatorSystemA.EntryAt(idx); + } else { + return generatorSystemB.EntryAt(idx - nTermsA); + } + }; + return new GeneratorSystem { NumEntries = nTermsA + nTermsB, EntryAt = generatorIndexFunction }; +} + +/// # Summary +/// Represents a dynamical generator as a set of simulatable gates and an +/// expansion in the JordanWigner basis. +/// +/// # Input +/// ## generatorIndex +/// A generator index to be represented as unitary evolution in the JordanWigner. +/// ## stepSize +/// A multiplier on the duration of time-evolution by the term referenced +/// in `generatorIndex`. +/// ## qubits +/// Register acted upon by time-evolution operator. +operation JordanWignerFermionImpl(generatorIndex : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + let ((idxTermType, idxDoubles), idxFermions) = generatorIndex!; + let termType = idxTermType[0]; + + if (termType == 0) { + ApplyJordanWignerZTerm(generatorIndex, stepSize, qubits); + } elif (termType == 1) { + ApplyJordanWignerZZTerm(generatorIndex, stepSize, qubits); + } elif (termType == 2) { + ApplyJordanWignerPQandPQQRTerm(generatorIndex, stepSize, qubits); + } elif (termType == 3) { + ApplyJordanWigner0123Term(generatorIndex, stepSize, qubits); + } +} + +/// # Summary +/// Represents a dynamical generator as a set of simulatable gates and an +/// expansion in the JordanWigner basis. +/// +/// # Output +/// An evolution set function that maps a `GeneratorIndex` for the JordanWigner basis to +/// an evolution unitary operation. +function JordanWignerFermionEvolutionSet() : GeneratorIndex -> (Double, Qubit[]) => Unit is Adj + Ctl { + generatorIndex -> (stepSize, qubits) => JordanWignerFermionImpl(generatorIndex, stepSize, qubits) +} diff --git a/library/chemistry/src/JordanWigner/JordanWignerOptimizedBlockEncoding.qs b/library/chemistry/src/JordanWigner/JordanWignerOptimizedBlockEncoding.qs new file mode 100644 index 0000000000..545e599878 --- /dev/null +++ b/library/chemistry/src/JordanWigner/JordanWignerOptimizedBlockEncoding.qs @@ -0,0 +1,1128 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. + +export + OptimizedBETermIndex, + OptimizedBEGeneratorSystem, + OptimizedBlockEncodingGeneratorSystem, + MixedStatePreparation, + MixedStatePreparationRequirements, + PurifiedMixedState, + ObliviousAmplitudeAmplificationFromStatePreparation, + ObliviousAmplitudeAmplificationFromPartialReflections, + ApplyObliviousAmplitudeAmplification, + PurifiedMixedStateRequirements, + BlockEncodingByLCU, + QuantumWalkByQubitization, + PauliBlockEncoding; + +import JordanWigner.OptimizedBEOperator.JordanWignerSelect; +import JordanWigner.OptimizedBEOperator.JordanWignerSelectQubitCount; +import JordanWigner.OptimizedBEOperator.JordanWignerSelectQubitManager; +import JordanWigner.Utils.JWOptimizedHTerms; +import JordanWigner.Utils.MultiplexOperationsFromGenerator; +import JordanWigner.Utils.RangeAsIntArray; +import Std.Arrays.*; +import Std.Math.*; +import Std.Convert.IntAsDouble; +import Std.Arithmetic.ApplyIfGreaterLE; +import Std.StatePreparation.PreparePureStateD; +import Std.Diagnostics.Fact; +import Utils.GeneratorIndex; +import Utils.GeneratorSystem; +import Utils.HTermToGenIdx; +import Utils.IsNotZero; + +/// # Summary +/// Term data in the optimized block-encoding algorithm. +struct OptimizedBETermIndex { + Coefficient : Double, + UseSignQubit : Bool, + ZControlRegisterMask : Bool[], + OptimizedControlRegisterMask : Bool[], + PauliBases : Int[], + RegisterIndices : Int[], +} + +/// # Summary +/// Function that returns `OptimizedBETermIndex` data for term `n` given an +/// integer `n`, together with the number of terms in the first `Int` and +/// the sum of absolute-values of all term coefficients in the `Double`. +struct OptimizedBEGeneratorSystem { + NumTerms : Int, + Norm : Double, + SelectTerm : (Int -> OptimizedBETermIndex) +} + +// Get OptimizedBETermIndex coefficients +function GetOptimizedBETermIndexCoeff(term : OptimizedBETermIndex) : Double { + let (a, b, c, d, e, f) = term!; + return a; +} + + +// Get OptimizedBEGeneratorSystem coefficients +function OptimizedBEGeneratorSystemCoeff(optimizedBEGeneratorSystem : OptimizedBEGeneratorSystem) : Double[] { + let (nTerms, oneNorm, intToGenIdx) = optimizedBEGeneratorSystem!; + mutable coefficients = [0.0, size = nTerms]; + + for idx in 0..nTerms - 1 { + coefficients w/= idx <- GetOptimizedBETermIndexCoeff(intToGenIdx(idx)); + } + + return coefficients; +} + + +/// # Summary +/// Converts a `GeneratorIndex` describing a Z term to +/// an expression `GeneratorIndex[]` in terms of Paulis. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a Z term. +/// +/// # Output +/// `GeneratorIndex[]` expressing Z term as Pauli terms. +function ZTermToPauliMajIdx(term : GeneratorIndex) : OptimizedBETermIndex { + let ((idxTermType, coeff), idxFermions) = term!; + let signQubit = coeff[0] < 0.0; + let selectZControlRegisters = [true]; + let optimizedBEControlRegisters = []; + let pauliBases = []; + let indexRegisters = idxFermions; + return new OptimizedBETermIndex { + Coefficient = coeff[0], + UseSignQubit = signQubit, + ZControlRegisterMask = selectZControlRegisters, + OptimizedControlRegisterMask = optimizedBEControlRegisters, + PauliBases = pauliBases, + RegisterIndices = indexRegisters + }; +} + + +/// # Summary +/// Converts a GeneratorIndex describing a ZZ term to +/// an expression `GeneratorIndex[]` in terms of Paulis. +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a ZZ term. +/// +/// # Output +/// `GeneratorIndex[]` expressing ZZ term as Pauli terms. +function ZZTermToPauliMajIdx(term : GeneratorIndex) : OptimizedBETermIndex { + let ((idxTermType, coeff), idxFermions) = term!; + let signQubit = coeff[0] < 0.0; + let selectZControlRegisters = [true, true]; + let optimizedBEControlRegisters = []; + let pauliBases = []; + let indexRegisters = idxFermions; + return new OptimizedBETermIndex { + Coefficient = 2.0 * coeff[0], + UseSignQubit = signQubit, + ZControlRegisterMask = selectZControlRegisters, + OptimizedControlRegisterMask = optimizedBEControlRegisters, + PauliBases = pauliBases, + RegisterIndices = indexRegisters + }; +} + + +/// # Summary +/// Converts a `GeneratorIndex` describing a PQ term to +/// an expression `GeneratorIndex[]` in terms of Paulis +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a PQ term. +/// +/// # Output +/// `GeneratorIndex[]` expressing PQ term as Pauli terms. +function PQTermToPauliMajIdx(term : GeneratorIndex) : OptimizedBETermIndex { + let ((idxTermType, coeff), idxFermions) = term!; + let sign = coeff[0] < 0.0; + let selectZControlRegisters = []; + let optimizedBEControlRegisters = [true, true]; + let pauliBases = [1, 2]; + let indexRegisters = idxFermions; + return new OptimizedBETermIndex { + Coefficient = 2.0 * coeff[0], + UseSignQubit = sign, + ZControlRegisterMask = selectZControlRegisters, + OptimizedControlRegisterMask = optimizedBEControlRegisters, + PauliBases = pauliBases, + RegisterIndices = indexRegisters + }; +} + + +/// # Summary +/// Converts a `GeneratorIndex` describing a PQ or PQQR term to +/// an expression `GeneratorIndex[]` in terms of Paulis +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a PQ or PQQR term. +/// +/// # Output +/// `GeneratorIndex[]` expressing PQ or PQQR term as Pauli terms. +function PQandPQQRTermToPauliMajIdx(term : GeneratorIndex) : OptimizedBETermIndex { + let ((idxTermType, coeff), idxFermions) = term!; + let sign = coeff[0] < 0.0; + + if Length(idxFermions) == 2 { + return PQTermToPauliMajIdx(term); + } else { + let qubitPidx = idxFermions[0]; + let qubitQidx = idxFermions[1]; + let qubitRidx = idxFermions[3]; + let selectZControlRegisters = [false, true]; + let optimizedBEControlRegisters = [true, false, true]; + let pauliBases = [1, 2]; + let indexRegisters = [qubitPidx, qubitQidx, qubitRidx]; + return new OptimizedBETermIndex { + Coefficient = 2.0 * coeff[0], + UseSignQubit = sign, + ZControlRegisterMask = selectZControlRegisters, + OptimizedControlRegisterMask = optimizedBEControlRegisters, + PauliBases = pauliBases, + RegisterIndices = indexRegisters + }; + } +} + + +/// # Summary +/// Converts a `GeneratorIndex` describing a PQRS term to +/// an expression `GeneratorIndex[]` in terms of Paulis +/// +/// # Input +/// ## term +/// `GeneratorIndex` representing a PQRS term. +/// +/// # Output +/// `GeneratorIndex[]` expressing PQRS term as Pauli terms. +function V0123TermToPauliMajIdx(term : GeneratorIndex) : OptimizedBETermIndex[] { + let ((idxTermType, v0123), idxFermions) = term!; + let qubitsPQ = idxFermions[0..1]; + let qubitsRS = idxFermions[2..3]; + let qubitsPQJW = RangeAsIntArray(qubitsPQ[0] + 1..qubitsPQ[1] - 1); + let qubitsRSJW = RangeAsIntArray(qubitsRS[0] + 1..qubitsRS[1] - 1); + let ops = [[1, 1, 1, 1], [1, 1, 2, 2], [1, 2, 1, 2], [1, 2, 2, 1], [2, 2, 2, 2], [2, 2, 1, 1], [2, 1, 2, 1], [2, 1, 1, 2]]; + mutable majIdxes = Repeated( + new OptimizedBETermIndex { + Coefficient = 0.0, + UseSignQubit = false, + ZControlRegisterMask = [], + OptimizedControlRegisterMask = [], + PauliBases = [], + RegisterIndices = [] + }, + 4 + ); + mutable nonZero = 0; + let selectZControlRegisters = []; + let optimizedBEControlRegisters = [true, true, true, true]; + let indexRegisters = idxFermions; + + for idxOp in 0..3 { + if IsNotZero(v0123[idxOp]) { + let newCoeff = (2.0 * 0.25) * v0123[idxOp]; + majIdxes w/= nonZero <- new OptimizedBETermIndex { + Coefficient = newCoeff, + UseSignQubit = v0123[idxOp] < 0.0, + ZControlRegisterMask = selectZControlRegisters, + OptimizedControlRegisterMask = optimizedBEControlRegisters, + PauliBases = ops[idxOp], + RegisterIndices = indexRegisters + }; + nonZero = nonZero + 1; + } + } + + return majIdxes[0..nonZero - 1]; +} + + +/// # Summary +/// Converts a Hamiltonian described by `JWOptimizedHTerms` +/// to a `GeneratorSystem` expressed in terms of the Pauli +/// `GeneratorIndex`. +/// +/// # Input +/// ## data +/// Description of Hamiltonian in `JWOptimizedHTerms` format. +/// +/// # Output +/// Representation of Hamiltonian as `GeneratorSystem`. +function OptimizedBlockEncodingGeneratorSystem(data : JWOptimizedHTerms) : OptimizedBEGeneratorSystem { + let (ZData, ZZData, PQandPQQRData, h0123Data) = data!; + mutable majIdxes = Repeated( + new OptimizedBETermIndex { + Coefficient = 0.0, + UseSignQubit = false, + ZControlRegisterMask = [], + OptimizedControlRegisterMask = [], + PauliBases = [], + RegisterIndices = [] + }, + ((Length(ZData) + Length(ZZData)) + Length(PQandPQQRData)) + 4 * Length(h0123Data) + ); + mutable startIdx = 0; + + for idx in IndexRange(ZData) { + // Array of Arrays of Length 1 + majIdxes w/= idx <- ZTermToPauliMajIdx(HTermToGenIdx(ZData[idx], [0])); + } + + startIdx = Length(ZData); + + for idx in IndexRange(ZZData) { + // Array of Arrays of Length 1 + majIdxes w/= startIdx + idx <- ZZTermToPauliMajIdx(HTermToGenIdx(ZZData[idx], [1])); + } + + startIdx = startIdx + Length(ZZData); + + for idx in IndexRange(PQandPQQRData) { + + // Array of Arrays of Length 1 + majIdxes w/= startIdx + idx <- PQandPQQRTermToPauliMajIdx(HTermToGenIdx(PQandPQQRData[idx], [2])); + } + + startIdx = startIdx + Length(PQandPQQRData); + mutable finalIdx = startIdx; + + for idx in 0..Length(h0123Data) - 1 { + + // Array of Arrays of Length up to 4 + let genArr = V0123TermToPauliMajIdx(HTermToGenIdx(h0123Data[idx], [3])); + + for idx0123 in IndexRange(genArr) { + majIdxes w/= finalIdx <- genArr[idx0123]; + finalIdx = finalIdx + 1; + } + } + + mutable oneNorm = 0.0; + + for idx in 0..finalIdx - 1 { + oneNorm = oneNorm + AbsD(GetOptimizedBETermIndexCoeff(majIdxes[idx])); + } + + let majIdxes = majIdxes[0..finalIdx - 1]; + return new OptimizedBEGeneratorSystem { NumTerms = finalIdx, Norm = oneNorm, SelectTerm = idx -> majIdxes[idx] }; +} + + +operation ToJordanWignerSelectInput( + idx : Int, + optimizedBEGeneratorSystem : OptimizedBEGeneratorSystem, + signQubit : Qubit, + selectZControlRegisters : Qubit[], + optimizedBEControlRegisters : Qubit[], + pauliBasesIdx : Qubit[], + indexRegisters : Qubit[][] +) : Unit is Adj + Ctl { + let (nTerms, oneNorm, intToGenIdx) = optimizedBEGeneratorSystem!; + let (coeff, signQubitSet, selectZControlRegistersSet, OptimizedBEControlRegistersSet, pauliBasesSet, indexRegistersSet) = (intToGenIdx(idx))!; + + // Write bit to apply - signQubit + if (signQubitSet == true) { + X(signQubit); + } + + // Write bit to activate selectZ operator + for i in IndexRange(selectZControlRegistersSet) { + if selectZControlRegistersSet[i] { + X(selectZControlRegisters[i]); + } + } + + // Write bit to activate OptimizedBEXY operator + for i in IndexRange(OptimizedBEControlRegistersSet) { + if OptimizedBEControlRegistersSet[i] { + X(optimizedBEControlRegisters[i]); + } + } + + // Write bitstring to apply desired XZ... or YZ... Pauli string + for i in IndexRange(indexRegistersSet) { + ApplyXorInPlace(indexRegistersSet[i], indexRegisters[i]); + } + + // Crete state to select uniform superposition of X and Y operators. + if Length(pauliBasesSet) == 2 { + // for PQ or PQQR terms, create |00> + |11> + ApplyXorInPlace(0, pauliBasesIdx); + } elif Length(pauliBasesSet) == 4 { + // for PQRS terms, create |abcd> + |a^ b^ c^ d^> + if pauliBasesSet[2] == 1 and pauliBasesSet[3] == 1 { + ApplyXorInPlace(1, pauliBasesIdx); + } elif pauliBasesSet[2] == 2 and pauliBasesSet[3] == 2 { + ApplyXorInPlace(2, pauliBasesIdx); + } elif pauliBasesSet[2] == 1 and pauliBasesSet[3] == 2 { + ApplyXorInPlace(3, pauliBasesIdx); + } elif pauliBasesSet[2] == 2 and pauliBasesSet[3] == 1 { + ApplyXorInPlace(4, pauliBasesIdx); + } + } +} + +operation ToPauliBases(idx : Int, pauliBases : Qubit[]) : Unit is Adj + Ctl { + let pauliBasesSet = [[1, 1, 1, 1], [1, 1, 2, 2], [1, 2, 1, 2], [1, 2, 2, 1]]; + H(pauliBases[0]); + + if idx > 0 { + for idxSet in 1..Length(pauliBasesSet[0]) - 1 { + if (pauliBasesSet[idx - 1])[idxSet] == 2 { + X(pauliBases[idxSet]); + } + + CNOT(pauliBases[0], pauliBases[idxSet]); + } + } +} + +// This prepares the state that selects _JordanWignerSelect_; +operation JordanWignerOptimizedBlockEncodingStatePrep(targetError : Double, optimizedBEGeneratorSystem : OptimizedBEGeneratorSystem, qROMIdxRegister : Qubit[], qROMGarbage : Qubit[], signQubit : Qubit, selectZControlRegisters : Qubit[], optimizedBEControlRegisters : Qubit[], pauliBases : Qubit[], pauliBasesIdx : Qubit[], indexRegisters : Qubit[][]) : Unit is Adj + Ctl { + let (nTerms, _, _) = optimizedBEGeneratorSystem!; + let coefficients = OptimizedBEGeneratorSystemCoeff(optimizedBEGeneratorSystem); + let purifiedState = PurifiedMixedState(targetError, coefficients); + let unitaryGenerator = (nTerms, idx -> ToJordanWignerSelectInput(idx, optimizedBEGeneratorSystem, _, _, _, _, _)); + let pauliBasesUnitaryGenerator = (5, idx -> (qs => ToPauliBases(idx, qs))); + + purifiedState.Prepare(qROMIdxRegister, [], qROMGarbage); + MultiplexOperationsFromGenerator(unitaryGenerator, qROMIdxRegister, (signQubit, selectZControlRegisters, optimizedBEControlRegisters, pauliBasesIdx, indexRegisters)); + MultiplexOperationsFromGenerator(pauliBasesUnitaryGenerator, pauliBasesIdx, pauliBases); +} + +function JordanWignerOptimizedBlockEncodingQubitManager(targetError : Double, nCoeffs : Int, nZ : Int, nMaj : Int, nIdxRegQubits : Int, ctrlRegister : Qubit[]) : ((Qubit[], Qubit[], Qubit, Qubit[], Qubit[], Qubit[], Qubit[], Qubit[][]), (Qubit, Qubit[], Qubit[], Qubit[], Qubit[][]), Qubit[]) { + let requirements = PurifiedMixedStateRequirements(targetError, nCoeffs); + let parts = Partitioned([requirements.NumIndexQubits, requirements.NumGarbageQubits], ctrlRegister); + let ((qROMIdx, qROMGarbage), rest0) = ((parts[0], parts[1]), parts[2]); + let ((signQubit, selectZControlRegisters, optimizedBEControlRegisters, pauliBases, indexRegisters, tmp), rest1) = JordanWignerSelectQubitManager(nZ, nMaj, nIdxRegQubits, rest0, []); + let registers = Partitioned([3], rest1); + let pauliBasesIdx = registers[0]; + return ((qROMIdx, qROMGarbage, signQubit, selectZControlRegisters, optimizedBEControlRegisters, pauliBases, pauliBasesIdx, indexRegisters), (signQubit, selectZControlRegisters, optimizedBEControlRegisters, pauliBases, indexRegisters), registers[1]); +} + +function JordanWignerOptimizedBlockEncodingQubitCount(targetError : Double, nCoeffs : Int, nZ : Int, nMaj : Int, nIdxRegQubits : Int, nTarget : Int) : ((Int, Int), (Int, Int, Int, Int, Int, Int, Int, Int[], Int)) { + let (nSelectTotal, (a0, a1, a2, a3, a4)) = JordanWignerSelectQubitCount(nZ, nMaj, nIdxRegQubits); + let (nQROMTotal, b0, b1) = PurifiedMixedStateRequirements(targetError, nCoeffs)!; + let pauliBasesIdx = 3; + return (((nSelectTotal + nQROMTotal) + pauliBasesIdx, nTarget), (b0, b1, a0, a1, a2, a3, pauliBasesIdx, a4, nTarget)); +} + + +operation JordanWignerOptimizedBlockEncodingStatePrepWrapper(targetError : Double, nCoeffs : Int, optimizedBEGeneratorSystem : OptimizedBEGeneratorSystem, nZ : Int, nMaj : Int, nIdxRegQubits : Int, ctrlRegister : Qubit[]) : Unit is Adj + Ctl { + let (statePrepRegister, selectRegister, rest) = JordanWignerOptimizedBlockEncodingQubitManager(targetError, nCoeffs, nZ, nMaj, nIdxRegQubits, ctrlRegister); + let statePrepOp = JordanWignerOptimizedBlockEncodingStatePrep(targetError, optimizedBEGeneratorSystem, _, _, _, _, _, _, _, _); + statePrepOp(statePrepRegister); +} + + +operation JordanWignerOptimizedBlockEncodingSelect(targetError : Double, nCoeffs : Int, optimizedBEGeneratorSystem : OptimizedBEGeneratorSystem, nZ : Int, nMaj : Int, nIdxRegQubits : Int, ctrlRegister : Qubit[], targetRegister : Qubit[]) : Unit is Adj + Ctl { + let (statePrepRegister, selectRegister, rest) = JordanWignerOptimizedBlockEncodingQubitManager(targetError, nCoeffs, nZ, nMaj, nIdxRegQubits, ctrlRegister); + let selectOp = JordanWignerSelect(_, _, _, _, _, targetRegister); + selectOp(selectRegister); +} + + +function JordanWignerOptimizedBlockEncoding(targetError : Double, data : JWOptimizedHTerms, nSpinOrbitals : Int) : ((Int, Int), (Double, (Qubit[], Qubit[]) => Unit is Adj + Ctl)) { + let nZ = 2; + let nMaj = 4; + let optimizedBEGeneratorSystem = OptimizedBlockEncodingGeneratorSystem(data); + let (nCoeffs, oneNorm, tmp) = optimizedBEGeneratorSystem!; + let nIdxRegQubits = Ceiling(Lg(IntAsDouble(nSpinOrbitals))); + let ((nCtrlRegisterQubits, nTargetRegisterQubits), rest) = JordanWignerOptimizedBlockEncodingQubitCount(targetError, nCoeffs, nZ, nMaj, nIdxRegQubits, nSpinOrbitals); + let statePrepOp = JordanWignerOptimizedBlockEncodingStatePrepWrapper(targetError, nCoeffs, optimizedBEGeneratorSystem, nZ, nMaj, nIdxRegQubits, _); + let selectOp = JordanWignerOptimizedBlockEncodingSelect(targetError, nCoeffs, optimizedBEGeneratorSystem, nZ, nMaj, nIdxRegQubits, _, _); + let blockEncodingReflection = BlockEncodingByLCU(statePrepOp, selectOp); + return ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, blockEncodingReflection)); +} + +function JordanWignerOptimizedQuantumWalkByQubitization(targetError : Double, data : JWOptimizedHTerms, nSpinOrbitals : Int) : ((Int, Int), (Double, ((Qubit[], Qubit[]) => Unit is Adj + Ctl))) { + let ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, blockEncodingReflection)) = JordanWignerOptimizedBlockEncoding(targetError, data, nSpinOrbitals); + return ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, QuantumWalkByQubitization(blockEncodingReflection))); +} + +/// # Summary +/// Represents a particular mixed state that can be prepared on an index +/// and a garbage register. +/// +/// # Input +/// ## Requirements +/// Specifies the size of the qubit registers required to prepare the +/// mixed state represented by this UDT value. +/// ## Norm +/// Specifies the 1-norm of the coefficients used to define this mixed +/// state. +/// ## Prepare +/// An operation that, given an index register, a data register, and a +/// garbage register initially in the $\ket{0}$, $\let{00\dots 0}$, and +/// $\ket{00\dots 0}$ states (respectively), +/// prepares the state represented by this UDT value on those registers. +struct MixedStatePreparation { + Requirements : MixedStatePreparationRequirements, + Norm : Double, + Prepare : ((Qubit[], Qubit[], Qubit[]) => Unit is Adj + Ctl), +} + +/// # Summary +/// Represents the number of qubits required in order to prepare a given +/// mixed state. +/// +/// # Input +/// ## NTotalQubits +/// The total number of qubits required by the represented state preparation +/// operation. +/// ## NIndexQubits +/// The number of qubits required for the index register used by the +/// represented state preparation operation. +/// ## NGarbageQubits +/// The number of qubits required for the garbage register used by the +/// represented state preparation operation. +struct MixedStatePreparationRequirements { + NumTotalQubits : Int, + NumIndexQubits : Int, + NumGarbageQubits : Int, +} + +/// # Summary +/// Returns an operation that prepares a a purification of a given mixed state. +/// A "purified mixed state" refers to states of the form |ψ⟩ = Σᵢ √𝑝ᵢ |𝑖⟩ |garbageᵢ⟩ specified by a vector of +/// coefficients {𝑝ᵢ}. States of this form can be reduced to mixed states ρ ≔ 𝑝ᵢ |𝑖⟩⟨𝑖| by tracing over the "garbage" +/// register (that is, a mixed state that is diagonal in the computational basis). +/// +/// See https://arxiv.org/pdf/1805.03662.pdf?page=15 for further discussion. +/// +/// # Description +/// Uses the Quantum ROM technique to represent a given density matrix, +/// returning that representation as a state preparation operation. +/// +/// In particular, given a list of $N$ coefficients $\alpha_j$, this +/// function returns an operation that uses the Quantum ROM technique to +/// prepare an approximation +/// $$ +/// \begin{align} +/// \tilde\rho = \sum_{j = 0}^{N - 1} p_j \ket{j}\bra{j} +/// \end{align} +/// $$ +/// of the mixed state +/// $$ +/// \begin{align} +/// \rho = \sum_{j = 0}^{N-1} \frac{|\alpha_j|}{\sum_k |\alpha_k|} \ket{j}\bra{j}, +/// \end{align} +/// $$ +/// where each $p_j$ is an approximation to the given coefficient $\alpha_j$ +/// such that +/// $$ +/// \begin{align} +/// \left| p_j - \frac{ |\alpha_j| }{ \sum_k |\alpha_k| } \right| \le \frac{\epsilon}{N} +/// \end{align} +/// $$ +/// for each $j$. +/// +/// When passed an index register and a register of garbage qubits, +/// initially in the state $\ket{0} \ket{00\cdots 0}$, the returned operation +/// prepares both registers into the purification of $\tilde \rho$, +/// $$ +/// \begin{align} +/// \sum_{j=0}^{N-1} \sqrt{p_j} \ket{j}\ket{\text{garbage}_j}, +/// \end{align} +/// $$ +/// such that resetting and deallocating the garbage register enacts the +/// desired preparation to within the target error $\epsilon$. +/// +/// # Input +/// ## targetError +/// The target error $\epsilon$. +/// ## coefficients +/// Array of $N$ coefficients specifying the probability of basis states. +/// Negative numbers $-\alpha_j$ will be treated as positive $|\alpha_j|$. +/// +/// # Output +/// An operation that prepares $\tilde \rho$ as a purification onto a joint +/// index and garbage register. +/// +/// # Remarks +/// The coefficients provided to this operation are normalized following the +/// 1-norm, such that the coefficients are always considered to describe a +/// valid categorical probability distribution. +/// +/// # Example +/// The following code snippet prepares an purification of the $3$-qubit state +/// $\rho=\sum_{j=0}^{4}\frac{|\alpha_j|}{\sum_k |\alpha_k|}\ket{j}\bra{j}$, where +/// $\vec\alpha=(1.0, 2.0, 3.0, 4.0, 5.0)$, and the target error is +/// $10^{-3}$: +/// ```qsharp +/// let coefficients = [1.0, 2.0, 3.0, 4.0, 5.0]; +/// let targetError = 1e-3; +/// let purifiedState = PurifiedMixedState(targetError, coefficients); +/// using (indexRegister = Qubit[purifiedState.Requirements.NIndexQubits]) { +/// using (garbageRegister = Qubit[purifiedState.Requirements.NGarbageQubits]) { +/// purifiedState.Prepare(LittleEndian(indexRegister), [], garbageRegister); +/// } +/// } +/// ``` +/// +/// # References +/// - [Encoding Electronic Spectra in Quantum Circuits with Linear T Complexity](https://arxiv.org/abs/1805.03662) +/// Ryan Babbush, Craig Gidney, Dominic W. Berry, Nathan Wiebe, Jarrod McClean, Alexandru Paler, Austin Fowler, Hartmut Neven +function PurifiedMixedState(targetError : Double, coefficients : Double[]) : MixedStatePreparation { + let nBitsPrecision = -Ceiling(Lg(0.5 * targetError)) + 1; + let positiveCoefficients = Mapped(AbsD, coefficients); + let (oneNorm, keepCoeff, altIndex) = QuantumROMDiscretization(nBitsPrecision, positiveCoefficients); + let nCoeffs = Length(positiveCoefficients); + let nBitsIndices = Ceiling(Lg(IntAsDouble(nCoeffs))); + + let op = PrepareQuantumROMState(nBitsPrecision, nCoeffs, nBitsIndices, keepCoeff, altIndex, [], _, _, _); + let qubitCounts = PurifiedMixedStateRequirements(targetError, nCoeffs); + return new MixedStatePreparation { Requirements = qubitCounts, Norm = oneNorm, Prepare = op }; +} + +operation PrepareQuantumROMState( + nBitsPrecision : Int, + nCoeffs : Int, + nBitsIndices : Int, + keepCoeff : Int[], + altIndex : Int[], + data : Bool[][], + indexRegister : Qubit[], + dataQubits : Qubit[], + garbageRegister : Qubit[] +) : Unit is Adj + Ctl { + let garbageIdx0 = nBitsIndices; + let garbageIdx1 = garbageIdx0 + nBitsPrecision; + let garbageIdx2 = garbageIdx1 + nBitsPrecision; + let garbageIdx3 = garbageIdx2 + 1; + + let altIndexRegister = garbageRegister[0..garbageIdx0 - 1]; + let keepCoeffRegister = garbageRegister[garbageIdx0..garbageIdx1 - 1]; + let uniformKeepCoeffRegister = garbageRegister[garbageIdx1..garbageIdx2 - 1]; + let flagQubit = garbageRegister[garbageIdx3 - 1]; + let dataRegister = dataQubits; + let altDataRegister = garbageRegister[garbageIdx3...]; + + // Create uniform superposition over index and alt coeff register. + PrepareUniformSuperposition(nCoeffs, indexRegister); + ApplyToEachCA(H, uniformKeepCoeffRegister); + + // Write bitstrings to altIndex and keepCoeff register. + let unitaryGenerator = (nCoeffs, idx -> WriteQuantumROMBitString(idx, keepCoeff, altIndex, data, _, _, _, _)); + MultiplexOperationsFromGenerator(unitaryGenerator, indexRegister, (keepCoeffRegister, altIndexRegister, dataRegister, altDataRegister)); + + // Perform comparison + ApplyIfGreaterLE(X, uniformKeepCoeffRegister, keepCoeffRegister, flagQubit); + + let indexRegisterSize = Length(indexRegister); + + // Swap in register based on comparison + let register = indexRegister + dataRegister; + let altRegister = altIndexRegister + altDataRegister; + for i in IndexRange(indexRegister) { + Controlled SWAP([flagQubit], (register[i], altRegister[i])); + } +} + +/// # Summary +/// Creates a uniform superposition over states that encode 0 through `nIndices - 1`. +/// +/// # Description +/// +/// This operation can be described by a unitary matrix $U$ that creates +/// a uniform superposition over all number states +/// $0$ to $M-1$, given an input state $\ket{0\cdots 0}$. In other words, +/// $$ +/// \begin{align} +/// U \ket{0} = \frac{1}{\sqrt{M}} \sum_{j=0}^{M - 1} \ket{j}. +/// \end{align} +/// $$. +/// +/// # Input +/// ## nIndices +/// The desired number of states $M$ in the uniform superposition. +/// ## indexRegister +/// The qubit register that stores the number states in `LittleEndian` format. +/// This register must be able to store the number $M-1$, and is assumed to be +/// initialized in the state $\ket{0\cdots 0}$. +/// +/// # Remarks +/// The operation is adjointable, but requires that `indexRegister` is in a uniform +/// superposition over the first `nIndices` basis states in that case. +/// +/// # Example +/// The following example prepares the state $\frac{1}{\sqrt{6}}\sum_{j=0}^{5}\ket{j}$ +/// on $3$ qubits. +/// ``` Q# +/// let nIndices = 6; +/// using(indexRegister = Qubit[3]) { +/// PrepareUniformSuperposition(nIndices, LittleEndian(indexRegister)); +/// // ... +/// } +/// ``` +operation PrepareUniformSuperposition(nIndices : Int, indexRegister : Qubit[]) : Unit is Adj + Ctl { + if nIndices == 0 { + fail "Cannot prepare uniform superposition over 0 basis states."; + } elif nIndices == 1 { + // Superposition over one state, so do nothing. + } elif nIndices == 2 { + H(indexRegister[0]); + } else { + let nQubits = BitSizeI(nIndices - 1); + if nQubits > Length(indexRegister) { + fail $"Cannot prepare uniform superposition over {nIndices} states as it is larger than the qubit register."; + } + + use flagQubit = Qubit[3]; + let targetQubits = indexRegister[0..nQubits - 1]; + let qubits = flagQubit + targetQubits; + let stateOracle = PrepareUniformSuperpositionOracle(nIndices, nQubits, 0, _); + + let phases = ([0.0, PI()], [PI(), 0.0]); + + ObliviousAmplitudeAmplificationFromStatePreparation( + phases, + stateOracle, + (qs0, qs1) => I(qs0[0]), + 0 + )(qubits, []); + + + ApplyToEachCA(X, flagQubit); + } +} + +function ObliviousAmplitudeAmplificationFromStatePreparation( + phases : (Double[], Double[]), + startStateOracle : Qubit[] => Unit is Adj + Ctl, + signalOracle : (Qubit[], Qubit[]) => Unit is Adj + Ctl, + idxFlagQubit : Int +) : (Qubit[], Qubit[]) => Unit is Adj + Ctl { + let startStateReflection = (phase, qs) => { + within { + ApplyToEachCA(X, qs); + } apply { + Controlled R1(Rest(qs), (phase, qs[0])); + } + }; + let targetStateReflection = (phase, qs) => R1(phase, qs[idxFlagQubit]); + let obliviousSignalOracle = (qs0, qs1) => { startStateOracle(qs0); signalOracle(qs0, qs1); }; + return ObliviousAmplitudeAmplificationFromPartialReflections( + phases, + startStateReflection, + targetStateReflection, + obliviousSignalOracle + ); +} + +function ObliviousAmplitudeAmplificationFromPartialReflections( + phases : (Double[], Double[]), + startStateReflection : (Double, Qubit[]) => Unit is Adj + Ctl, + targetStateReflection : (Double, Qubit[]) => Unit is Adj + Ctl, + signalOracle : (Qubit[], Qubit[]) => Unit is Adj + Ctl +) : (Qubit[], Qubit[]) => Unit is Adj + Ctl { + return ApplyObliviousAmplitudeAmplification( + phases, + startStateReflection, + targetStateReflection, + signalOracle, + _, + _ + ); +} + +/// # Summary +/// Oblivious amplitude amplification by specifying partial reflections. +/// +/// # Input +/// ## phases +/// Phases of partial reflections +/// ## startStateReflection +/// Reflection operator about start state of auxiliary register +/// ## targetStateReflection +/// Reflection operator about target state of auxiliary register +/// ## signalOracle +/// Unitary oracle $O$ of type `ObliviousOracle` that acts jointly on the +/// auxiliary and system registers. +/// ## auxiliaryRegister +/// Auxiliary register +/// ## systemRegister +/// System register +/// +/// # Remarks +/// Given a particular auxiliary start state $\ket{\text{start}}\_a$, a +/// particular auxiliary target state $\ket{\text{target}}\_a$, and any +/// system state $\ket{\psi}\_s$, suppose that +/// \begin{align} +/// O\ket{\text{start}}\_a\ket{\psi}\_s= \lambda\ket{\text{target}}\_a U \ket{\psi}\_s + \sqrt{1-|\lambda|^2}\ket{\text{target}^\perp}\_a\cdots +/// \end{align} +/// for some unitary $U$. +/// By a sequence of reflections about the start and target states on the +/// auxiliary register interleaved by applications of `signalOracle` and its +/// adjoint, the success probability of applying $U$ may be altered. +/// +/// In most cases, `auxiliaryRegister` is initialized in the state $\ket{\text{start}}\_a$. +/// +/// # References +/// - See [*D.W. Berry, A.M. Childs, R. Cleve, R. Kothari, R.D. Somma*](https://arxiv.org/abs/1312.1414) +/// for the standard version. +/// - See [*G.H. Low, I.L. Chuang*](https://arxiv.org/abs/1610.06546) +/// for a generalization to partial reflections. +operation ApplyObliviousAmplitudeAmplification( + phases : (Double[], Double[]), + startStateReflection : (Double, Qubit[]) => Unit is Adj + Ctl, + targetStateReflection : (Double, Qubit[]) => Unit is Adj + Ctl, + signalOracle : (Qubit[], Qubit[]) => Unit is Adj + Ctl, + auxiliaryRegister : Qubit[], + systemRegister : Qubit[] +) : Unit is Adj + Ctl { + let (aboutStart, aboutTarget) = phases; + Fact(Length(aboutStart) == Length(aboutTarget), "number of phases about start and target state must be equal"); + let numPhases = Length(aboutStart); + + for idx in 0..numPhases-1 { + if aboutStart[idx] != 0.0 { + startStateReflection(aboutStart[idx], auxiliaryRegister); + } + + if aboutTarget[idx] != 0.0 { + // In the last iteration we do not need to apply `Adjoint signalOracle` + if idx == numPhases - 1 { + signalOracle(auxiliaryRegister, systemRegister); + targetStateReflection(aboutTarget[idx], auxiliaryRegister); + } else { + within { + signalOracle(auxiliaryRegister, systemRegister); + } apply { + targetStateReflection(aboutTarget[idx], auxiliaryRegister); + } + } + } + } + + // We do need one more `signalOracle` call, if the last phase about the target state was 0.0 + if numPhases == 0 or aboutTarget[numPhases - 1] == 0.0 { + signalOracle(auxiliaryRegister, systemRegister); + } +} + +operation PrepareUniformSuperpositionOracle(nIndices : Int, nQubits : Int, idxFlag : Int, qubits : Qubit[]) : Unit is Adj + Ctl { + let targetQubits = qubits[3...]; + let flagQubit = qubits[0]; + let auxillaryQubits = qubits[1..2]; + let theta = ArcSin(Sqrt(IntAsDouble(2^nQubits) / IntAsDouble(nIndices)) * Sin(PI() / 6.0)); + + ApplyToEachCA(H, targetQubits); + use compareQubits = Qubit[nQubits] { + within { + ApplyXorInPlace(nIndices - 1, compareQubits); + } apply { + ApplyIfGreaterLE(X, targetQubits, compareQubits, auxillaryQubits[0]); + X(auxillaryQubits[0]); + } + } + Ry(2.0 * theta, auxillaryQubits[1]); + (Controlled X)(auxillaryQubits, flagQubit); +} + +// Classical processing +// This discretizes the coefficients such that +// |coefficient[i] * oneNorm - discretizedCoefficient[i] * discretizedOneNorm| * nCoeffs <= 2^{1-bitsPrecision}. +function QuantumROMDiscretization(bitsPrecision : Int, coefficients : Double[]) : (Double, Int[], Int[]) { + let oneNorm = PNorm(1.0, coefficients); + let nCoefficients = Length(coefficients); + Fact(bitsPrecision <= 31, $"Bits of precision {bitsPrecision} unsupported. Max is 31."); + Fact(nCoefficients > 1, "Cannot prepare state with less than 2 coefficients."); + Fact(oneNorm >= 0.0, "State must have at least one coefficient > 0"); + + let barHeight = 2^bitsPrecision - 1; + + mutable altIndex = RangeAsIntArray(0..nCoefficients - 1); + mutable keepCoeff = Mapped( + coefficient -> Round((AbsD(coefficient) / oneNorm) * IntAsDouble(nCoefficients) * IntAsDouble(barHeight)), + coefficients + ); + + // Calculate difference between number of discretized bars vs. maximum + mutable bars = 0; + for idxCoeff in IndexRange(keepCoeff) { + bars += keepCoeff[idxCoeff] - barHeight; + } + + // Uniformly distribute excess bars across coefficients. + for idx in 0..AbsI(bars) - 1 { + keepCoeff w/= idx <- keepCoeff[idx] + (bars > 0 ? -1 | + 1); + } + + mutable barSink = []; + mutable barSource = []; + + for idxCoeff in IndexRange(keepCoeff) { + if keepCoeff[idxCoeff] > barHeight { + barSource += [idxCoeff]; + } elif keepCoeff[idxCoeff] < barHeight { + barSink += [idxCoeff]; + } + } + + for rep in 0..nCoefficients * 10 { + if Length(barSink) > 0 and Length(barSource) > 0 { + let idxSink = Tail(barSink); + let idxSource = Tail(barSource); + barSink = Most(barSink); + barSource = Most(barSource); + + keepCoeff w/= idxSource <- keepCoeff[idxSource] - barHeight + keepCoeff[idxSink]; + altIndex w/= idxSink <- idxSource; + + if keepCoeff[idxSource] < barHeight { + barSink += [idxSource]; + } elif keepCoeff[idxSource] > barHeight { + barSource += [idxSource]; + } + } elif Length(barSource) > 0 { + let idxSource = Tail(barSource); + barSource = Most(barSource); + keepCoeff w/= idxSource <- barHeight; + } else { + return (oneNorm, keepCoeff, altIndex); + } + } + + return (oneNorm, keepCoeff, altIndex); +} + +/// # Summary +/// Returns the total number of qubits that must be allocated +/// in order to apply the operation returned by +/// `PurifiedMixedState`. +/// +/// # Input +/// ## targetError +/// The target error $\epsilon$. +/// ## nCoefficients +/// The number of coefficients to be specified in preparing a mixed state. +/// +/// # Output +/// A description of how many qubits are required in total, and for each of +/// the index and garbage registers used by the +/// `PurifiedMixedState` function. +function PurifiedMixedStateRequirements(targetError : Double, nCoefficients : Int) : MixedStatePreparationRequirements { + Fact(targetError > 0.0, "targetError must be positive"); + Fact(nCoefficients > 0, "nCoefficients must be positive"); + + let nBitsPrecision = -Ceiling(Lg(0.5 * targetError)) + 1; + let nIndexQubits = Ceiling(Lg(IntAsDouble(nCoefficients))); + let nGarbageQubits = nIndexQubits + 2 * nBitsPrecision + 1; + let nTotal = nGarbageQubits + nIndexQubits; + return new MixedStatePreparationRequirements { NumTotalQubits = nTotal, NumIndexQubits = nIndexQubits, NumGarbageQubits = nGarbageQubits }; +} + +/// # Summary +/// Encodes an operator of interest into a `BlockEncoding`. +/// +/// This constructs a `BlockEncoding` unitary $U=P\cdot V\cdot P^\dagger$ that encodes some +/// operator $H = \sum_{j}|\alpha_j|U_j$ of interest that is a linear combination of +/// unitaries. Typically, $P$ is a state preparation unitary such that +/// $P\ket{0}\_a=\sum_j\sqrt{\alpha_j/\|\vec\alpha\|\_2}\ket{j}\_a$, +/// and $V=\sum_{j}\ket{j}\bra{j}\_a\otimes U_j$. +/// +/// # Input +/// ## statePreparation +/// A unitary $P$ that prepares some target state. +/// ## selector +/// A unitary $V$ that encodes the component unitaries of $H$. +/// +/// # Output +/// A unitary $U$ acting jointly on registers `a` and `s` that block- +/// encodes $H$, and satisfies $U^\dagger = U$. +/// +/// # Remarks +/// This `BlockEncoding` implementation gives it the properties of a +/// `BlockEncodingReflection`. +function BlockEncodingByLCU<'T, 'S>( + statePreparation : ('T => Unit is Adj + Ctl), + selector : (('T, 'S) => Unit is Adj + Ctl) +) : (('T, 'S) => Unit is Adj + Ctl) { + return ApplyBlockEncodingByLCU(statePreparation, selector, _, _); +} + +/// # Summary +/// Implementation of `BlockEncodingByLCU`. +operation ApplyBlockEncodingByLCU<'T, 'S>( + statePreparation : ('T => Unit is Adj + Ctl), + selector : (('T, 'S) => Unit is Adj + Ctl), + auxiliary : 'T, + system : 'S +) : Unit is Adj + Ctl { + within { + statePreparation(auxiliary); + } apply { + selector(auxiliary, system); + } +} + +/// # Summary +/// Converts a block-encoding reflection into a quantum walk. +/// +/// # Description +/// Given a block encoding represented by a unitary $U$ +/// that encodes some operator $H$ of interest, converts it into a quantum +/// walk $W$ containing the spectrum of $\pm e^{\pm i\sin^{-1}(H)}$. +/// +/// # Input +/// ## blockEncoding +/// A unitary $U$ to be converted into a Quantum +/// walk. +/// +/// # Output +/// A quantum walk $W$ acting jointly on registers `a` and `s` that block- +/// encodes $H$, and contains the spectrum of $\pm e^{\pm i\sin^{-1}(H)}$. +/// +/// # References +/// - [Hamiltonian Simulation by Qubitization](https://arxiv.org/abs/1610.06546) +/// Guang Hao Low, Isaac L. Chuang +function QuantumWalkByQubitization(blockEncoding : (Qubit[], Qubit[]) => Unit is Adj + Ctl) : ((Qubit[], Qubit[]) => Unit is Adj + Ctl) { + return ApplyQuantumWalkByQubitization(blockEncoding, _, _); +} + +/// # Summary +/// Implementation of `Qubitization`. +operation ApplyQuantumWalkByQubitization( + blockEncoding : (Qubit[], Qubit[]) => Unit is Adj + Ctl, + auxiliary : Qubit[], + system : Qubit[] +) : Unit is Adj + Ctl { + Exp([PauliI], -0.5 * PI(), [Head(system)]); + within { + ApplyToEachCA(X, auxiliary); + } apply { + Controlled R1(Rest(auxiliary), (PI(), Head(system))); + } + blockEncoding(auxiliary, system); +} + +operation WriteQuantumROMBitString(idx : Int, keepCoeff : Int[], altIndex : Int[], data : Bool[][], keepCoeffRegister : Qubit[], altIndexRegister : Qubit[], dataRegister : Qubit[], altDataRegister : Qubit[]) : Unit is Adj + Ctl { + if keepCoeff[idx] >= 0 { + ApplyXorInPlace(keepCoeff[idx], keepCoeffRegister); + } + ApplyXorInPlace(altIndex[idx], altIndexRegister); + if Length(dataRegister) > 0 { + for i in IndexRange(data[idx]) { + if data[idx][i] { + X(dataRegister[i]); + } + } + for i in IndexRange(data[altIndex[idx]]) { + if data[altIndex[idx]][i] { + X(altDataRegister[i]); + } + } + } +} + +/// # Summary +/// Creates a block-encoding unitary for a Hamiltonian. +/// +/// The Hamiltonian $H=\sum_{j}\alpha_j P_j$ is described by a +/// sum of Pauli terms $P_j$, each with real coefficient $\alpha_j$. +/// +/// # Input +/// ## generatorSystem +/// A `GeneratorSystem` that describes $H$ as a sum of Pauli terms +/// +/// # Output +/// ## First parameter +/// The one-norm of coefficients $\alpha=\sum_{j}|\alpha_j|$. +/// ## Second parameter +/// A block encoding unitary $U$ of the Hamiltonian $H$. As this unitary +/// satisfies $U^2 = I$, it is also a reflection. +/// +/// # Remarks +/// This is obtained by preparing and unpreparing the state $\sum_{j}\sqrt{\alpha_j/\alpha}\ket{j}$, +/// and constructing a multiply-controlled unitary `PrepareArbitraryStateD` and `MultiplexOperationsFromGenerator`. +function PauliBlockEncoding(generatorSystem : GeneratorSystem) : (Double, (Qubit[], Qubit[]) => Unit is Adj + Ctl) { + let multiplexer = unitaryGenerator -> MultiplexOperationsFromGenerator(unitaryGenerator, _, _); + return PauliBlockEncodingInner(generatorSystem, coeff -> (qs => PreparePureStateD(coeff, Reversed(qs))), multiplexer); +} + +/// # Summary +/// Creates a block-encoding unitary for a Hamiltonian. +/// +/// The Hamiltonian $H=\sum_{j}\alpha_j P_j$ is described by a +/// sum of Pauli terms $P_j$, each with real coefficient $\alpha_j$. +/// +/// # Input +/// ## generatorSystem +/// A `GeneratorSystem` that describes $H$ as a sum of Pauli terms +/// ## statePrepUnitary +/// A unitary operation $P$ that prepares $P\ket{0}=\sum_{j}\sqrt{\alpha_j}\ket{j}$ given +/// an array of coefficients $\{\sqrt{\alpha}_j\}$. +/// ## statePrepUnitary +/// A unitary operation $V$ that applies the unitary $V_j$ controlled on index $\ket{j}$, +/// given a function $f: j\rightarrow V_j$. +/// +/// # Output +/// ## First parameter +/// The one-norm of coefficients $\alpha=\sum_{j}|\alpha_j|$. +/// ## Second parameter +/// A block encoding unitary $U$ of the Hamiltonian $U$. As this unitary +/// satisfies $U^2 = I$, it is also a reflection. +/// +/// # Remarks +/// Example operations the prepare and unpreparing the state $\sum_{j}\sqrt{\alpha_j/\alpha}\ket{j}$, +/// and construct a multiply-controlled unitary are +/// `PrepareArbitraryStateD` and `MultiplexOperationsFromGenerator`. +function PauliBlockEncodingInner( + generatorSystem : GeneratorSystem, + statePrepUnitary : (Double[] -> (Qubit[] => Unit is Adj + Ctl)), + multiplexer : ((Int, (Int -> (Qubit[] => Unit is Adj + Ctl))) -> ((Qubit[], Qubit[]) => Unit is Adj + Ctl)) +) : (Double, (Qubit[], Qubit[]) => Unit is Adj + Ctl) { + let (nTerms, intToGenIdx) = generatorSystem!; + let op = idx -> Sqrt(AbsD({ + let ((idxPaulis, coeff), idxQubits) = intToGenIdx(idx)!; + coeff[0] + })); + let coefficients = Mapped(op, RangeAsIntArray(0..nTerms-1)); + let oneNorm = PNorm(2.0, coefficients)^2.0; + let unitaryGenerator = (nTerms, idx -> PauliLCUUnitary(intToGenIdx(idx))); + let statePreparation = statePrepUnitary(coefficients); + let selector = multiplexer(unitaryGenerator); + let blockEncoding = (qs0, qs1) => BlockEncodingByLCU(statePreparation, selector)(qs0, qs1); + return (oneNorm, blockEncoding); +} + +/// # Summary +/// Used in implementation of `PauliBlockEncoding` +function PauliLCUUnitary(generatorIndex : GeneratorIndex) : (Qubit[] => Unit is Adj + Ctl) { + return ApplyPauliLCUUnitary(generatorIndex, _); +} + +/// # Summary +/// Used in implementation of `PauliBlockEncoding` +operation ApplyPauliLCUUnitary(generatorIndex : GeneratorIndex, qubits : Qubit[]) : Unit is Adj + Ctl { + let ((idxPaulis, coeff), idxQubits) = generatorIndex!; + let paulis = [PauliI, PauliX, PauliY, PauliZ]; + let pauliString = IntArrayAsPauliArray(idxPaulis); + let pauliQubits = Subarray(idxQubits, qubits); + + ApplyPauli(pauliString, pauliQubits); + + if (coeff[0] < 0.0) { + // -1 phase + Exp([PauliI], PI(), [Head(pauliQubits)]); + } +} + +function IntArrayAsPauliArray(arr : Int[]) : Pauli[] { + let paulis = [PauliI, PauliX, PauliY, PauliZ]; + mutable pauliString = []; + for idxP in arr { + pauliString += [paulis[idxP]]; + } + pauliString +} diff --git a/library/chemistry/src/JordanWigner/JordanWignerVQE.qs b/library/chemistry/src/JordanWigner/JordanWignerVQE.qs new file mode 100644 index 0000000000..e2c596a0a2 --- /dev/null +++ b/library/chemistry/src/JordanWigner/JordanWignerVQE.qs @@ -0,0 +1,222 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +export + EstimateTermExpectation, + EstimateFrequency, + MeasurementOperators, + ExpandedCoefficients; + +import Std.Arrays.IndexRange; +import Std.Convert.IntAsDouble; +import Std.Math.AbsD; + +/// # Summary +/// Computes the energy associated to a given Jordan-Wigner Hamiltonian term +/// +/// # Description +/// This operation estimates the expectation value associated to each measurement operator and +/// multiplies it by the corresponding coefficient, using sampling. +/// The results are aggregated into a variable containing the energy of the Jordan-Wigner term. +/// +/// # Input +/// ## inputStateUnitary +/// The unitary used for state preparation. +/// ## ops +/// The measurement operators of the Jordan-Wigner term. +/// ## coeffs +/// The coefficients of the Jordan-Wigner term. +/// ## nQubits +/// The number of qubits required to simulate the molecular system. +/// ## nSamples +/// The number of samples to use for the estimation of the term expectation. +/// +/// # Output +/// The energy associated to the Jordan-Wigner term. +operation EstimateTermExpectation(inputStateUnitary : (Qubit[] => Unit), ops : Pauli[][], coeffs : Double[], nQubits : Int, nSamples : Int) : Double { + mutable jwTermEnergy = 0.; + + for i in IndexRange(coeffs) { + // Only perform computation if the coefficient is significant enough + if (AbsD(coeffs[i]) >= 1e-10) { + // Compute expectation value using the fast frequency estimator, add contribution to Jordan-Wigner term energy + let termExpectation = EstimateFrequency(inputStateUnitary, Measure(ops[i], _), nQubits, nSamples); + jwTermEnergy += (2. * termExpectation - 1.) * coeffs[i]; + } + } + + return jwTermEnergy; +} + +/// # Summary +/// Given a preparation and measurement, estimates the frequency +/// with which that measurement succeeds (returns `Zero`) by +/// performing a given number of trials. +/// +/// # Input +/// ## preparation +/// An operation $P$ that prepares a given state $\rho$ on +/// its input register. +/// ## measurement +/// An operation $M$ representing the measurement of interest. +/// ## nQubits +/// The number of qubits on which the preparation and measurement +/// each act. +/// ## nMeasurements +/// The number of times that the measurement should be performed +/// in order to estimate the frequency of interest. +/// +/// # Output +/// An estimate $\hat{p}$ of the frequency with which +/// $M(P(\ket{00 \cdots 0}\bra{00 \cdots 0}))$ returns `Zero`, +/// obtained using the unbiased binomial estimator $\hat{p} = +/// n\_{\uparrow} / n\_{\text{measurements}}$, where $n\_{\uparrow}$ is +/// the number of `Zero` results observed. +/// +/// This is particularly important on target machines which respect +/// physical limitations, such that probabilities cannot be measured. +operation EstimateFrequency(preparation : (Qubit[] => Unit), measurement : (Qubit[] => Result), nQubits : Int, nMeasurements : Int) : Double { + mutable nUp = 0; + + for idxMeasurement in 0..nMeasurements - 1 { + use register = Qubit[nQubits]; + preparation(register); + let result = measurement(register); + + if (result == Zero) { + // NB!!!!! This reverses Zero and One to use conventions + // common in the QCVV community. That is confusing + // but is confusing with an actual purpose. + nUp = nUp + 1; + } + + // NB: We absolutely must reset here, since preparation() + // and measurement() can each use randomness internally. + ResetAll(register); + } + + return IntAsDouble(nUp) / IntAsDouble(nMeasurements); +} + + + +/// # Summary +/// Computes all the measurement operators required to compute the expectation of a Jordan-Wigner term. +/// +/// # Input +/// ## nQubits +/// The number of qubits required to simulate the molecular system. +/// ## indices +/// An array containing the indices of the qubit each Pauli operator is applied to. +/// ## termType +/// The type of the Jordan-Wigner term. +/// +/// # Output +/// An array of measurement operators (each being an array of Pauli). +function MeasurementOperators(nQubits : Int, indices : Int[], termType : Int) : Pauli[][] { + // Compute the size and initialize the array of operators to be returned + mutable nOps = 0; + if (termType == 2) { + nOps = 2; + } elif (termType == 3) { + nOps = 8; + } else { + nOps = 1; + } + + mutable ops = [[], size = nOps]; + + // Z and ZZ terms + if termType == 0 or termType == 1 { + mutable op = Repeated(PauliI, nQubits); + for idx in indices { + op w/= idx <- PauliZ; + } + ops w/= 0 <- op; + } + + // PQRS terms set operators between indices P and Q (resp R and S) to PauliZ + elif termType == 3 { + let compactOps = [[PauliX, PauliX, PauliX, PauliX], [PauliY, PauliY, PauliY, PauliY], [PauliX, PauliX, PauliY, PauliY], [PauliY, PauliY, PauliX, PauliX], [PauliX, PauliY, PauliX, PauliY], [PauliY, PauliX, PauliY, PauliX], [PauliY, PauliX, PauliX, PauliY], [PauliX, PauliY, PauliY, PauliX]]; + + for iOp in 0..7 { + mutable compactOp = compactOps[iOp]; + + mutable op = Repeated(PauliI, nQubits); + for idx in IndexRange(indices) { + op w/= indices[idx] <- compactOp[idx]; + } + for i in indices[0] + 1..indices[1] - 1 { + op w/= i <- PauliZ; + } + for i in indices[2] + 1..indices[3] - 1 { + op w/= i <- PauliZ; + } + ops w/= iOp <- op; + } + } + + // Case of PQ and PQQR terms + elif termType == 2 { + let compactOps = [[PauliX, PauliX], [PauliY, PauliY]]; + + for iOp in 0..1 { + mutable compactOp = compactOps[iOp]; + + mutable op = Repeated(PauliI, nQubits); + + let nIndices = Length(indices); + op w/= indices[0] <- compactOp[0]; + op w/= indices[nIndices-1] <- compactOp[1]; + for i in indices[0] + 1..indices[nIndices - 1] - 1 { + op w/= i <- PauliZ; + } + + // Case of PQQR term + if nIndices == 4 { + op w/= indices[1] <- (indices[0] < indices[1] and indices[1] < indices[3]) ? PauliI | PauliZ; + } + ops w/= iOp <- op; + } + } + + return ops; +} + + +/// # Summary +/// Expands the compact representation of the Jordan-Wigner coefficients in order +/// to obtain a one-to-one mapping between these and Pauli terms. +/// +/// # Input +/// ## coeff +/// An array of coefficients, as read from the Jordan-Wigner Hamiltonian data structure. +/// ## termType +/// The type of the Jordan-Wigner term. +/// +/// # Output +/// Expanded arrays of coefficients, one per Pauli term. +function ExpandedCoefficients(coeff : Double[], termType : Int) : Double[] { + // Compute the numbers of coefficients to return + mutable nCoeffs = 0; + if (termType == 2) { + nCoeffs = 2; + } elif (termType == 3) { + nCoeffs = 8; + } else { + nCoeffs = 1; + } + + mutable coeffs = [0.0, size = nCoeffs]; + + // Return the expanded array of coefficients + if termType == 0 or termType == 1 { + coeffs w/= 0 <- coeff[0]; + } elif termType == 2 or termType == 3 { + for i in 0..nCoeffs - 1 { + coeffs w/= i <- coeff[i / 2]; + } + } + + return coeffs; +} diff --git a/library/chemistry/src/JordanWigner/OptimizedBEOperator.qs b/library/chemistry/src/JordanWigner/OptimizedBEOperator.qs new file mode 100644 index 0000000000..4c562ec7cf --- /dev/null +++ b/library/chemistry/src/JordanWigner/OptimizedBEOperator.qs @@ -0,0 +1,149 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. + +export OptimizedBEXY; + +import JordanWigner.Utils.MultiplexOperationsFromGenerator; +import Std.Arrays.IndexRange; +import Std.Arrays.Partitioned; +import Std.Math.Max; + +/// # Summary +/// Applies a sequence of Z operations and either an X or Y operation to +/// a register of qubits, where the selection of target qubits and basis +/// are conditioned on the state of a control register. +/// +/// # Description +/// This operation can be described by a unitary matrix $U$ that applies +/// the Pauli string on $(X^{z+1}\_pY^{z}\_p)Z\_{p-1}...Z_0$ on +/// qubits $0..p$ conditioned on an index $z\in\{0,1\}$ and $p$. +/// +/// That is, +/// $$ +/// \begin{align} +/// U\ket{z}\ket{p}\ket{\psi} = \ket{z}\ket{p}(X^{z+1}\_pY^{z}\_p)Z\_{p-1}...Z_0\ket{\psi} +/// \end{align} +/// $$ +/// +/// # Input +/// ## pauliBasis +/// When this qubit is in state $\ket{0}$, an `X` operation is applied. When it is in state $\ket{1}$, `Y` is applied. +/// ## indexRegister +/// The state $\ket{p}$ of this register determines the qubit on which `X` or `Y` is applied. +/// ## targetRegister +/// Register of qubits on which the Pauli operators are applied. +/// +/// # References +/// - [Encoding Electronic Spectra in Quantum Circuits with Linear T Complexity](https://arxiv.org/abs/1805.03662) +/// Ryan Babbush, Craig Gidney, Dominic W. Berry, Nathan Wiebe, Jarrod McClean, Alexandru Paler, Austin Fowler, Hartmut Neven +operation OptimizedBEXY(pauliBasis : Qubit, indexRegister : Qubit[], targetRegister : Qubit[]) : Unit is Adj + Ctl { + let unitaryGenerator = (Length(targetRegister), idx -> OptimizedBEXYImpl(idx, _, _, _)); + + use accumulator = Qubit(); + + // This assumes that MultiplexOperationsFromGenerator applies unitaries indexed in unitaryGenerator in ascending order. + X(accumulator); + MultiplexOperationsFromGenerator(unitaryGenerator, indexRegister, (pauliBasis, accumulator, targetRegister)); + // If indexRegister encodes an integer that is larger than Length(targetRegister), + // _OptimizedBEXY_ will fail due to an out of range error. In this situation, + // releasing the accumulator qubit will throw an error as it will be in the One state. +} + +// Subroutine of OptimizedBEXY. +operation OptimizedBEXYImpl(targetIndex : Int, pauliBasis : Qubit, accumulator : Qubit, targetRegister : Qubit[]) : Unit is Adj + Ctl { + + body (...) { + // This should always be called as a controlled operation. + fail "_OptimizedBEXY should always be called as a controlled operation."; + } + + controlled (ctrl, ...) { + if Length(targetRegister) <= targetIndex { + fail "targetIndex out of range."; + } + + Controlled X(ctrl, accumulator); + within { + Controlled Adjoint S([pauliBasis], targetRegister[targetIndex]); + } apply { + Controlled X(ctrl, targetRegister[targetIndex]); + } + Controlled Z([accumulator], targetRegister[targetIndex]); + } + +} + +/// # Summary +/// Applies a Z operation to a qubit indicated by the state of another +/// register. +/// +/// # Description +/// The operation can be represented by a unitary matrix $U$ that applies +/// the `Std.Intrinsic.Z` operation on a qubit $p$ +/// conditioned on an index state $\ket{p}$. That is, +/// $$ +/// \begin{align} +/// U\ket{p}\ket{\psi} = \ket{p}Z\_p\ket{\psi}. +/// \end{align} +/// $$ +/// +/// # Input +/// ## indexRegister +/// A register in the state $\ket{p}$, determining the qubit on which $Z$ is applied. +/// ## targetRegister +/// Register of qubits on which the Pauli operators are applied. +operation SelectZ(indexRegister : Qubit[], targetRegister : Qubit[]) : Unit is Adj + Ctl { + let unitaryGenerator = (Length(targetRegister), idx -> (qs => Z(qs[idx]))); + MultiplexOperationsFromGenerator(unitaryGenerator, indexRegister, targetRegister); + // If indexRegister encodes an integer that is larger than Length(targetRegister), + // _SelectZ_ will fail due to an out of range error. In this situation, + // releasing the accumulator qubit will throw an error as it will be in the One state. +} + +operation JordanWignerSelect( + signQubit : Qubit, + selectZControlRegisters : Qubit[], + OptimizedBEControlRegisters : Qubit[], + pauliBases : Qubit[], + indexRegisters : Qubit[][], + targetRegister : Qubit[] +) : Unit is Adj + Ctl { + Z(signQubit); + + for idxRegister in IndexRange(OptimizedBEControlRegisters) { + Controlled OptimizedBEXY([OptimizedBEControlRegisters[idxRegister]], (pauliBases[idxRegister], indexRegisters[idxRegister], targetRegister)); + } + + for idxRegister in IndexRange(selectZControlRegisters) { + Controlled SelectZ([selectZControlRegisters[idxRegister]], (indexRegisters[idxRegister], targetRegister)); + } +} + +function JordanWignerSelectQubitCount(nZ : Int, nMaj : Int, nIdxRegQubits : Int) : (Int, (Int, Int, Int, Int, Int[])) { + let signQubit = 1; + let selectZControlRegisters = nZ; + let OptimizedBEControlRegisters = nMaj; + let pauliBases = nMaj; + let indexRegisters = Repeated(nIdxRegQubits, Max([nZ, nMaj])); + let nTotal = ((1 + nZ) + 2 * nMaj) + Max([nZ, nMaj]) * nIdxRegQubits; + return (nTotal, (signQubit, selectZControlRegisters, OptimizedBEControlRegisters, pauliBases, indexRegisters)); +} + +function JordanWignerSelectQubitManager(nZ : Int, nMaj : Int, nIdxRegQubits : Int, ctrlRegister : Qubit[], targetRegister : Qubit[]) : ((Qubit, Qubit[], Qubit[], Qubit[], Qubit[][], Qubit[]), Qubit[]) { + let (nTotal, (a, b, c, d, e)) = JordanWignerSelectQubitCount(nZ, nMaj, nIdxRegQubits); + let split = [a, b, c, d] + e; + let registers = Partitioned(split, ctrlRegister); + let signQubit = registers[0]; + let selectZControlRegisters = registers[1]; + let OptimizedBEControlRegisters = registers[2]; + let pauliBases = registers[3]; + let indexRegistersTmp = registers[4..(4 + Length(e)) - 1]; + let rest = registers[Length(registers) - 1]; + mutable indexRegisters = []; + + for idx in IndexRange(e) { + indexRegisters += [indexRegistersTmp[idx]]; + } + + return ((signQubit[0], selectZControlRegisters, OptimizedBEControlRegisters, pauliBases, indexRegisters, targetRegister), rest); +} diff --git a/library/chemistry/src/JordanWigner/StatePreparation.qs b/library/chemistry/src/JordanWigner/StatePreparation.qs new file mode 100644 index 0000000000..894f8f5f68 --- /dev/null +++ b/library/chemistry/src/JordanWigner/StatePreparation.qs @@ -0,0 +1,206 @@ +import Std.StatePreparation.ApproximatelyPreparePureStateCP; +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +export + PrepareSparseMultiConfigurationalState, + PrepareUnitaryCoupledClusterState; + +import JordanWigner.JordanWignerClusterOperatorEvolutionSet.JordanWignerClusterOperatorEvolutionSet; +import JordanWigner.JordanWignerClusterOperatorEvolutionSet.JordanWignerClusterOperatorGeneratorSystem; +import JordanWigner.Utils.JordanWignerInputState; +import JordanWigner.Utils.TrotterSimulationAlgorithm; +import Std.Arrays.*; +import Std.Convert.ComplexAsComplexPolar; +import Std.Convert.IntAsDouble; +import Std.StatePreparation.PreparePureStateD; +import Std.Math.*; +import Utils.EvolutionGenerator; + +operation PrepareTrialState(stateData : (Int, JordanWignerInputState[]), qubits : Qubit[]) : Unit { + let (stateType, terms) = stateData; + + // State type indexing from FermionHamiltonianStatePrep + // public enum StateType + //{ + // Default = 0, Single_Configurational = 1, Sparse_Multi_Configurational = 2, Unitary_Coupled_Cluster = 3 + //} + + if stateType == 2 { + // Sparse_Multi_Configurational + if IsEmpty(terms) { + // Do nothing, as there are no terms to prepare. + } elif Length(terms) == 1 { + let (_, qubitIndices) = terms[0]!; + PrepareSingleConfigurationalStateSingleSiteOccupation(qubitIndices, qubits); + } else { + PrepareSparseMultiConfigurationalState(qs => I(qs[0]), terms, qubits); + } + } elif stateType == 3 { + // Unitary_Coupled_Cluster + let nTerms = Length(terms); + let trotterStepSize = 1.0; + + // The last term is the reference state. + let referenceState = PrepareTrialState((2, [terms[nTerms - 1]]), _); + + PrepareUnitaryCoupledClusterState(referenceState, terms[...nTerms - 2], trotterStepSize, qubits); + } else { + fail ("Unsupported input state."); + } +} + +/// # Summary +/// Simple state preparation of trial state by occupying +/// spin-orbitals +/// +/// # Input +/// ## qubitIndices +/// Indices of qubits to be occupied by electrons. +/// ## qubits +/// Qubits of Hamiltonian. +operation PrepareSingleConfigurationalStateSingleSiteOccupation(qubitIndices : Int[], qubits : Qubit[]) : Unit is Adj + Ctl { + ApplyToEachCA(X, Subarray(qubitIndices, qubits)); +} + +function PrepareSingleConfigurationalStateSingleSiteOccupationWrapper(qubitIndices : Int[]) : (Qubit[] => Unit is Adj + Ctl) { + return PrepareSingleConfigurationalStateSingleSiteOccupation(qubitIndices, _); +} + +/// # Summary +/// Sparse multi-configurational state preparation of trial state by adding excitations +/// to initial trial state. +/// +/// # Input +/// ## initialStatePreparation +/// Unitary to prepare initial trial state. +/// ## excitations +/// Excitations of initial trial state represented by +/// the amplitude of the excitation and the qubit indices +/// the excitation acts on. +/// ## qubits +/// Qubits of Hamiltonian. +operation PrepareSparseMultiConfigurationalState( + initialStatePreparation : (Qubit[] => Unit), + excitations : JordanWignerInputState[], + qubits : Qubit[] +) : Unit { + let nExcitations = Length(excitations); + + mutable coefficientsSqrtAbs = [0.0, size = nExcitations]; + mutable coefficientsNewComplexPolar = Repeated(new ComplexPolar { Magnitude = 0.0, Argument = 0.0 }, nExcitations); + mutable applyFlips = [[], size = nExcitations]; + + for idx in 0..nExcitations - 1 { + let ((r, i), excitation) = excitations[idx]!; + coefficientsSqrtAbs w/= idx <- Sqrt(AbsComplexPolar(ComplexAsComplexPolar(new Complex { Real = r, Imag = i }))); + coefficientsNewComplexPolar w/= idx <- new ComplexPolar { + Magnitude = coefficientsSqrtAbs[idx], + Argument = ArgComplexPolar(ComplexAsComplexPolar(new Complex { Real = r, Imag = i })) + }; + applyFlips w/= idx <- excitation; + } + + let nBitsIndices = Ceiling(Lg(IntAsDouble(nExcitations))); + + mutable success = false; + repeat { + use auxillary = Qubit[nBitsIndices + 1]; + use flag = Qubit(); + + let arr = Mapped(PrepareSingleConfigurationalStateSingleSiteOccupationWrapper, applyFlips); + let multiplexer = MultiplexerBruteForceFromGenerator(nExcitations, idx -> arr[idx]); + ApproximatelyPreparePureStateCP(0.0, coefficientsNewComplexPolar, Reversed(auxillary)); + multiplexer(auxillary, qubits); + Adjoint PreparePureStateD(coefficientsSqrtAbs, Reversed(auxillary)); + ApplyControlledOnInt(0, X, auxillary, flag); + + // if measurement outcome one we prepared required state + let outcome = M(flag); + success = outcome == One; + ResetAll(auxillary); + Reset(flag); + } until success + fixup { + ResetAll(qubits); + } +} + +/// # Summary +/// Unitary coupled-cluster state preparation of trial state +/// +/// # Input +/// ## initialStatePreparation +/// Unitary to prepare initial trial state. +/// ## qubits +/// Qubits of Hamiltonian. +operation PrepareUnitaryCoupledClusterState( + initialStatePreparation : (Qubit[] => Unit), + clusterOperator : JordanWignerInputState[], + trotterStepSize : Double, + qubits : Qubit[] +) : Unit { + let clusterOperatorGeneratorSystem = JordanWignerClusterOperatorGeneratorSystem(clusterOperator); + let evolutionGenerator = new EvolutionGenerator { + EvolutionSet = JordanWignerClusterOperatorEvolutionSet(), + System = clusterOperatorGeneratorSystem + }; + let trotterOrder = 1; + let simulationAlgorithm = TrotterSimulationAlgorithm(trotterStepSize, trotterOrder); + let oracle = simulationAlgorithm(1.0, evolutionGenerator, _); + initialStatePreparation(qubits); + oracle(qubits); +} + +/// # Summary +/// Returns a multiply-controlled unitary operation $U$ that applies a +/// unitary $V_j$ when controlled by n-qubit number state $\ket{j}$. +/// +/// $U = \sum^{2^n-1}_{j=0}\ket{j}\bra{j}\otimes V_j$. +/// +/// # Input +/// ## unitaryGenerator +/// A tuple where the first element `Int` is the number of unitaries $N$, +/// and the second element `(Int -> ('T => () is Adj + Ctl))` +/// is a function that takes an integer $j$ in $[0,N-1]$ and outputs the unitary +/// operation $V_j$. +/// +/// # Output +/// A multiply-controlled unitary operation $U$ that applies unitaries +/// described by `unitaryGenerator`. +function MultiplexerBruteForceFromGenerator(unitaryGenerator : (Int, (Int -> (Qubit[] => Unit is Adj + Ctl)))) : ((Qubit[], Qubit[]) => Unit is Adj + Ctl) { + return MultiplexOperationsBruteForceFromGenerator(unitaryGenerator, _, _); +} + +/// # Summary +/// Applies multiply-controlled unitary operation $U$ that applies a +/// unitary $V_j$ when controlled by n-qubit number state $\ket{j}$. +/// +/// $U = \sum^{N-1}_{j=0}\ket{j}\bra{j}\otimes V_j$. +/// +/// # Input +/// ## unitaryGenerator +/// A tuple where the first element `Int` is the number of unitaries $N$, +/// and the second element `(Int -> ('T => () is Adj + Ctl))` +/// is a function that takes an integer $j$ in $[0,N-1]$ and outputs the unitary +/// operation $V_j$. +/// +/// ## index +/// $n$-qubit control register that encodes number states $\ket{j}$ in +/// little-endian format. +/// +/// ## target +/// Generic qubit register that $V_j$ acts on. +/// +/// # Remarks +/// `coefficients` will be padded with identity elements if +/// fewer than $2^n$ are specified. This version is implemented +/// directly by looping through n-controlled unitary operators. +operation MultiplexOperationsBruteForceFromGenerator<'T>(unitaryGenerator : (Int, (Int -> ('T => Unit is Adj + Ctl))), index : Qubit[], target : 'T) : Unit is Adj + Ctl { + let nIndex = Length(index); + let nStates = 2^nIndex; + let (nUnitaries, unitaryFunction) = unitaryGenerator; + for idxOp in 0..MinI(nStates, nUnitaries) - 1 { + ApplyControlledOnInt(idxOp, unitaryFunction(idxOp), index, target); + } +} diff --git a/library/chemistry/src/JordanWigner/Utils.qs b/library/chemistry/src/JordanWigner/Utils.qs new file mode 100644 index 0000000000..3f5345742d --- /dev/null +++ b/library/chemistry/src/JordanWigner/Utils.qs @@ -0,0 +1,453 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. + +export + JWOptimizedHTerms, + JordanWignerInputState, + JordanWignerEncodingData, + DecomposedIntoTimeStepsCA, + TrotterStep, + TrotterSimulationAlgorithm, + MultiplexOperationsFromGenerator; + +import Std.Arrays.*; +import Std.Convert.IntAsDouble; +import Std.Math.*; +import Utils.EvolutionGenerator; + +/// # Summary +/// Format for a four term Jordan-Wigner optimized Hamiltonian. +/// The meaning of the data represented is determined by the algorithm that receives it. +struct JWOptimizedHTerms { + HTerm0 : (Int[], Double[])[], + HTerm1 : (Int[], Double[])[], + HTerm2 : (Int[], Double[])[], + HTerm3 : (Int[], Double[])[], +} + +/// # Summary +/// Represents preparation of the initial state +/// The meaning of the data represented is determined by the algorithm that receives it. +struct JordanWignerInputState { + Amplitude : (Double, Double), + FermionIndices : Int[], +} + +/// # Summary +/// Format of data to represent all information for Hamiltonian simulation. +/// The meaning of the data represented is determined by the algorithm that receives it. +struct JordanWignerEncodingData { + NumQubits : Int, + Terms : JWOptimizedHTerms, + InputState : (Int, JordanWignerInputState[]), + EnergyOffset : Double, +} + +/// # Summary +/// Implementation of the first-order Trotter–Suzuki integrator. +/// +/// # Type Parameters +/// ## 'T +/// The type which each time step should act upon; typically, either +/// `Qubit[]` or `Qubit`. +/// +/// # Input +/// ## nSteps +/// The number of operations to be decomposed into time steps. +/// ## op +/// An operation which accepts an index input (type `Int`) and a time +/// input (type `Double`) and a quantum register (type `'T`) for decomposition. +/// ## stepSize +/// Multiplier on size of each step of the simulation. +/// ## target +/// A quantum register on which the operations act. +/// +/// # Example +/// The following are equivalent: +/// ```qsharp +/// op(0, deltaT, target); +/// op(1, deltaT, target); +/// ``` +/// and +/// ```qsharp +/// Trotter1ImplCA((2, op), deltaT, target); +/// ``` +operation Trotter1ImplCA<'T>((nSteps : Int, op : ((Int, Double, 'T) => Unit is Adj + Ctl)), stepSize : Double, target : 'T) : Unit is Adj + Ctl { + for idx in 0..nSteps - 1 { + op(idx, stepSize, target); + } +} + +/// # Summary +/// Implementation of the second-order Trotter–Suzuki integrator. +/// +/// # Type Parameters +/// ## 'T +/// The type which each time step should act upon; typically, either +/// `Qubit[]` or `Qubit`. +/// +/// # Input +/// ## nSteps +/// The number of operations to be decomposed into time steps. +/// ## op +/// An operation which accepts an index input (type `Int`) and a time +/// input (type `Double`) and a quantum register (type `'T`) for decomposition. +/// ## stepSize +/// Multiplier on size of each step of the simulation. +/// ## target +/// A quantum register on which the operations act. +/// +/// # Example +/// The following are equivalent: +/// ```qsharp +/// op(0, deltaT / 2.0, target); +/// op(1, deltaT / 2.0, target); +/// op(1, deltaT / 2.0, target); +/// op(0, deltaT / 2.0, target); +/// ``` +/// and +/// ```qsharp +/// Trotter2ImplCA((2, op), deltaT, target); +/// ``` +operation Trotter2ImplCA<'T>( + (nSteps : Int, op : ((Int, Double, 'T) => Unit is Adj + Ctl)), + stepSize : Double, + target : 'T +) : Unit is Adj + Ctl { + for idx in 0..nSteps - 1 { + op(idx, stepSize * 0.5, target); + } + + for idx in nSteps - 1.. -1..0 { + op(idx, stepSize * 0.5, target); + } +} + +/// # Summary +/// Computes Trotter step size in recursive implementation of +/// Trotter simulation algorithm. +/// +/// # Remarks +/// This operation uses a different indexing convention than that of +/// [quant-ph/0508139](https://arxiv.org/abs/quant-ph/0508139). In +/// particular, `DecomposedIntoTimeStepsCA(_, 4)` corresponds to the +/// scalar $p_2(\lambda)$ in quant-ph/0508139. +function TrotterStepSize(order : Int) : Double { + return 1.0 / (4.0 - 4.0^(1.0 / (IntAsDouble(order) - 1.0))); +} + + +/// # Summary +/// Recursive implementation of even-order Trotter–Suzuki integrator. +/// +/// # Type Parameters +/// ## 'T +/// The type which each time step should act upon; typically, either +/// `Qubit[]` or `Qubit`. +/// +/// # Input +/// ## order +/// Order of Trotter-Suzuki integrator. +/// ## nSteps +/// The number of operations to be decomposed into time steps. +/// ## op +/// An operation which accepts an index input (type `Int`) and a time +/// input (type `Double`) and a quantum register (type `'T`) for decomposition. +/// ## stepSize +/// Multiplier on size of each step of the simulation. +/// ## target +/// A quantum register on which the operations act. +operation TrotterArbitraryImplCA<'T>( + order : Int, + (nSteps : Int, op : ((Int, Double, 'T) => Unit is Adj + Ctl)), + stepSize : Double, + target : 'T +) : Unit is Adj + Ctl { + if (order > 2) { + let stepSizeOuter = TrotterStepSize(order); + let stepSizeInner = 1.0 - 4.0 * stepSizeOuter; + TrotterArbitraryImplCA(order - 2, (nSteps, op), stepSizeOuter * stepSize, target); + TrotterArbitraryImplCA(order - 2, (nSteps, op), stepSizeOuter * stepSize, target); + TrotterArbitraryImplCA(order - 2, (nSteps, op), stepSizeInner * stepSize, target); + TrotterArbitraryImplCA(order - 2, (nSteps, op), stepSizeOuter * stepSize, target); + TrotterArbitraryImplCA(order - 2, (nSteps, op), stepSizeOuter * stepSize, target); + } elif (order == 2) { + Trotter2ImplCA((nSteps, op), stepSize, target); + } else { + Trotter1ImplCA((nSteps, op), stepSize, target); + } +} + +/// # Summary +/// Returns an operation implementing the Trotter–Suzuki integrator for +/// a given operation. +/// +/// # Type Parameters +/// ## 'T +/// The type which each time step should act upon; typically, either +/// `Qubit[]` or `Qubit`. +/// +/// # Input +/// ## nSteps +/// The number of operations to be decomposed into time steps. +/// ## op +/// An operation which accepts an index input (type `Int`) and a time +/// input (type `Double`) for decomposition. +/// ## trotterOrder +/// Selects the order of the Trotter–Suzuki integrator to be used. +/// Order 1 and even orders 2, 4, 6,... are currently supported. +/// +/// # Output +/// Returns a unitary implementing the Trotter–Suzuki integrator, where +/// the first parameter `Double` is the integration step size, and the +/// second parameter is the target acted upon. +/// +/// # Remarks +/// When called with `order` equal to `1`, this function returns an operation +/// that can be simulated by the lowest-order Trotter–Suzuki integrator +/// $$ +/// \begin{align} +/// S_1(\lambda) = \prod_{j = 1}^{m} e^{H_j \lambda}, +/// \end{align} +/// $$ +/// where we have followed the notation of +/// [quant-ph/0508139](https://arxiv.org/abs/quant-ph/0508139) +/// and let $\lambda$ be the evolution time (represented by the first input +/// of the returned operation), and have let $\{H_j\}_{j = 1}^{m}$ be the +/// set of (skew-Hermitian) dynamical generators being integrated such that +/// `op(j, lambda, _)` is simulated by the unitary operator +/// $e^{H_j \lambda}$. +/// +/// Similarly, an `order` of `2` returns the second-order symmetric +/// Trotter–Suzuki integrator +/// $$ +/// \begin{align} +/// S_2(\lambda) = \prod_{j = 1}^{m} e^{H_k \lambda / 2} +/// \prod_{j' = m}^{1} e^{H_{j'} \lambda / 2}. +/// \end{align} +/// $$ +/// +/// Higher even values of `order` are implemented using the recursive +/// construction of [quant-ph/0508139](https://arxiv.org/abs/quant-ph/0508139). +/// +/// # References +/// - [ *D. W. Berry, G. Ahokas, R. Cleve, B. C. Sanders* ](https://arxiv.org/abs/quant-ph/0508139) +function DecomposedIntoTimeStepsCA<'T>( + (nSteps : Int, op : ((Int, Double, 'T) => Unit is Adj + Ctl)), + trotterOrder : Int +) : ((Double, 'T) => Unit is Adj + Ctl) { + if (trotterOrder == 1) { + return Trotter1ImplCA((nSteps, op), _, _); + } elif (trotterOrder == 2) { + return Trotter2ImplCA((nSteps, op), _, _); + } elif (trotterOrder % 2 == 0) { + return TrotterArbitraryImplCA(trotterOrder, (nSteps, op), _, _); + } else { + fail $"Odd order {trotterOrder} not yet supported."; + } +} + +/// # Summary +/// Implements time-evolution by a term contained in a `GeneratorSystem`. +/// +/// # Input +/// ## evolutionGenerator +/// A complete description of the system to be simulated. +/// ## idx +/// Integer index to a term in the described system. +/// ## stepsize +/// Multiplier on duration of time-evolution by term indexed by `idx`. +/// ## qubits +/// Qubits acted on by simulation. +operation TrotterStepImpl(evolutionGenerator : EvolutionGenerator, idx : Int, stepsize : Double, qubits : Qubit[]) : Unit is Adj + Ctl { + let (evolutionSet, generatorSystem) = evolutionGenerator!; + let (nTerms, generatorSystemFunction) = generatorSystem!; + let generatorIndex = generatorSystemFunction(idx); + (evolutionSet(generatorIndex))(stepsize, qubits); +} + +/// # Summary +/// Implements a single time-step of time-evolution by the system +/// described in an `EvolutionGenerator` using a Trotter–Suzuki +/// decomposition. +/// +/// # Input +/// ## evolutionGenerator +/// A complete description of the system to be simulated. +/// ## trotterOrder +/// Order of Trotter integrator. This must be either 1 or an even number. +/// ## trotterStepSize +/// Duration of simulated time-evolution in single Trotter step. +/// +/// # Output +/// Unitary operation that approximates a single step of time-evolution +/// for duration `trotterStepSize`. +function TrotterStep(evolutionGenerator : EvolutionGenerator, trotterOrder : Int, trotterStepSize : Double) : (Qubit[] => Unit is Adj + Ctl) { + let (evolutionSet, generatorSystem) = evolutionGenerator!; + let (nTerms, generatorSystemFunction) = generatorSystem!; + + // The input to DecomposeIntoTimeStepsCA has signature + // (Int, ((Int, Double, Qubit[]) => () is Adj + Ctl)) + let trotterForm = (nTerms, TrotterStepImpl(evolutionGenerator, _, _, _)); + return (DecomposedIntoTimeStepsCA(trotterForm, trotterOrder))(trotterStepSize, _); +} + +/// # Summary +/// Makes repeated calls to `TrotterStep` to approximate the +/// time-evolution operator exp(_-iHt_). +/// +/// # Input +/// ## trotterStepSize +/// Duration of simulated time-evolution in single Trotter step. +/// ## trotterOrder +/// Order of Trotter integrator. This must be either 1 or an even number. +/// ## maxTime +/// Total duration of simulation $t$. +/// ## evolutionGenerator +/// A complete description of the system to be simulated. +/// ## qubits +/// Qubits acted on by simulation. +operation TrotterSimulationAlgorithmImpl(trotterStepSize : Double, trotterOrder : Int, maxTime : Double, evolutionGenerator : EvolutionGenerator, qubits : Qubit[]) : Unit is Adj + Ctl { + let nTimeSlices = Ceiling(maxTime / trotterStepSize); + let resizedTrotterStepSize = maxTime / IntAsDouble(nTimeSlices); + + for idxTimeSlice in 0..nTimeSlices - 1 { + (TrotterStep(evolutionGenerator, trotterOrder, resizedTrotterStepSize))(qubits); + } +} + +/// # Summary +/// `SimulationAlgorithm` function that uses a Trotter–Suzuki +/// decomposition to approximate the time-evolution operator _exp(-iHt)_. +/// +/// # Input +/// ## trotterStepSize +/// Duration of simulated time-evolution in single Trotter step. +/// ## trotterOrder +/// Order of Trotter integrator. This must be either 1 or an even number. +/// +/// # Output +/// A `SimulationAlgorithm` type. +function TrotterSimulationAlgorithm(trotterStepSize : Double, trotterOrder : Int) : (Double, EvolutionGenerator, Qubit[]) => Unit is Adj + Ctl { + return TrotterSimulationAlgorithmImpl(trotterStepSize, trotterOrder, _, _, _); +} + +/// # Summary +/// Applies a multiply-controlled unitary operation $U$ that applies a +/// unitary $V_j$ when controlled by n-qubit number state $\ket{j}$. +/// +/// $U = \sum^{N-1}_{j=0}\ket{j}\bra{j}\otimes V_j$. +/// +/// # Input +/// ## unitaryGenerator +/// A tuple where the first element `Int` is the number of unitaries $N$, +/// and the second element `(Int -> ('T => () is Adj + Ctl))` +/// is a function that takes an integer $j$ in $[0,N-1]$ and outputs the unitary +/// operation $V_j$. +/// +/// ## index +/// $n$-qubit control register that encodes number states $\ket{j}$ in +/// little-endian format. +/// +/// ## target +/// Generic qubit register that $V_j$ acts on. +/// +/// # Remarks +/// `coefficients` will be padded with identity elements if +/// fewer than $2^n$ are specified. This implementation uses +/// $n-1$ auxiliary qubits. +/// +/// # References +/// - [ *Andrew M. Childs, Dmitri Maslov, Yunseong Nam, Neil J. Ross, Yuan Su*, +/// arXiv:1711.10980](https://arxiv.org/abs/1711.10980) +operation MultiplexOperationsFromGenerator<'T>(unitaryGenerator : (Int, (Int -> ('T => Unit is Adj + Ctl))), index : Qubit[], target : 'T) : Unit is Ctl + Adj { + let (nUnitaries, unitaryFunction) = unitaryGenerator; + let unitaryGeneratorWithOffset = (nUnitaries, 0, unitaryFunction); + if Length(index) == 0 { + fail "MultiplexOperations failed. Number of index qubits must be greater than 0."; + } + if nUnitaries > 0 { + let auxiliary = []; + Adjoint MultiplexOperationsFromGeneratorImpl(unitaryGeneratorWithOffset, auxiliary, index, target); + } +} + +/// # Summary +/// Implementation step of `MultiplexOperationsFromGenerator`. +operation MultiplexOperationsFromGeneratorImpl<'T>(unitaryGenerator : (Int, Int, (Int -> ('T => Unit is Adj + Ctl))), auxiliary : Qubit[], index : Qubit[], target : 'T) : Unit { + body (...) { + let nIndex = Length(index); + let nStates = 2^nIndex; + + let (nUnitaries, unitaryOffset, unitaryFunction) = unitaryGenerator; + + let nUnitariesLeft = MinI(nUnitaries, nStates / 2); + let nUnitariesRight = MinI(nUnitaries, nStates); + + let leftUnitaries = (nUnitariesLeft, unitaryOffset, unitaryFunction); + let rightUnitaries = (nUnitariesRight - nUnitariesLeft, unitaryOffset + nUnitariesLeft, unitaryFunction); + + let newControls = Most(index); + + if nUnitaries > 0 { + if Length(auxiliary) == 1 and nIndex == 0 { + // Termination case + + (Controlled Adjoint (unitaryFunction(unitaryOffset)))(auxiliary, target); + } elif Length(auxiliary) == 0 and nIndex >= 1 { + // Start case + let newauxiliary = Tail(index); + if nUnitariesRight > 0 { + MultiplexOperationsFromGeneratorImpl(rightUnitaries, [newauxiliary], newControls, target); + } + within { + X(newauxiliary); + } apply { + MultiplexOperationsFromGeneratorImpl(leftUnitaries, [newauxiliary], newControls, target); + } + } else { + // Recursion that reduces nIndex by 1 and sets Length(auxiliary) to 1. + let controls = [Tail(index)] + auxiliary; + use newauxiliary = Qubit(); + use andauxiliary = Qubit[MaxI(0, Length(controls) - 2)]; + within { + if Length(controls) == 0 { + X(newauxiliary); + } elif Length(controls) == 1 { + CNOT(Head(controls), newauxiliary); + } else { + let controls1 = controls[0..0] + andauxiliary; + let controls2 = Rest(controls); + let targets = andauxiliary + [newauxiliary]; + for i in IndexRange(controls1) { + AND(controls1[i], controls2[i], targets[i]); + } + } + + } apply { + if nUnitariesRight > 0 { + MultiplexOperationsFromGeneratorImpl(rightUnitaries, [newauxiliary], newControls, target); + } + within { + (Controlled X)(auxiliary, newauxiliary); + } apply { + MultiplexOperationsFromGeneratorImpl(leftUnitaries, [newauxiliary], newControls, target); + } + } + } + } + } + adjoint auto; + controlled (controlRegister, ...) { + MultiplexOperationsFromGeneratorImpl(unitaryGenerator, auxiliary + controlRegister, index, target); + } + controlled adjoint auto; +} + +function RangeAsIntArray(range : Range) : Int[] { + mutable arr = []; + for i in range { + arr += [i]; + } + return arr; +} diff --git a/library/chemistry/src/Tests.qs b/library/chemistry/src/Tests.qs new file mode 100644 index 0000000000..2b44b29239 --- /dev/null +++ b/library/chemistry/src/Tests.qs @@ -0,0 +1,382 @@ +import Std.StatePreparation.ApproximatelyPreparePureStateCP; +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. + +import JordanWigner.JordanWignerOptimizedBlockEncoding.PrepareUniformSuperposition; +import JordanWigner.JordanWignerClusterOperatorEvolutionSet.JordanWignerClusterOperatorPQRSTermSigns; +import JordanWigner.OptimizedBEOperator.OptimizedBEXY; +import JordanWigner.OptimizedBEOperator.SelectZ; +import JordanWigner.StatePreparation.PrepareSparseMultiConfigurationalState; +import JordanWigner.StatePreparation.PrepareUnitaryCoupledClusterState; +import JordanWigner.Utils.JordanWignerInputState; +import Std.Arrays.IndexRange; +import Std.Convert.ComplexAsComplexPolar; +import Std.Convert.IntAsDouble; +import Std.Diagnostics.CheckAllZero; +import Std.Diagnostics.CheckZero; +import Std.Diagnostics.DumpRegister; +import Std.Diagnostics.Fact; +import Std.Math.Ceiling; +import Std.Math.Complex; +import Std.Math.ComplexPolar; +import Std.Math.Lg; +import Std.Arrays.Reversed; +import Std.Math.Sqrt; + +@Config(Unrestricted) +@Test() +operation PrepareUniformSuperpositionTest() : Unit { + let NBits = 4; + use qs = Qubit[NBits]; + for NStates in 1..2^NBits-1 { + PrepareUniformSuperposition(NStates, qs); + let t = Std.Math.Sqrt(1.0 / IntAsDouble(NStates)); + Adjoint Std.StatePreparation.PreparePureStateD(Repeated(t, NStates), Reversed(qs)); + Fact(CheckAllZero(qs), "All Z"); + ResetAll(qs); + } +} + +@Config(Unrestricted) +@Test() +operation PrepareSparseMultiConfigurationalState0Test() : Unit { + let nQubits = 6; + let expectedResult = 39; + let excitations = [new JordanWignerInputState { Amplitude = (0.1, 0.0), FermionIndices = [0, 1, 2, 5] }]; + + use qubits = Qubit[nQubits]; + PrepareSparseMultiConfigurationalState(qs => I(qs[0]), excitations, qubits); + + Fact(MeasureInteger(qubits) == expectedResult, "PrepareSparseMultiConfigurationalState0Test failed."); +} + +@Config(Unrestricted) +operation OptimizedBEOperatorZeroTestHelper(pauliBasis : Pauli, targetRegisterSize : Int, targetIndex : Int) : Unit { + let indexRegisterSize = Ceiling(Lg(IntAsDouble(targetRegisterSize))); + use pauliBasisQubit = Qubit(); + use indexRegister = Qubit[indexRegisterSize]; + use targetRegister = Qubit[targetRegisterSize]; + + // Choose X or Y operator. + if pauliBasis == PauliX { + // no op + } elif pauliBasis == PauliY { + X(pauliBasisQubit); + } + + // Create indexRegister state. + ApplyXorInPlace(targetIndex, indexRegister); + + // Initialize targetRegister states in |0> + OptimizedBEXY(pauliBasisQubit, indexRegister, targetRegister); + + for idxTest in 0..targetRegisterSize - 1 { + let testQubit = targetRegister[idxTest]; + + if targetIndex > idxTest { + + // Test Z Pauli + // |0> -> |0> + // |+> -> |-> + Message($"Test Z Pauli on qubit {idxTest}"); + Fact(CheckZero(testQubit), $"Error: Test {idxTest} {idxTest} Z Pauli |0>"); + } elif targetIndex == idxTest { + + // Test X Pauli + // |0> -> |1> + // |+> -> |+> + + // Test Y Pauli + // |0> -> i|1> + // |+> -> -i|-> + Message($"Test X or Y Pauli on qubit {idxTest}"); + within { + X(testQubit); + } apply { + Fact(CheckZero(testQubit), $"Error: Test {idxTest} X or Y Pauli |0>"); + } + } else { + + // Test Identitfy Pauli + // |0> -> |0> + // |+> -> |+> + Message($"Test ZI Pauli on qubit {idxTest}"); + Fact(CheckZero(testQubit), $"Error: Test {idxTest} I Pauli |0>"); + } + } + + OptimizedBEXY(pauliBasisQubit, indexRegister, targetRegister); + Adjoint ApplyXorInPlace(targetIndex, indexRegister); + + // Choose X or Y operator. + if pauliBasis == PauliX { + // no op + } elif pauliBasis == PauliY { + X(pauliBasisQubit); + } +} + +@Config(Unrestricted) +@Test() +operation OptimizedBEOperatorZeroTest() : Unit { + let paulis = [PauliX, PauliY]; + let targetRegisterSize = 7; + + for idxPauli in 0..1 { + let pauliBasis = paulis[idxPauli]; + + for targetIndex in 0..targetRegisterSize - 1 { + Message($"pauliBasis {pauliBasis}; targetIndex {targetIndex}"); + OptimizedBEOperatorZeroTestHelper(pauliBasis, targetRegisterSize, targetIndex); + } + } +} + +@Config(Unrestricted) +operation OptimizedBEOperatorPlusTestHelper(pauliBasis : Pauli, targetRegisterSize : Int, targetIndex : Int) : Unit { + let indexRegisterSize = Ceiling(Lg(IntAsDouble(targetRegisterSize))); + use pauliBasisQubit = Qubit(); + use indexRegister = Qubit[indexRegisterSize]; + use targetRegister = Qubit[targetRegisterSize]; + // Choose X or Y operator. + if (pauliBasis == PauliX) { + // no op + } elif pauliBasis == PauliY { + X(pauliBasisQubit); + } + + // Create indexRegister state. + ApplyXorInPlace(targetIndex, indexRegister); + + // Initialize targetRegister states in |+> + ApplyToEachCA(H, targetRegister); + OptimizedBEXY(pauliBasisQubit, indexRegister, targetRegister); + for idxTest in 0..targetRegisterSize - 1 { + let testQubit = targetRegister[idxTest]; + if (targetIndex > idxTest) { + // Test Z Pauli + // |0> -> |0> + // |+> -> |-> + Message($"Test Z Pauli on qubit {idxTest}"); + within { + H(testQubit); + X(testQubit); + } apply { + Fact(CheckZero(testQubit), $"Error: Test {idxTest} Z Pauli |->"); + } + } elif (targetIndex == idxTest) { + // Test X Pauli + // |0> -> |1> + // |+> -> |+> + if (pauliBasis == PauliX) { + Message($"Test X Pauli on qubit {idxTest}"); + within { + H(testQubit); + } apply { + Fact(CheckZero(testQubit), $"Error: Test {idxTest} X Pauli |+>"); + } + } + + // Test Y Pauli + // |0> -> i|1> + // |+> -> -i|-> + if (pauliBasis == PauliY) { + Message($"Test Y Pauli on qubit {idxTest}"); + within { + H(testQubit); + X(testQubit); + } apply { + Fact(CheckZero(testQubit), $"Error: Test {idxTest} Y Pauli |+>"); + } + } + + } else { + // Test Identitfy Pauli + // |0> -> |0> + // |+> -> |+> + Message($"Test I Pauli on qubit {idxTest}"); + within { + H(testQubit); + } apply { + Fact(CheckZero(testQubit), $"Error: Test {idxTest} I Pauli |+>"); + } + } + } + OptimizedBEXY(pauliBasisQubit, indexRegister, targetRegister); + ApplyToEachCA(H, targetRegister); + + (Adjoint ApplyXorInPlace)(targetIndex, indexRegister); + + // Choose X or Y operator. + if pauliBasis == PauliX { + // no op + } elif pauliBasis == PauliY { + X(pauliBasisQubit); + } +} + +@Config(Unrestricted) +@Test() +operation OptimizedBEOperatorPlusTest() : Unit { + + let paulis = [PauliX, PauliY]; + let targetRegisterSize = 7; + + for idxPauli in 0..1 { + let pauliBasis = paulis[idxPauli]; + + for targetIndex in 0..targetRegisterSize - 1 { + Message($"pauliBasis {pauliBasis}; targetIndex {targetIndex}"); + OptimizedBEOperatorPlusTestHelper(pauliBasis, targetRegisterSize, targetIndex); + } + } +} + +@Config(Unrestricted) +@Test() +operation SelectZTest() : Unit { + let targetRegisterSize = 7; + let indexRegisterSize = Ceiling(Lg(IntAsDouble(targetRegisterSize))); + + use targetRegister = Qubit[targetRegisterSize]; + use indexRegister = Qubit[indexRegisterSize]; + for idxTest in 0..targetRegisterSize - 1 { + H(targetRegister[idxTest]); + ApplyXorInPlace(idxTest, indexRegister); + SelectZ(indexRegister, targetRegister); + within { + H(targetRegister[idxTest]); + X(targetRegister[idxTest]); + } apply { + Fact(CheckZero(targetRegister[idxTest]), $"Error: Test {idxTest} X Pauli |+>"); + } + Z(targetRegister[idxTest]); + Adjoint ApplyXorInPlace(idxTest, indexRegister); + H(targetRegister[idxTest]); + } +} + +@Config(Unrestricted) +function JordanWignerClusterOperatorPQRSTermSignsTestHelper(idx : Int) : (Int[], Double[], Double) { + let cases = [ + ([1, 2, 3, 4], [1.0,-1.0,-1.0,-1.0, 1.0, 1.0, 1.0,-1.0], 1.0), + ([2, 1, 4, 3], [1.0,-1.0,-1.0,-1.0, 1.0, 1.0, 1.0,-1.0], 1.0), + ([3, 4, 1, 2], [1.0,-1.0,-1.0,-1.0, 1.0, 1.0, 1.0,-1.0],-1.0), + ([2, 1, 3, 4], [1.0,-1.0,-1.0,-1.0, 1.0, 1.0, 1.0,-1.0],-1.0), + ([1, 3, 2, 4], [-1.0,-1.0,-1.0, 1.0,-1.0, 1.0, 1.0, 1.0], 1.0), + ([4, 2, 3, 1], [-1.0,-1.0,-1.0, 1.0,-1.0, 1.0, 1.0, 1.0],-1.0), + ([1, 4, 2, 3], [1.0, 1.0,-1.0, 1.0,-1.0, 1.0,-1.0,-1.0], 1.0), + ([2, 3, 4, 1], [1.0, 1.0,-1.0, 1.0,-1.0, 1.0,-1.0,-1.0], 1.0) + ]; + return cases[idx]; +} + +@Config(Unrestricted) +@Test() +function JordanWignerClusterOperatorPQRSTermSignsTest() : Unit { + for idx in 0..7 { + let (testCase, expectedSigns, expectedGlobalSign) = JordanWignerClusterOperatorPQRSTermSignsTestHelper(idx); + let (sortedIndices, signs, globalSign) = JordanWignerClusterOperatorPQRSTermSigns(testCase); + + let p = sortedIndices[0]; + let q = sortedIndices[1]; + let r = sortedIndices[2]; + let s = sortedIndices[3]; + + Fact(p < q and q < r and r < s, "Expected p < q < r < s."); + NearEqualityFactD(globalSign, expectedGlobalSign); + for i in IndexRange(signs) { + NearEqualityFactD(signs[i], expectedSigns[i]); + } + } +} + +@Config(Unrestricted) +function DoublesToComplexPolar(input : Double[]) : ComplexPolar[] { + mutable arr = [new ComplexPolar { Magnitude = 0.0, Argument = 0.0 }, size = Length(input)]; + for idx in 0..Length(input)-1 { + arr w/= idx <- ComplexAsComplexPolar(new Complex { Real = input[idx], Imag = 0. }); + } + return arr; +} + +@Config(Unrestricted) +operation JordanWignerUCCTermTestHelper(nQubits : Int, excitations : Int[], term : JordanWignerInputState[], result : Double[]) : Unit { + use qubits = Qubit[nQubits]; + for idx in excitations { + X(qubits[idx]); + } + PrepareUnitaryCoupledClusterState(qs => I(qs[0]), term, 1.0, qubits); + DumpRegister(qubits); + Adjoint ApproximatelyPreparePureStateCP(0.0, DoublesToComplexPolar(result), Reversed(qubits)); + Fact(CheckAllZero(qubits), "Expected qubits to all be in ground state."); + ResetAll(qubits); +} + +@Config(Unrestricted) +@Test() +operation JordanWignerUCCSTermTest() : Unit { + // test using Exp(2.0 (a^\dag_1 a_3 - h.c.)) + let term0 = [new JordanWignerInputState { Amplitude = (2.0, 0.0), FermionIndices = [1, 3] }]; + let state0 = [0., 0.,-0.416147, 0., 0., 0., 0., 0.,-0.909297, 0., 0., 0., 0., 0., 0., 0.]; + JordanWignerUCCTermTestHelper(4, [1], term0, state0); + + // test using Exp(2.0 (a^\dag_3 a_1 - h.c.)) + let term1 = [new JordanWignerInputState { Amplitude = (2.0, 0.0), FermionIndices = [3, 1] }]; + let state1 = [0., 0.,-0.416147, 0., 0., 0., 0., 0., 0.909297, 0., 0., 0., 0., 0., 0., 0.]; + JordanWignerUCCTermTestHelper(4, [1], term1, state1); +} + +@Config(Unrestricted) +@Test() +operation JordanWignerUCCDTermPQRSTest() : Unit { + // test using Exp(2.0 (a^\dag_0 a^\dag_1 a_3 a_4 - h.c.)) + let term0 = [new JordanWignerInputState { Amplitude = (2.0, 0.0), FermionIndices = [0, 1, 2, 4] }]; + let state0 = [0., 0., 0.,-0.416147, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.909297, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]; + JordanWignerUCCTermTestHelper(5, [0, 1], term0, state0); + + // test using Exp(2.0 (a^\dag_0 a^\dag_1 a_3 a_4 - h.c.)) + let term1 = [new JordanWignerInputState { Amplitude = (2.0, 0.0), FermionIndices = [0, 1, 2, 4] }]; + let state1 = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,-0.416147, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,-0.909297, 0., 0., 0.]; + JordanWignerUCCTermTestHelper(5, [0, 1, 3], term1, state1); + + // test using Exp(2.0 (a^\dag_1 a^\dag_0 a_2 a_4 - h.c.)) + let term2 = [new JordanWignerInputState { Amplitude = (2.0, 0.0), FermionIndices = [1, 0, 2, 4] }]; + let state2 = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,-0.416147, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.909297, 0., 0., 0.]; + JordanWignerUCCTermTestHelper(5, [0, 1, 3], term2, state2); + + // test using Exp(2.0 (a^\dag_1 a^\dag_0 a_2 a_4 - h.c.)) + let term3 = [new JordanWignerInputState { Amplitude = (-2.0, 0.0), FermionIndices = [4, 2, 1, 0] }]; + let state3 = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,-0.416147, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.909297, 0., 0., 0.]; + JordanWignerUCCTermTestHelper(5, [0, 1, 3], term2, state2); +} + +@Config(Unrestricted) +// @Test() +operation JordanWignerUCCDTermPRQSTest() : Unit { + let term0 = [new JordanWignerInputState { Amplitude = (2.0, 0.0), FermionIndices = [2, 0, 4, 1] }]; + let state0 = [0., 0., 0., 0., 0.,-0.416147, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.909297, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]; + JordanWignerUCCTermTestHelper(5, [0, 2], term0, state0); + + let term1 = [new JordanWignerInputState { Amplitude = (2.0, 0.0), FermionIndices = [2, 0, 4, 1] }]; + let state1 = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,-0.416147, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,-0.909297, 0., 0., 0., 0., 0.]; + JordanWignerUCCTermTestHelper(5, [0, 2, 3], term1, state1); +} + +@Config(Unrestricted) +@Test() +operation JordanWignerUCCDTermPRSQTest() : Unit { + let term3 = [new JordanWignerInputState { Amplitude = (2.0, 0.0), FermionIndices = [0, 4, 2, 3] }]; + let state3 = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.909297, 0., 0., 0., 0.,-0.416147, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]; + JordanWignerUCCTermTestHelper(5, [0, 4], term3, state3); +} + + +@Config(Unrestricted) +function NearEqualityFactD(actual : Double, expected : Double) : Unit { + let tolerance = 1e-10; + let delta = actual - expected; + if (delta > tolerance or delta < -tolerance) { + fail $"Values were not equal within tolerance\nActual: {actual}, Expected: {expected}, Tolerance: {tolerance}"; + } +} diff --git a/library/chemistry/src/Utils.qs b/library/chemistry/src/Utils.qs new file mode 100644 index 0000000000..fedce81151 --- /dev/null +++ b/library/chemistry/src/Utils.qs @@ -0,0 +1,124 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. + +export + GeneratorIndex, + GeneratorSystem, + EvolutionGenerator, + HTermToGenIdx, + HTermsToGenSys, + HTermsToGenIdx; + +import Std.Math.AbsD; + +/// # Summary +/// Represents a single primitive term in the set of all dynamical generators, e.g. +/// Hermitian operators, for which there exists a map from that generator +/// to time-evolution by that generator, through an evolution set function. +/// +/// The first element +/// (Int[], Double[]) is indexes that single term -- For instance, the Pauli string +/// XXY with coefficient 0.5 would be indexed by ([1,1,2], [0.5]). Alternatively, +/// Hamiltonians parameterized by a continuous variable, such as X cos φ + Y sin φ, +/// might for instance be represented by ([], [φ]). The second +/// element indexes the subsystem on which the generator acts on. +/// +/// # Remarks +/// > [!WARNING] +/// > The interpretation of an `GeneratorIndex` is not defined except +/// > with reference to a particular set of generators. +/// +/// # Example +/// Using Pauli evolution set, the operator +/// $\pi X_2 X_5 Y_9$ is represented as: +/// ```qsharp +/// let index = new GeneratorIndex { +/// Term = ([1, 1, 2], [PI()]), +/// Subsystem = [2, 5, 9] +/// }; +/// ``` +struct GeneratorIndex { + Term : (Int[], Double[]), + Subsystem : Int[], +} + +/// # Summary +/// Represents a collection of `GeneratorIndex`es. +/// +/// We iterate over this +/// collection using a single-index integer, and the size of the +/// collection is assumed to be known. +struct GeneratorSystem { NumEntries : Int, EntryAt : (Int -> GeneratorIndex) } + +/// # Summary +/// Represents a dynamical generator as a set of simulatable gates and +/// an expansion in terms of that basis. +/// +/// Last parameter for number of terms. +struct EvolutionGenerator { EvolutionSet : GeneratorIndex -> (Double, Qubit[]) => Unit is Adj + Ctl, System : GeneratorSystem } + +/// # Summary +/// Converts a Hamiltonian term to a GeneratorIndex. +/// +/// # Input +/// ## term +/// Input data in `(Int[], Double[])` format. +/// ## termType +/// Additional information added to GeneratorIndex. +/// +/// # Output +/// A GeneratorIndex representing a Hamiltonian term represented by `term`, +/// together with additional information added by `termType`. +function HTermToGenIdx(term : (Int[], Double[]), termType : Int[]) : GeneratorIndex { + let (idxFermions, coeff) = term; + new GeneratorIndex { Term = (termType, coeff), Subsystem = idxFermions } +} + + +/// # Summary +/// Converts an index to a Hamiltonian term in `(Int[], Double[])[]` data format to a GeneratorIndex. +/// +/// # Input +/// ## data +/// Input data in `(Int[], Double[])[]` format. +/// ## termType +/// Additional information added to GeneratorIndex. +/// ## idx +/// Index to a term of the Hamiltonian +/// +/// # Output +/// A GeneratorIndex representing a Hamiltonian term represented by `data[idx]`, +/// together with additional information added by `termType`. +function HTermsToGenIdx(data : (Int[], Double[])[], termType : Int[], idx : Int) : GeneratorIndex { + return HTermToGenIdx(data[idx], termType); +} + + +/// # Summary +/// Converts a Hamiltonian in `(Int[], Double[])[]` data format to a GeneratorSystem. +/// +/// # Input +/// ## data +/// Input data in `(Int[], Double[])[]` format. +/// ## termType +/// Additional information added to GeneratorIndex. +/// +/// # Output +/// A GeneratorSystem representing a Hamiltonian represented by the input `data`. +function HTermsToGenSys(data : (Int[], Double[])[], termType : Int[]) : GeneratorSystem { + new GeneratorSystem { NumEntries = Length(data), EntryAt = HTermsToGenIdx(data, termType, _) } +} + + +/// # Summary +/// Checks whether a `Double` number is not approximately zero. +/// +/// # Input +/// ## number +/// Number to be checked +/// +/// # Output +/// Returns true if `number` has an absolute value greater than `1e-15`. +function IsNotZero(number : Double) : Bool { + AbsD(number) > 1e-15 +} diff --git a/samples/chemistry/SPSA/qsharp.json b/samples/chemistry/SPSA/qsharp.json new file mode 100644 index 0000000000..295f1ad88d --- /dev/null +++ b/samples/chemistry/SPSA/qsharp.json @@ -0,0 +1,7 @@ +{ + "dependencies": { + "Chemistry": { + "path": "../../../library/chemistry" + } + } +} diff --git a/samples/chemistry/SPSA/src/SPSA.qs b/samples/chemistry/SPSA/src/SPSA.qs new file mode 100644 index 0000000000..ee454914d9 --- /dev/null +++ b/samples/chemistry/SPSA/src/SPSA.qs @@ -0,0 +1,165 @@ +// Copyright (c) Microsoft Corporation. +// Licensed under the MIT License. + +export + SpsaOptions, + DefaultSpsaOptions, + FindMinimumWithSpsa; + +import Std.Random.DrawRandomInt; +import Std.Arrays.Zipped; +import Std.Arrays.DrawMany; +import Std.Arrays.Mapped; +import Std.Convert.IntAsDouble; + + +/// # Summary +/// Options for use with optimizing objectives via the simultaneous +/// perturbative stochastic approximation (SPSA) algorithm. +/// +/// # Named Items +/// ## StepScale +/// The coefficient by which steps along gradient vectors should be scaled. +/// ## StepPower +/// The power to which the iteration number should be raised when computing +/// how far to step along the gradient vector. +/// ## StepOffset +/// A number to be added to the number of iterations when computing +/// how far to step along the gradient vector. +/// ## SearchScale +/// The coefficient by which searches should be scaled when estimating +/// gradient vectors. +/// ## SearchPower +/// The power to which the iteration number should be raised when computing +/// how far to search in order to estimate gradient vectors. +/// ## NIterations +/// The number of iterations of SPSA to run before stopping. +/// ## MaximumSetback +/// Whether the maximum setback rule is enabled (requiring an additional +/// objective evaluation at each iteration), and if so, the maximum +/// allowed increase in objective values at each iteration. +struct SpsaOptions { + StepScale : Double, + StepPower : Double, + StepOffset : Int, + SearchScale : Double, + SearchPower : Double, + NIterations : Int, + MaximumSetback : (Bool, Double), +} + +/// # Summary +/// Returns a default set of options for use with SPSA optimization. +function DefaultSpsaOptions() : SpsaOptions { + new SpsaOptions { + SearchScale = 0.1, + SearchPower = 0.101, + StepScale = 1.0, + StepPower = 0.602, + StepOffset = 0, + MaximumSetback = (false, 0.1), + NIterations = 30, + } +} + +/// # Summary +/// Given an operation that evaluates an objective at a given point, +/// attempts to find the minimum value of the objective by using the +/// simulntaneous perturbative stochastic approximation (SPSA). +/// +/// # Input +/// ## oracle +/// An operation that evaluates the objective function at a given point. +/// ## startingPoint +/// An initial guess to be used in optimizing the objective function +/// provided. +/// ## options +/// Options used to control the optimization algorithm. +/// +/// # Output +/// The coordinates and final objective value found by the SPSA algorithm. +operation FindMinimumWithSpsa(oracle : (Double[] => Double), startingPoint : Double[], options : SpsaOptions) : (Double[], Double) { + let nParameters = Length(startingPoint); + // The SPSA algorithm relies on projecting gradients onto random vectors + // where each element is either +1 or −1. We can implement that in Q# + // by choosing an element out of [-1.0, +1.0] uniformly at random. + let drawDelta = () => [-1.0, 1.0][DrawRandomInt(0, 1)]; + + mutable currentPoint = startingPoint; + + // Depending on what options are enabled, we may reject certain + // updates, so we keep a counter as to how many iterations have been + // accepted. + mutable nAcceptedUpdates = 0; + mutable lastObjective = 0.0; + + // The SPSA algorithm proceeds by estimating the gradient of the + // objective, projected onto a random vector Δ of ±1 elements. At each + // iteration, the step size used to evaluate the gradient and the + // step taken along the estimated gradient decay to zero, + // such that the algorithm converges to a local optimum by follow + // a directed random walk that is biased by gradients of the objective. + for idxStep in 1..options.NIterations { + Message($"Iteration {idxStep}:"); + + // Following this strategy, we'll start by using the options + // passed into this operation to set αₖ, the amount that we look + // along Δ when using the midpoint formula to evaluate the gradient + // of the objective function 𝑜, and βₖ, the amount that we step + // along the gradient to find the next evaluation point. + let searchSize = options.SearchScale / IntAsDouble(1 + nAcceptedUpdates)^options.SearchPower; + let stepSize = options.StepScale / IntAsDouble(1 + nAcceptedUpdates + options.StepOffset)^options.StepPower; + + // We next draw Δ itself, then use it to find 𝑥ₖ + αₖ Δ and + // 𝑥ₖ − αₖ Δ. + let delta = DrawMany(drawDelta, nParameters, ()); + let search = Mapped(d -> searchSize * d, delta); + let fwd = Mapped((a, b) -> a + b, Zipped(currentPoint, search)); + let bwd = Mapped((a, b) -> a + b, Zipped(currentPoint, Mapped(d -> -d, search))); + + // We then evaluate 𝑜 at each of these two points to find the + // negative gradient 𝑔ₖ = 𝑜(𝑥ₖ − αₖ Δ) − 𝑜(𝑥ₖ + αₖ Δ). + let valueAtForward = oracle(fwd); + let valueAtBackward = oracle(bwd); + let negGradient = (oracle(bwd) - oracle(fwd)) / (2.0 * searchSize); + Message($" obj({fwd}) = {valueAtForward}"); + Message($" obj({bwd}) = {valueAtBackward}"); + + // We can step along 𝑔ₖ to find 𝑥ₖ₊₁. Depending on whether options + // such as the maximum setback rule are enabled, we may reject + // the update. Either way, we report out to the caller at this + // point. + let step = Mapped(d -> negGradient * stepSize * d, delta); + let proposal = Mapped((a, b) -> a + b, Zipped(step, currentPoint)); + if Fst(options.MaximumSetback) { + // Is this our first update? If so, accept and set the + // lastObjective. + if nAcceptedUpdates == 0 { + Message($" First update; accepting."); + lastObjective = oracle(proposal); + nAcceptedUpdates += 1; + currentPoint = proposal; + } else { + // How much did our objective get worse (increase) by? + let thisObjective = oracle(proposal); + if thisObjective - lastObjective <= Snd(options.MaximumSetback) { + Message($" Proposed update gave objective of {thisObjective}, which is within maximum allowable setback of previous objective {lastObjective}; accepting."); + // Within the limit, so we're good. + lastObjective = thisObjective; + nAcceptedUpdates += 1; + currentPoint = proposal; + } else { + Message($" Proposed update gave objective of {thisObjective}, which exceeds maximum allowable setback from previous objective {lastObjective}; rejecting."); + } + } + } else { + // No maximum setback rule, so always accept the proposed + // update. + nAcceptedUpdates += 1; + currentPoint = proposal; + } + + } + + return (currentPoint, oracle(currentPoint)); +} diff --git a/samples/chemistry/Variational Quantum Algorithms.ipynb b/samples/chemistry/Variational Quantum Algorithms.ipynb new file mode 100644 index 0000000000..bc6235cceb --- /dev/null +++ b/samples/chemistry/Variational Quantum Algorithms.ipynb @@ -0,0 +1,2083 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a56dec3d", + "metadata": {}, + "source": [ + "# Variational Quantum Algorithms in Q\\#" + ] + }, + { + "cell_type": "markdown", + "id": "96818c17", + "metadata": {}, + "source": [ + "## Abstract" + ] + }, + { + "cell_type": "markdown", + "id": "c55ddd55", + "metadata": {}, + "source": [ + "In this sample, we will explore how the rich classical control provided by Q# can be used to efficiently write out variational quantum algorithms. In particular, we'll focus on the example of the _variational quantum eigensolver_, also known as VQE." + ] + }, + { + "cell_type": "markdown", + "id": "4fc6c47b", + "metadata": {}, + "source": [ + "## Preamble\n", + "\n", + "$\n", + " \\renewcommand{\\ket}[1]{\\left|#1\\right\\rangle}\n", + "$" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "11fa63fc", + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": "// Copyright (c) Microsoft Corporation.\n// Licensed under the MIT License.\n\n// This file provides CodeMirror syntax highlighting for Q# magic cells\n// in classic Jupyter Notebooks. It does nothing in other (Jupyter Notebook 7,\n// VS Code, Azure Notebooks, etc.) environments.\n\n// Detect the prerequisites and do nothing if they don't exist.\nif (window.require && window.CodeMirror && window.Jupyter) {\n // The simple mode plugin for CodeMirror is not loaded by default, so require it.\n window.require([\"codemirror/addon/mode/simple\"], function defineMode() {\n let rules = [\n {\n token: \"comment\",\n regex: /(\\/\\/).*/,\n beginWord: false,\n },\n {\n token: \"string\",\n regex: String.raw`^\\\"(?:[^\\\"\\\\]|\\\\[\\s\\S])*(?:\\\"|$)`,\n beginWord: false,\n },\n {\n token: \"keyword\",\n regex: String.raw`(namespace|open|as|operation|function|body|adjoint|newtype|controlled|internal)\\b`,\n beginWord: true,\n },\n {\n token: \"keyword\",\n regex: String.raw`(if|elif|else|repeat|until|fixup|for|in|return|fail|within|apply)\\b`,\n beginWord: true,\n },\n {\n token: \"keyword\",\n regex: String.raw`(Adjoint|Controlled|Adj|Ctl|is|self|auto|distribute|invert|intrinsic)\\b`,\n beginWord: true,\n },\n {\n token: \"keyword\",\n regex: String.raw`(let|set|use|borrow|mutable)\\b`,\n beginWord: true,\n },\n {\n token: \"operatorKeyword\",\n regex: String.raw`(not|and|or)\\b|(w/)`,\n beginWord: true,\n },\n {\n token: \"operatorKeyword\",\n regex: String.raw`(=)|(!)|(<)|(>)|(\\+)|(-)|(\\*)|(/)|(\\^)|(%)|(\\|)|(&&&)|(~~~)|(\\.\\.\\.)|(\\.\\.)|(\\?)`,\n beginWord: false,\n },\n {\n token: \"meta\",\n regex: String.raw`(Int|BigInt|Double|Bool|Qubit|Pauli|Result|Range|String|Unit)\\b`,\n beginWord: true,\n },\n {\n token: \"atom\",\n regex: String.raw`(true|false|Pauli(I|X|Y|Z)|One|Zero)\\b`,\n beginWord: true,\n },\n ];\n let simpleRules = [];\n for (let rule of rules) {\n simpleRules.push({\n token: rule.token,\n regex: new RegExp(rule.regex, \"g\"),\n sol: rule.beginWord,\n });\n if (rule.beginWord) {\n // Need an additional rule due to the fact that CodeMirror simple mode doesn't work with ^ token\n simpleRules.push({\n token: rule.token,\n regex: new RegExp(String.raw`\\W` + rule.regex, \"g\"),\n sol: false,\n });\n }\n }\n\n // Register the mode defined above with CodeMirror\n window.CodeMirror.defineSimpleMode(\"qsharp\", { start: simpleRules });\n window.CodeMirror.defineMIME(\"text/x-qsharp\", \"qsharp\");\n\n // Tell Jupyter to associate %%qsharp magic cells with the qsharp mode\n window.Jupyter.CodeCell.options_default.highlight_modes[\"qsharp\"] = {\n reg: [/^%%qsharp/],\n };\n\n // Force re-highlighting of all cells the first time this code runs\n for (const cell of window.Jupyter.notebook.get_cells()) {\n cell.auto_highlight();\n }\n });\n}\n", + "text/plain": [] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/x.qsharp-config": "{\"targetProfile\":\"unrestricted\",\"languageFeatures\":null,\"manifest\":\"{\\n \\\"dependencies\\\": {\\n \\\"Chemistry\\\": {\\n \\\"path\\\": \\\"../../../library/chemistry\\\"\\n }\\n }\\n}\\n\",\"projectRoot\":\"file:///Users/swernli/Programming/qsharp/samples/chemistry/SPSA\"}", + "text/plain": [ + "Q# initialized with configuration: {'targetProfile': 'unrestricted', 'languageFeatures': None, 'manifest': '{\\n \"dependencies\": {\\n \"Chemistry\": {\\n \"path\": \"../../../library/chemistry\"\\n }\\n }\\n}\\n', 'projectRoot': 'file:///Users/swernli/Programming/qsharp/samples/chemistry/SPSA'}" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from itertools import product\n", + "import qutip as qt\n", + "import qsharp\n", + "\n", + "qsharp.init(project_root=\"SPSA\")" + ] + }, + { + "cell_type": "markdown", + "id": "9a80f41f", + "metadata": {}, + "source": [ + "For Q# code embedded in this notebook, it will be helpful to open a few namespaces before we proceed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "373cb753", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "import Std.Arrays.*;\n", + "import Std.Convert.*;\n", + "import Std.Diagnostics.*;\n", + "import Std.Random.*;\n", + "import Std.Math.*;" + ] + }, + { + "cell_type": "markdown", + "id": "dce084a9", + "metadata": {}, + "source": [ + "## Introducing the Variational Quantum Eigensolver" + ] + }, + { + "cell_type": "markdown", + "id": "983fbce7", + "metadata": {}, + "source": [ + "In variational quantum algorithms, rather than performing all of our computation on the quantum device itself, we use a quantum program to estimate some quantity, then use a classical optimizer to find inputs to our quantum program that minimize or maximize that quantity. That is, rather than thinking of classical computation purely as a pre-processing or post-processing step, our computation uses both classical and quantum computation together." + ] + }, + { + "cell_type": "markdown", + "id": "150b8540", + "metadata": {}, + "source": [ + "> **💡 TIP:** We could also consider using classical computation while qubits are still alive, rather than returning an estimate to a classical optimizer. This approach is indeed very useful, for instance in iterative phase estimation. For this sample, however, we'll focus on the variational case." + ] + }, + { + "cell_type": "markdown", + "id": "8f8c35bf", + "metadata": {}, + "source": [ + "\n", + "For example, consider the problem of finding the _ground state energy_ of a given operator $H$, often called a _Hamiltonian_. The ground state energy of $H$ is defined as its smallest eigenvalue $E_0$, as the Hamiltonian represents the possible energy configurations of a given physical system.\n", + "\n", + "As stated, even though this is a minimization, we need to know the eigenvalues of $H$ to make any progress. It's pretty straightforward to rephrase the problem in terms of arbitrary states, however. In particular, the expectation value $\\left\\langle H \\right\\rangle H = \\left\\langle \\psi | H | \\psi \\right\\rangle$ must be at least as large as $E_0$. To see this, we can expand the expectation value in terms of the eigenvectors $\\{\\ket{\\phi_i}\\}$ of $H$, such that $H\\ket{\\phi_i} = E_i$.\n", + "\n", + "Since the decomposition of $H$ into eigenvectors gives us a basis, we know that $\\left\\langle \\phi_i | \\phi_j \\right\\rangle$ is zero whenever $i \\ne j$, allowing us to expand $\\ket{\\psi}$ into the eigenbasis of $H$ as $\\ket{\\psi} = \\sum_i \\alpha_i \\ket{\\phi_i}$ for some complex coefficients $\\{\\alpha_i\\}$. Using this decomposition, we can expand the expectation $\\left\\langle \\psi | H | \\psi \\right\\rangle$ in terms of the eigenvalues of $H$:\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + " \\left\\langle \\psi | H | \\psi \\right\\rangle\n", + " & = \\sum_i |\\alpha_i|^2 \\left\\langle \\phi_i | H | \\phi_i \\right\\rangle \\\\\n", + " & = \\sum_i |\\alpha_i|^2 \\left\\langle \\phi_i | E_i \\phi_i \\right\\rangle \\\\\n", + " & = \\sum_i |\\alpha_i|^2 E_i.\n", + "\\end{aligned}\n", + "$$\n", + "\n", + "Since each $|\\alpha_i|^2$ is nonnegative, and since $\\sum_i |\\alpha_i|^2 = 1$, the above gives us that\n", + "$$\n", + " E_0 = \\min_i E_i \\le \\left\\langle \\psi | H | \\psi \\right\\rangle \\le \\max_i E_i.\n", + "$$\n", + "\n", + "Thus, we can rephrase the original problem as a minimization not just over eigenstates, but of all arbitrary states,\n", + "$$\n", + " E_0 = \\min_{\\ket{\\psi}} \\left\\langle \\psi | H | \\psi \\right\\rangle,\n", + "$$\n", + "where the minimum is achieved when $\\ket{\\psi} = \\ket{\\phi_0}$.\n", + "\n", + "Using this rephrasing, the _variational quantum eigensolver_ algorithm is just the variational algorithm that we get by using a classical optimizer to find $\\min_{\\ket{\\psi}} \\left\\langle \\psi | H | \\psi \\right\\rangle$. In pseudocode:\n", + "\n", + "```\n", + "operation FindMinimumEigenvalue {\n", + " Pick an initial guess |ψ⟩.\n", + " until target accuracy reached {\n", + " Estimate E = ⟨ψ | H | ψ⟩ using a quantum operation.\n", + " Use a classical optimizer to pick the next |ψ⟩.\n", + " }\n", + " Return the minimum E that we found and the state |ψ⟩ that achieved that minimum.\n", + "}\n", + "```" + ] + }, + { + "attachments": { + "image.png": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAocAAAExCAIAAACBDZ9AAAAAAXNSR0IArs4c6QAAAERlWElmTU0AKgAAAAgAAYdpAAQAAAABAAAAGgAAAAAAA6ABAAMAAAABAAEAAKACAAQAAAABAAACh6ADAAQAAAABAAABMQAAAACLM+G+AAA5IElEQVR4Ae2dC9RcVZXnkwGXNjrqaFpEbEVISIBAQIG15KFgq0RBeYgKCBiFVhMhPERACYoY1CACRkwUoYmICQoSVFrjY+Sta3hIAgQCBFQUUZue0Z5ubZc63/yr9lf7O7n31v3qcd/1q/Wtyqlz9tlnn9+p3P/d596qmjo2NjaFBwQgAAEIQAACFSDw3yoQAyFAAAIQgAAEINAigCrzPoAABCAAAQhUhQCqXJWVIA4IQAACEIAAqsx7AAIQgAAEIFAVAqhyVVaCOCAAAQhAAAKoMu8BCEAAAhCAQFUIoMpVWQnigAAEIAABCKDKvAcgAAEIQAACVSGAKldlJYgDAhCAAAQggCrzHoAABCAAAQhUhQCqXJWVIA4IQAACEIAAqsx7AAIQgAAEIFAVAqhyVVaCOCAAAQhAAAKoMu8BCEAAAhCAQFUIoMpVWQnigAAEIAABCKDKvAcgAAEIQAACVSGAKldlJYgDAhCAAAQggCrzHoAABCAAAQhUhQCqXJWVIA4IQAACEIAAqsx7AAIQgAAEIFAVAqhyVVaCOCAAAQhAAAKoMu8BCEAAAhCAQFUIoMpVWQnigAAEIAABCKDKvAcgAAEIQAACVSGAKldlJYgDAhCAAAQggCrzHoAABCAAAQhUhQCqXJWVIA4IQAACEIAAqsx7AAIQgAAEIFAVAqhyVVaCOCAAAQhAAAKoMu8BCEAAAhCAQFUIoMpVWQnigAAEIAABCKDKvAcgAAEIQAACVSGAKldlJYgDAhCAAAQgsDkIRorA0Zc9MFLzZbIQ6J3AVcfv2LsxlhDIiQCqnBPY6rp9667TqhsckUGgJALXrH2qpJEZFgKbEECVN8ExEi/GRmKWTBICEIBAHQlwXbmOq0bMEIAABCDQTALkys1c15RZkSqnwKEJAhCAQLkEUOVy+ZcyOrpcCnYGhQAEIDA5AVR5ckYNsxhDlBu2okwHAhBoEAGuKzdoMZkKBCAAAQjUnACqXPMFJHwIQAACEGgQAXawG7SYvU2FHezeOGEFAQhAoAQC5MolQGdICEAAAhCAQCIBVDkRC5UQgAAEIACBEgiwg10C9HKHZAe7XP6MDgEIQCCFALlyChyaIAABCEAAAoUSIFcuFHcVBhubwgeWq7AOxAABCEAggQC5cgIUqiAAAQhAAAKlECBXLgV7qYOSKpeKn8EhAAEIpBBAlVPgNLMJUW7mujIrCECgEQTYwW7EMjIJCEAAAhBoBAFy5UYsYz+TIFfuhxa2EIAABAolgCoXirsSgyHLlVgGgoAABCCQQIAd7AQoVEEAAhCAAARKIUCuXAr2Mgfl88pl0mdsCEAAAqkEyJVT8dAIAQhAAAIQKJAAuXKBsKsxFN+DXY11IAoIQAACCQTIlROgUAUBCEAAAhAohQCqXAp2BoUABCAAAQgkEGAHOwFKs6vYwW72+jI7CECg1gTIlWu9fAQPAQhAAAKNIkCu3KjlZDKDEfj2ymXr777t3jtu8u5Hzl/0pqMW+MshC5/6wFFyvsue+535mZUDu8rEycCj0xECECiGAKpcDOcKjcIOdrgYD917x/VXXhzqsbWuWr5YOn3GBf2J6JLTWuorD1+95dc+ioZw/wPDz8SJh0QBAhCoLAFUubJLk19gfOXmOFtJ3bknHGovjnzfWQe1k+O2Tn9WOqq/G1Z+3ir7WgzlxFOmhJDHy4cce9Km9X15zcRJXyNiDAEIlEAAVS4BerlDhnJRbiSlj379lZ+1GM6+ZPXMXfY0MtvvsufBx55k2e2qL5x3YD/72BM5cTA3ObzqliesYmD4mTgJgqIIAQhUlACqXNGFIay8CdywcpmJqFJbSXI4XOSlmpRAf7ydVUtf1fHqL5ynSnWUfpvx+ae9wyVZhaNftbUMTIzd3l6mu3Lj0Llceb05Mf+qDx/WZDVhPEe876wZs3f3SVmTKmVpE7GTktAVZQhAoCwCqHJZ5Msbd+B8rbyQ8xj5gZ/ebm5P//RXN9lvjgzWxvXI/XdZdah2Ul/9XXVzKw92SfberX3sdl9TvpYKbko+4so6up/QuZq6OYkP99B94ycQ3jQhvTu3Tj5sCE3fx5qp+k1j874UIACBggnwyaiCgZc/nA6//ImAadLOe+4XpyFhs3XyVl829VJm+ZWbn1CTVd6wapk8hDUq6++Dn/6q6t2VjG0gF/iIK72M1KiLuic6sSHCQZW1m6Xl9ArPbBStxalxLQB7qbHcxupH/Nmw8AyB0gmgyqUvAQGUQOBfVi1LGdWF020s3dRLidz27YyzfeuWt08UXK0nqtol7SFbzfpOju6u3NJrdnr53l4ZFtyJVZ7/wXfc177l2zv6lfLWBkD74TX20ieuON3GmniGAASqQIAd7CqsQqExKCXi4RB2fPneXnYsLpyR1rfrAu3O4zeFea/ps3e3sgmknHiTyg93tr5VGTFzV97RazwAq4k7kWfpq3VcdMlq7+iujnl168K2P6TBbzxygcegekvl3YACBCBQEQKockUWosAwQtEocNhqDeUQQqXqhGjatvMe+x14xAJp6cOdDe0ZO+3ukvvIfeNXmls1gc2Ou+3tNvL3wN3jV6+3n926duuu3v7e6GXmsMYDMFcRJ3IrP19r33GmIM2zVXZmMP6vWlXSucWBR7Ymooe5atU7gVY1DwhAoCoEUOWqrERhcXA0FmoluAa8LambsL/gg++w1zt00ujEVPVrXzzPzCLprKfO1nrfnTepIBU07O7KzVyn4zUeQNTJfXecd2LrY9Zye1r76rWN5St71ufGt9mtXs/eZK7CGrehAAEIVIEA15WrsArEUBqBr3/xPNdFBSFJdgls5ZebPr7V+XyzK7f0z0we7FwtDnu4Z+lrWK+yXZxWIX4N2y3tKnLEiV6GkuzGYcHjVKU2uj1at4nH400UIACBcgmQK5fLv4TRPW0qYezKDKkE98OfW/2JdsZpIheGphz0A0EOKuW2Vgn2O/ebuF4rM7+gK50zOTdvci7dddSeB7srb/Jx3ZXrtGzsz2zMibeGwbztvWfpsrE8dIvThvtOcI9bPACPhAIEIFAiAVS5RPglDc3xuA1el2M/vHT1xvV3uVJKZdXypmNOaiWyMUpSPiXEJr0yU9/Q7I1HLPBW+bFrvRs7t3ptoq5TpsiV+7fRW0N3RvS0O9GJt7YnMf40vXPBO3FGHzg/+oFsRevDhX4oQwACpROYOjbw9+WXHjsB9E/g6MseeO305/bfb3R7PHLfHZ9Y2LqIK8HzZHR0cTR35j/c+Purjt+xufNjZrUhwHXl2iwVgZZC4JH1nXutSxmeQSEAgREjwA72iC14cDvuyM18uAlPD64TD+eJ3hCAAAS6EiBX7oqGBgiIwDXtW71mty85AwQCEIBA3gTIlfMmXDn/3EjQ15Jc/qMnzB5ufXHDGAIQGIwAufJg3OgFAQhAAAIQyJ4AuXL2TCvusfMBnIqHSXgQgAAERpEAufIorjpzhgAEIACBahIgV67muuQYFddHc4SLawhAAALDESBXHo4fvSEAAQhAAALZESBXzo5lTTyN8V2LNVkpwoQABEaQALnyCC46U4YABCAAgYoSQJUrujCEBQEIQAACI0iAHeyRW3Tu9hq5JWfCEIBAfQiQK9dnrYgUAhCAAASaToBcuekrHJsf3yISQ0IFBCAAgaoQIFeuykoQBwQgAAEIQIBceeTeA1xXHrklZ8IQgEB9CJAr12etiBQCEIAABJpOgFy56Sscmx/XlWNIqIAABCBQFQLkylVZCeKAAAQgAAEIkCuP3HuAXHnklpwJQwAC9SGAKtdnrbKKFFnOiiR+IAABCGRNAFXOmmjl/SHKlV8iAoQABEaXANeVR3ftmTkEIAABCFSNALly1VYk93jIlXNHzAAQgAAEBiWAKg9Krr79kOX6rh2RQwACTSeAKjd9hWPzQ5RjSKiAAAQgUBUCqHJVVqKwOP7X438obCwGggAEIACBvghMHeNrkfsChjEEIAABCEAgNwLcg50bWhxDAAIQgAAE+iSAKvcJDHMIQAACEIBAbgRQ5dzQ4hgCEIAABCDQJwFUuU9gmEMAAhCAAARyI4Aq54YWxxCAAAQgAIE+CaDKfQLDHAIQgAAEIJAbAVQ5N7Q4hgAEIAABCPRJAFXuExjmEIAABAYlcNTpH9h/3jGD9qZfCQROXfLJWQfN1Z8KxQzPd3sVw5lR+iDw+//4i6yf+6yn9dEH06YQ+I8//fVZf1f749KXv7n6k1/6oq/Jh/7pve88+FB/SaEuBLSO37n15qUfWvT6vfcpLObav/sLI8VAxRD49//8y0/W/2+Ntffs5z37mQhzMdSzGUWCes4VG9Y9OvhXus7Z7jnnvGtWNtGU50UJ8U8fWP+Niz+30/QZFoUyLRUQ5gzX5Pu337bwk4vzPt35zVNPKeYiJVnDocoZvk9wNSwBHdZ/ePe/Lv/mz+ToT39+2ev3eEED0qZhodC/VgSWXP4lSbKyK5dkhb/hhjW1mgTBjhP4bVuVC8aBKhcMnOG6Evjjf/3tm7c9uWLN42Yhbf7Tn/926L4v2uIZm3XtQwMEKkZgza23vHzHndKzq/UbH3nLySda4JHdUUsBfU69t4Y+33XoW8447p/cia5kP/nUv9rLMIN3AxV00fQ3//bUyvM/Y2m9at6476svPONDbiMD7eXay7BJZyFXrP6G3NqMbOjILPykxIJUgrtuwwbzttW0v79xxVfC4CMRmn8b1+dluxGq1GUC/ZkTs/H49dLHVVn16i6VtXHVFA7qns2JnsNW82mBhfUyC5N1IbrnwQdUoyQ+0uRueylwt1cvlLDJnYA2rlf+8JcuyTaeXqpSTbkPzwAQyIKADtnSvzkz0zbhZSABkyroT/KmI7h62eC6imm7stYqqdBLKVwvrfIpe+u47qEN7lOKstULXuAOZeZNkRkrxZexzgNkrGeplxTRbKQ3W06bZk4iTWYgtxItGehsQAFLKc1Yz5LMyA1uajVv6iIaGtSB6IRGZQ9M45rky4+MVbaQdPagMGQmCVSTdF1lzUuuhNSGVkEv3ZUK6q5na1VBAyUSsy7a6pClnHgX1WhqYS+Nrrlo1ayLnjUdLZkNMfAFC1TZeVIojcAf//y3L3/vl1+78Yl4BKpUkwziTdRAoGoEnvjtbxXSC6dNSw9MR20zOO4th6tw1/r77aUO8ZIBP5pL4SRpKzoH/ZRWE9rdZu1gfiRatn/uGmb15vDyb1xrL+PPkjrL8vUsgVTebzZKmj35VpOikvCH3SVvvmMvA5NJM5i776ukVeGpgOZo3tRFo8jMgRzQvqnKTkTURWcGUj7zrGeNYsoaDu1lzUuBeX5vBT+xMDNv7UbMvSUWtAQK2FFopfRyxerrQmMFHL4coIwqDwCNLlkS+Ovfxi657rFv3f5kN6dqkoHMuhlQD4EaEZByeLSmN3ZLkUnRnFmb5Nm77bCjUljZp7fKj9yGibUNIe002fMRlTdrp9pfRgomyVapjD+ipm4sJ162wkGv3i9S4y/j5yhKlCdanz8tBGJljSsDO1nZfafZE8btjqHAe5MK2j0WrrBG3sILw5b4mkE3YmH3SFnjKrDIRkickp9URbr3/jLL68q+12/D++lP79FgOYIEFn/lodvv+7f0if/grt8pXT5n3iYHrPQutEKgjgRCiVL8oYDpZUqr0lNtFEuYZaac0vK5J3/3OwlJZCM3otM9Uooc3tOdRD4Y1uMQETM7WQk3tCMGkZea6ZO33myXjb1pNy/FConEYlbRishJRuRl1Hqg19mosk4ixE7r5EqsVdFbwV8OFFvGnexuhUqFlPEM6+ZO6e+iyx64++Hf9xK4lPvML65ffPyOm282tRd7bCBQPAHLNXUr05SDBxzc0kTvHKZ6qkxvtX1juxYrtbCkTYdlbWi7w34Lls3bhWE/eEqhU/xo01j7zL4ZPqRC2z1WKcOFTcqGfY86rO9WTiTWzdjq7VzBbSIvvX6YQjY72CcsPlcnceHa6w3hSzhMfPRtKgHlvqd+/r4eJdkgyFhduMbc1LdEM+YlYVC61m2XNWWOE4oeGGlX1rLS9NagxxSTJVML7egqXQ5b08u2T2422v221Nx2bnV5OL2vt1rHcDPcm/oq2GVyu1Qf77j1lltGKgUqZXM+Yhy+DImF9ZGybXpHLqjbZO3cJWI/8MsMVFkLqTO4eYcelhKETt+UOttfeD+e1luVOplyA2u1erMP399mbIm4tYZvI9WE1/bNTAbmzXY2rJebha4UQ49TUC+3VC/FHAYctsosbEqcu048FZU32UuL057lUK2R81Mj5mHUq6Dbqk9aeu+Dv/i//YatLurIXdn9csO+MAI6xEvMtHcYHrj0/zdyWEiMRzvPOkz5MU3/x3VoPes97zPjlFY598OXDWQ7q7qbTB68SX50GEmJRBvgFrZsdD3bjuomOa5G8maXuhOnoMoXPn+aBjU/movukOpmmV4vXZfQ2p68WSoqPwyOR6Vtic5j3sGHKjA/tqta2B1mx2r8X7lyLAakl71oAdEQZi9HGksvh7+9KxJbBjvY92x4UE7Da/KRMcRRoXvqLFL6s60Ds9Sy6Q0nAy2k3s2SItWbvfqqxvuqXsZ601uNsGrNJt3i0PrJXsZ6x4euxDTcabFxEzdAwilome2N4lf19RZUkBaGfFqEdqpoxn6ZR34S5+7bDDLQua3PTgGbW535KlTx8ZMyNYU3L0SYV/nlE0/91xlfuP+3/+fPgwX589/8ccFF65a8b/bW054xmAd6QSBXAjq46WijY4KPogO3Hy68Ml6wi8GhDoXHq5RWOT91QyvzMZ8+nB36VK/DRaQpPrpqdLTxsN2J15t/HXb0l5KV6hCqVvejrexwRonjdqvUgVEkfV6R3XhFqIOtpmb1OuRa/H6ftu+ix/13Ixa3DGtsETWo/qzejs+hzfDlqWNjw97aGle7MCydVmgCIR0TKltyk2GtsWthqH/yE+mu5ZEku6Jbd9c8tXo53jceZ6J9nLIFHE5BrrSz5NckXDtt4qFbTUeVLroWcDj3+HT8P0Pc2GcXD8mGrv7zY0/+8YPL7/v3//zrkKE++5mbf3r+zttutcWQfuieIYGsvnGTL3TLcFF6dBU/PPbYEbPMCWSwg50eU+vGh02/R9SSSKu3vuF9htr9kFC5TysrGfWa8N53nQnKwLdW3KaXgoRNZv4JP5VtrPhlDNsMsLDNswIOQ1Klp7Dmx+/R0CZBeCe9mYU3CPR4tUYddT7oHx9cc9utijYMyQKr+LP2n0+4eN3wkqxpyolcDbAHXnFEhAcBCIw4gQx2sCclaGoXmsVrwtZiyiarvWytmMQqA+43MCW76qLtFN9RSfdg0vu922+zfRL7sL+V1VGfr9eug3zKrHbb1/oZKO08f3D5+LclpHPosfUvf/1/C5fe++n5s7d54Rb8wFSP0DCDAAQqTiADVdbH3iUSSj27pW6RtFJEVJPyMbJ+kSm97reL24f70l6ZWAgv8CQadKv0beduBmG9cmu99DMAbad7q+RZqnzDzTdZNj93n329qeIF3Zyln4G68Osb84hTSn/q26bzA1N5sMXn6BDQNUS/jDg6s67mTDPYwbZkzr8WLjJP+6oa2y62JitHvsIm0ivlpS7oeqsSRwm8b4ArBfetY9mEm+R66WbW3W5Pi58xuHMvWKiW+HplLwUltZGQ0nsZGSmxzgDsTx7CLroAr03sem1f289A5STJBkfO9UtTGihkRRkCEIBAHQlkoMqatjJOJXl2Z5NRkMBYwifNljL5HWtq1aaxLpH6xmy/1KSjPtB5l35B3f1bSXXJWVm7yafuhVY5dG43vvv5gQRPImd7wmamJt0gHXaxsk1Bn8n2JjnXzRH+MqWgO+kVht9JL0th6Sbw2mwQK927KBv/C411fqDpy2GPV6NTAiumyX4Gyn6ZMdcRNYR+b0rD5ToKziEAAQjkTSCDHWyFKDlRbmdC4hH7lq/uVZbaqdWawjuu3bj3gror63VvPoo8hHfky0znCuFlY4mrLtlajd3nLHu58jv4pYh+d3ckHtXrVMAH7X0Kpuga1M9L4vd4+1gSYImu34Otemm/wvMu8mZ+Ur511r2VXtDG9ddvfCLxNyfyiE0/MKVffnzb/ls/+5lPy8M/PiEAAQgUQCCDT0YVEKUPIV3sXRG9V10K4QeuLObIB8NUqfMbfTW8f9SqslPTN3Bd/i+/SPnNiZwif/PeWx134Eu3eDo/yZwT4DS3fDIqjU6sLfJ9CbH2iQqdr+vsPDxfn2ij1DgC2exgNw5LORNS4q5c2ffYFYR+I0wZvN9GpyYZ2I+dlRNib6NO+jNQvbkZxIofmBqEGn3KIGAf6fT/3SkhxH89KcWYproTyGYHu+4UKhJ//Lt7It9lY7fUaR+7IgF3C0M/IHH6kTP0183gtGX3r3v0D91aJ62fs91zLlgw8ftuk9pjAIEiCfSY2va+42U3rkbu/UyZkU7fddWs9w+Y2J6cXylL8UxTAQRqpsrhVeQC6BQ/hITZb16Lj977f+N43+rUZHKztJzwDVDVWVMiCQlkntrqCyx1gh4OkV6Of/FRur1UX3tyvat+ujdahyRQM1UecrZ0hwAEIJAfAcuSzb/dRmoXg3XLiJRVP59gd5sqK5WNDMIvM5CNfWxEAqlesvTkVZ9w0f00iWHrRhNd1VKT3XBjWbJZ2t2pnjH7zaoaVJ8gVTw6y7csObT3QXV/q319glq5pG2IinlGlYvhzCgQgEDzCSjd1H6efXQz3NmSBOpXZ/QBCt/ts9tH/Bt/TQKtVXeBmXhb8iqlF7jEL3iwez/tkyPyIJ/2cRjV62Oi/q0gdq7g98makJvM63KY/iTY4fmBhlONsnOPR5HrCx5Ipot5B3O3VzGcGQUCEBgVAkoxw2+/17RVo4w2/OBluMmshFUGSmoNkK5hSRF9y7rbfrg0OLz3UycBduOYNFj1oYrrex2Uf7tIX7LoIxrIDSLnB2qSwMvezyrsmx7iPxBg0fKcOQFUOXOkOIQABEaXgIlc+GO9VqNkNISiG7ClfFajj1pIg8ObsUNd73ar19Zbbqnu/kUI7txU3J1rdHmzX0o2m4jM2/mBeZOBRF32/j1Feqkh5C0Mz8eikAcBVDkPqviEAARGlICJXPh783bdN/LNP9rQtp+/G09tZ85yXqpR2XW9261e2k+29Fq7zdqU9u6m4i6i8XjsN+t8O9rOD/ylvmlfrvSDOnKrP7v4HWb5PhCFnAhwXTknsLiFAARGkYBupFJm6SInBPFk15R47rTWz+rYzrBrsGpMF13Xlbl2u9VL0qtLv3bHlm4Wsz3qiIrbTwOE8ei79H17XMP5+YHKepi9X/+2Sp6LJECuPEH7tttumzt37tSpU/Wsshq8MGGUQ6mscXOYCi4hMOoE9PM5+va9kEJEJtVkSuy3eoXGEmz74VfTUcub/RpwaOll3a7l+9WqlIqn/IyebiVT7u4Gdn4Q+eUe90yhFAKo8gT2xYsXL1q0aGxs7DWvec2+++4redbLffbZZ8Iin1JZ4+YzG7xCAAKbEJBMRm7+Cm/1sq1my6fVTTvGkljPZSPXgN2vxFV/9lK5soS220+7mqLbtW1ZKlFWrxSZtyZZmnPJdrg97gFQyI8AqjzBds2aNabBp59+urRZjwIkWcOXNe7EzClBAAIZEdB9VZJhXZG1z0fFb/7SOOGtXnqpTwPrk8p2Hde2jl3F47vfFqbu05YT66K7sfQhY7+QrNvKzJv2tGWsTFoar49ayVjeNJYqfXtcGbm2x+0qssm87OVBPs25zhKsi43LcwEEKvHrFMpKE6e6ZMkSCWRiU36V559//hlnnCFJjgyRWG+VZnnrrbeaivt0DjjgACmutaoyfTqJ/iMxNOMlv2HQjHWMz4KVjTOhBgL9EqhEriw9U9wSrXaC2nqymuIlWWH86Ec/kpr2yFERmrFLsjq6orsk9+Ktr3F7cYgNBCAAAQjUjkAlVDlOLe+tYyWmdj9XfOjvfe97uq4cr+9WI3s1hQGbZ51kdOuSWJ84rlwp1ER7KiEAAQhAoHkEKvHJqB//+MciG8mMPePMA3riHrUPtNdee3m5l0JEgG06/TrRQPEuEnvdd6b6UPV7CQkbCEAAAhCoI4FK5MoFb94q+zQd9Y8ked5sialJoMq6GOxNiatrrRE1leTLOEVH+xpXoZrMJwZAJQQgMFIEzjnnnGrO96abbqpmYLWLqhKqbJvAxs4UK1eOUk3TUX0kyS79qmAj9nt+YHppH6OShNtj0uD7GlehmsxP6hYDCJRLIJPf1szESbkc8h69msIsVd5///3znvso+K/EDrZAS5glaUY8siEcXwa3jDe5h8h+uFtadmuJrN+NlXghOZIBu4ewYCruftQk/xLp9Cm4fS/jWqhym5J8hyFRhkB+BP76t7ELv77xB3f9Lqch1j36h9d94PZuzl+3+wtOfdv0zTdL/shGt141rZf03nzzzYnBS//2az8SW/OuTAlMQ0uYb7zxxrxjaLb/8lXZNo0lY6ajUtxuguorMcwl54hkmki7Z50cRAQ1XQvj9j1eVO5rXIWksIeZtU+QAgSGISBFPOGwbf/u6Zt96/Ynh/EzQN83773VcQe+tO6SrP/42irTcUOf3bAvKdIXCCZ+W1G3hNjqJcoDMMykS7fAdK7wsY997KMf/Wgmo4yyk/J3sJVuagFciXPVHtPCMAk2EbXRrdUjUVP6R6Ti3vydlK7lMutrXAvYhvMhKECgFAJbPH2zdx7wD2/ff+siR9dwGlRDFzloHmNJkof5AkG7dlvNZFSxSZJLPF3IY71K8Vm+Kmva6eIX52KXb1OeLf+OdzQtDCUzvJBsrd5LTYk7zG4Qsbf6Hq8B9zWuBZw4nAdDAQKFEXj2M5921Gv/Yd7clxQzogbScBq0mOFyHUVXr+y/s87+lYHoER6OJh1amtctVZ20b94GCgxJzgRy+aqszZyI+GlLp5us2pzt3Zzy7PluhFFi0qkAzCxUSgugm5/Qbfw/VWQPPDQOy/2OG6b4oR/KECiewBbP2Ozgfbaaf/DL8h5aQ2ggDZf3QFn516Fj0htfUmzUNH369DAYHQ89Awk3zFRW/WOPPebGqpGxv+ylYE7Mvx301MvCs0p3kjicBxYZV/VXX321940UUqYfsdTLASYVd1KvmpJV2d8HTk01kqte5NC79F6Qgiov1yaSd9FukspaeD37+YHeYVJovyfLjcOCusTTYpuO+oaWieW+xlU8Cjsu/4meqYRAMQR0s/RrX/H3uv0qv+HkXEM0767sMAGI0NNR5T3veY9XmtpZBqKDgO4vsYOVDHRA2G677S644AI3Vo2f63tlekFdLIuQcz/qqqCX6mjfsWgedNhU5bbbbmsvFYaCkYFi07PGDYV5/vz5dohLHz1slYq7xofnJQNMKnRbx3KZqqx1NQHTG9HXQ2V7Q+REU++V8I2rJddbSm84BaARNbreW8rd0yVZlnpHWoT+XnSdlv/42UZkOn2NK4f9vsUjw/ESAnkQ0K7yK3d63qfnz87DudzKeTM2riN89D86skFoBpZfHn744fZSh5Tw/74dlMIrWdLv5cuXh8518PQjUlg/aTnxCKPDlHVURq5I5s2b535MpM1Az5J2GXjrUUcd9eijjyp+r0kv6IB55JFHmsZL5tU3FOaBJ5U+aGVby1RlraXeZ/GN6EkVcRiaGlRrHKqmhWEnjApGo/s5Y8pAHrZHKz9e2YuHHsdVqApYxinB0ASBsgg891lP23X6c5Yu3OVpm2d2MJErOZRbOS9rXnmPm3hB6pZbblH66/moCXD4f1+HAmUOHpvpdyh+EvtQHd1y0kI4iozlJLwSd8cdd6hyzz33dD8yCE8sbDp+XJU3TSQ8gfCOiQVNSsN5DKtWrZIw+x74wJNKHKv6lZn9R6r+VD1C6Wj4zvZ6veO9XGQhfVxtJySexhYZIWNBIJ3ADi/975ecPOfZz9w83ayXVjmRKznsxbj6NpJMJa/ainPtNOkyBVI5bFLi+/rXv94nFd/oDoVQZqbfK1eu9C6WD/SVLutgGDkEWajhecOKFSvC04W4gQuqR6Jk99JLL/WXKQVT33A4k//HH3/ceg0wqZThqt80iqqsVdG70P+T2CLZXnrxC5Y+rgUZf8cXHycjQiCdwLZbbbF04Zwt/8fT083SW9VdTuQq3axGrdrpte00FSzsuNZavd239apXvaqv2elQFt7wpb6RzeR0b3aEiYh9vMvGjRvD04W4QbxGPpXvxuvjNaa+L3rRi7zJzjbCY2Nfk3I/NS2MqCrr/0lE6iIbMoUtZ/q4CtJ3yAsLiYEgMBiBrac9Y9kpc7Z54YCaqo7qLieDjV7NXv7/N1H5wgTx17/+tabw4he/OJyIjg9+z40K8U0+CZgkM+ximaVvJofdvRxJpsMw5Cq+cy593WabbcJRVNa9Ne5QhUjrS17S+uBc5IwhYpPyUql52BqZVNjUvPKIqnJ8IXVJ2BY+3mQ1apVNt9ZJ67v571Y/qUMMIFBBAro567MLdxlg/1ld1LGR93ZFtuXiJ+KWIfzqV7+KL6hSYR0i/BFe63XjeEqqXq7f3jcs+LlCXIDlVn3Dgbopq9+cZZ49nrBgpxphzcDlcFIDO6lFR1S5FstEkBCoDQF9A9eF79/5Fds/t/eIZawuDfjqrsQpm/LZSb8ptCcAapLYJPYaptLy8sjZQKLDcJc40aDESr/rzWLofVIlxpzJ0KhyJhhxAgEITBDQt1V/6r077b3z8yequpdkJuO6f8F19/lNCS8km0K7sZp8Zzuydy0bNSmxdmMVQldeH9nsVb0lypaCh5vMXvYdbPmPnBaYlod72hF1lH/zHM4l3svCC68We8CRgu11h1m1svP4nnk4qYiHhr1ElRu2oEwHAlUhsOiYmfqVp/RoZCCzdJsGtLq4hrJql349bzYBC/exTRrDlDe++y0BCz/aK1Zm71vQ4ca1l30HO842/CCyt0r4f/7zn/tLFaTlmovXmEKbWlul3cMVV3Tv4gW74zrUePsgVnhmEJmU921kAVVu5LIyKQiUT8B+YEq/9dQtFDXpF6ganCXbxO2TjaYrLqvKVqVqoTqagOkjy47LRE4yaX0twXUVNzM5jCifyVvEzH1GChJXeTD/atIQehn/KKaE//vf/37Y1/J4O7FQ98ilaFlqdvEkPvTgZcWv7wKTB7+ArW8UUWChxvc1Kfdc0wKqXNOFI2wI1IBAyg9MNeZnoCZdBqlLj18gKCmKiJ86Sv/sbmcNpGQ3HM7UNPJhqrhAhl0iZZ0WaFC/m1pjacRQDs1e9dpSdtVUpVRf6bikV7vi6q5ekfMAqXvvH6ZatmyZhFkqbnvsKofnKxqur0lF5li7lxl86r92cyZgCECgMAK6rVq/+KSfZF6x5nEfVD8Ddei+L6rRb0545IMVpHOSGWWWUpeIsoYO9ZWWShMlfp7+qqMeEcHzLpZBHnHEEV5jOh3u/XpTt0JE/xLN9CViilwbyx6YzBRVt8AsDH3vZqK3xEoJsx6JTQNMKtFPXSrJleuyUsQJgboSiPzAVO1+BipD7kpMU7yZvl577bUpNmGTvjxLaWVYI51WCishDyuHL0uMFbm+4atHV/q6MSW+WYWR06R6nEvxZuTKxTNnRAiMHAH7gSllzJr53rOf17yfgeplRbXfO6mZNFVa2y0HDbsrg9Su8mmnnRZWKqPVZnJYk1VZF5u1Ux3m8Sme9dWh+i7rFIO+mvKbVF9hFGaMKheGmoEmCFywIJdfGZoYgFL1CGgrW78Bpbga+VUhvfDWpVa/NbqbfXtX+PRurWG9MtHIZrh2yJXRZpWhhmOpHB8uYhC+jAQWNvVbznVS/QZTjP3UDPEVEzGjQAACEIAABJpKgOvKTV1Z5gUBCEAAAvUjgCrXb82IGAIQgAAEmkoAVW7qyjIvCEAAAhCoHwFUuX5rRsQQgAAEINBUAqhyU1eWeUEAAhCAQP0IoMr1WzMihgAEIACBphJAlZu6sswLAhCAAATqRwBVrt+aETEEIAABCDSVAKrc1JVlXhCAAAQgUD8CqHL91oyIIQABCECgqQRQ5aauLPOCAAQgAIH6EUCV67dmRAwBCEAAAk0lwG9GNXVlx+d1ynVzGz5DpgeBfAhcdNiafBzjFQJpBFDlNDrNaHvdrHnNmAizgEBhBH6wYUVhYzEQBEICqHJIo5nlsSljzZwYs4IABCDQOAJcV27ckjIhCEAAAhCoLQFy5douXe+Bkyr3zgpLCEAAAqUSQJVLxV/Q4MhyQaAZBgIQgMCQBFDlIQHWoDuaXINFIkQIQAACbQJcV+aNAAEIQAACEKgKAXLlqqxEnnGQLedJF98QgAAEsiNArpwdSzxBAAIQgAAEhiNArjwcvzr0HiNVrsMyESMEIAABEUCVR+FtgCyPwiozRwhAoAkE2MFuwioyBwhAAAIQaAYBVLkZ68gsIAABCECgCQTYwW7CKqbPge/BTudDKwQgAIHqECBXrs5aEAkEIAABCIw6AXLlEXgHcLPXCCwyU4QABJpBAFVuxjqmzwJZTudDKwQgAIGqEECVq7IS+cWBJufHFs8QgAAEsiXAdeVseeINAhCAAAQgMDgBcuXB2dWnJ9lyfdaKSCEAgdEmQK482uvP7CEAAQhAoEoEyJWrtBr5xML3YOfDFa8QgAAEsieAKmfPtHoe2cGefE0evOeRay79luzu+fH9Zn3syW+dteuMHXabMXlnLCAAAQhkRABVzggkbupM4Nz5n3Ex9nlcefE1Kn9ixYdLEWYPafW6KzwkChCAQOMJcF258UvMBCch4Pq3216zJYH2p0TZun143icm6Z9ns0LK0z2+IQCByhEgV67ckmQeEN+DnYJ09RXftSx5171mn738VGd1yLveYLmy+j5wz8PFp8sWla49eEgps6AJAhBoDAFUuTFLyUQGIfCV9ja1en5k+amR/sec/FZr3bB2o6nyufMvXNu+6nzdun82Y4m62XiN6lV53x0PmqWcqCa00QXss+Z9UpXq4t11TvDW97wpMops5OSwOe82Y+8on4e+6w2q1MNavcZtEp0nDmd+eIYABCpCAFWuyELkGQY3e3Whu3rFd61FoqicNO3RbjWhDY2lvuO9Ot3PXTCu3Fbvqj/Ryy07Gi9LedbfdWtbYm+jjLtt/2N9dXJglbPmTLdopcETZua2i3Mzc8/hcBMeKEEAAhUgwHXlCixC7iHoUM1fAgHX1J33nBVH5K2z5mzXaW0t1abGrZqWarYJr17xHVM+1Vy39vLzrjiz3bRJrw3rxqVUljKQmds8eM/D8hPWqKy/jyw7JQxgh92md16Oi7BH2M25yXDicIGrBEQj3NpaNR4QKJ4Aqlw886JH5FjbjYCvxCHz3hC38dZZu81QqyfWqndj02DptNV85eJr1SqVPXvZKapRx7aEj3sym/vu2GCvF19xpnl2G3drBvLjNSqYczV55YZ1ney5HaHq487Hx54yJWU4d0jBCTg3ChAomACqXDBwhqsQAdPUbgFZqyeybiYJt/L1nQ1we+n7ya6yqneZ9F4+qN9B5jZe4zY+qBeOOflwL7tOe413jLvymvhw3p0CBCBQOgGuK5e+BAUEoASARzqBKKKPL7jIOrQlttXqG9rtZLXV6DWHzJuryoc6mevM8R1vc+CeWwVX7ra4jjd15H8n8+w2PnTYse002ZV3DJ1bEGFNZDgz4BkCEKgIAVS5IguRYxguCzmOUXPX0jNtJvskNtzzyNofr9fLXffa6eB5cw2g19hLt5GZ1YScrXz9ijXRXusetVGk3O7HambvucN4TcxGBhtildd+6QbrqF7dOirISYczA54hAIGKEGAHuyILkWcYdtjmOUbg6JPGd4OlcBt++ohdsL3+ijWL3rXE1uPw4w+auIrbrmqp7NgUGbuNlDti436ual9mVr/Ze+wQsdF3eUZrdtnOau4P7+sOY24HsGHtozJbvOAi03vVzep0bHVvP9x5y9geMT+b9ApbKRuBcXD8A4GiCZArF02c8apDQDvP99+pDxavt79IYJLbMIH21sN3O97LKrQUt/1o3wjdKi1697ioy4Nppzdd9dnW7WDhw5Ngr5RD62V+Fv/zGWEY8hBx4q2ReneogtvEhwvNKEMAAqUTIFcufQkKCID0pyuBRZ8/+eiT3hKugaS0lf62Pje8fvH7L2pvAre6Sx3dTAbX3vMle+mfSpq123S3kYHczt5DH7hqPdTU3maWn9ajPeJ4SDotsEq3OWTeARaA6lWw+rBS9RrIwm5b+uxankLnptOhTXy4TmDuhIIRaMHkAYHiCUwd43f+iqde4IinXDd3720n7totcOR6D3Xe+y+2hFXTuKYjwPWeEtH3Q+D2x6696LA1/fTAFgLZEGAHOxuO1faic38e/RE46/Mnbbhn49nvPr+TaPbXHWsIQAACgxFAlQfjRq/mE9DW8TX3XNr8eTJDCECgSgRQ5SqtRk6xkCrnBBa3EIAABLImgCpnTbR6/hDl6q0JEUEAAhBIJoAqJ3NpVi263Kz1ZDYQgEBzCfDJqOauLTODAAQgAIG6EUCV67ZixAsBCEAAAs0lwA52c9e2M7Mx/ybGTg3/QgACEIBANQmQK1dzXYgKAhCAAARGkQC58iisOnd7jcIqM0cIQKAJBFDlJqxi+hzQ5HQ+tEIAAhCoDgFUuTprkVskfNV5bmhxDAEIQCBbAlxXzpYn3iAAAQhAAAKDEyBXHpxdXXpyD3ZdVoo4IQABCJAr8x6AAAQgAAEIVIUAuXJVViLPOLjfK5nuQ2sfPee4i9S26u5Lvv3lH65cer3Ku7xyh8OOnztz1+2sj9ucc/kp11225t6fPCiDD13yfrWq6eF1P7Ne8Y7e3W3UcfYeM9/0ztfGm+LdFc/9dz6k4dR01MJDtp/zsnhI8V4KyYK0ISJzUeUnT/i8fMqhyha55uWerRfPEIBAiQRQ5RLhFzQ0mtwN9EPrfmZNplVWlmjpb+Xdl9hLs5G8mX6rcqc9Zgqp9O9jbUU3Mz1bxyMXHhLqbmhjBjPa+tqtu437qbZ2umfJpwI4s30qILVe1T57sFb5VMGauvn8aKC7Zu96r77b77od7xBHTQECpRNgB7v0JSCA8glIqyRdUkSJn0Uj8QvDMjEzG4mu65/srVLPZr/+zoesELeRYMteiWliUziuDSdjhWRR6VTADEyS401xnx5SeGZgTuTfPVgNzxCAQEUIkCtXZCFyDYNcKBmvJ50fvfzkmbtuO2XK2KHHH2CKqHL7b0rcRr5WX7bGPMreOupZOqe+7e4t4HGbN73zH/Untw+ve8y6n3nJAhUeWrvRtXz7OdvopbVak1TczNRR0utNFp43+XBe4yG1u7RC8lMNhdo2440R4KQIgWoQQJWrsQ55RsGhN53ukQsP9l1cZ6V9Zi+re2gjaTTlDitl05Hz8a8dT7SxSFYt/aYVjnrFiVawZ4mlIlHZBV5ONIpteqterd6kviqf0dZ1NdlwqgnDNrd6jlSqV6TGLSlAAALlEmAHu1z+xYyuIzB/UQKed0rznM8jnSy2vTATuWloE6xZxGerRboob+7cPXshaJrwdOTCN3/k8pPPuGS+mamgGmuWhJ973MXePWySEi85YZma3OdOe8xwSxVcqq1y/Z0Pe4ShGeUkAhOrQwkCRRIgVy6SdkljSTt4xAhMCLBpa9tg1dJvmeHMOdvqTCbRZpPEs8P2hivHr0MfctwB46dAPmLHZryi8/Ijl9m2udu1paHz6qBjX6u/JScuN2UNBw2bWq1B/K3eHf8eUkuq25XjrgKbzmj8CwEIVIUAuXJVViK/OPQtIvwlERhHfv3la6x1yYnKO1uPsy87qWNvFdo61ob2BEarVepplRvWPmpyvssrZ7ll3Eb+zd6afFx1//aVP7RWlVXQc8dyXGP1MqVJrZHh5NBDOvDYf+x4a1ntuMcMe8lzCgHjyTMEiidArlw8c0asBIGrO2nxvT/ZcPTuJ3lMUtb2DVytCrfxVhXad1HNUq9IRzUdrES5/Ui0OaK9KZ3YpE7W+si9PzPP5seedZagQkpTN5+ai4d0w5X/M/RJGQIQqCYBVLma60JUBRGQFj5w58MSQhtP+ieFi4xtehlWnv65+eefuFw13lE2Bx2r+6snHiaHocGMXXQBu/Ww7t6kGnW3VnPiZwMmqxaSNXm01sWjlU/pbthxxz22j4TUHnxKYqU18QwBCJROYOoYPyhU+iLkGcAp1819xUvfkOcItfT98NrHPn78UoV+9mULt4/JcC2nRNCZErj7F9+96LDxz79l6hhnEJiEANeVJwFEcyMJaDe4kfNiUhCAQN0JsINd9xXsJf7xW4F6MR01G92cNXHX8qhNnvlCAALVI4AqV29Nso4ITY4TvXrpt1W58ytnAScOhxoIQKBEAqhyifALGxrpiaK+8q4LO1XA6ZDgXwhAoAIEuK5cgUUgBAhAAAIQgECbALly898I+qqE5k+SGUIAAhBoBAFUuRHLmD4JRDmdD60QgAAEKkMAVa7MUuQYCLKcI1xcQwACEMiQAKqcIcyKukKTK7owhAUBCEAgRgBVjiFpYAW63MBFZUoQgEAjCXAPdiOXlUlBAAIQgEAtCZAr13LZ+gqaTLkvXBhDAAIQKJEAqlwi/MKGRpcLQ81AEIAABIYiwA72UPjoDAEIQAACEMiQALlyhjAr64pcubJLQ2AQgAAENiGAKm+Co5Ev0ORGLiuTggAEGkmAHexGLiuTggAEIACBWhIgV67lsvUV9NgY2XJfwDCGAAQgUBoBcuXS0DMwBCAAAQhAIEKAXDkCpJEvyZUbuaxMCgIQaCABcuUGLipTggAEIACBmhIgV67pwvURNr+v3AcsTCEAAQiUSgBVLhV/IYPf/8ubCxmHQSAAAQhAYFgCU7lBd1iE9IcABCAAAQhkRIDryhmBxA0EIAABCEBgaAKo8tAIcQABCEAAAhDIiACqnBFI3EAAAhCAAASGJoAqD40QBxCAAAQgAIGMCKDKGYHEDQQgAAEIQGBoAqjy0AhxAAEIQAACEMiIAKqcEUjcQAACEIAABIYmgCoPjRAHEIAABCAAgYwIoMoZgcQNBCAAAQhAYGgCqPLQCHEAAQhAAAIQyIgAqpwRSNxAAAIQgAAEhiaAKg+NEAcQgAAEIACBjAigyhmBxA0EIAABCEBgaAKo8tAIcQABCEAAAhDIiACqnBFI3EAAAhCAAASGJoAqD40QBxCAAAQgAIGMCKDKGYHEDQQgAAEIQGBoAqjy0AhxAAEIQAACEMiIAKqcEUjcQAACEIAABIYmgCoPjRAHEIAABCAAgYwIoMoZgcQNBCAAAQhAYGgCqPLQCHEAAQhAAAIQyIgAqpwRSNxAAAIQgAAEhiaAKg+NEAcQgAAEIACBjAigyhmBxA0EIAABCEBgaAKo8tAIcQABCEAAAhDIiACqnBFI3EAAAhCAAASGJoAqD40QBxCAAAQgAIGMCKDKGYHEDQQgAAEIQGBoAqjy0AhxAAEIQAACEMiIwP8H4sV0TKzmO+EAAAAASUVORK5CYII=" + } + }, + "cell_type": "markdown", + "id": "41037e71", + "metadata": {}, + "source": [ + "![image.png](attachment:image.png)" + ] + }, + { + "cell_type": "markdown", + "id": "68e94097", + "metadata": {}, + "source": [ + "To turn the above pseudocode into a real quantum program, we still need to figure out two things:\n", + "\n", + "- How to estimate $\\left\\langle \\psi | H | \\psi \\right\\rangle$ for a given state $\\ket{\\psi}$.\n", + "- How to optimize over all quantum states $\\ket{\\psi}$.\n", + "\n", + "In the rest of this sample, we'll see how to do each of these, and how you can use Q# to write a VQE implementation." + ] + }, + { + "cell_type": "markdown", + "id": "2f7b662b", + "metadata": {}, + "source": [ + "## Estimating expectation values of $H$" + ] + }, + { + "cell_type": "markdown", + "id": "c3bb1ba6", + "metadata": {}, + "source": [ + "Before proceeding to see how to estimate $\\left\\langle \\psi | H | \\psi \\right\\rangle$, however, it helps to consider what that expectation value _is_. To do so, let's consider a concrete example of a two-qubit Hamiltonian, using QuTiP to construct a Python object representing $H$." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "a992fe41", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=True$$\\left(\\begin{array}{cc}-2 & 0 & 0 & 3\\\\0 & 7 & 0 & 1\\\\0 & 0 & -4 & 0\\\\3 & 1 & 0 & -1\\end{array}\\right)$$" + ], + "text/plain": [ + "Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=True\n", + "Qobj data =\n", + "[[-2. 0. 0. 3.]\n", + " [ 0. 7. 0. 1.]\n", + " [ 0. 0. -4. 0.]\n", + " [ 3. 1. 0. -1.]]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "H = qt.Qobj([\n", + " [-2, 0, 0, 3],\n", + " [0, 7, 0, 1],\n", + " [0, 0, -4, 0],\n", + " [3, 1, 0, -1]\n", + "], dims=[[2, 2]] * 2)\n", + "H" + ] + }, + { + "cell_type": "markdown", + "id": "4b09b667", + "metadata": {}, + "source": [ + "Here, `H` is a Python object representing the 4 × 4 matrix $H$. We can use the `eigenstates` method provided by QuTiP to quickly find the eigenvalues and eigenvectors of $H$:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "33551a8c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(array([-4.57776672, -4. , 1.43800535, 7.13976137]),\n", + " array([Quantum object: dims=[[2, 2], [1, 1]], shape=(4, 1), type='ket', dtype=Dense\n", + " Qobj data =\n", + " [[ 0.75726547]\n", + " [ 0.05620122]\n", + " [ 0. ]\n", + " [-0.65068458]] ,\n", + " Quantum object: dims=[[2, 2], [1, 1]], shape=(4, 1), type='ket', dtype=Dense\n", + " Qobj data =\n", + " [[0.]\n", + " [0.]\n", + " [1.]\n", + " [0.]] ,\n", + " Quantum object: dims=[[2, 2], [1, 1]], shape=(4, 1), type='ket', dtype=Dense\n", + " Qobj data =\n", + " [[ 0.65152827]\n", + " [-0.13424187]\n", + " [-0. ]\n", + " [ 0.74665256]] ,\n", + " Quantum object: dims=[[2, 2], [1, 1]], shape=(4, 1), type='ket', dtype=Dense\n", + " Qobj data =\n", + " [[0.04538633]\n", + " [0.9893536 ]\n", + " [0. ]\n", + " [0.13827341]] ],\n", + " dtype=object))" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "H.eigenstates()" + ] + }, + { + "cell_type": "markdown", + "id": "ba4c191e", + "metadata": {}, + "source": [ + "In particular, we can use the `min` function provided by Python to minimize over the eigenvalues of $H$ and find its ground state $\\ket{\\phi_0}$." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "3199b955", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(-4.57776672447103,\n", + " Quantum object: dims=[[2, 2], [1, 1]], shape=(4, 1), type='ket', dtype=Dense\n", + " Qobj data =\n", + " [[ 0.75726547]\n", + " [ 0.05620122]\n", + " [ 0. ]\n", + " [-0.65068458]])" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "min_energy, ground_state = min(zip(*H.eigenstates()), key=lambda eig: eig[0])\n", + "min_energy, ground_state" + ] + }, + { + "cell_type": "markdown", + "id": "2ce8638c", + "metadata": {}, + "source": [ + "Here, `ground_state` represents the eigenvector $\\ket{\\phi_0}$ corresponding to the smallest eigenvalue $E_0$ of $H$. By the above argument, we would expect that $\\left\\langle \\phi_0 | H | \\phi_0 \\right\\rangle = E_0$. We can check that using QuTiP as well:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "a2140642", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(-4.57776672447103+0j)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(ground_state.dag() * H * ground_state)" + ] + }, + { + "cell_type": "markdown", + "id": "853df02f", + "metadata": {}, + "source": [ + "What's going on here? Effectively, for any operator $O$ such that $O^{\\dagger} = O$, expressions like $\\left\\langle \\psi | O | \\psi \\right\\rangle$ represent another way of thinking about quantum measurement. In particular, if we think of the eigenvalues of $O$ as labels for its corresponding various eigenvectors, then $\\left\\langle \\psi | O | \\psi \\right\\rangle$ is the average label that we get when we measure $O$ in the basis of its eigenvectors." + ] + }, + { + "cell_type": "markdown", + "id": "8537e34c", + "metadata": {}, + "source": [ + "> **💡 TIP:** This way of thinking about quantum measurement is sometimes called the _observable framework_, with $O$ being called an _observable_. We avoid using this terminology here to avoid confusion, however, as the expectation value of $O$ cannot be observed directly, but only inferred from repeated measurements." + ] + }, + { + "cell_type": "markdown", + "id": "817f94ec", + "metadata": {}, + "source": [ + "For example, if $O = Z$, then the eigenstate $\\ket{0}$ is labeled by its eigenvalue $+1$, while $\\ket{1}$ is labeled by $-1$. The expectation value $\\left\\langle \\psi | Z | \\psi \\right\\rangle$ is then the probability of getting a `Zero` minus the probability of getting a `One`.\n", + "\n", + "More generally, finding the expectation value of an arbitrary Pauli operator for an arbitrary input state is straightforward using the `Measure` and `EstimateFrequencyA` operations." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "d6558958", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "import Chemistry.JordanWigner.JordanWignerVQE.EstimateFrequency;\n", + "\n", + "operation EstimatePauliExpectation(pauli : Pauli[], preparation : (Qubit[] => Unit is Adj), nShots : Int) : Double {\n", + " return 2.0 * EstimateFrequency(\n", + " preparation,\n", + " Measure(pauli, _),\n", + " Length(pauli),\n", + " nShots\n", + " ) - 1.0;\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "2d07f00c", + "metadata": {}, + "source": [ + "Here, we've represented the state $\\ket{\\psi}$ by an operation `preparation` that prepares that state. Since each Pauli operator other than the identity operator has exactly two eigenvalues, $+1$ and $-1$, corresponding to `Zero` and `One` respectively, we can turn the estimate of the probability $p_0$ with which `Measure(pauli, _)` returns `Zero` into an expectation value $p_0 - p_1 = p_0 - (1 - p_0) = 2 p_0 - 1$.\n", + "\n", + "For example, doing nothing prepares the state $\\ket{0}$, so we get an expectation value of $1$:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "22a23913", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "operation EstimateExpectationOfZero() : Double {\n", + " return EstimatePauliExpectation(\n", + " [PauliZ],\n", + " qs => I(qs[0]),\n", + " 100\n", + " );\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "1135713f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.0" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "qsharp.code.EstimateExpectationOfZero()" + ] + }, + { + "cell_type": "markdown", + "id": "aef4f547", + "metadata": {}, + "source": [ + "Similarly, using `X` to prepare $\\ket{1}$ gives us an expectation of $-1$, while using `H` to prepare $\\ket{+} = \\frac{1}{\\sqrt{2}} (\\ket{0} + \\ket{1})$ gives an expectation value of $0$." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "8dbf0842", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "operation EstimateExpectationOfOne() : Double {\n", + " return EstimatePauliExpectation(\n", + " [PauliZ],\n", + " ApplyToEachCA(X, _),\n", + " 100\n", + " );\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "2cab1fca", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "-1.0" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "qsharp.code.EstimateExpectationOfOne()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "408d1da3", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "operation EstimateExpectationOfPlus() : Double {\n", + " return EstimatePauliExpectation(\n", + " [PauliZ],\n", + " ApplyToEachCA(H, _),\n", + " 100\n", + " );\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "fc6d8273", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "-0.040000000000000036" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "qsharp.code.EstimateExpectationOfPlus()" + ] + }, + { + "cell_type": "markdown", + "id": "6d023ef6", + "metadata": {}, + "source": [ + "Note that in practice, we won't always get exactly 0 due to using a finite number of measurements." + ] + }, + { + "cell_type": "markdown", + "id": "b2d9316d", + "metadata": {}, + "source": [ + "In any case, to recap, it's easy to use a quantum program to find $\\left\\langle \\psi | H | \\psi \\right\\rangle$ in the special case that $H$ is a multi-qubit Pauli operator. What about the more general case? The linearity of quantum mechanics saves us again! It turns out that we can expand the expectation of any other operator in terms of expectations of Pauli operators. To see how that works, suppose that $H = 2 Z - X$." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "ae9236c7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=CSR, isherm=True$$\\left(\\begin{array}{cc}2 & -1\\\\-1 & -2\\end{array}\\right)$$" + ], + "text/plain": [ + "Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=CSR, isherm=True\n", + "Qobj data =\n", + "[[ 2. -1.]\n", + " [-1. -2.]]" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "2 * qt.sigmaz() - qt.sigmax()" + ] + }, + { + "cell_type": "markdown", + "id": "8d7db2d3", + "metadata": {}, + "source": [ + "We can expand $\\left\\langle \\psi | (2Z - X) | \\psi \\right\\rangle$ by using linearity:\n", + "\n", + "$$\n", + " \\left\\langle \\psi | (2Z - X) | \\psi \\right\\rangle = 2 \\left\\langle \\psi | Z | \\psi \\right\\rangle - \\left\\langle \\psi | X | \\psi \\right\\rangle.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "797bd006", + "metadata": {}, + "source": [ + "Each of the terms in this expansion is something that we can estimate easily using our `EstimatePauliExpectation` operation from above." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "d10b40ab", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "operation EstimateExpectation(terms : (Double, Pauli[])[], preparation : (Qubit[] => Unit is Adj), nShots : Int) : Double {\n", + " mutable sum = 0.0;\n", + " for (coefficient, pauli) in terms {\n", + " sum += coefficient * EstimatePauliExpectation(\n", + " pauli, preparation, nShots\n", + " );\n", + " }\n", + " return sum;\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "49bf9fcc", + "metadata": {}, + "source": [ + "With this, we almost have everything we need to estimate expectations $\\left\\langle \\psi | H | \\psi \\right\\rangle$. We just need a way of finding the decomposition of $H$ into Pauli operators. Thankfully, QuTiP can help here as well." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "5496baa1", + "metadata": {}, + "outputs": [], + "source": [ + "def pauli_to_qobj(pauli):\n", + " if pauli == qsharp.Pauli.I:\n", + " return qt.Qobj(qt.identity(2))\n", + " elif pauli == qsharp.Pauli.X:\n", + " return qt.Qobj(qt.sigmax())\n", + " elif pauli == qsharp.Pauli.Y:\n", + " return qt.Qobj(qt.sigmay())\n", + " elif pauli == qsharp.Pauli.Z:\n", + " return qt.Qobj(qt.sigmaz())" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "731556db", + "metadata": {}, + "outputs": [], + "source": [ + "def pauli_basis(n_qubits):\n", + " scale = 2 ** n_qubits\n", + " return [\n", + " tuple([P, qt.tensor(*(pauli_to_qobj(p) for p in P)) / scale])\n", + " for P in product([qsharp.Pauli.I, qsharp.Pauli.X, qsharp.Pauli.Y, qsharp.Pauli.Z], repeat=n_qubits)\n", + " ]" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "c8e0b922", + "metadata": {}, + "outputs": [], + "source": [ + "def expand_in_pauli_basis(op):\n", + " return [\n", + " (coeff.real, list(label))\n", + " for label, coeff in [\n", + " tuple([label, (P.dag() * op).tr()])\n", + " for label, P in pauli_basis(n_qubits=len(op.dims[0]))\n", + " ]\n", + " if abs(coeff) >= 1e-10\n", + " ]" + ] + }, + { + "cell_type": "markdown", + "id": "ff5c27bc", + "metadata": {}, + "source": [ + "For example, we can use the two functions above to find that $H = -3 𝟙 \\otimes Z + 0.5 X \\otimes 𝟙 + 1.5 X \\otimes X - 0.5 X \\otimes Z - 1.5 Y \\otimes Y + 2.5 Z \\otimes 𝟙 - 1.5 Z \\otimes Z$:" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "cd8e5438", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[(-3.0, [Pauli.I, Pauli.Z]),\n", + " (0.5, [Pauli.X, Pauli.I]),\n", + " (1.5, [Pauli.X, Pauli.X]),\n", + " (-0.5, [Pauli.X, Pauli.Z]),\n", + " (-1.5, [Pauli.Y, Pauli.Y]),\n", + " (2.5, [Pauli.Z, Pauli.I]),\n", + " (-1.5, [Pauli.Z, Pauli.Z])]" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "H_decomposition = expand_in_pauli_basis(H)\n", + "H_decomposition" + ] + }, + { + "cell_type": "markdown", + "id": "8879b674", + "metadata": {}, + "source": [ + "This decomposition is exactly what we need to pass as `terms` above. For example, to estimate the expectation of the $\\ket{++}$ state, $\\left\\langle ++ | H | ++ \\right\\rangle$, we can pass `terms` to `EstimateExpectation`. In doing so, we'll use the name \"energy\" for our new operation, pointing to that expectation values of Hamiltonian operators $H$ represent the average energy of a system given the quantum state of that system." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "fad87b56", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "operation EstimateEnergyOfPlus(terms : (Double, Pauli[])[]) : Double {\n", + " return EstimateExpectation(\n", + " terms,\n", + " ApplyToEachCA(H, _),\n", + " 100\n", + " );\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "1f0e22df", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.9399999999999995" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "qsharp.code.EstimateEnergyOfPlus(H_decomposition)" + ] + }, + { + "cell_type": "markdown", + "id": "29a48a9f", + "metadata": {}, + "source": [ + "This is pretty far from the minimum $E_0$ we found from the eigenvalue decomposition above, so in the next part we'll see one more trick we can use to write our VQE implementation." + ] + }, + { + "cell_type": "markdown", + "id": "9bfad16a", + "metadata": {}, + "source": [ + "## Preparing ansatz states" + ] + }, + { + "cell_type": "markdown", + "id": "3ae1ddfe", + "metadata": {}, + "source": [ + "Recall that our goal is to use a classical optimizer to find the minimum energy $E_0$ of some Hamiltonian operator $H$:\n", + "$$\n", + " E_0 = \\min_{\\ket{\\psi}} \\left\\langle \\psi | H | \\psi \\right\\rangle.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "fe83f982", + "metadata": {}, + "source": [ + "In practice, we can't reasonably optimize over all possible quantum states $\\ket{\\psi}$, so that the above optimization problem is intractable for all but the smallest systems. Instead, we note that we can find an upper-bound for $E_0$ by solving a simpler problem. Suppose that there is a set $A$ of states that are easier to prepare. Then,\n", + "$$\n", + " E_0 = \\min_{\\ket{\\psi}} \\left\\langle \\psi | H | \\psi \\right\\rangle \\le \\min_{\\ket{\\psi} \\in A} \\left\\langle \\psi | H | \\psi \\right\\rangle.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "3926dbbf", + "metadata": {}, + "source": [ + "If we choose $A$ carefully so that $\\ket{\\phi_0} \\in A$, this inequality will be tight. More often, however, $\\ket{\\phi_0}$ won't actually be in $A$, but will be close to some other state in $A$, giving us a reasonable approximation to $E_0$." + ] + }, + { + "cell_type": "markdown", + "id": "c6bce67e", + "metadata": {}, + "source": [ + "We say that $A$ is our _ansatz_ for running VQE; in this way, $A$ plays a similar role to a model in an ML training loop, representing what we believe to be a reasonable set of guesses for $\\ket{\\phi_0}$. One can get pretty involved with their choice of ansatz, but for this sample, we'll choose a really simple one, and set $A$ to be those states that can be prepared by a small number of Pauli rotations." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "950a5184", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "operation PrepareAnsatz(axes : Pauli[][], angles : Double[], register : Qubit[]) : Unit is Adj {\n", + " for i in IndexRange(angles) {\n", + " Exp(axes[i], angles[i], register);\n", + " }\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "b377850f", + "metadata": {}, + "source": [ + "Our minimization $\\min_{\\ket{\\psi} \\in A} \\left\\langle \\psi | H | \\psi \\right\\rangle$ is now a minimization over `angles` for some fixed list of Pauli rotation axes. Using `DumpMachine` and `DumpRegister`, we can visualize how this ansatz works for a few different choices of parameters `angles`, convincing ourselves that $A$ contains enough interesting states to find $E_0$." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "ec65a772", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "operation DumpAnsatz(axes : Pauli[][], angles : Double[]) : Unit {\n", + " use register = Qubit[Length(angles)];\n", + " within {\n", + " PrepareAnsatz(axes, angles, register);\n", + " } apply {\n", + " DumpRegister(register);\n", + " }\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "b99c1d78", + "metadata": {}, + "outputs": [], + "source": [ + "ansatz_axes = [[qsharp.Pauli.Z, qsharp.Pauli.Z], [qsharp.Pauli.X, qsharp.Pauli.Y]]" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "7df99a83", + "metadata": {}, + "outputs": [ + { + "data": { + "text/markdown": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + "
Basis State
(|𝜓₁…𝜓ₙ⟩)		AmplitudeMeasurement ProbabilityPhase
\n", + " |00⟩\n", + " \n", + " −0.1171−0.3013𝑖\n", + " \n", + " \n", + " 10.4516%\n", + " ↑\n", + " -1.9416\n", + "
\n", + " |11⟩\n", + " \n", + " −0.3429−0.8820𝑖\n", + " \n", + " \n", + " 89.5484%\n", + " ↑\n", + " -1.9416\n", + "
\n" + ], + "text/plain": [ + "STATE:\n", + "|00⟩: −0.1171−0.3013𝑖\n", + "|11⟩: −0.3429−0.8820𝑖" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "qsharp.code.DumpAnsatz(ansatz_axes, [1.2, 1.9])" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "fc79ef91", + "metadata": {}, + "outputs": [ + { + "data": { + "text/markdown": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + "
Basis State
(|𝜓₁…𝜓ₙ⟩)		AmplitudeMeasurement ProbabilityPhase
\n", + " |00⟩\n", + " \n", + " 0.2771+0.7129𝑖\n", + " \n", + " \n", + " 58.4984%\n", + " ↑\n", + " 1.2000\n", + "
\n", + " |11⟩\n", + " \n", + " 0.2334+0.6004𝑖\n", + " \n", + " \n", + " 41.5016%\n", + " ↑\n", + " 1.2000\n", + "
\n" + ], + "text/plain": [ + "STATE:\n", + "|00⟩: 0.2771+0.7129𝑖\n", + "|11⟩: 0.2334+0.6004𝑖" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "qsharp.code.DumpAnsatz(ansatz_axes, [1.2, -0.7])" + ] + }, + { + "cell_type": "markdown", + "id": "6d2f1e81", + "metadata": {}, + "source": [ + "At this point, it's helpful to write some new operations that directly estimate the energy of a state given a parameterization of our ansatz, rather than an opaque operation that prepares the state." + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "11583999", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "operation EstimateEnergyAtAnsatz(terms : (Double, Pauli[])[], axes : Pauli[][], angles : Double[], nShots : Int) : Double {\n", + " return EstimateExpectation(\n", + " terms,\n", + " PrepareAnsatz(axes, angles, _),\n", + " nShots\n", + " );\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "id": "418ace1a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.7390000000000003" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "qsharp.code.EstimateEnergyAtAnsatz(\n", + " H_decomposition,\n", + " ansatz_axes,\n", + " [1.2, 1.9],\n", + " 1000\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "6f07737a", + "metadata": {}, + "source": [ + "## Running VQE in Q\\#" + ] + }, + { + "cell_type": "markdown", + "id": "b0970ee4", + "metadata": {}, + "source": [ + "Using this new operation, we have something that looks much more like classical optimization problems that we may be used to! Indeed, we can simply use our favorite optimization algorithm to pick different values for `angles` until we find a good approximation for $E_0$. In this sample, we'll use the _SPSA algorithm_ to optimize `angles`, as this algorithm works especially well for objective functions that have some amount of noise, such as that added by taking a finite number of shots above.\n", + "\n", + "We won't go through the details of SPSA here, but if you're interested to learn more, check out [`SPSA.qs`](./SPSA/src/SPSA.qs) in this folder to see how we implemented SPSA in Q#." + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "7f3049df", + "metadata": { + "vscode": { + "languageId": "qsharp" + } + }, + "outputs": [], + "source": [ + "%%qsharp\n", + "\n", + "import SPSA.*;\n", + "\n", + "operation FindMinimumEnergy(terms : (Double, Pauli[])[], axes : Pauli[][], initialGuess : Double[], nShots : Int) : (Double[], Double) {\n", + " let oracle = EstimateEnergyAtAnsatz(terms, axes, _, nShots);\n", + " let options = DefaultSpsaOptions();\n", + " return FindMinimumWithSpsa(\n", + " oracle,\n", + " initialGuess,\n", + " options\n", + " );\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "id": "337c5318", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Iteration 1:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([1.0999999999999999, 2.0]) = 1.1939999999999995" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([1.3, 1.7999999999999998]) = 0.25900000000000034" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 2:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([5.043238648643683, -1.9432386486436832]) = -3.138" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([4.856761351356317, -1.756761351356317]) = -2.083" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 3:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([9.084878774549992, -5.80588383667864]) = -4.171000000000001" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([8.90588383667864, -5.984878774549993]) = -3.568" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 4:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([10.535640921734934, -4.355121689493698]) = 0.757000000000001" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([10.36177200971693, -4.528990601511703]) = -0.07900000000000018" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 5:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([8.553949974759709, -6.336812636468924]) = -1.702" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([8.383955805601582, -6.5068068056270505]) = -0.6140000000000003" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 6:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([11.097417268504492, -3.9602378039359896]) = 1.4259999999999997" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([10.930524807292644, -3.7933453427241397]) = 1.1429999999999998" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 7:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([10.432605046730536, -3.295425582162034]) = -1.021" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([10.596919245566806, -3.459739780998304]) = -0.22099999999999986" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 8:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([8.721082781473008, -1.7460163408815679]) = -2.064000000000001" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([8.883195805450072, -1.583903316904506]) = -1.0920000000000005" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 9:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([7.222224052110989, -3.4050710091614342]) = -0.4760000000000002" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([7.062028113193143, -3.244875070243587]) = -1.3249999999999997" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 10:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([5.311727845447455, -1.653075068593994]) = -1.5130000000000003" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([5.470228111543549, -1.4945748024978993]) = -0.4889999999999999" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 11:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.744449775267648, -3.0633713279524746]) = -2.4690000000000003" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.901431586088974, -3.2203531387738003]) = -1.5310000000000001" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 12:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.4885394395173024, -1.8074609922021285]) = -2.3770000000000002" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.3329311639951906, -1.6518527166800168]) = -1.5099999999999998" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 13:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([4.123518366011016, -3.442439918695842]) = -0.09000000000000008" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.969163005265829, -3.2880845579506555]) = -0.9939999999999993" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 14:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.797839071673697, -2.1167606243585233]) = -3.904" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.644634733149064, -1.9635562858338904]) = -3.2169999999999996" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 15:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.635425375809424, -2.9543469284942505]) = -3.1829999999999994" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.7875658545441793, -3.1064874072290056]) = -2.2909999999999995" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 16:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.4604292880354066, -1.779350840720233]) = -2.1570000000000005" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.6115812808670293, -1.9305028335518557]) = -3.0879999999999996" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 17:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.838960914575619, -3.1578824672604453]) = -1.8079999999999998" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.688731610781285, -3.0076531634661112]) = -2.878" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 18:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.729021712127196, -2.047943264812022]) = -3.569" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.8783862407562566, -2.197307793441083]) = -4.287999999999999" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 19:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.289955839527876, -2.6088773922127024]) = -4.449000000000001" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.438506942923287, -2.757428495608113]) = -4.083" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 20:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.62494808322752, -1.9438696359123462]) = -3.3240000000000003" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.4771645770228794, -1.7960861297077058]) = -2.428" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 21:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.417129210089135, -2.5889937137135655]) = -4.404" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.270072161028739, -2.7360507627739614]) = -3.9350000000000005" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 22:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.5668971716627156, -2.439225752139985]) = -4.494" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.420529452636529, -2.5855934711661717]) = -4.356" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 23:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.8513765864502267, -2.154746337352474]) = -3.9999999999999996" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.705664531092227, -2.3004583927104734]) = -4.415" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 24:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.0719990233208585, -2.7890368472072855]) = -3.735" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([3.217086076595415, -2.934123900481842]) = -3.1500000000000004" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 25:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.3814949330501607, -2.0985327569365877]) = -3.7260000000000004" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.52598502078897, -2.243022844675397]) = -4.365" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 26:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.954317511810305, -2.671355335696732]) = -4.126" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.81039865934241, -2.527436483228837]) = -4.388" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 27:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.560155148540727, -2.277192972427154]) = -4.227" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.416783837871512, -2.133821661757939]) = -4.069" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 28:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.79396677547435, -2.511004599360777]) = -4.646" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.936812430378211, -2.653850254264638]) = -4.301" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 29:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.539000396034041, -2.398378494045976]) = -4.606" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.681340670159549, -2.256038219920468]) = -4.342" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Iteration 30:" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.53489122838635, -2.402487661693667]) = -4.338" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " obj([2.3930375019740335, -2.544341388105983]) = -4.519" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "([2.621354808572941, -2.316024081507076], -4.361)" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "qsharp.code.FindMinimumEnergy(H_decomposition, ansatz_axes, [1.2, 1.9], 1000)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "d6276b06", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "-4.57776672447103" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "min_energy" + ] + } + ], + "metadata": { + "interpreter": { + "hash": "464ea5c5cae2889b0c59ecc40909802f0909ef10120e6598fbb774b74247ba9f" + }, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}