Chemistry library, SPSA sample, and notebook (#2105)

Converting the old chemistry library and Variational Quantum Algorithms
sample to modern QDK. As part of this, I tried where possible to update
to our newer Q# programming patterns. I've tried to make sure we export
on the reasonable part of the APIs as public facing, but if there is
something that is not exported and should be or something that is
exported and shouldn't be we can adjust.

For the original library code, see
https://github.com/microsoft/QuantumLibraries/tree/main/Chemistry/src/Runtime.
For the original sample notebook, see
https://github.com/microsoft/Quantum/blob/main/samples/algorithms/variational-algorithms/Variational%20Quantum%20Algorithms.ipynb

---------

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

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)

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"
+  ]
+}

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, _)));
+}

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;
+}