Some cleanup of chemistry library (#2158)

Redundant code removed (Especially state preparation and LE comparison),
some renames, some cleanup of old syntax.
Added test for superposition preparation.

---------

Co-authored-by: Dmitry Vasilevsky <[email protected]>

Author: DmitryVasilevsky

library/chemistry/src/JordanWigner/Convenience.qs

@@ -12,23 +12,21 @@ import JordanWigner.JordanWignerEvolutionSet.JordanWignerGeneratorSystem;
 import JordanWigner.JordanWignerOptimizedBlockEncoding.JordanWignerOptimizedBlockEncoding;
 import JordanWigner.JordanWignerOptimizedBlockEncoding.PauliBlockEncoding;
 import JordanWigner.JordanWignerOptimizedBlockEncoding.QuantumWalkByQubitization;
-import JordanWigner.StatePreparation.PrepareArbitraryStateD;
 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;
-import Utils.GeneratorSystem;
 
 // Convenience functions for performing simulation.
 
 /// # Summary
 /// Returns Trotter step operation and the parameters necessary to run it.
 ///
 /// # Input
-/// ## qSharpData
+/// ## jwHamiltonian
 /// Hamiltonian described by `JordanWignerEncodingData` format.
 /// ## trotterStepSize
 /// Step size of Trotter integrator.
@@ -39,24 +37,22 @@ import Utils.GeneratorSystem;
 /// A tuple where: `Int` is the number of qubits allocated,
 /// `Double` is `1.0/trotterStepSize`, and the operation
 /// is the Trotter step.
-function TrotterStepOracle(qSharpData : JordanWignerEncodingData, trotterStepSize : Double, trotterOrder : Int) : (Int, (Double, (Qubit[] => Unit is Adj + Ctl))) {
-    let (nSpinOrbitals, data, statePrepData, energyShift) = qSharpData!;
-    let generatorSystem = JordanWignerGeneratorSystem(data);
+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 = nSpinOrbitals;
+    let nTargetRegisterQubits = jwHamiltonian.NumQubits;
     let rescaleFactor = 1.0 / trotterStepSize;
     return (nTargetRegisterQubits, (rescaleFactor, oracle));
 }
 
 
-function QubitizationOracleSeperatedRegisters(qSharpData : JordanWignerEncodingData) : ((Int, Int), (Double, ((Qubit[], Qubit[]) => Unit is Adj + Ctl))) {
-    let (nSpinOrbitals, data, statePrepData, energyShift) = qSharpData!;
-    let generatorSystem = JordanWignerBlockEncodingGeneratorSystem(data);
+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 = nSpinOrbitals;
+    let nTargetRegisterQubits = jwHamiltonian.NumQubits;
     let nCtrlRegisterQubits = Ceiling(Lg(IntAsDouble(nTerms)));
     return ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, QuantumWalkByQubitization(blockEncodingReflection)));
 }
@@ -65,15 +61,15 @@ function QubitizationOracleSeperatedRegisters(qSharpData : JordanWignerEncodingD
 /// Returns Qubitization operation and the parameters necessary to run it.
 ///
 /// # Input
-/// ## qSharpData
+/// ## 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(qSharpData : JordanWignerEncodingData) : (Int, (Double, (Qubit[] => Unit is Adj + Ctl))) {
-    let ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, oracle)) = QubitizationOracleSeperatedRegisters(qSharpData);
+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, _)));
 }
@@ -84,9 +80,8 @@ operation MergeTwoRegisters(oracle : ((Qubit[], Qubit[]) => Unit is Adj + Ctl),
 }
 
 
-function OptimizedQubitizationOracleSeperatedRegisters(qSharpData : JordanWignerEncodingData, targetError : Double) : ((Int, Int), (Double, ((Qubit[], Qubit[]) => Unit is Adj + Ctl))) {
-    let (nSpinOrbitals, data, statePrepData, energyShift) = qSharpData!;
-    let ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, blockEncodingReflection)) = JordanWignerOptimizedBlockEncoding(targetError, data, nSpinOrbitals);
+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)));
 }
 
@@ -96,7 +91,7 @@ function OptimizedQubitizationOracleSeperatedRegisters(qSharpData : JordanWigner
 /// and the parameters necessary to run it.
 ///
 /// # Input
-/// ## qSharpData
+/// ## jwHamiltonian
 /// Hamiltonian described by `JordanWignerEncodingData` format.
 /// ## targetError
 /// Error of auxillary state preparation step.
@@ -105,8 +100,8 @@ function OptimizedQubitizationOracleSeperatedRegisters(qSharpData : JordanWigner
 /// 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(qSharpData : JordanWignerEncodingData, targetError : Double) : (Int, (Double, (Qubit[] => Unit is Adj + Ctl))) {
-    let ((nCtrlRegisterQubits, nTargetRegisterQubits), (oneNorm, oracle)) = OptimizedQubitizationOracleSeperatedRegisters(qSharpData, targetError);
+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, _)));
 }

library/chemistry/src/JordanWigner/ConvenienceVQE.qs

@@ -13,7 +13,6 @@ import JordanWigner.StatePreparation.PrepareTrialState;
 import JordanWigner.Utils.JordanWignerEncodingData;
 import JordanWigner.Utils.JordanWignerInputState;
 import JordanWigner.Utils.JWOptimizedHTerms;
-import Std.Arrays.ForEach;
 
 /// # Summary
 /// Estimates the energy of the molecule by summing the energy contributed by the individual Jordan-Wigner terms.

library/chemistry/src/JordanWigner/JordanWignerBlockEncoding.qs

@@ -36,8 +36,8 @@ import Utils.IsNotZero;
 /// # Output
 /// `GeneratorIndex[]` expressing Z term as Pauli terms.
 function ZTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] {
-    let ((idxTermType, coeff), idxFermions) = term!;
-    return [new GeneratorIndex { Term = ([3], coeff), Subsystem = idxFermions }];
+    let (_, coeff) = term.Term;
+    return [new GeneratorIndex { Term = ([3], coeff), Subsystem = term.Subsystem }];
 }
 
 /// # Summary
@@ -51,8 +51,8 @@ function ZTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] {
 /// # Output
 /// `GeneratorIndex[]` expressing ZZ term as Pauli terms.
 function ZZTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] {
-    let ((idxTermType, coeff), idxFermions) = term!;
-    return [new GeneratorIndex { Term = ([3, 3], coeff), Subsystem = idxFermions }];
+    let (_, coeff) = term.Term;
+    return [new GeneratorIndex { Term = ([3, 3], coeff), Subsystem = term.Subsystem }];
 }
 
 /// # Summary
@@ -66,10 +66,10 @@ function ZZTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] {
 /// # Output
 /// `GeneratorIndex[]` expressing PQ term as Pauli terms.
 function PQTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] {
-    let ((idxTermType, coeff), idxFermions) = term!;
+    let (_, coeff) = term.Term;
     let newCoeff = [coeff[0]];
-    let qubitPidx = idxFermions[0];
-    let qubitQidx = idxFermions[1];
+    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 },
@@ -88,15 +88,15 @@ function PQTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] {
 /// # Output
 /// `GeneratorIndex[]` expressing PQ or PQQR term as Pauli terms.
 function PQandPQQRTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] {
-    let ((idxTermType, coeff), idxFermions) = term!;
+    let (_, coeff) = term.Term;
     let newCoeff = [coeff[0]];
 
-    if Length(idxFermions) == 2 {
+    if Length(term.Subsystem) == 2 {
         return PQTermToPauliGenIdx(term);
     } else {
-        let qubitPidx = idxFermions[0];
-        let qubitQidx = idxFermions[1];
-        let qubitRidx = idxFermions[3];
+        let qubitPidx = term.Subsystem[0];
+        let qubitQidx = term.Subsystem[1];
+        let qubitRidx = term.Subsystem[3];
 
         if (qubitPidx < qubitQidx and qubitQidx < qubitRidx) {
 
@@ -129,9 +129,9 @@ function PQandPQQRTermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] {
 /// # Output
 /// `GeneratorIndex[]` expressing PQRS term as Pauli terms.
 function V0123TermToPauliGenIdx(term : GeneratorIndex) : GeneratorIndex[] {
-    let ((idxTermType, v0123), idxFermions) = term!;
-    let qubitsPQ = idxFermions[0..1];
-    let qubitsRS = idxFermions[2..3];
+    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]];

library/chemistry/src/JordanWigner/JordanWignerEvolutionSet.qs

@@ -6,7 +6,6 @@ export
     JordanWignerFermionEvolutionSet;
 
 import JordanWigner.Utils.JWOptimizedHTerms;
-import Std.Arrays.Fold;
 import Std.Arrays.IndexRange;
 import Std.Arrays.Subarray;
 import Utils.GeneratorIndex;
@@ -40,9 +39,9 @@ import Utils.IsNotZero;
 /// ## qubits
 /// Qubits of Hamiltonian.
 operation ApplyJordanWignerZTerm(term : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl {
-    let ((idxTermType, coeff), idxFermions) = term!;
+    let (_, coeff) = term.Term;
     let angle = (1.0 * coeff[0]) * stepSize;
-    let qubit = qubits[idxFermions[0]];
+    let qubit = qubits[term.Subsystem[0]];
     Exp([PauliZ], angle, [qubit]);
 }
 
@@ -58,9 +57,9 @@ operation ApplyJordanWignerZTerm(term : GeneratorIndex, stepSize : Double, qubit
 /// ## qubits
 /// Qubits of Hamiltonian.
 operation ApplyJordanWignerZZTerm(term : GeneratorIndex, stepSize : Double, qubits : Qubit[]) : Unit is Adj + Ctl {
-    let ((idxTermType, coeff), idxFermions) = term!;
+    let (_, coeff) = term.Term;
     let angle = (1.0 * coeff[0]) * stepSize;
-    let qubitsZZ = Subarray(idxFermions[0..1], qubits);
+    let qubitsZZ = Subarray(term.Subsystem[0..1], qubits);
     Exp([PauliZ, PauliZ], angle, qubitsZZ);
 }
 
@@ -190,48 +189,18 @@ function JordanWignerGeneratorSystem(data : JWOptimizedHTerms) : GeneratorSystem
 /// A `GeneratorSystem` representing a system that is the sum of the
 /// input generator systems.
 function AddGeneratorSystems(generatorSystemA : GeneratorSystem, generatorSystemB : GeneratorSystem) : GeneratorSystem {
-    let nTermsA = GetGeneratorSystemNTerms(generatorSystemA);
-    let nTermsB = GetGeneratorSystemNTerms(generatorSystemB);
-    let generatorIndexFunctionA = GetGeneratorSystemFunction(generatorSystemA);
-    let generatorIndexFunctionB = GetGeneratorSystemFunction(generatorSystemB);
+    let nTermsA = generatorSystemA.NumEntries;
+    let nTermsB = generatorSystemB.NumEntries;
     let generatorIndexFunction = idx -> {
         if idx < nTermsA {
-            return generatorIndexFunctionA(idx);
+            return generatorSystemA.EntryAt(idx);
         } else {
-            return generatorIndexFunctionB(idx - nTermsA);
+            return generatorSystemB.EntryAt(idx - nTermsA);
         }
     };
     return new GeneratorSystem { NumEntries = nTermsA + nTermsB, EntryAt = generatorIndexFunction };
 }
 
-/// # Summary
-/// Retrieves the number of terms in a `GeneratorSystem`.
-///
-/// # Input
-/// ## generatorSystem
-/// The `GeneratorSystem` of interest.
-///
-/// # Output
-/// The number of terms in a `GeneratorSystem`.
-function GetGeneratorSystemNTerms(generatorSystem : GeneratorSystem) : Int {
-    let (nTerms, generatorIndexFunction) = generatorSystem!;
-    return nTerms;
-}
-
-/// # Summary
-/// Retrieves the `GeneratorIndex` function in a `GeneratorSystem`.
-///
-/// # Input
-/// ## generatorSystem
-/// The `GeneratorSystem` of interest.
-///
-/// # Output
-/// An function that indexes each `GeneratorIndex` term in a Hamiltonian.
-function GetGeneratorSystemFunction(generatorSystem : GeneratorSystem) : (Int -> GeneratorIndex) {
-    let (nTerms, generatorIndexFunction) = generatorSystem!;
-    return generatorIndexFunction;
-}
-
 /// # Summary
 /// Represents a dynamical generator as a set of simulatable gates and an
 /// expansion in the JordanWigner basis.

library/chemistry/src/JordanWigner/JordanWignerOptimizedBlockEncoding.qs

@@ -19,22 +19,15 @@ export
 import JordanWigner.OptimizedBEOperator.JordanWignerSelect;
 import JordanWigner.OptimizedBEOperator.JordanWignerSelectQubitCount;
 import JordanWigner.OptimizedBEOperator.JordanWignerSelectQubitManager;
-import JordanWigner.StatePreparation.PrepareArbitraryStateD;
 import JordanWigner.Utils.JWOptimizedHTerms;
 import JordanWigner.Utils.MultiplexOperationsFromGenerator;
 import JordanWigner.Utils.RangeAsIntArray;
-import Std.Arithmetic.MAJ;
-import Std.Arrays.Head;
-import Std.Arrays.IndexRange;
-import Std.Arrays.Mapped;
-import Std.Arrays.Most;
-import Std.Arrays.Partitioned;
-import Std.Arrays.Rest;
-import Std.Arrays.Subarray;
-import Std.Arrays.Tail;
+import Std.Arrays.*;
+import Std.Math.*;
 import Std.Convert.IntAsDouble;
+import Std.Arithmetic.ApplyIfGreaterLE;
+import Std.StatePreparation.PreparePureStateD;
 import Std.Diagnostics.Fact;
-import Std.Math.*;
 import Utils.GeneratorIndex;
 import Utils.GeneratorSystem;
 import Utils.HTermToGenIdx;
@@ -623,7 +616,7 @@ operation PrepareQuantumROMState(
     MultiplexOperationsFromGenerator(unitaryGenerator, indexRegister, (keepCoeffRegister, altIndexRegister, dataRegister, altDataRegister));
 
     // Perform comparison
-    CompareUsingRippleCarry(uniformKeepCoeffRegister, keepCoeffRegister, flagQubit);
+    ApplyIfGreaterLE(X, uniformKeepCoeffRegister, keepCoeffRegister, flagQubit);
 
     let indexRegisterSize = Length(indexRegister);
 
@@ -828,63 +821,14 @@ operation PrepareUniformSuperpositionOracle(nIndices : Int, nQubits : Int, idxFl
         within {
             ApplyXorInPlace(nIndices - 1, compareQubits);
         } apply {
-            CompareUsingRippleCarry(targetQubits, compareQubits, auxillaryQubits[0]);
+            ApplyIfGreaterLE(X, targetQubits, compareQubits, auxillaryQubits[0]);
             X(auxillaryQubits[0]);
         }
     }
     Ry(2.0 * theta, auxillaryQubits[1]);
     (Controlled X)(auxillaryQubits, flagQubit);
 }
 
-/// # Summary
-/// This operation tests if an integer represented by a register of qubits
-/// is greater than another integer, applying an XOR of the result onto an
-/// output qubit.
-///
-/// # Description
-/// Given two integers `x` and `y` stored in equal-size qubit registers,
-/// this operation checks if they satisfy `x > y`. If true, 1 is
-/// XORed into an output qubit. Otherwise, 0 is XORed into an output qubit.
-/// In other words, this operation can be represented by the unitary
-/// $$
-/// \begin{align}
-///     U\ket{x}\ket{y}\ket{z} = \ket{x}\ket{y}\ket{z\oplus (x>y)}.
-/// \end{align}
-/// $$
-///
-/// # Input
-/// ## x
-/// First number to be compared stored in `LittleEndian` format in a qubit register.
-/// ## y
-/// Second number to be compared stored in `LittleEndian` format in a qubit register.
-/// ## output
-/// Qubit that stores the result of the comparison $x>y$.
-///
-/// # References
-/// - [A new quantum ripple-carry addition circuit](https://arxiv.org/abs/quant-ph/0410184)
-///   Steven A. Cuccaro, Thomas G. Draper, Samuel A. Kutin, David Petrie Moulton
-operation CompareUsingRippleCarry(x : Qubit[], y : Qubit[], output : Qubit) : Unit is Adj + Ctl {
-    if (Length(x) != Length(y)) {
-        fail "Size of integer registers must be equal.";
-    }
-
-    use auxiliary = Qubit();
-    within {
-        let nQubitsX = Length(x);
-
-        // Take 2's complement
-        ApplyToEachCA(X, x + [auxiliary]);
-
-        MAJ(y[0], auxiliary, x[0]);
-        for i in IndexRange(Most(x)) {
-            MAJ(Most(x)[i], Rest(y)[i], Rest(x)[i])
-        }
-    } apply {
-        X(output);
-        CNOT(Tail(x), output);
-    }
-}
-
 // Classical processing
 // This discretizes the coefficients such that
 // |coefficient[i] * oneNorm - discretizedCoefficient[i] * discretizedOneNorm| * nCoeffs <= 2^{1-bitsPrecision}.
@@ -893,7 +837,7 @@ function QuantumROMDiscretization(bitsPrecision : Int, coefficients : Double[])
     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");
+    Fact(oneNorm >= 0.0, "State must have at least one coefficient > 0");
 
     let barHeight = 2^bitsPrecision - 1;
 
@@ -1103,7 +1047,7 @@ operation WriteQuantumROMBitString(idx : Int, keepCoeff : Int[], altIndex : Int[
 /// 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 => PrepareArbitraryStateD(coeff, qs)), multiplexer);
+    return PauliBlockEncodingInner(generatorSystem, coeff -> (qs => PreparePureStateD(coeff, Reversed(qs))), multiplexer);
 }
 
 /// # Summary

library/chemistry/src/JordanWigner/JordanWignerVQE.qs

@@ -92,7 +92,7 @@ operation EstimateFrequency(preparation : (Qubit[] => Unit), measurement : (Qubi
 
         // NB: We absolutely must reset here, since preparation()
         //     and measurement() can each use randomness internally.
-        ApplyToEach(Reset, register);
+        ResetAll(register);
     }
 
     return IntAsDouble(nUp) / IntAsDouble(nMeasurements);