From 1dd53d6c93c2913f924647590ab7d686b61d00ba Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Tue, 29 Oct 2024 17:27:46 +0100 Subject: [PATCH 01/13] Updated transformfermionic to pauli to work with reduced and non reduced Hamiltonians --- src/qrisp/algorithms/vqe/problems/electronic_structure.py | 4 +--- src/qrisp/operators/fermionic/fermionic_hamiltonian.py | 8 +++++++- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/qrisp/algorithms/vqe/problems/electronic_structure.py b/src/qrisp/algorithms/vqe/problems/electronic_structure.py index ccd47fc4..e1484ba8 100644 --- a/src/qrisp/algorithms/vqe/problems/electronic_structure.py +++ b/src/qrisp/algorithms/vqe/problems/electronic_structure.py @@ -506,15 +506,13 @@ def create_electronic_hamiltonian_fermionic(arg, active_orb=None, active_elec=No for i in range(K): for j in range(K): if F[I+i][I+j]!=0: - H += F[I+i][I+j]*c(i)*a(j) # cre(i,K,mapping_type)*ann(j,K,mapping_type) + H += F[I+i][I+j]*c(i)*a(j) for i in range(K): for j in range(K): for k in range(K): for l in range(K): if two_int[I+i][I+j][I+k][I+l]!=0 and i!=j and k!=l: - #h = #cre2(i,j,K,mapping_type)*ann2(k,l,K,mapping_type) - #h *= (0.5*two_int[I+i][I+j][I+k][I+l]) H += (0.5*two_int[I+i][I+j][I+k][I+l])*c(i)*c(j)*a(k)*a(l) # apply threshold diff --git a/src/qrisp/operators/fermionic/fermionic_hamiltonian.py b/src/qrisp/operators/fermionic/fermionic_hamiltonian.py index 82a7923e..fb09364b 100644 --- a/src/qrisp/operators/fermionic/fermionic_hamiltonian.py +++ b/src/qrisp/operators/fermionic/fermionic_hamiltonian.py @@ -475,11 +475,17 @@ def to_pauli_hamiltonian(self, mapping='jordan_wigner'): H = 0 for term,coeff in self.terms_dict.items(): - h = coeff + + h = coeff/2 for ladder in term.ladder_list[::-1]: h *= jordan_wigner(ladder) H += h + h = coeff/2 + for ladder in term.dagger().ladder_list[::-1]: + h *= jordan_wigner(ladder) + H += h + jordan_wigner.cache_clear() return H From e90df3821739e5c751d1a4646b11f593c2bfc74d Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Tue, 29 Oct 2024 17:30:48 +0100 Subject: [PATCH 02/13] Minor update: printing of Pauli terms --- src/qrisp/operators/pauli/bound_pauli_term.py | 8 ++++++++ src/qrisp/operators/pauli/pauli_term.py | 3 +++ 2 files changed, 11 insertions(+) diff --git a/src/qrisp/operators/pauli/bound_pauli_term.py b/src/qrisp/operators/pauli/bound_pauli_term.py index dfad0709..10373d31 100644 --- a/src/qrisp/operators/pauli/bound_pauli_term.py +++ b/src/qrisp/operators/pauli/bound_pauli_term.py @@ -52,6 +52,14 @@ def copy(self): # Printing # + def __str__(self): + # Convert the sympy expression to a string and return it + expr = self.to_expr() + return str(expr) + + def __repr__(self): + return str(self) + def to_expr(self): """ Returns a SymPy expression representing the BoundPauliTerm. diff --git a/src/qrisp/operators/pauli/pauli_term.py b/src/qrisp/operators/pauli/pauli_term.py index a19ee24d..f3108701 100644 --- a/src/qrisp/operators/pauli/pauli_term.py +++ b/src/qrisp/operators/pauli/pauli_term.py @@ -57,6 +57,9 @@ def __str__(self): # Convert the sympy expression to a string and return it expr = self.to_expr() return str(expr) + + def __repr__(self): + return str(self) def to_expr(self): """ From b93cf986e5b463bb65af310dc77f6d079a1ba4d0 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Tue, 29 Oct 2024 17:57:04 +0100 Subject: [PATCH 03/13] Added documentation and test for fermionic to pauli transform --- .../fermionic/fermionic_hamiltonian.py | 15 +++++++ .../test_fermionic_to_pauli.py | 41 +++++++++++++++++++ 2 files changed, 56 insertions(+) create mode 100644 tests/hamiltonian_tests/test_fermionic_to_pauli.py diff --git a/src/qrisp/operators/fermionic/fermionic_hamiltonian.py b/src/qrisp/operators/fermionic/fermionic_hamiltonian.py index fb09364b..ccb4d26c 100644 --- a/src/qrisp/operators/fermionic/fermionic_hamiltonian.py +++ b/src/qrisp/operators/fermionic/fermionic_hamiltonian.py @@ -472,6 +472,21 @@ def ground_state_energy(self): # def to_pauli_hamiltonian(self, mapping='jordan_wigner'): + """ + Transforms the fermionic Hamiltonian to a :ref:`PauliHamiltonian`. + + Parameters + ---------- + mapping : str, optional + The mapping to transform the Hamiltonian. Available is ``jordan_wigner``. + The default is ``jordan_wigner``. + + Returns + ------- + H : :ref:`PauliHamiltonian`` + The resulting Pauli Hamiltonian. + + """ H = 0 for term,coeff in self.terms_dict.items(): diff --git a/tests/hamiltonian_tests/test_fermionic_to_pauli.py b/tests/hamiltonian_tests/test_fermionic_to_pauli.py new file mode 100644 index 00000000..fad612ed --- /dev/null +++ b/tests/hamiltonian_tests/test_fermionic_to_pauli.py @@ -0,0 +1,41 @@ +""" +\******************************************************************************** +* Copyright (c) 2024 the Qrisp authors +* +* This program and the accompanying materials are made available under the +* terms of the Eclipse Public License 2.0 which is available at +* http://www.eclipse.org/legal/epl-2.0. +* +* This Source Code may also be made available under the following Secondary +* Licenses when the conditions for such availability set forth in the Eclipse +* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 +* with the GNU Classpath Exception which is +* available at https://www.gnu.org/software/classpath/license.html. +* +* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 +********************************************************************************/ +""" + +from qrisp.operators.fermionic import a, c +from qrisp.operators.pauli import X,Y,Z + +def test_fermionic_to_pauli(): + + # Check if transformation works for both, reduced and non-reduced FermionicHamiltonians + + H = c(0)*c(1)*a(3)*a(2) + c(2)*c(3)*a(1)*a(0) + + G1 = H.to_pauli_hamiltonian() + + H.reduce() + + G2 = H.to_pauli_hamiltonian() + + assert str(G1-G2)=='0' + + K = (1/8)*(X(0)*X(1)*X(2)*X(3) - X(0)*X(1)*Y(2)*Y(3) \ + + X(0)*Y(1)*X(2)*Y(3) + X(0)*Y(1)*Y(2)*X(3) \ + + Y(0)*X(1)*X(2)*Y(3) + Y(0)*X(1)*Y(2)*X(3) \ + - Y(0)*Y(1)*X(2)*X(3) + Y(0)*Y(1)*Y(2)*Y(3)) + + assert str(G1-K)=='0' \ No newline at end of file From d6cc10070a75ce8b5a8eae63499d6be213b9bd5e Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Tue, 29 Oct 2024 19:01:56 +0100 Subject: [PATCH 04/13] Updates electronic structure: create_electronic_hamiltonian returns FermionicHamiltonian --- .../vqe/problems/electronic_structure.py | 219 +----------------- .../fermionic/fermionic_hamiltonian.py | 4 +- .../operators/fermionic/transformations.py | 3 +- 3 files changed, 17 insertions(+), 209 deletions(-) diff --git a/src/qrisp/algorithms/vqe/problems/electronic_structure.py b/src/qrisp/algorithms/vqe/problems/electronic_structure.py index e1484ba8..b05502d3 100644 --- a/src/qrisp/algorithms/vqe/problems/electronic_structure.py +++ b/src/qrisp/algorithms/vqe/problems/electronic_structure.py @@ -204,139 +204,7 @@ def apply_threshold(matrix, threshold): return data -# -# Fermion to qubit mappings -# - -# Jordan-Wigner annihilation operaror -@cache -def a_jw(j): - #return PauliHamiltonian({tuple([(i,"Z") for i in range(j)]+[(j,"X")]):0.5,tuple([(i,"Z") for i in range(j)]+[(j,"Y")]):0.5j}) - d1={i:'Z' for i in range(j)} - d1[j]='X' - d2={i:'Z' for i in range(j)} - d2[j]='Y' - return PauliHamiltonian({PauliTerm(d1):0.5,PauliTerm(d2):0.5j}) - -# Jordan-Wigner creation operator -@cache -def c_jw(j): - #return PauliHamiltonian({tuple([(i,"Z") for i in range(j)]+[(j,"X")]):0.5,tuple([(i,"Z") for i in range(j)]+[(j,"Y")]):-0.5j}) - d1={i:'Z' for i in range(j)} - d1[j]='X' - d2={i:'Z' for i in range(j)} - d2[j]='Y' - return PauliHamiltonian({PauliTerm(d1):0.5,PauliTerm(d2):-0.5j}) - -# Parity annihilation operator -@cache -def a_par(j,M): - if j>0: - d1={i:'X' for i in range(j,M)} - d1[j-1]='Z' - d2={i:'X' for i in range(j+1,M)} - d2[j]='Y' - return PauliHamiltonian({PauliTerm(d1):0.5,PauliTerm(d2):0.5j}) - #return PauliHamiltonian({tuple([(j-1,"Z"),(j,"X")]+[(i,"X") for i in range(j+1,M)]):0.5,tuple([(j,"Y")]+[(i,"X") for i in range(j+1,M)]):0.5j}) - else: - d1={i:'X' for i in range(j,M)} - #d1[j-1]='Z' - d2={i:'X' for i in range(j+1,M)} - d2[j]='Y' - return PauliHamiltonian({PauliTerm(d1):0.5,PauliTerm(d2):0.5j}) - #return PauliHamiltonian({tuple([(j,"X")]+[(i,"X") for i in range(j+1,M)]):0.5,tuple([(j,"Y")]+[(i,"X") for i in range(j+1,M)]):0.5j}) - -# Parity creation operator -@cache -def c_par(j,M): - if j>0: - d1={i:'X' for i in range(j,M)} - d1[j-1]='Z' - d2={i:'X' for i in range(j+1,M)} - d2[j]='Y' - return PauliHamiltonian({PauliTerm(d1):0.5,PauliTerm(d2):-0.5j}) - #return PauliHamiltonian({tuple([(j-1,"Z"),(j,"X")]+[(i,"X") for i in range(j+1,M)]):0.5,tuple([(j,"Y")]+[(i,"X") for i in range(j+1,M)]):-0.5j},0) - else: - d1={i:'X' for i in range(j,M)} - #d1[j-1]='Z' - d2={i:'X' for i in range(j+1,M)} - d2[j]='Y' - return PauliHamiltonian({PauliTerm(d1):0.5,PauliTerm(d2):-0.5j}) - #return PauliHamiltonian({tuple([(j,"X")]+[(i,"X") for i in range(j+1,M)]):0.5,tuple([(j,"Y")]+[(i,"X") for i in range(j+1,M)]):-0.5j},0) - -@cache -def ann(i,M,mapping_type): - """ - Returns the qubit operator for the fermionic annihilation operator $a_i$. - - Parameters - ---------- - i : int - The index of the annihilation operator $a_i$. - M: int - The number of fermions. - mapping_type : str - The mapping type. Available are ``jordan_wigner``, ``parity``. - - Returns - ------- - PauliHamiltonian - The qubit PauliHamiltonian for the annihilation operator $a_i$. - - """ - - if mapping_type=='jordan_wigner': - return a_jw(i) - if mapping_type=='parity': - return a_par(i,M) - -@cache -def cre(i,M,mapping_type): - """ - Returns the qubit operator for the fermionic creation operator $a_i^{\dagger}$. - - Parameters - ---------- - i : int - The index of the annihilation operator $a_i^{\dagger}$. - M: int - The number of fermions. - mapping_type : str - The mapping type. Available are ``jordan_wigner``, ``parity``. - - Returns - ------- - PauliHamiltonian - The qubit PauliHamiltonian for the annihilation operator $a_i$. - - """ - - if mapping_type=='jordan_wigner': - return c_jw(i) - if mapping_type=='parity': - return c_par(i,M) - -@cache -def cre2(i,j,M,mapping_type): - - if mapping_type=='jordan_wigner': - return c_jw(i)*c_jw(j) - if mapping_type=='parity': - return c_par(i,M)*c_par(j,M) - -@cache -def ann2(i,j,M,mapping_type): - - if mapping_type=='jordan_wigner': - return a_jw(i)*a_jw(j) - if mapping_type=='parity': - return a_par(i,M)*a_par(j,M) - -# -# Hamiltonian -# - -def create_electronic_hamiltonian(arg, active_orb=None, active_elec=None, mapping_type='jordan_wigner', threshold=1e-4): +def create_electronic_hamiltonian(arg, active_orb=None, active_elec=None): """ Creates the qubit Hamiltonian for an electronic structure problem. If an Active Space (AS) is specified, the Hamiltonian is calculated following this `paper `_. @@ -360,20 +228,17 @@ def create_electronic_hamiltonian(arg, active_orb=None, active_elec=None, mappin The number of active spin orbitals. active_elec : int, optional The number of active electrons. - mapping_type : string, optinal - The mapping from the fermionic Hamiltonian to the qubit Hamiltonian. Available are ``jordan_wigner``, ``parity``. - The default is ``jordan_wigner``. - threshold : float, optional - The threshold for the absolute value of the coefficients of Pauli products in the quantum Hamiltonian. The default is 1e-4. Returns ------- - H : Hamiltonian - The quantum :ref:`Hamiltonian`. + H : :ref:`FermionicHamiltonian` + The fermionic Hamiltonian. Examples -------- + We calucalte the fermionic Hamiltonian for the Hydrogen molecule, and transform it to a Pauli Hamiltonian via Jordan-Wigner transform. + :: from pyscf import gto @@ -384,7 +249,7 @@ def create_electronic_hamiltonian(arg, active_orb=None, active_elec=None, mappin basis = 'sto-3g') H = create_electronic_hamiltonian(mol) - H + H.to_pauli_hamiltonian() Yields: @@ -397,71 +262,9 @@ def create_electronic_hamiltonian(arg, active_orb=None, active_elec=None, mappin &+0.165927850337703 Z_1Z_2 + 0.120625234833904 Z_1Z_3 - 0.223431536908133 Z_2 + 0.174412876122615Z Z_2Z_3 - 0.223431536908133 Z_3 """ - import pyscf - if isinstance(arg,pyscf.gto.Mole): - data = electronic_data(arg) - elif isinstance(arg,dict): - data = arg - if not verify_symmetries(data['two_int']): - raise Warning("Failed to verify symmetries for two-electron integrals") - else: - raise TypeError("Cannot create electronic Hamiltonian from type "+str(type(arg))) - - one_int = data['one_int'] - two_int = data['two_int'] - M = data['num_orb'] - N = data['num_elec'] - K = active_orb - L = active_elec - - if K is None or L is None: - K = M - L = N - - if L>N or K>M or KM: - raise Exception("Invalid number of active electrons or orbitals") - - # number of inactive electrons - I = N-L - - # inactive Fock operator - F = one_int.copy() - for p in range(M): - for q in range(M): - for i in range(I): - #F[p][q] += (two_int[i][p][i][q]-two_int[i][q][p][i]) - F[p][q] += (two_int[i][p][q][i]-two_int[i][q][i][p]) - - # inactive energy - E = 0 - for j in range(I): - E += (one_int[j][j]+F[j][j])/2 - - # Hamiltonian - H=E - for i in range(K): - for j in range(K): - if F[I+i][I+j]!=0: - H += F[I+i][I+j]*cre(i,K,mapping_type)*ann(j,K,mapping_type) - - for i in range(K): - for j in range(K): - for k in range(K): - for l in range(K): - if two_int[I+i][I+j][I+k][I+l]!=0 and i!=j and k!=l: - h = cre2(i,j,K,mapping_type)*ann2(k,l,K,mapping_type) - h *= (0.5*two_int[I+i][I+j][I+k][I+l]) - H += h - # apply threshold - H.apply_threshold(threshold) - return H - -def create_electronic_hamiltonian_fermionic(arg, active_orb=None, active_elec=None, mapping_type='jordan_wigner', threshold=1e-4): - """ - - """ import pyscf + if isinstance(arg,pyscf.gto.Mole): data = electronic_data(arg) elif isinstance(arg,dict): @@ -515,8 +318,6 @@ def create_electronic_hamiltonian_fermionic(arg, active_orb=None, active_elec=No if two_int[I+i][I+j][I+k][I+l]!=0 and i!=j and k!=l: H += (0.5*two_int[I+i][I+j][I+k][I+l])*c(i)*c(j)*a(k)*a(l) - # apply threshold - #H.apply_threshold(threshold) H.reduce() return H @@ -771,4 +572,8 @@ def electronic_structure_problem(arg, active_orb=None, active_elec=None, mapping ansatz, num_params = create_QCCSD_ansatz(K,L) - return VQEProblem(create_electronic_hamiltonian(data,K,L,mapping_type=mapping_type,threshold=threshold), ansatz, num_params, init_function=create_hartree_fock_init_function(K,L,mapping_type)) \ No newline at end of file + fermionic_hamiltonian = create_electronic_hamiltonian(data,K,L) + hamiltonian = fermionic_hamiltonian.to_pauli_hamiltonian(mapping_type=mapping_type,num_qubits=K) + hamiltonian.apply_threshold(threshold) + + return VQEProblem(hamiltonian, ansatz, num_params, init_function=create_hartree_fock_init_function(K,L,mapping_type)) \ No newline at end of file diff --git a/src/qrisp/operators/fermionic/fermionic_hamiltonian.py b/src/qrisp/operators/fermionic/fermionic_hamiltonian.py index ccb4d26c..6377ed74 100644 --- a/src/qrisp/operators/fermionic/fermionic_hamiltonian.py +++ b/src/qrisp/operators/fermionic/fermionic_hamiltonian.py @@ -471,7 +471,7 @@ def ground_state_energy(self): # Transformations # - def to_pauli_hamiltonian(self, mapping='jordan_wigner'): + def to_pauli_hamiltonian(self, mapping_type='jordan_wigner', num_qubits=None): """ Transforms the fermionic Hamiltonian to a :ref:`PauliHamiltonian`. @@ -480,6 +480,8 @@ def to_pauli_hamiltonian(self, mapping='jordan_wigner'): mapping : str, optional The mapping to transform the Hamiltonian. Available is ``jordan_wigner``. The default is ``jordan_wigner``. + num_qubits : int, optional + The number of qubits. This information is necessary for, e.g., parity transform. Returns ------- diff --git a/src/qrisp/operators/fermionic/transformations.py b/src/qrisp/operators/fermionic/transformations.py index 8225f30e..88236d2d 100644 --- a/src/qrisp/operators/fermionic/transformations.py +++ b/src/qrisp/operators/fermionic/transformations.py @@ -44,4 +44,5 @@ def jordan_wigner(ladder): if ladder[1]: return c_jw(ladder[0]) else: - return a_jw(ladder[0]) \ No newline at end of file + return a_jw(ladder[0]) + \ No newline at end of file From 5d9a876497766911311d69097793461d693f665e Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Wed, 30 Oct 2024 13:00:50 +0100 Subject: [PATCH 05/13] Examples Molecular Potential Energy Curves, Ground State Energy with QPE --- documentation/source/_static/BeH2_PEC_4_2.png | Bin 0 -> 27186 bytes documentation/source/_static/BeH2_PEC_6_2.png | Bin 0 -> 27582 bytes .../Examples/GroundStateEnergyQPE.rst | 54 +++++++++++ .../MolecularPotentialEnergyCurve.rst | 91 ++++++++++++++++-- .../source/reference/Examples/index.rst | 9 +- 5 files changed, 144 insertions(+), 10 deletions(-) create mode 100644 documentation/source/_static/BeH2_PEC_4_2.png create mode 100644 documentation/source/_static/BeH2_PEC_6_2.png create mode 100644 documentation/source/reference/Examples/GroundStateEnergyQPE.rst diff --git a/documentation/source/_static/BeH2_PEC_4_2.png b/documentation/source/_static/BeH2_PEC_4_2.png new file mode 100644 index 0000000000000000000000000000000000000000..f50f830cae1e950c1b9dfed38f162655aa6fa815 GIT binary patch literal 27186 zcmbSz2Rzq%`~J7IM0REoi9!h>l7j@%@bVxbN$}uIqmLpHMxvX(Q`K5{a})Q9)LNL|Sb~BCRUi zunzyyxkVua|B-N%JLRZpYwqZL_M#a{<*cLK1zX1p=Z&{rGP~$--quEB&;C981h$=X zbhLAj6cn`n>lgOeUbGP8*s;bJ37hWcCEnI2&K$!rD>}uaobO~3KGpy znKLc#-}g%?7z{YKCurM-6oDt3_Ld!5sIS~7vX`Qc^5xi^%j64*={b*vx#GE?Vu|=^ ziDTPtzc$^*#l@8%Z(%mNj)s<&R+LG(*bje#uK6?JuVG<(r|@$wVLBc|S^SZ&s0_qU zUn$&m!B1&f6e%>t?@OeorXIO;Nh~HdmXl_Ka?cigpy`S-kBq&&kcWrI8MCC*9E|jW z^tBBQz9)9>a}vQ76sYO^*Nk`NHJ|6;G35F1&YRaIrC(56(afS*5IWaZ^6 zs;jG8>V?>X0?mVh-TwMITa}!Wvf=216I1T`cwS`%FZaf(sw7=rLn>MmwxPnP6jUQX9fq>4Dxz#69bo*2gB)^qkZdG)g+91mVP*vMEz zpLH*VxF>f_<(mrAh&DsXwOo|8GS8@!_1im17e< zMR&z**=4zCe=icd;dc+*ez^afmRPe;e8%+cEV`l{oCCD+8PyI%D6vX{8c zR^Wr1q{#3FPaGL8Ck|{RWs9`OsYhbi zVCvJ26)l_gM)VXpZx>>3MTn9KT?wsV6^Ifo{(^K zst{|Pv@Jr3GeXJfT*=4KcxQ)$-YU#Y^~`Lue|;Sd4CWNFy5p~-tIN{sHnw}=^P%M< zzrK$C>zdVRiQsvXz`0Ra487o_tCM+xr zJGHT3_npIh5nP;{tBiA)#l*ytZEkj+(cfosp~d9Sujztf-wtTc?|oCrFCoFg6{e7G zovmMR<^KI`nIWEIZM!|bz2^Nn7%15P;X|Us`#6_D@-H1ZWZ9+jeCsx`t|3i+|6Xko z@9Eoa=$_eN($(2n-q4_$)LmVj+-WoubooD`%I^phZET`Yj5po5abv}u-KTTAo>{+t zCb`ueD{XRiR?fvGe~{5cRh5dgO;C{fxlPZN_;`+tjD3fk_j)b3NKMy?C>%d7Zi)Z#U&*|)?K;s5|_-)d9sYkZc=TPOeq$3_@qpFadkcQ zyR?((k|8EmR)JV?e1;tR?MS4(7djsMS_u%VjsaI-ZEi|B9u^;;t#xT?P)^{~W75~+ zX{vRU%%o$-j*&vvzcn*AuXy>AjGdjGHNllF}b9jjk>GHF2P1 zloR{)Q@!yWN+Na%JOA$@v4Eu{Iq~}v{<2rDTv?r8R3x*LA-Yh4v=WfQqMCUx2ZV&^ zWEhf385Zye3qQ$;P7z(nYo^a;6z_A>`8P?ELT}x=(%*mP_74yHAR+^bO5$pm#Uor{ z!?_`aYtl+h7#J9awQQ8X5ScMSKIyA9S~%ljN8QfcB2Pq(a@@x0YbO81uxABmA`nVNE6xq5X@Yddk* zgY#6U#AlwxbCQQ6W!}4X{N{5azDGaE$Ph^SI@{>{!cC_8JN?G%!;KGbbnZ}}KNEtd zRTrV>;0S4H(VF?vNq_peP1TFbWMd^}=M>&X6!d-)a&j`?c#Uc8{*@sW zJ>r{jjG2YSzajR(_T9TTH7BdcU`es%Bx$j&Z>Ps=9P!OMHas?#J`)>sl>ij#^rj7C zw71MxCbZ5c%W7NnwA-~rmKtKrqV7sx{jf^)>pKP^VrYdqW!AkWi|KukIW6zwD|5$FdE{F?#qSF+uYWDhvGNP7p;=bEwf3PN(=Q0bl5%B##I{~M_b4_ zIXU-V`f@7v{d22!syC5>I(*YpQ;lg^o7lwqTXl@89Nr$hA>r_;|5ey_m+5+Obv?a+ z0*B9|o>X`=M|-^uO}3mA^k@lS58S!JHhfP~sd%yel zjOxmCB(%+&HwU!K$;pvOsahLTwJcwa?9=@=ZIksmoBXfbQkwI+wzgvN*AFDqbw)-; znLnn#fA@X<{JEhMD}6Rg{a=Z#=H^yTI)55Ex+}***g3_-7>`Hpk&~BScWI(q_j9T9 z*WT-Sd7_9b54CKptSFqFolPw)th$}&zI{CP#MEEO)x$$-%a$!O2UsY=LPO7t^4#6= z7dZ4Z*&e14dgo5ti%cPoV)K2?4ELIvG)h9wBcH0MY~thR&+BSyXrNks`S7D5+pb*2 zky#4CVxQsRa~HbvX_x2TdcQ04_747aCnBP}txbE}efQ~Svbwra)t;Z<#)n|4cN#Kn z-h5;y1DCKc-JU&rBxA28g(@W}C@A=x63BY~JV)Z}pUJLSR7_V_x1-Q0U1=`&e6!=| zoUk*V=)cZ^SI2_<{mZL^+3fpYT$ZF4ba!_TRy>Od{xUu;I$>pDvHrk;0}&V5*zTEo zH;Swq71|RmBu0Ea&vx$yHAg~nvb?n`v1uQeC0AL=i>TxXaT}UPxnNnGeRF5G9ubPNoF-G5B%>aoel$dEm9tHTT%U$HnZ%DUz@F2a< z!EE9rQe0e|lTG`>W7MPEkK`Lt)-=vc=3%I z&+W4#O-Ur5ftNSjFA57?8%a)zn zW67~v8Hr0(CC~}FesdgY+Ch2&z_Zyq=XY)sv1mBcl&E-j-+8*#WMt~v+P**~TMoK0 zV80BvrbiCC{A;o!bCIeZdr5t$lGqi{oAjsr3g2)?Jd@wvUMH&`1NZ@`PKE?_w1pPkdS!%ouSj`v1cRO!9^%cw5yV9|i8FA}4`O}QP0!ciH$)~i>qer(mgC-Wu~ zJA2mDT@mv-V!Jz9a-bwI($>B+ShPbrLO~@kB9#qWzU9>5(9mXwi`G4bn*r*Qt!YHv ziGcWuOOcx{iI;xWwjDdx()m})?&dNk${fbD`}Yq&y2xd6X60}G^bc-H{QiNWG_9p| z%!UZ%`h+-3ZZ590-B`S_2b|^g^jMQs?o~fxef^<>2&ato?d>`v2f}u(MD2es75C?( zPKjRgS5lsSEo=FN$RA;yfk?;E*QEsnC;pNKPSXlG6?R)uI%E}=e-<$-IY7QD$vZ`2<@e=Z;@g6i zI!Soro<*5TB4GqwJHAY_@s4eR^CCK__kr?vHzni z{_~dqMd`zt_ulu$4Vi5;DIUShLPygJ+cHnz_zkR%Y@<1YAhyYw3mCyPN5I@}@qdw7 zP!!$I&u3nopM8nKV%A24fDnQxm}U`@v6Jl_9AqwB5J*T!813$DZQZPytlXHM_3{Qa zzgrO!!HncOT6N8usw{s;3q<>;PoFmY(mHwaL2oYMXLmzG zb^vF#5-)?PX!Q{L#*ESqW;&6}$~?D6N}RW+f_yMDOVWN1Xphm`8S%QtWyg*koWPUm zL$b$?nG972BONo(|M;5ibUHRBh5+gpkRp4nnFI}kb{xNx68;qh0}3y_c14VOvL+Jp z>)X5!A3vrcg7vu7or8w2U%#gDm({a&|5&?u&Dymr`WDK{lwvmBeyFU5!(B#^P!93Rw=g&A)rMoSU2Ex9RENE`mdNqHNwdl=LsrHAJdxrE>VcMpyn? z@41gmiAhPTTs$~_F3e2_6vn%LVv(Avq!0GW!W^{Bpf1QeQ@%Q>{hqk3_QL5nm$_fR zOeXuw{Q?7JP<|o#&}`WfG&|P5T~bmn-R=9A&a2sFOL|lDz#~8)0zNgzX@3VKCjq{x zCUf&3A0N_}uDs*vL)UNINI&=qdoC|pn&3!Mo|%hRuh3drE3o!ILPk!C5U_e8o4;&R zBF~W>mZQYaMKy_Pd6zS3+keEZk3{%LFjwFT+quJls3BJS%Ni#z(uz{A+hK_Kzp0f9 zyCVTrNo=iu%30=c7nNAa&C{RX)`V?8cEZ4*=Y%y;OMIG`$O73TN_n-~W$x$C{+W^H z{SKcf^$iSYfwF10%ilLBy0ocPnv2VF&Cv%mXV0GHGb!HxA}+(Wgm);UKHL0mJ@;A5mqqxHo+^2Qg-!Ekqa1+=5{dd(w=~%+b)kJ-LeYsGl zVb|=p>n9;?-pI|}7jwv)T~$?;Ca}t7H1$#c*ROpv<4N+>$KLoyGd(!*>aAoLt!RfM z9e#yb#GDe4@aoKrt4-??+c&6_uYxQRz5^`o%O_m~qdUHEofYI$KAS5-FAQgOR=NN>+>J2@wRvW!@z zu;10$k|%zRAn!U2nY3ziC z#+6-~$@@%aO?@y6;W1yYeG8%rcycG!_Wu3*J}b2e+4<1QF(9dqB4+_`ue~M8#dQjgllzFZO22$YCH`E`>pO4X~c=Lt=?AB3Z;{?u=nq>6h1*BK4 zCseWUHPP%?;sM;o zBF=T)L;2vNhY$Ox{sP$MN7wzYI#%Otk30G8xa z^LopaOQ9^xY|ZQkeG7(q$AxpJo?A4@wOlhbGuu`A+_tYbH<*(vm{Ko2OhD~$ipbKG zG8>P9tcn6~^fQ~DjU-@Heqari_uIIdUfpNUe%|MGx=2Yc*9UJZdFhMamjWJh4UU^`#CSmWB)~sG-Djd~iZCMrxmJJg(NH9qG_Cu7Z(SClD zE>zMgACYHH9upagU6>B|y z^f-4sK1;T+u<-7YTKqZ*axe+CJv{~hIl{X!Ui>;){WTokmB%9wD=Car-i1z$9Q{xpEN!rs z=Q^JxC1gB#LWgvPwLpZ*4~QnSsP{WNJ5N0|-&FhN&B@;El#~};c~)Vuv6+qsg!%a? zfW!a#^?mtX^D`qvkn~!d)#2W;^+1b(7{mv6xZr!@2^=zc`DQ~aGv0M?;t z2sZwOQQ8^pO{Vh`$4pgm0iPN&i9<_2O;DZh)5A}IquoBXOB_Eu(~?eHh0mlE*OhpuqqJcqA$QWy|lFS`ppzA0wqEE8yu(@7O8TCZZp61lM2E#uAiT`zL}YTR4wio|bs>5}i^Z7&|k+-R|3V`rzjeE9x1QR7NK zx3Mq&<}6fT12gGB0Q}9P=T9(N)vr$)Z*gLeJyVrW*L?Q5b#d(k?IC=^70&$ zn$hl3WPwwq9t9dvV$YSl=SO%mPUWmNsSYY{N!3(8dGa~{!w~^Slj%?O0guG}{bo4B z6i%EtL2`EfA#e2kSf9P-JkQZp4;pCHew<7b`DWJVX3I2M6klx7&^3`Q+?&9yDYA1a z8kJI0b2BYH{ne4?M-RQ1m(Dd-^iR#Zjr0_WgBgXJxaIxla@FJxF;e)59%@C ziH?!y=7VlCFDkB6M7hw0^&(Orb|<)OWn^qS_hy+{*hZ;Ob^{5MdJNMfNnP7+QC8^o zDYr9TO zezd!IvB@porl&A99J}D`hgV^z@+>z81P0DpQqYRp()r8o8d>cV9?r04{d#eww~-GN z{Y%LiCimMn^M*A#9g?wC|F%oXk`hy2P*A{UM$Im@_d}@s?Xa*_U{$rWwTXC56kd3n zvnD1F+ovvGyx8C0f8w)-0yqwkGbf7o@bRrCfeOd}N=B_mKYZ@de``b!is5*$3IJ?I zhSrUqPai*)zkgrUyXX`gSoCgP+a*!UH+cu;i?Vh%W?3DN9TiQ!WN`ZE!V__-uWm7i zb>F|HnMx4vaZfE&vWK;L#$-G^WwY{ezuK| zFRj6_1Pt$QE#SbVFa2E+h%73y=g;rrQHu+zs5pv^Gn|{P%VQ`iDw>j!(PDWTw8Dlp ztCOD0(+N5yJbajAKP129z$HpBIqij+2=7G7B*I4^lVF>~VIzZ($P+tdy}U{r4(&e< zDFY-7Dd$v%v+ctj#`^o7e(2!l;o&>Tyzk5l66ID&t#r@X;5XxeJ6^pJxcVr*L~iy| zOS07kVvkCZ>%RX|-?N)f{?3A~u$E`eETR5hY}@jtI&g zzs8^IEuU19)dJS)7;aWLu}JCa;-e;`VL+)`T>5j6VqszWZRl&$$IqW1J9%=mWnYP- zvvUo03Q0{(P3lLtgWlHW_wU)j$cp5imn`@L^=qb&_j&!Y#=Y0Zdwyr(QHgH3tFoPo zi~il`k>k>v6laCbeVRwTLgc*RmQ?z-EQ7#!uX#nYdk`v6HPF|DaEK?St0{GuA zAwdv5=UUSus>dZAhkbzaSv4ZGKfG>EIvz)vmkIb zP*Q?fV*$9dn;BL|dG}B!+c02gajJ#`Vfs|=Ia)?RLrM~=gqH`t7GDKfOQiEYkLh3y zwn}3mWG~iVTK;~o>X;?PSy#|GG@J^0H7RuqU^u{0`KUgt@+)hVvc+wCGk!jon)>$R zM@@hERaz#d8xZnN>FN>}9Ba=$=IH2H;zZOiq%4CXUnWDBBOezN6S)aex;~VX{+Y(~ z*-7T?GPRHQRY8+=J(^>U785a9vu4d16li)|3BlvzUF;6+882|(f$y1)GxY7A9Xa4& z@k&Ff#RFM=Zhn3cSv(~zttnsL1OvRZFs=6U*Dp4W&9t=VX?~3C=JwiN{w!5Y0cSILg?DgqjoTl7K1aU6srpdd`G=a{ z7;qgA;Ss6QIOk3A%90^SN1}#guM;Upr&z^pwqWM3fh*H5^z-J;^P4io~$=w6qzKpN#SH8jKf?)fCSP-{))dojf(@8#f z=fT0Mnm4;=Mq4S6q=qu`L8SSR72cEPd0QV?es%10p5=A0Yy${amB>C}w{Me0DX#+! zI0HTS$en4(L$e>q3hB$jcGm8{Uv&QLx9;M79eTlA4_uN3Ps}!Ihv*e8%O!NKUd|vW z_^ySXeq!oaUXKE0XGJ z`cg**nnUZ4UF9|e!mC1bR7VO;)k?d>2DzFjC1jmsWzCQORq+w3QKnnz&*QmsgI%HZ zSC9BdYVlP;(F?tMcfm4RilicB7?PtFVPoSH|#u%GNdf-Qj; zTMg`IXJ>aV1h&NM_J4XWW&??AK{{H?*cNM5<@})UxX-N-G_z}XLkm2m? zrN93kHHWlpytLaN4{?|z(D$ohI-688a+jrZcex+ry$h#_G}}aU%l-{-oa5K7Nfl@x2O1)b2e{8SnSJf#m`f-#a2H} zn_MeB-O|JK&M!HRPTdL_VHGmIR_l&MEDyzw(Ao*%qv*z=kNd(7RUS96B$E&{Jrk~bG5mjm_erZAC10ke|0=+r%MGom7lN}v- ziXW+IL{EPXB~(dEXr^OYdn|_Wvqa_iUSWr&Eu%3mJGi;GadI9>(%(jN&LhMxC6)HK zqRBFJvfWz##p5XfETXKH%bfUJc^tgl$*{xQ-dP>u;?ndF3X%o$O(4dJ2}9ea9gn`c zFmhS$$0N5`tv-o0b#y0#WvpFm@|e-EY|CWbl*226*ZrH8LRN3{V^N0iV%m+x>;2Cr^Zx0?u!9&Y))T%N2FC)k>JF?>eO$+d3 zqko6Pr)1^iS_?{W`()*)1c!q*Ha7dc7V9i5W?zOVALy$Z~Cp6X7F9!r{tQrH?NAl6) zvwk-8@;_jzN!6$CRdmQFm-Bjz%>DwIo0d&7A529Pq~+h3*FF~Y`nwC-CK?)aD7#Pg zBF;#n;QpI9hut7TcwbNnTAr-=(dQk41is7g$4K*fBzk!T1;1)fWa@vVUt5o3RFhBs z8~yq(Dx9jB_QJno&TQvj-QD{7`F5)i_bMU0gx$I0i?G#NkX^_wh#0xB;1TO>cjxGz zW|LE6z15$5>>KM&u{%OMe{SaiZgDc*#r}M@@I%|-3O6UZi7PsDw ziJ=TrNNjxu+s=^_Cqm1esw*qc5d5-TzmK+mB<2DDgG3st38?|SuqNkJkX3uu`7J|c z{~2K`dg%fnWi|44L|!kxv&ofYTMfnYx2SRo5`URZ_C0KK9je|$@_F>=k>Rh%h=^6BjT97Y87i_Zd;eCH z`M!FkkojG-Ca>2oW#q}-sIk10|BCTQ%+enD3B|d&3?m~Wuin1Z3IF){vrNz9d>Cyu zz8E_4UDi|@4hn(^+qG*K=nNi9toV=_|d+$T1|&4NPr$PKEkC-n7$ z`n;F7l9Q7&A9ULXZP2n(gePGh?u>H+?ceWgr&X89NvZt5fH^qGIFT0mJb!FO>S;`n zqkxdO78|yL&9=DlF|VxbItORxc9`>0gmkqAlj`f3SFKvrk!KaDT|+@lts`NMjgZ?L z@29K=tye%m*(`_&G|?*tFwk6F8oUtD5SAqmAt$uF1C^})bQ{CQC_CD)*P=x$_B}u4 z>;JG_-kbuhJYcRT1VxPNHuS1ZgDw?Lun77Xov@l=S`T}RtkR%#Nm6ddA0&ym!BxMM zPH1@RPXA}_LvCl0%u_@)d`(4nSg0=cppmSpBP6kP5^;7V+D)?^v$ChRsO8&$$J|xq zT%-)zECH)Lkn&d4(uz=MHDsF&j#}M>xW#XGVpz^^vQZ~h;G3CIuQ@l@?S5uKVFhpd zZ5lU}!fnJa1*Ci+xeN1~m~=g?+3PHg0pO{@HN`A~{nSKnby93_G6xbN{*6@@pPCO; z>Dg#W_2kJDhCUPmwq4aUv?G2d#eONT{X#-&K(wqP!QLY0;Zb~^4Ls;}fbJG<1#UVb zM@A$iCLZPryGuAK3kii!!{~~VuV5;A+Nj^3+Hn0S=6zSh_3PK6&_}O) zYDhC!w#5_%)C4$xVCZ@YCi5qZIieCU{Uq>Itq7g~kAy@S9S-CkUP@0Zx&t&J9`M%; zqxSL`7wODw*&eQm=F`w%X~S+3(|sEmo*9H}crBfn874KyUymShS028*$RHWR8n&Xo zsvKqWc85~nS4Rs+jwRjdqp_TESqt1IAKsL1*x)3l4z9wW+i`x&Y;=KnUa8lQ+Y@|xShO-F9%HD{J#SkiUh1jc1NR=ZLmdr{oE zz#;PXMx!m~60{-z-l*!fbOe2O@7_Jt6;(4aQ8s%X3|k1Y5(kIJ&YU7x;U2@i(RI&u+c?_H%7&{EQYc zb3hlTgA0c{a<;=W2YV_}a8*`U8>7mF2JiwNsj;4AXk7$5{e7Z`!ph1jbL^iw=xF$i zF5R4~c>8pF<*lKn1Kk`ol~MU&7ar~$+bQ zmoH_8yuOVu*`Sm{CGsvUmH<0?;B@Xvxe6nD3=9q?G?yZ@N30~oBHoy+--gV;lPxf# zWFz}?A)rEMZ7H}pH~rZJmHdh$h-Tw;zo4L&x%C#i5;b3Cums7*Lho*fxw*m#}9P92^7=DlrMFx2youC75uh84iH znuOH>-;Sv2@w5L~?)`0n`NIek5X?Cx{HF5FAAE^rthh}#iR9wu)>aTJ%T*2``>{z4 zRm2ZiI&IgrW&puM&f&3QkQ)QAql%_+J=TsyZ@JJx7>$EieA)C~)A#v-`pn zEBcLva06?asryGtxx%9RZ*I*0tENguf(G|cBfI#!&DumM`RN{;{+ z*V8J5(-|0iAA4KxsTEn$)AI6e+q7o=@s4TKHmhLq3p31$$-e>NPRR^%o88CCgY%8{0qi zFIg>We8e;qAK}qQ;Wqeh_W4(dsysa>=5Zs_*LPjQUS~u{T>JV-I_|=BQNfSLmfE-c zPh(_S@?m>XCUUjDTJH5cgW!ZOZV9{L`9J1uf1muce9KR8XKCr_pzB-D)PyS~ooh^3 z3+YB3TJV{gyfrIkJoeUIVQyE%y@>*FesF{A6(Ft*#XmYKFUFOk`{V3AeHRs z3340u#WA`s4C=4n2bty&FLUbD2GEc2YVG?uWo2o3LRZ%x9_nwA!+m=R&)KL1f@xbD zBZTIZtSsN)V3n)v@%~Ivb%^_k)kboLou0~ijqAv^7u_vXtOtLC9lr|d5GTV%m7Xnu z+u?z~{knhnm0YENUF_AXWKd+;9JeIHwJMxinn2Z2B)-*CZmV)nxcP2rX?Dd!28VyE z#++_qYhW-+wz(rWc#4&4ns+&-*7V7`Lt^xT3FfL?VP}@R^usM?bi&8q7GCZ@L+)M* zrh(aZlupz^N3~nV)%RkEiL3m2eF^)w#M*hHk0lmekY+qhN2Yk*i1+kkehs_0};G7yK;8kW{KeII=Z*6V`ek3^N;yP?7yN=~jMdVBqrE zm%dSjEbCI?%?o7Qy`PEkJb@@zxpYrgT?>p3fPs?`NEMWov2xgv_|}FgT!SBvuP1>a zu!_LZa847}0kC{8Yiiab&260jKA!Gahx!Sona#YaF)xt3`-g_Eq5v_e^rH>1Zfa_p zY&&f`K(M_KaVLL$EgOPY5DgR@@Z?{^QZqrw&PNIF;pr*lG^#~Z(@2CrA)-U#nXI4? zybmGwXa7}dtp>398)3y=UiwkSChBD>Nd%abaw)f&4FG6-AES8^P$6q{_JUZ1Ii}^* zF{C`kJ^FrdgA%&W?% zs0|&-C2*+r@HaY+wW+A@L9n!J<#U|qW`Xk-UhwsXt{*o;N-?gxD@5|aNIWboB(lC1 zw+yDO`ia?_$JXz~f>=d2Lh}eP!G~~{(i42|bI;SmpifA<?T|y0t|qBMi9Yz*WxyWIf9s3&=N0NHwnbRWKA%;#BdH)k`rtR;*x@r(oyty z)DqOvevhAzmF8eC$btIlHr6!r0HUh|XsAaH@1d&h2vZ2ZwTFb2b;Qo@v7$Ena8puX z*9bS5C+zh^jfv6G0K*cu5v^lrR}vNxX_ekYPD0GqavL8RNqrIw=mmvx^Luz7h(PcX z-s|L6AuHaCB)hM@L*zF+Jv_9&&rVI1!{iiE3-O{YL-#634yf#$5V;VSRhEvBrN!et zpn+6DC4QHx$xDcR=tgp~>B}-WAYfFgpQ_xmUY;IOThr!V4DUFbG7Rt=H)?%8b35|% z=q)vbc#M|IY^zRd>GBWRooQ@pJ4)AcKb0AO!Qn76cE;Ib#-+oo3D6(@Z`e=MqoNOb;bLoaJI1LXQ~wJ zIQRb9Zkz7>%W((Y#a>xUfEe4qfB$_rkw5o}7ohozM0)b{>B(MAolk>E}Q$*pg@x>#YNB7y_9%E{j4U)J(Y*aq*g6}RlSJDj=%BUE0^8Xtt!0XR^8 zetRQ$7G}ds=O}gA14sUHLaOWPk_!q7V$YFLFzvY(92~6u%NEfR?9X-dL4yD6^)wc5Hrn^V|9o0| zkQN_t9qgp8BcFVk#V#-tOujp$dsWSz(cW|`cNK6U0YudO!=q7% zwjI=7D37^$d50c`NFT+7a?sZryWfc)mDs$y{T=#|uaXYR}d2gcI8;(t@#TQah#n`^qUB`b=6 zO#t%*%*e>NFga9ISXcqqE2)Lb*!IFV!#Fe$q-K`EIsh?S=_`SJ($bW&zb~O@t{l72 z7h^YsSPCl&(ZUwzHo6P)_s*@F?6z}f++2Im2N*C8?2&V!Q<$*W5x#K2vW3+IGkyIi z;av+D|DeWu=gyt9+#)Q2yIAIz;1a?TBT7jC`${5!Ap_)nmVn;90FjeGC_@*x1~zlti${+^er#O*SbRlt@&EE zeDiDRuPpVfz{0zddp^OJNeO}w?VV(ZgIc|UE8_fm6pckv7Jqv#ZM(agGQ_ z_VY96Bc>QHB2ih!U(3nKIj>Q6n3c@re%u~=bZc&C-JJGP?5qCniWdrR0@s9 zl)5@)1oFvhc+Ak1hKeB}Mp&uft#N%(ouc8KuSKQOS0^y|+M!1Bz;DRCAC}wS#SXsM1uv0zo zHTl=NH6Ur<2M*u=u4GlZ{rR6FJ<+r|ITUip7+yI7w*lVdxK5j!nwuM=AqBf{ZP>0^ zru7MB9W9ygpAkM@y~k%RH;rK<@h|OGD_&S!B)V=ZtExD8dDjEWY@nhd8c=DsZZ)vo zdUf5Vqn4I@rSqfw;Dlc^HqsLQc%*MM=f3kj$k_z9F1(*4(1DBx)w$u?WJGslZWPP* zZQDp_b!-}S0SMr$0WhWs?2J951lzB|{K)FG<6F**Yi(?N%E{WUgDfe`iyAxZ<81}* zjxsT(k1XCQb2*2`)?K-}cE26=;pf%jrGgMIEUsLf zqL|8esfFnGeigHKDuHK5eq^@lN0pQ4>}YaEtwBlr3dDL}m%@v;HqSi~1Q)2%lHf}~ zu59f^66vR9D>X<_aeFMnlm~Wp2Sp|RpjkmF`id=mXN+)s&A}epJ3%4aeQY4Cv?Ja% zkeZrCG!$WX5FA(yyjraS?Y#QWF0fcyT52U5aSeedRM9w`@>Fhz{>5)M)@O`cGMM4N zsCNr61*-@^3Ve57|BEa7s}xjKtHEwyKYvcx1dI0>sf-tgV2d_vix#_J@{NdH=}Oy=d2tP%I%jxdvz3jcVe(m;8Nv zd?bhDg3-1Z#^SmVdbq&%BG)tn`ud^n5ZS@<_CrKVlrjlw@zMc_W6+#jKF2zTMn(Ap zDoOqN@)XIZ687#%cvFVp%)1NTNa;#Tx{hF;)55{m>Yp+8$DEw|qp4L{>}iKsz=8oO0ZM*f@Uaprz+gJ|U@%4Tjd{y(L{NSjel$)}_p> z9%;U#COx;tdufiBsOT`D{9A3|86Zf%!CCaG=eW&U61HgcrD=|C-U4$9A?88|SVg+Y zelP&}$OLT|eB{PvW;+r!Qfplg&V9HId*aDrFVs_ZfUF)1Q#C&)Us6xaeU4ARP=->d zE?P9bG&Cv#01044bS`Nwx^I&Wmq^l5K7JU&T;XriUDM#-khm&mg_pg1_SYU4pv}bVj`TwN?-9336CG2_;}nv?SwVV(Ig*C$gYU7<~t@e zH(LwTtcfv#%^75TxXs*ou09^P7}CuqN4aeK`sf9{y}fby)JpAGqV1K#dqD{N#MEed zwlP*L@i~IE(a=c9wZ=gvIJ_VPR;d069mC7$);eW1Nv}hzV>mT2Ssun z)fPYW37qMCV#YTcf54fF07{{pGJDbIPZV$1bG9Y3jD$JTpPrXs81}nNZUZ5!)nMec z*m?<0Kyb{esq5SBf>?S*IroPn$$ylduzX#(^>*&g+9^jP7@KTc#0aA~dEdWLPy8o9?nV|2X8;WHJ{5dzm6e=XsA3=a=KgH!+? z;2IKP(S*ZM)6yF)X8j{0{y0M+!F##vPSkNFB^*FQh<~z$uwd)+kGKh()ns&|bPyli zZ3*QC0Ty>$kJg-kHM}Sz z!}xH7|KcnNmFTe)Y3_$w6Z|O~j)^C`#EEMd?oKvEPl)a1!p?VTdz8luCSHMJ33D!SA~!_tXBE zp1yiwr2C|t%WR#!JSGrNuoV-Sva zSa?$hEr`%Z1Wu&C%;^O-3VK0HZGww7pC|jdHc#H${>!yE{zC zAh8Ev6(hdW5?KXwAz}9uGO6B~tVCK;I$B#>+y5c4R7WT7ajs>nQffD3BeR%0oHKt6 zjknsivBBYCKNb=5Gf*a8;cOQoQ6UKs?ox<=Ep8WXNewiGZ%ez!_3#;r*;l5*5pmoF zjJv}%H^6fdVInCZK^!!pR8v=1SNU^lsFi&hy*^6W%e^w>B$%#bz_XpKC3q5o9zi2` zs>`4@FgQGusJWhX%?DrhXA)&gH{UVfmargL^~WJR_0M^t^*mdo@}_i_kDt_ZH5u(L z5JAr#669XjZ>GV|9l=cz_EPLjGW1Vf>oR%mtt+mTeGQZiii`+1=^5)l57iPlS-7~l z6C2#{2xqLVGs8WG38(mW`Jfj76$F3Bh3q(1d3-4`F*|!y6NeO4S-!EtaY;nb2NxGk zXHftqk*`Gi`6g0Qj^*FD7@IMWL^ECQc$pU$7aKwBfR7*{GJ}5szp%-dlEwI){igEr z*%krF^@jZ9Zhr_2&;2lu!q0in0-D<4cogC#;29(iFTs`ye-;+eNgmW@E4cq;_{Z9hShbp?LARsAg(a&(LdozcY+BABM_VoSffYx8b4r1kk}7Z zojQ$UPSB(|8XMjTsu<6X~vHM5SdCa=M*uRrTtM@dk-^I91QE`^GlA$T2>ns=e^8T ztEgSy=#6;v6Z}If)tQN|Jbr?B!j7xbyq-riz&THJmtP~N-}gOXI1So8Y00gyG4vU+ z_K8zX;P@uyHNq33qu7Ox`!0*%o)1I01AAs!bK1<*bZ49CB60DXB{+--2*#EHt_Z2JxELD3PKyiChS|nNsEg4a9k)ajL@Q^``C5s+MSXS4e0m|p! zh4b_8UTSy_@-zX2f`W?5T2h4K;g3Y=3F?pyRR0xpjwmQ>z>A^z4`&&miF^a$1LynI zQ{Zr{@wVn^9CLKi?ULp6w{J%Qf$cxlu>>&-RrEhAo^Cim9L{t{K=%r0`-c{9W{-9vy@(1$iL4MgkZ+=)WfV9+M&sD_~iM zKzSADiy)*Cs0qh0apDD}-r!&9A^?vBA)T(1txcJjjV%SKmgoc|m?HM3vU_{at|kzS z>-47q6u5G@1A(Do2d=evn&9WZ_hgE<6x8Tgh0DfeBn3L7Q!UJ|g& z$*v`!HF@k;;_wd`^1*^tU=&d6xQt~!XDHbV))EAF70l806WtO(L@zr#Squst%Mm6w zt1Y*FM~OxKVrc3FeZtj^psxcdI;>H5-aDWr%!Gi}I zQEFA~)-iy|{Va}!>Ue5FO|@me57Jly+}KbUwTnuVByo68IpOC?jpzJF@i1bmI`Zm0icjE^W1Y;y9jLl<_qT(8<3`8}9Gf<%W zHxJK_bU8Rl4M;w#*B&o$v(Jlc`!zPL^)P3sDSTgRM|fo95E8+L4I56dgm=@RN(IBd z2I6JJ?52AGkfv~O0$Xcq>(9RBrAMvT9(IM+RD;?DYf^=?5+++T)uIx2CtpIpOI?fv zJ$8^BvRm|msOEd@zbW_!1~6~h1#+<{W8?xE&$gzWoRy|oCFQkowl;!)7#-_fQ)he^ zmTcz3c$$vDpu|($s58d>>RparVcW5m7*`jUeu7s;-O=m$!+?QF{ENT>44QYb2nps7 z#l`JsRa(eU{zO$%f4P-m(Q{WS0s#@#30ddB9+3Ve2(ppyzk-?pfnqp>ZW;{h%k`26_>^1>M)J7Q5_S!L%Mz?(^M z?C4}XgK}`T!N43;g5T(i9C(gCJu)iu!iwEM_)uVHN3~a3S(&1I=v(744r03zDGYL% zNkSVRK$KR9gq3W}2jJ$> zbq`Tm?&_kP0Zz4nW{LSfT^Xai$)9kcs`G%BM3);z58EOEmapqL)J zdHE{zFB0=U1t#iTURjc%h;L<8Q+8Atjxl3F$8w?b*R0m>IHZTz+=POTFhsCDXj}Ge zn7#PP&5{9nBq1>O<96^15mxmvPuL-d!Wk%k0J<}d|5#5nSAi*tk@jSbm2|Q!_7xYF z8RYs=zMf;-wrz(~Gl-^(=EL;5ImT-Vb^~dTI3EGLR(Z6j)jMR%Er+~{8K3G+U=D*S%MeqXwI7w{y#zyt@ zag7yYTE^C>xf!Or_pO@D;+=< zaGJ1s-5$tkpIYCP#iW8T09Oh)$d{7OpMJ2C-)cRiaSsoQph57Hpxu&ZrE$;U_s<0h z1I12?jT^C5*TDLfq`9bE{OSv#u;B~<4c(~{1Y^YQR=4huB($nROS=qn3o&c@MVIPZ zR4f?)F434eIXRi)Q1t0i$vHMAiQ{F4><8VhA@?l+gsjXJ5H4|UUV~@!OVNjF`X2=5 zqEoZ7pzRf3CU_hiOQtujiGT{8@gV>y-&Y9B)CAjx&7WsOR^t-ZSpepVIB+218=&TD zLWa#-dU1*lt_kIPVk%Sdg8&=^*GRc#KaE90449C8PV{Uoz@ULB;5;VF*WuOg9^2wEiT}=1@$Vhc+Y3+?ByuP#t?ce%eWa;3Ehx9IRwp(I%RK zBPfcO=Es=#o*^Zx#_D~3RzKp_xMNTKN4Gf~NJ;cxZiluA?(ryU7x)9mx*tjCi_EWV zaLCsJ`uQ}3(-U6#wqLKaZ}=9PgMq2LSyEEcCBz6eDd5LQ4m=LTaRKi}vc(-Dfl$G_ zl|x;`WJ6=ha`^S>V#zQolas)d8#&I9Bn}k-DjI|?B!|v=`4IL^x?A7Sh~Kc1kIzpr zQEab+T?WnG;)FF;a-6q}opm%C^t(AxKY<9$QYt8}tEs8!P!YDmkst8>i#v(1%AP&T zNic#^9zSWIO0|z#^lZsu@*GU?-(es)8pLsk9h)=5Ax7IF-+W2+Uj8{_^7c9awn-eY zFjz0%CkvJb>-su4Mx64q1=$XpsmbLgh@{6)pISY(@`Wuc77MH&s!~d7syg2xm|k$) z*1G}+3;Mh|sh)g$nMXDb4krx_Z=$z&I~r2Tmge~(HK(-FXl zD7qoBJ6`oE4^{N{?Et_m$~FTFH)3NXoXzam}TW+gV#{QNpOKqC%y9r7$FtT)Mi{Nf#uP_I$I) z<2=qE=dts6_`}0Pzkc8E_w)I@->=L2_2r`pnA70zOK%+@G!?#X34^q|p9yCwr*cJw zdt-(e&twjWDO2W$wf*4XCqS~U>ZpFp$sHX|8ZExck<*s~7L=)sh7FtMmyM=X17oVo z$d*M+c@@>Y_0sBBTDWrZnc1>{!cMQWI?<@we$uV&X%FEPTuh>!6K)W=1BB6j+dn9K zzIwzKHa1P#gxRxyEDCwr+5iq|d1_+yUr}wWu-7sv& z9`S9!YZ|cfBx;r;NZd6_E*n#griY!?ckb#gE7;-)wnNX*&hsAuEwNj={xH*qfqEWg z*^6S$poGBLavtHD5XnKt?ipK4K`~2PTse^c#4OO>w=b2Li+4(>@~wJy9t)Nv-y$`}z{|b&`1lU0Ew1 zCz8`t%w|0N#N!e0@;E2RdAw47*n5ssUED4#DtS>@&m>Fq(+5|^-D_wt{Vg?h2PL#v zydLgTr$#QAdO)r3py7_bll{9#UgLD(yKxlrxnX$Oj^oWJX6$C(rygnB(*tgM0PM{x zxRe8FbKX3?1nNgpOs3bX%$K7h91Y&3Y+g{QobJ`t8GiTT;XyxGq^p@zUuR_T`Fv=o0b0`CSi+|!+mS?Rj9B`2%IGt?p#%rcFNh*;xh zHlC~e;--^hVLqSJu)Hc)O_igirDYYcch8>Q@as`yUfs~W81F&s07hvt8GIWOHD=g( z?8{PfoOgz5DM(!jXP%PD%C!uAy>D%OBdzaexq#@)@#|N1m5&U0G@CA(e2RSOz{eNT zew9AJ$n&w8nIlWu{M#r0P)h8)r!c^+PiXh#+v-LDy)~HrD~#4Lm)BupY_JVL-ACa< z9hh@zK`0aS8rqazzH8??zM?LpI7&EZRh6SwwO+?77iU>gXxQ^Mzwa?WOSX^mq)C$7 z0pas?TV1iKoqf)qMoFfp)zAp@`MyiQvl2eW(IipM1P?{QREbWA{(3T`Sx@4X18!MM z&t-bwaNw5e`ug<%q4D*np=Jb)Ltv=jA*vrd=%cT%FR{>YV4Tg@O#=uE>wuG>81=<2 zE}=(&nevlign9k`^Oz{wCmb4{Gh$<6w%}~gNxvV|e8*0FC2VPNX})P;wxN-Ea~%+;u8Te)7=w@p3Vt<$DESlN zCr*rH5>NayXrf;p*{5Zt4ceP{6st=$#B@q(Y90zi6@DskPfc6>d2zxK8yt5d0~n6B zrG|9Ui`yl(!d5R1NbPQWgeCjuby{T@WOy1aPo5M|vZXJU0aF<*WOX=uvM9|LUN=h% zxII`BySax;+v>icbs_I`Q7$}(jN(yB;#l`GOGEeS?d!W1Wz;|T4Ci7Y&NZA+(wvLb z#u0upxN+d?7^tUZE~dp1HaVSb?{KDd6Ymzsx1|{^B3(d!DU!?FBN?Pi?AhD4wg8d{|6YN8|}_ED&!dEH%iyFXwa0~!oc;WHdE+p#}!$7FMtXVq8-Hw<7DJn@0wp(d`l1EoSyaJ}Y1(yH0Ms?Yv9 zzqb%PZlDN}Pj-iC_3v-Kw8pWsac(yRB+*p)G`E_FupxI=z6B8GNjpL5C<^^u9<9GP zKaciLCIzSpuf}-fE|J2J7fkPf0j$GbSh77D?WR(SCEgP670}1t0m4+xR_~J&B(Xd7 ziyuK@+gMZMBb#=g)xtfcJjUDOP^uO0%W-3cpo)BNRF7Y#Cyi*&t&J_BH8FWQ0(s8| zsiyLEcoS=mAS6lxUrUh`>LP{8Q$$wsKuGCcF%;4f!f;D@|$@Imf(h0U4VYFDX_MU>Pip>aG-G!L>2^g2POg6zYO5&&vGr2orVvbxJ(05E{OHX>yR%j{{qd9F^NXh$W)}*iQ zUd1vZLs5PLlA#wO@1d)UeyrI52;xp? zZQq2)`sA=M#T>;QH&KP~SVBqN2a;?7P|8n=I}lvmQ=vIZ#rZ9s3c`hB<;tLK(8u55 zR^1>I4O(Hj%`AhULW-2@t?ZC~$OGZg(b1(%7GPl6dRVH&G&{9um$7TQ)~HdQo}SX7 zJ6B#CIC${k;Vq}QFN&VA71+-7<+-~|A+0*Scr*!jig#JkjHdyfIx8!wP&4lY;U|O- zOSk_DcY>lIv=q|2vp*M1uT;dxDMx|;Xh|7iSr&54z$2~&grvX3xZr-J*4>jfU-`Ib z)?RBb&LU(t=VG_Tn8_sA1jJ*Nu`S6&#yh?gFDp!ux}L3D!a-(6dQc@$WaF;a9#q62 zq*T?{C$xt>k81;M?SAHeH+E-q)C7M-(3@kM+wK@0a`gzTLGK^pr^!)FbSC9Kk4cKowknh$mF3*FI>FXtHs*= zNS8yn4PTfkDGx>iZV88`-B!7kbiWSHq@zz_J9VBh@x7h(<1+l0XdY;b#aFa!+40g9 zg@*dMXF{62)<%99n4Q%cDE{!U7inXQLa4E}6)O?XoHA|O;g$WLF-51=GV8bt)owV$ zc4M1%Mi;dLdyj-vV$#XSTt-D8@@K}{vT0=lTwF(COd52ZAQ17)G6X z=FY86;h%4Y34ntqtH$m`xmg!~DlxIQ49C-P(C5$FS+G$0mh^NfW!jSri;E9NU!K3@ zB|Voc^@0LNl|xWfYf}Y9c|TFDjU&iusDdUptrsvb`6CzSc!JzyZs8V;&EhsA{XRKJ zAPE*;#QPQEU`Pm#X?V=3Bk+|UDS&8>+j)7;neW}NWY5Q=0#_Etz_*kW;Lps+y@dN% z4jQ(@s7kht9Q}H>q2W*(Cb5TIiV(axDf56QOK4!lJp}J6Mo_Sw)M_Ash~hnatLIXj zuAp`ZYfl?dAi>EWP2WAy5MzNA7w}ENZ(Wnvy|1TAWOu_KwSZcM7qeL7PfvtdqT#16 z<2g=VzwUz-bK&PP`Y;^w2>?VP!4SZX-fwfp+l>ycPHTg+IG5=o01xp@8^M`$n~68B z7Z&!Gjmh~U-Va%qAPEQQiRs3~VgRxy#lcz(EwYP|k3{zpmj`ICA^6=V%`mGegp~?? z0Kywb`la}@30sgy{s#8B*s6HQj<{=u@~7*Zz&Tyr)U@sJVN;n@;gO7YaBw)iwmISl z$=pBizP{_+S%G_hoM?EKTPp1ymZ4sZbh#0|F3GcW9HMm7SFC}r@VPTI)lIHawR0@H z9>O%A>UeGFpm~!Sr(x)iF?I2JOE#&ou`$MRF`oE;6Floa@)(h*Z!wqu9S$QI(Bek! zvk%L14)mB4``x2ap~xE;Zsa6kM6LSKqey_8cV|bdm=sFefsA=n&Im*=aSPCY2^!>} zlXjC4s4CJg&f?fx`tKVS5cIxSx;Aazd=~3b@Y^$gLcd5|)%C5|_9T4}$0jNJhmw-p zwW3g0#Yxn=pWaph@|^MVb{4uO@p`bw(_Mlzto(YolY=Jj^@+E6^U8+5dTG)B7>53f eHy*vx4^A@wVxq7B literal 0 HcmV?d00001 diff --git a/documentation/source/_static/BeH2_PEC_6_2.png b/documentation/source/_static/BeH2_PEC_6_2.png new file mode 100644 index 0000000000000000000000000000000000000000..0d139db27ede43f090d8f01bf262e600c8c97864 GIT binary patch literal 27582 zcmbTe2|ShS+CRP;6-g*lqRbLPB84)`OooifP!=LlG9*dJlq4c#NQR6>rpPQRWgapo zgrcaV^#4BgKIiPS-}9c&`~U6FIcwXN^{nT)hwHk&!*vH~YbsMyvQUypBx+R^1zi$p zwK<8js&vD8{HBLhB^>`c;Hqfks^?_udco4ghNNNX>g?d;>R@+r=Xo2Kb9PQ=#DpY; z_U+ku+SS$hoRqMz<3ImE$jRlDFgrJ=8!ocZS>^aS5{br=`2VVO*#~wclD&Yc!jWU{ zcW3%98r{jQ?HSQ~e0t0v)h zYTvOC195F_xpvdDj2iy5Yu8387Oi^!EuivF!&%!C3^B~tzD<=DJHL9{Tp~NLhW*un zZwFQfDR8W!#XsT)#pDey2L=YF7oHeVR8UYz(#?M_qw_f?@-ePt+lbSnVp;Sy>rLb(MY{~ z6Q8g!ja<^}*Bv`?CBZ-ij>q<0nkn8f?#wq*QaCSD{64FxskuHU)H1N{Bo{aLE*_rs zw{PFRmJ{~ogy}IWhg(+vc(vcN|Ni~^TPXbt21jF0jMyA*{iw*x`@kUHy{D(A>iP2^ z_MZWJg)&9tUI!phzhlz0~#$s|Ky)f~t(w6_Xktk%*xTV4+clNA^ zshQcXt(4SeJMih%^um$1Zr!>%KkwQzH9PBETU*PNr-Y$fo0s+G9UeaQ!2@~I>5xnZ zr@-L&h=hdA5h~UbnfU%f4~IEXk-EA%KG`-}{o^w>Hv8jql8TFC%);>T_4N$G^m=Z6 zM}#IZ-567sFxJ`8mzKB1_7!eWi$J%=4yM}_L|I@%z&)M18%ZGO!mXo_;JieQc z@7bdk87dkYn(p4-WW!wZ?vis4!VC}Nj?-RQ)8#ft3tKII9&3F0Qt|ZZ^xh>wqdd~K zZQICB7RJWcF*nn_g;c`jv#Sjr8m&@RR#v|&NjvEC`{{mtUJj0AgKAs>-)v zZ*Olu?YR**V?CnIGv|U`5*!??wDjVR7{zFPXw6WSU&p|BZsqcova*BMHt*OiASNbO zQyWAdk(e0nM?phAMuGK2mc~|FCErtWj#pGv)Ol`1@on*0`I%ptSy`6$_F39zPMxA4 zNuC{9o1BuO|Ix(EOwiNJ*qD;ECAbd%V`(37Rmm+YG_Twp5=;&Wb@}IOXq{4YOkq=P zRzU%sQ^g{Wxz8`RwQJX&zUA)r@h$Vx*BSZq=lB2o^~+|==hp%wX{UsQ;aElUZNZw6 zy5Mb1kvDHz2btJM-M#BDv73pzUQxVW>A*i<_4D6YJA6+G5uMA-%{BCVc&o5beC=vq zrH0>?zn5gO%gdKOma%lrVAvXdEKG@<9Z~rCt@5}3g9p6lr@MV?YZXW|H>LF7zl`IR zSbq5U@m9H{Eg~W!&Mq#k#X@O^_b_lCIKceo&6~W(kN3H{g{&VmhUUTbMSy}hxn0))#Klkg~;AkDA zMeV&Kek78)x%v0`Hm$iV)R&{PR^q-@!E~ zT8GEHs;a8H{Nf|yJ`Yb%d)Mi8R1D9$x`-HBAXNYJrx)qjOy85UqxIC^7ZdENX``B()(!juAqVGxW^g@OA5{XoEjoQAni4K85RRq&}K=63M zT17=g7n|`{@ijjezs!B@bFzANc~#z%Cls?MW=%CUH@|)RcIf^4w9heyIi^*ZkhxlV zt-LxS&>s`csI`B7lbdw;%XrHgQjufdH3T-&v!zd!G^7nQG=|+;y`>#boGe0i95gUk z{x;~tDJn`w#VAyptQc%PW>)HY#KB?D&!0b?fBkTKX!2yUkdTm3P|@9V!Y6?UZBh-#2+Xg%qvnj7#p*;@g_6 z9L}Z`tzhRTNF?-ZbzNOL#Ff{NPH8#x6s%iAvDx8qvv%}t0jrs>Vrz?t#8}ADyZL-q z;aJlx<80R|;%1W#6r!s>I+V1ojB_%BdfS+*>yGI&A!|BK=vKx~@YepegDPd4#6(5k zeYFz3&%RPs81FRe9P4r7j$?_8hSyz}Dawu6uN zO%p>Ry3>r}xnA>1!P0nHQl-LeC9hi_X5|_zbq>MKjPN~s!)vOF$WG}S?MlxlZ>4zt zWb?|M#s(I($!Skji4cRNak%60mJgbNN3;E&bD2Lm(_f7&K+DX$5jkVS+SRO*&W9(9 zFi=}ItfjQzR0`2v^4czOMgeo<;Nb8&-%8nQd2w9Bl^fGRch-1=ZoqR_aq|i<+mWZ;mvlqTvAhHV)mD;pqLnQ{>*ARi&+;c|!*s)_n z@fQZ({9wV2kBvOulDaukwcO*|$Iy|Hlf}iwJ9q6`m7kwaedg2jbV!EY1F}~_Q`6D% zZv&o8(jNQh=;%nKN478JyuGDme$L9**3~6B_7xeMt-(~Aw%YAt+c8?V>dX9IT$6Rj z4)K8pMMVty`ufh!&f1=d?d`|5GBR3yo$hfN6R!LzTCw!;;Lo`Rc4KZ=Jl8SG&vb`{E@9?)s=OK*AC8=dJi;Lk$%iiavBblU2uxIb7#0cLvJ^I z?uQTLV}^#!=FzuruO`)=6sbL#FBIK&QuS}?tGAvxih>{4} z#_@WI?k(w^MhM8qv((F_TqiC!HmY2`c5P(&*OH2g%Cno?YB!>y{1J#Hia%wWJUMaZ z%w`l1B`yL20>!1J9Kym{|DNMcoI3TOzb+>5-8*Y}!z#XiR+V(G~%um7fU*H=dD7@m-_MT zQ?AP(rp4&^aenWX`ug=zQBk|U#Y;G`4Y-d|_e{06>Td9tL-qLV`E%72`@P#v<}*Ei z{yZ=&?9AiBd~wTv)WmK3_Lr2DbQjnJ%Pod5i|r#Xik>e&bqu9)E1Q5Llg&sJN@)y3 zVbN3fIgTrSer__Syu1$Wcj&(Cz4TQ)#=WPvH&dN5D>L&1$_W>X|2^Mn^mo54yFjzB zuyC|_&Fa;$Nw(MOJe9i6reP_1h}Dgd(SG;CUzN919Kf6-@A8N zi%UwXUcRIz?RA~d#(p@_U+U%()K%&>h4~D4lcGE|c@#UR#2JfKh$qTjHS)nfvK!A* z@ED2#3ew#}7wIB-wd&tJkEGqY^(rQci-%|Ab&`BwaIk!?`BUqLCFDQ^fgpb>#*=T7 z6-jI%b3bNlXy!gQ3v=$@&ww4>VBz!2->LE!gIctZKT1Eg+?c??*kElw>uN=fSE+Ql zi(kf@pAI~>i<};8t!*A0WJlGfbn+y}=PzF@ZEdgJxUuD-MJ3zWvuDZir?BI&eQKMc z1m>o4E9LmMj#U2sCH(%yoiELh&%5Vh%LY61~}J&e@V zH=&-vK2d)WBjTqeezPi{s4CQBb8~awUf_^F`JaK7GlA9$H*iF3K;m?Ea~q2cYCpO1 z9VyYbZgELadAqwyo;-PSlMyG!OkP$5RufC~@-DrVQ2$?ltXjZ>&Eun6=;*2+uw>L( z5qBVXWo1eP*~(agkfHwv^bI4P`{_ulYvsmi850*nt5 zs{8v{c_9i2T-ojRo+FQuP#ZpR^a?atn@p`m*nB7uR%QbOl7QV;ws zO|bB4EiW&#BnAgI-g3^j=a`9k!6d8!Mxbh7fGsR6Z1(f#YE+@S7=o$I1UWcT8tQU_ z9-k5wG>@HlowSC7laok?OW!_mB8g#RA2BdsC~_Jca2U_d$x(81E1Y(-v*RZL13c+X z)$-`iUihK^M=WeRWqM|?f?7WjUKSZ;sh(xf&a! zb)qPkL9luS0agB9_U!KI$?(?1;t;)i{W_)k9np2%+}s38Ic<^RFn3fAhhvF|-W0$%K$9`9KHiFu5Pv4bTz^l`({?f)SP%S2!;#q%6U5Jfu5 zD=ISJdx9jqAR%pa?+?Iqw_-%p>$|h-Wb1uJ>5fLv6&UxAm=1Dh(p$IHzAI-H>1}`i zf`>UO&FjHt^maKN(j7*|#=9k+?I3Z{ZM?nb`0Al5i?h0b*Tt@rbl0z6pEp(}Mn*%a zXXL)1-DlV3pBG8b03fWx4`J+Q&>bU+*Hdl0CnDe0eFGSiqlS zj*W=>EPIu7#qD>z`O75MeQ{~YYIdkfMO|It(j^%rY;6213sYU3LPDC^UkLm^==CLZ zdI43E^M`hR63uOnKY~&dE_MIa3#NYyrXDW@O^VXRz1L!;ys0tnU#UVlu#1;Lr(K zhwIUDqJnG-+`LeBlgP=ny@mE&#}+s_l5KZxrKDwJ%h_8Wykq5l`&XyuGgByXakk~cD9YHDuIcN&zL z?s-IU=gu9;bMH3+R&Uv{BQ%0jxq5K0lK+F+%aHKf^o2rlNjDPRIHbu1lU(?dz2~Qm z02&#hh0Lly_LrTV9Xf&x%ON5{3(#1Vu9LFYx^WANBk>NqjE!?|T!0SVckKB>lo>UM zB&&T%o?oX|k?zX+urUgmJ~OoN_CJ+yalKVTm}8h;O*^m z#zl|TBjEX(T3i+M4LpzE(Rg6bCnRJ`B#(B5SgyI!^D`{xKMjQ5y0!Ifr4RcvzkCa= zLvw4_Q?sl+ckUeNa7%);$i>grSqNG z-^fV%`f&De)GS$db*tOj=%0El$Q?a;UFbm?4^O)5YVCW;1nU=j1vB{|@e(CzqV+yi zg8RtEL(AVs$<%id2&EAfvzs$Af-ZjVBOlYmK8bDPpY!&cHTfYDD{ zS-Bk$7taUvr_eI5__^Nn`0*ECDmH97x9>zYzCt}1Um=iN-g<0dC+A)v8Itt1lH|Q* zbZpB->ZO6kH{QNA?fG{3GO5PkTHvu(yKE&4FsP=XQHP{9`}0eS)z9y9rFUZk17}ZcmeC(WA`S`14I|#n z0zM5t8N5_88ZIl;73;S2)e(q|M9lKUn-uJ1U(hKd0qmEaMO$3jC@U)q?EXMi=n5!t z7dJP8N8P0MpmqvDe}(3Iqs*Rk;nRT2r$VPews7gMOxLKlUjwwF*|lp|%sx9hRM-T= zGKz~4aUqp6Njnf*-rCx_v8lSIrl(iYI@OiMv=X&3kjZW?Gn%RXG6~eC*REb&OFA6R zCjGj?+takvRW(|Vp7zaeMtXWb;A9dBMRHZ$Oi)OSW{{z#W|;BI9)hn3QYh%EsHiBo z?~lsk0)BEC3~H2|yyf20-Q9SbN>`C1VyA!r)&Bka3kqUmS&*?NUC4;b$qdR5YO_wY zBs4ZR*Tjn3YwPRZpS|SfCf=E27F=F_sN`cKX0*7hZ0etXUbTRjTZ>3Ji^`W%Tznhw zLlpP1Cr*R)bE6G9=H{DGC>&avenh%FJ?-$e{9>fqdJL~lR*K%p7(k$w?60j^hPm|^ z-iX-P%sdS&^SNI8+$t<<($LIIC~z6JwmI?KsQ)A}-_QJf@3Q;+6#4x5vorEXkM$)y zGnRO8^~o*42Y~+m{rf)sr&zgx>S3CGQ%E^Zx1yrN z_03j(6{|@gkLcD1%+$l3*(NauF9;#i3{OmK@I8DPEC^r!esKA#eGd=Y*9OvtCMa{! zSj!d3NQi=*-0F{SOV%(lMyF*)8VFfb zgp4&uryweApr)=Ht_j%ww!*CJyit+E9n(8fX7?+9ex@R7v3#qBlYqqO+1Yit8(M@< zjkx{K69!O5b9mgsLIw^wJK4%z+t6Tzntx)jQg+WCSGL^Ts6QEM@M=J~r(H`a|E$(V zzkvDV;N_(Ni&jv;XyHW*Q1$G}Ix4X65t<4dw18s6A3iiKFvO{REph4jb}zw^+H^a4 z#=buCQahL{9CUz@MY^IMKkKzrQ7hkLsF`cS5=Pe@3J5w~Yx*keReQ&YntT!2(fW!w)R=eXSiVwvI$+s#h$+}QW{^+_u3 z8M3;D2Il%|^%C~M)x<=0qPYIO?9Ii+HPmI{QwQcg-+SqNXv>9#5Bk{K9Q*edRT(^0 z3}zr!4&@l61?jFF%qV%+?%lih+B9v2PPBsVy+c=+a%ByTP~$#h6BGHJTN#)H9XT<~3WJb%{I^=bgj)n`xxRtPXJ}xMY%KEa(QvK zwzh6zX3pp=ViY!4aCXjnH69Wk-hjxMhti;sKDV?~N1dpYgYCPL&BE9`~>aboD1FQrxa>eZ|FNrlObxJWJFyGLDZ?M~3Y6N%(KFLI6@AJ86g zlYXR5h;;8lOCYWlzPYMlGCHzmi&Q)F__@}-r}QLv%&eWIc+Ck3;?l1LL+ms{z5>fe z^Ub`pv~&cCoREYPF8#1844WJq^GE!jFcuFr{gWWc>HO;Sbb(>N(?AFjQ&H=kd6h}4 zBfznhGIam&&Ci7TXL_{ns)pfmKdZA0#k4tB%DuK13!sm_gHt3J$%kfT42V6!)%y{2 z1Mp>GtOt^tGh)^4>@wa;3Xo8f8WmV?B%J_wEcgs?gM_q^P!alkP&gRnpIkTK{ml$Y z)`LX>4Fzk3)itcN(t4NP5+y?@k5(^_gF>1D()-}0C2GV6;#!a@Qe4Wu&h(8!VJN+O zRIv5kB!D34GRVj0^=jq+Pnci>=|)UUAhP>*l&%EhTJh_vNr{w#vT{vZ+oKZUsK9hp z)0RE*mgl(-{d5#3yRk5`UuLT@74J2S&69rHpRQS*W^Wc05`+*Nx?SS+Q_pYfBZWxN zM!K4Rnu8UJ_IzsB#?w=IKPbAgQg%HJTL?hPyYcaw7cVv;<3r1k24G~7R<^a}Aq|0@ ziP`TMTfLN#k%1BP-ggEq);(gCF}W^z=^I zhCF+Ac)B}ZKAX@F_*l;&j}>$Z18FsFS<`=ipGNnlClx78Vw!B`)g<3k!EY z$v`NO{rz)|kB^U;CO)W!KC+VTLDC`;>!~_! z4pGtG&8`wxZyz&EzBX~FxB?l)sMM7o#5nt1$+N2$7Zr%yLN_%O+HV4R@$WgmgPy)^&-hKv!V}A{L!|GLeH-1yRad0uy9Y092FQY3 zdh`kHhP5$LZe+c?VS$aOgA_PQQmII-ZEdu4bgg>?Ev^~=&J!O}UuqyL)jp zQ*OO+nd_pBEDG@?`S`OkywvCRyR}VZDCwr!e-kX@%h-iA0BKV`QN}yQm6lc-;A%Y< zEQ@qU`joKe?=zbv6Tb`{oK%YbXmVah*UHscG;kv$ZG;Q7Clm_Wo`0VymmQUqhHotiJ za(HAUdEX-4KgtbG4$+II3=hUyi#l|>9`N<(b{0(Ba^JRfYva?L>}*Al%bMFzCy_Je z`xT$xF1U@#R)C*BBjdvn6lQeu9g5Ev>y3HhS4c^Y{?kiaZR7=Qtjmj4_7fVAU7O-_ zeidV$P~&p=0DKVvc6$q>V!-0n|1lsj=;hb*=#;(rK|=Q8*!zMuTzWF;yhKtFE@8w{ z=&st}N>FmcM4l&{E6~!@sy+Sv{fUQv{E!2sEPqP#RxadJo57M7x0Y9=CmpIZJ2mFn zA-hJDL2F14FMVXCGaS(N9GtpiBl*KVsb1`}eOY40Fc^m9^)QWucADNAd3)6QVXBsV ztjX_O_5KID3gauPwcE^@M#pbJlZ>DQkE^1k75x;FhO)Z4ddkw-Gu)!Ky|1>X&CGAp z5Sj7czwuW7rOiFi1yis7m2*h1`c9=Wj(9myU82t>nKjV zN;^Zkd?2pWa(er!(2psO)oIO-8GCuO>wg|-cMmaiaS^>48M&=#-D=;aZ{QSJ{(wvV zBIcpgd0pp|tU0{De0qQR%K@5-EVaH#fyez?eQTRyWCK`u&FO*^vgkLiU7h_EE96#8 z4Aa%CS5N;umLhacfR_27MptsXjjS0RpVco`vbgVVS{-|sbi`woq- zAR-z-du-XVg#;f6@Z!)ymqi1RKjBW8a@|3oj8gbWJ^(Hfwfo4(6(Gji`uhFFPai!J zOVfq0l3Ym0|dAU*2A%6ai=g*&y zwFulee4XN7zOO$#U>Vmn6cznD8@8%Ej=x^P-0>IQcbjcpblBX{KT3MzhM?eJ7$>$n zf1kCS|1?HA()sK3aYFF|8%?BljAO@8MU!#)Rk>Z@H;3N6OEr+z2sPTHghdSZ zO~)BW2WpV9^E=|E&y5ww>utOsYZ7|f;Ee>J#2@z~{*S^rDK%ALbZ`*d=zo>xGNRls zdUz1L4~ZyGz2@I2LjbIg5->crIp;}H(KSq|WHd&RE%xefHc%C+{HcH~%?cG5M zlfkfez4XvDn@jl9vGO0iDiY5R+*iAToUNSoi&s77GQnjKu+-dr*qnEm@N$r)z_+K= z0WyJ280pH z(WZ*!ANC}g^*!_vXI1}#B%$E2ylto%G)Mt2=mx5mTO;9BqAU%myT~zaI$#D>It@wk z;%6=_oRH^{yqu_Dps#`ku1%5;*aB2#{?y}n*tdW8tB^=2R|PwBEy4lElfiEA2?*3A z`YZ=Q{)NP{UOr%R$LKbC9V-+armgEplarGhHg2@b8xufUcq8ULR=V9P;hT*;Bkg4E z{~vWPG=cNbAy#H)=I`*3X=`hfsHmu1VZs88uFlG>=zYwnpdclHdhz&xQS@l>chWr5 z`|H)yN8P{oiB58G#NYW3RJEin5`X(ZV&JBhjt>2l&&unr9S|MnNanE;C0H4L&J`k9 z$dCbSa3{%=BRS1hmvDqML2(A2QAl|JEi85oWGwPFtao}L4pGU!YgUGw;z7t)Uhfj1 zrqZgHtDF)%yUP;cL;i#z=cGc1Ty>XN^QL1&KpCy84R&>zv(8Z(LlDYlrWbAv@*njM zZGzxtEXHG=8hNAP(22ImSkpTizR<@%e3enn#J6YBJ2DtJJ9HjsGv&tnbx1>THNg-)bDTMY%(9dS(Gpj z&sSR;m!Ih`^)ER$4!}SMwGajhN#~)uelr(p-t*y;wiE-PuN%w6E z*)^Gly7_LJt2RH(s_Hzj&S>Fi=E{gM^0bY(Eq=DyV^P6Mpu^Lr`8Ko5ydC;F(_2_m zVAK4r$;-=YpVTS%V9dQ2^uacbg7#wfa;{O_)C$Q0M*}lDR=qs`;5JRTs_|FenxjX5CR@mr;4;@lFg3}8dOI!_YF8)3( zZFntErA0k03$8_4aCIy*dL|~ppIv4;J&NMRB?io?H8f$TN;|=o*_D?VAs~-@0%K7uGq!dM0>+-kbQ+pLc256VbX1V+gN{s*09g3J2 z2wD?k_i?Z?90kW7*a=~W?FLDmvXpq>+ki`>vYs3CL|wv_!XFSD5WWvFUu?swf5 zaxj#;a~-_`g^-s z{$>#e%V!rJwI806|Ek!ssFy>C4hsp{1iYv@dJeA8n%-XHq_?fD1s#L)A-@0e*8btM z%{-Z+R*tKJAgzbW*3sM9*H@VZT8T3bUxtxwp7 z4+*X#IPlF-bb-Y!>qD45q{pG+p-B7qk@^fTmKot$MV>Qs&jnK1Kt)vpUOg`5I>Kvv z_5U%?0`uz($H*>!=OmsCo`Pg>(=WZ&5!cW~R{2X--Hy<=iu9@Sw+v|rTr?rt61f#6 z-!TINLbVInyh9H2t-J z?^S@$13E>ki1lADEr7%|$6lPmd+762fUjFRJ5S;Fyaf8Pglc5#T3 zLCvVxBP{N`Xq7YV25&q;LYqQC31Ni=L`49R@ZC)zlFC3#F!I}anDZm%2ocwq|3ywa zVXP-O^xIrTO^iYM`cTw}>#6PkIF>i*s!{Rqi+Fc&^!9WZHtNIurf@lwGR_m1)^%3+ zrwcssE0JHf$4qVL1B)f$`eznIT#qn7iiJ32IcdCl`)`vxSHUBo%dg5MC1oioSU)D> zrE%~dZqGki`_$l-j2o!J3ssfC?XfF%&z1I##(M=VQO$_79HKVr<&7r*ET4TCf3=Ea zAbVjmo|>sh4io$`B#4oDZFJl&3~S(DeuPWs$iDKoApw^J!S)%CiQ zMn*>W!1xo=^LbOfXWg(?@D=3cT`noP_$49bc(`Wj71t8xv*(}vZj!R_3oUy(yj!vT z1xh_YivX1K5{tr8K?`n+1h%2iImiU-le_Mi;Cu;_+K7qNv)Pgc2J z{xW2>*Zy31^UE8m%`Ghr0Ou&X*Yfi6;!gJU2$%o(;Q@`7fL^w?tcsuauDH+>pyxU| z3Su|w7#VF)R&KC8K=bESQDl)wRMEPdSR2P|H%%h_Q9cqF9(D?TmiKnMs_Klp?H-Pii;wLZBJ zgMY)bX`e@$Fu*WVjvSB59tL@ih-U){{gVpy$Q_&CzkjbGwOQ{vzYz1S6<<`VtmaIT z1uM)WKT{L+VCj0-0HwLT=`5KB>lul^y-EB~A@<5JSbn zLeWKE=+c@Y{m`50PyO4d^tLmZrnQl%Ww!|!oOl^4KJm8XQT5W?+#HzjIshMQ8`v@x zj~wxB9sLD6CxQ~W=m#A$CB1O6fu>ZUvWJHhRxOD%Jbbs;Z$1TS0iKvEawY!je551n zM&dOT;9eo@$p2iR^YyUc0aJ;&io3#Gvj#=&D}O~YQMi!U9(#!lC*KCCxnB0?>>Aiw zoFU6Twrg7tAQ&U#U24i`~Pf~4t}%p#rJXU-})~QhWnprGbl!qz4T*}7hKv9`Ys@2382te2R>Cj zQDz5L1z}tuxMPs<>T&z2!HPFfuzTHe)MHm$t01J2%Rhzj}cb0nc?+_Z9@1#NRZ#(YZy-OmY$fzSkOp^PDj0Z#T@IpSc9BYb72jlyF zZgebcvH!wgRdw}h+*?6i8YJWE_tl~~_U_#RVM^inab^PX8yXt+T!KqNM@I+5ECn`$ z%{awnFn1}yTKMJe+<0Rdg3;j_LUJm;_=N{95enR{U&koS#X5@#1s(jbk|gb|w+D6q z9sL>EpC$SE`6b;xY&FO+4Z`?6-{FHohj>&NAa)Bd1}PX7lDEO~w-yWrVWl9NYcRQl z*MsydO(Ra*(2&rx8gC|vSk<3^<%pOFrEoTCjA4(6J()ZUzjI@h3JQwRo>PfY0$;v9qN7*$Ym+NueMO z!3E$FWC)3Rcx;T@vmfOSwDVoEvM_h_JV*U00(%M;ONrIj$u!hK%%r!Wh7no_uqxdfiKlaml3<+( z0x1krI@%DPf*1oO0fM!5Hz7yEvv7=$UB*lD>&O1SiIQX9d~iQpL$ypcy#RFmZ*=v^ z8zP?Z2QL^%#b5RFTXR%x8}FFg`SVJ!uI`gkh9v0X;&MxddLvPFM0(DA{Fsv$ie@}e zKGuEE#o?7$eg@?|sZ)eq^W6$JOB#Z6FfSJdqa5?6yA2tjQH$Gm-b0@T1X>k%-5ebq zgMD6KW;rDq)CJKy-)qsub-J52DJe-31~b3DC-z;Nfv+M3^uGG)a7#%UC-IlVfbhXE z2KIZyRfUB7;P)6+}Eg=*W8 z>k0AY>!-nRe!a94zn6bav}Z&5r@pv-Xjgsk3!$>QX3diheu;Gh0|UjN zFg)O~Pgv3`2zfw_8!5D3UBHFR7@i9)WlW!XZ4`+(R|% zoL=)u(a11KDXHqZI>i)kvc0{IdCso_K2Q$nnVGE~KFRN$A!9202-bp2kvH_8KPbJY(uMuoV*@aIA?vP1YqM~92XR-E$;uO9^C!_o%+xs#-1f&S`P(n;C`EttE zcDFPSY18)oo6z<69iZcct!s9ucSs%mK^d{i=Tg2EiSlK-j^VSy&mvp#F!5{${^WrB8z2Eyv$WQyvrh z*bmKR4ktL;;0Z+H~3fc$SPa0eJs**J7e=n z{PEM>NEvlieiZp03*@^p-Yi6G5ReyvtXo@W30Dmg5INrB$B&Zn1bNVH1SSh%JwOLH zDZ#^Fdl3#jSoYM;c-l00eC;`vfBFpv(M$)I8(7Q+C|%cKK*9j*b{IrUlH`E{y95QP z!Q=Qpdi2Qr_h9>?tbKNBenxh-{*OmW(TTFzqH`bqo9T378NygCwkItUQy|uW*IXT= z`};emv~_h`Pw347LvDQ^y^8Hngf;pE+QbW8p|D=dxuVm$Psn!1Wy@eTp%TCTvKDch zs1PN!vQvBZ>=}ZZ`}wmAO8k9`bMf)uRZkw6* zl6;mIdDm}dJ$>W1#9@xafAFU#RN&kwWeW6>S59^{B-(m1j4|H@v?eZGHtb z`0}UU%UZ)vVeuSAUabWpe!b_LqWq>Ej8+D4n8JylgatwLvba|w{9>o^@$Th4Y)4`9@<7xb@h!~k?BU2)gwF?h0|wrGEXr^eu?gZ9J>{ri>ci~pD4C5` zc6L3eF>-R!$iI|aXBl@pOrM9j@Aml?wD%p^x#ii57Y~nj&PSf#rrL7k zcS>$I{4;X+HGEbEw_(A=)+B%)B{7pnjvl21AfRrYszQr}9S0dA*Zg>hkgVrl0cBBs z(AReE9f2O>VJCxqWxn#Z>Ni@kb@xq@G7EZ}n;fJpE^)#^ujJU+42>BXR&%lwN#Q z_#mvT{fEC&L1r7dRqF{YGd<>mi21N)9#K`@NVGrYTF4&y{qu`UzYd}&d{4|MCJ03h zlx*YWLLP51qDK)V_-@xlRW3S$T|k1Jx(_gR$Wnjn@h~|m|HJaXvBSI) z7yk!@CH6jH02Hx(!Gso0!Zg6FGH?bXJ0&Hk(X$iPV%hb6Wula%Y}HuEj~QN^FvgOP7CS8kcF11DUgk1|0;1CN>^CZ=@a#7=b5CUv47 z;@BF-j9iVvB>RJMDc~?qwCk$F@OOeCCF=@WtO%2Spz53*LE~o3F6xkVC?7=$uS?D> z=ROJ{Fe}k>@8l#9sQ=!6tfPvJUX?sq+T<1Fh-gH>0J&(3X6I7z_FS4F}s- zc9)3QQ44y2^u41RI>*sA{^0{OqGbVg-Sxv2lYjbjFujLp=&dwTEyuFnP#*us+ zfLHMPiqb6!&~akahw63?H@ZzF?o84X)Z@Nmm*QK)bi{cUhDeb0U`ryORlj&s8$DT` zLK)6%Z`U8Dkff2hR_6He2p6Sj#uS(6EmZU?5dq{;|L!hc_*wYTskJT5Vg*tO-ptUX zGd4BdM9OpKkMhpV_0c|ay-O|466_jL@DhXRl0RK0G;fW zkWfa~v5Ss#TFlyhvSIuggp$?%g<#N+qup?uWkQOdghd#xt0Rhc)!WjfdSX5L5cL(3@ z-PQ1b<-3fnG;`I^$l^}|9Qr_FxqfIzyURxRsq7~#q0Hh$2b7Y6)HOj)PI&{{qgLIG z`>2r-AvB`)(}~?`V{J|FAPA-f1*nvnh|pH-ybIZqPk>$G3~R-Y4+7I~OX%Tfeuhjj z7S@qtc92*WSspyx+?MFRrp!Lv+dhW&?<#nB{Ln0Z0&NN|4RJ8|!`HRP(|fpj-9z)I zM440tzuK7bI>M)e&K@)qZYN#Fqzt1SJVwGPLGFqXYQ21;ht6h%=jU#Gyv<6pA=cX{ zs)l+vn*%ZWrj|<}_(q}(&F%JBE8WM*9mv~)m+UG!yQ(y--$7B^3x6Pq=q1`hPjA&t zNh@9rlb+vt>T04n;@d|v|8IJFdY0&*TeSwz^hO{t&nUfIg2LLdLWifOf&mz4!Jug9 zHe9)~8nD;wjI_Dqlnc%*V3qdRN_57ml#-q`HBkZ=UYTBiZHo$S+lyb`^UqwLd0QfO zjfQ$>`RHk6L8!ziIcac zso%MREPpReX}O93a74RqccqVPg3t0r!p;hewuVHQ@0`l#)&s4ELQz~X0m6GrSapf+ zej2(DuQk0+}xFN)QKAk< z@TAKy6*NWJ>qW0L5VzD|n-Wzr+^0=$iy0jKJO>UO&$gLm;m96LS#hx8^4m^ky|^&w z9Y8ExIFa0u*hJ2}mAotEMoI8N!YGXMA3ofK9V(B2HV0+lFuDyF z4<8;vsKp4o^xI;?kuJlYmVU?__h{26^Y&=rg-_lG`L#3 zxyO{096;ZbO;^&c25_yqW0OnGGYireWZ*=L8J?Rtit10a4yFc z1k7aM>HyLigB#<7&q~aoT&BB$>(6lgW0X%s}h~*A7kCSK`LyuHLQjoV#dOjuz~#= zd+1=7()d8`fObOcT@gXCjto6si07$kX(urNyUSdOBRM8s?;@faI*trI_Dl*-rIDj> z+zyeZO`8LX+#;4~-Ee9LV#$7Yb2L=M#Kkppv&RzHGmw^Z@9Eb&e=ZIFRzf#YZ&38S zs;Z5ZVsldy%MN*C3kw>=b%f$3?pM0t&Rue&9M>E^^Dq|711eKA}Pir3C=VlP+7gplrLA`%^btT19e!I&Z9&0+lWIFXvXzZ(Gw5N5z9Nm`jcP#(SI|mO_~CQj z-Hy*STJ}s&;v;rsg#>?f$;L?I@-qG*uYjFWI`GIn9>EOIpQza1X{FBZqnog`kZ1Y!rc<}`y| z2?hOnlIRRVOAh9A_(Q+~Z8#OeCFbAPrB*qc&&visx9Z&Go+<(qnQE}a?y9fM` zvzzO6;&IXd9z_u$9xmh6h?;NF0Zx!FFc^a~PjhwzJDpcpsOadJ?eP8wOz9~odX9Ek zAHu0m1jHt)>~p57i{HNwfzqx+fk>Png#=fHvKr26N!N)jMA_2aZ4|X6Am-4$m-UeQ z7PMmBrcQ<86aX@SAm$+5`!lpEYv2?tq~``W!sbxfSK(+1%%X2V00mOU*7`ezP-*b} zipIt)XycsiwhFg9Hl9{r8*J>q{N_j_)9#y7XQr9-+a^e!jtwnW=Rend(w?Ts2&0l?Xb*(UBGa3YK2IPD*-+1L$iY%;%^*h#iS;`!^-$PN~Q2J@R#autsKST8uoPadbg_JXlH!4P3Z$57C29GujW*Lz07}Kq{TZ!XJa9$WB zRsRKQlj~sR2uc;o=66)E?Xk|scC>aQ@>`)+Cc3bc!dUHYlrr^pWEqN}x+9(e$l7V{ z+lv}vPr^&E2F=M?a*S?DO6wVfOalSd#D5%KJ3GYR9?U4544&y|gWu%mrLr;QSd}%u zQvCMLbySeGxOD)&!3zwrUR=A`yXDY#HV?(@RXJ#!neC+M7y0#DY9yPzbi3~=EnwcadAT9(=q6L{AjkXug?ms98Q4a#UoK0Gs^c=cppT=#vTa@1GV`HU0rZn z-A|m@(Q29IYC434H9RvTC}EeCgd-%aJVBnw@e9q|0bTuP}iQrKlqG(i*D6bK3BCrnA(2#WY!5z(9q$A9RK^&GX<#bk_VK_2b|_Apx9RH*TGdv;Eq^k`c?l8<-c-4ejt-*x7UU{%TZa z_g;&L5W*y~8mE?}`jQ30aScwOH55zL%o9O@A%q24B8}W#e|iyrl~r3rCE= zG`MoAlJtn*JB50!9|s)6L@Z%d2%LWgt^m>?L`(@+eHaOF&*@oNdPW<{`HAB+RA>kf zK6FP+EWErdIGr&sFYk0gW%|_0**1QenY$xnA~X0`Os$`h;0VP7N9x_YNrxI%qUWYQ z4Plo^0z-W2D;;*A(Qry&T13cbpFGCJa2;-j=!AvF`J>Z=BAM9>-H~;kGb%r7#Qf|V|wB=|2%K! zkMk&r<8=o6diSiVuWi1JxJbGjER3$z!oxV|!ukC9&eQC^R|Se#)S+;I^&Xm-=q8_{ z7QznJ=7aVAN@QeYQkK;|_;s1K6RzDqUU7;83M51#!k1=#CP>hHy){a26!*-fEJ`cq zTY7Q-^RNH>%=OSTK{H)eQ!RKp{(5*PGjZr1HsA(GtMZ0*>3=spopO*YhHn$cmHctM z*k3;jZSAiIf41)-n*4DH1F;NYUV|x*IB%$TX%OaU9M<4`;X(k4W3*YG#d!|zv+@!q z=%R5_0&!l+ucnq3r-ADdkDy@>gM2vBm171rlLkUFaoi#>EP2K$^INdcS&^}&btFPl zLz;Sw-fEC7@L&*K+hMTfYGnV?SV(xoEtGbtYd_Gds?>Z~XoHli2vM&QXKp zBFLy=DA_ED(PCYCj|j>KIz9Y zUX>GcjQo1#7NngA4{pa=fR{lEv#E`9i9i(&qfM$7fH?>UK}9XaewI#AGfh=rNVwPy zZ@?ThFr1NMo;`z8|A>Qyu$zc;5{XlPe10!kARDYjrUK2K^@BHhi#V|hpM%oDVR#BI z)7>H)NzYeYW*1OM(~)jM*#Xf)Ma7Mqk5ErIRB|9QOQN*a2CyFQVm8FBsHsN*jBeuG zt&#$?+Xf+nzPocO!Oq!2ufoDLD42gxI1?Jh(dz1IE?n-6^Fh6eQ@3;^gt8|S7P~5! zxtXP0)2$nC?7E2>XMgE03=cPn3x{?*c{q*iM+ZmX^YV*d_DM;gL_1FS^-G=OLT+6jB-<83a8n?WpEuy zWPfokZmqO-%rfnKW%@}thae|mLd1*X`3J?C(zizk-1#gWiF2>uSM>TZ`GEW>5+w$F zBgr@qp`ZZLjsiF(&^j#fiLgG!BPkKbFH#Cre8Yh_*!1M^B@um;g*eOQN__l|1h3Dk z=)0?fp@#u=@49vC2zQWkStq`JJLbF6f-pxSbrM|a_fI~*gNfrbcJI!t@;0uKPSL<5 zFD7JUW!1pVeI0cm8#Hn9j2_k@3>){f5@jeq0UFXkCZRzrY0p{H;2!a&1Uwj<&(E{| zE>2EPq{HUs3BBx0dMClZfrDtg#Fv1|g>Z;tO=H`s0^P|Q=Bs)mA>zE#IxPCvL+7(+ zE)ml;b}0z$xWw@WB%(3_)07V+yZ)OU^|z4yXYu+11p&dDqOQw$+Jd>(hJmz9mhJ1T zbQ?;d!vC+cGmon|U*q_ZEn%2SN@OjONm)jsD{4ZKn8FxKBvi6QCQ1&C%A|{wZL&nA zX0p^MX|XkRTFlsr5yw&)TN0y1r2Bqm?(6<@mp|@byySF#=lA=5zt8i0w)1z@&k@_5 zJv|d{)TPpKE`D{^fq6lgQO{G{#t4kCXUMV z^~qW=Vw;oCxdljRT12fQ(F5RDf`T15@UP3UtP4yej1h?GBKFWnC==>~Z@iC7t}$3Q zGU|G2db*o+0k(6;`oGz9rH4})1x`73;r{k{w_4tM)^KG0}ltjiN{jus=*)na+4*H}028hCqf&PcOP3j&s-F+*L zem6)xw3ueH6*tIdRqRZk%%w&N&U0NYw+8;sP%p%NE$Ih+2}d5a;@(9hpky$#g$>24 znUv!IHOGz_b;l-8Zao^GRXk**D5w6=raBS%()H7~d(vWLFa7gO8jcszPO!sSl#)kJo`hk@y%IZGM*&|L z!W9+09O`QsWc4l!+ZEa)MB2!q+h|W66SS@L%#vSvh@%I-MpAqQFvgFrp%ASB7A~&Y zppi*g!b)+XD*|duOmJp^4eb8j%^H%1NBQOZqE9^%eCJxPs}AS2N7=m;`~o21RE)y1 zQm7zcBqJWWJ^OuJ%eyYZVSt_$cCLvzVC_FxaVG2i7ke}Jy=m4^faTkfW+l74v8kya z;v&(2fo^>ZHBv$5qvBdqjlsN%Mb8_y52GkcHDlbYD*nW1tm!DDlc2O5G@|$tbrf+- zi^_g*qmvUU@JOFeOP6K2_3I40=%g@h&I)eMq5+$fIPI``!<^?pS5fz=ojlgY##m2JFA5aW5nZ9lI}n4a3?JWy?UeYsKnDZQ;kTW} zeY1nA6$K17iJCnChj`S$Lje=~9v63uyjmm?qqy+I{3w6^cdG+G{kLH9{}3|WISc9l zgO*_RjGP8p*4Bu~_MBLqv{D9dK&Ou%{+M>|4U_>A@&4n+J?K;Sm(tR{bNShohpQj3 zaR#3)xpKvx|0@*NM%M_VC}wQP0unh4vj*FjYIc+$_3g6|o_^BPN}-S>5~})Y$Pa9Z z`WcI?f6H9mf5?#SFD5=1^D0NGH_Wzyy!486V1uzJbmA;|TApskGHh*ozHN#zAy6Se zf^k3gRpsTi4h{}>*xOZm`pwGYvEi=eFcb=L6K`yKdr54nFbNf?AflYKDkgw>2)*l(Sw$cXRGT7E3BoXs$d33O{m1jz{wwmfXehtqv#~X*$?Ir%9C8k>J%`F z-u^n_9yzs_JJvQdOa$8!u{cqIx_Ww-)#_g8)Kp(ivq@Mr74%BLL1>*HYv+?5jYUGp z2oB2j_4W10H$~OVp6eo3fr>rQi9#5$i$h{!y3u4TenLzHMF3#euk2dK-F8W&T{*_zJ`}zjoXubOTqvYV%?LJA8_z|=F@+w zp*r=n>dF*HL{G@vg8)6R+3RqrJ%*ld>?dgoX#0tXxkPL;xe6o1YVvjiEmSI5Jw16k zB;8=R8_88!xL;HVjlsuC^EKA)fW9?#=~bh-xf{B{f1##9Vco{WIaU)iD1DK0=?oDj zCP}AuL&}^ZFL<%H?-CRUc7Ftb#vS(KcS=zTXLiMBLdNhyz|RYR{j~s;dlmaCh>a-HLOwRl zZx;?KX@GrhrEE?}Fh;4GaI4M;R0vEjZ1#{)(oG1j`W75W_(rK!7JgNewpQ8imqdrd z8|w`iwV0&aagGZ@cPi(NmPDG)oTqGGRN^!=cYxDcSy||_b#-;!5gFo?+YF@-6k-uj z>n_9Yvb!0dNet{52OV8Q>RN=(JLV}X2(pgj*{x4>Bg2Np@e!N1-Wqo+i z(%hTphQKre%#(b+R$^;m)zOlKqLRLdx#MLfJ_5MlcvsLa$Ty2s4{ec^wNh}hh4SD; z7(}utAk0JPPf={bTc*w26V^>nuL4x4tpby(KJ*y)$KF@c`5z8$>cBl#5y%sFdL_nD zlc@*31s42-p_kUuY-Wlrq|5!__r7>vB$1i+os=^$LP0Fg6FQi#eF|8zf?e!D>89lm zLlIDXfA^#Tk^-NSlG4M*8L|o+V;3aX;wlmlSJ~dL>(ZMHL8=gaV|+C>Mw@uQL);J(?2gOEp>;R5*^x3dTCCSLDhTN!QtVW zId!+bPSY1Le8-$^x49%vkPsX#Ut#;Cy6v+G%y3XI;kZuoSsx94 zQBYaRavunhrQ|G=atUTuNob5%+7#P?(U}q`!A6S253*A8Et(VA^ z6-o|=QCYrr@&3fcx&&N>i=(TZ`F=L#sf|MH*KpJ02`X~-Z>^Mx2ikLX_PWm`Izb3! zP~}h#6u&iVVrp(eI(-%KBu|9pwtqW+;OE!$y~yo$ITaW!@%4>lqa962nc}GqpVJTC zch?qnE8YlELTG(v0|i)ESa=ajUYQXaiGU7;0i88-($Q`EykzF6DSOC1xr38XupK&G|U=W}y&tBZbR zASx_QNIeBanYCd3`q9c3<1xY!hS>MJo;-iX4Xcmz2azl8vfRYPNhbhk=Nt8hPI?=0 zW+k3H>AXk(3%&>9o;mpou(5?jMUpG4Y>CXH4RvfVV17dwA@Wy*Dy}oX_)Uh<V1`r02o9F%A_(wR5_`lXBeSB3OD;6JwvLYB!5Np~*2L<8FEs=;PhajbeRe(sqSTqB zB-8XKZycvhQ>l7y>cjxe%P+GLVWvdWBiD=a=uQSvU{7xz?;@HriL0YfgF4&N9!fAL z@!maQa;*HSG$ox*h#lZHaqS?w3a6s_lyz=0g~i}Q^cL)Z=0u{QIN1ck=3ErCfqz^) z&C;1DGByc!Cy%vqJ68-A1{0j=_gD+o0kFdS`?F?rl%#s38HX^*wc*G{QswY0YJ8I}Ul{<| zEp{e?EXisK!UTS7$2Kl}Sown8qug>cMvL3;F!2eb9 zmMpX_C|*N58Ej#UIzmj7-`xvn;iOj+r~s|Q5tKUA5=nq%$~w*L8^k#y_Y5C&Cro(a zeZoToNSS34=j{TKA>jus`G-J6oLovAP1QkN#=Ty~dJ*+7BSryBZ8ffH>2;agzOZ&sY%a3{n4NBv?v=o!SaJD}= I!F|g=0P2A44gdfE literal 0 HcmV?d00001 diff --git a/documentation/source/reference/Examples/GroundStateEnergyQPE.rst b/documentation/source/reference/Examples/GroundStateEnergyQPE.rst new file mode 100644 index 00000000..6678ffe6 --- /dev/null +++ b/documentation/source/reference/Examples/GroundStateEnergyQPE.rst @@ -0,0 +1,54 @@ +.. _GroundStateEnergyQPE: + +Ground State Energy with QPE +============================ + +.. currentmodule:: qrisp + +Example Hydrogen +================ + +We caluclate the ground state energy of the Hydrogen molecule with :ref:`Quantum Phase Estimation `. + +Utilizing symmetries, one can find a reduced two qubit Hamiltonian for the Hydrogen molecule. Frist, we define the Hamiltonian, and compute the +ground state energy classically. + +:: + + from qrisp import QuantumVariable, x, QPE + from qrisp.operators.pauli.pauli import X,Y,Z + import numpy as np + + # Hydrogen (reduced 2 qubit Hamiltonian) + H = -1.05237325 + 0.39793742*Z(0) -0.39793742*Z(1) -0.0112801*Z(0)*Z(1) + 0.1809312*X(0)*X(1) + E0 = H.ground_state_energy() + # Yields: -1.85727502928823 + +In the following, we utilize the ``trotterization`` method of the :ref:`PauliHamiltonian` to obtain a function **U** that applies **Hamiltonian Simulation** +via Trotterization. If we start in a state that is close to the ground state and apply :ref:`QPE`, we get an estimate of the ground state energy. + +:: + + # ansatz state + qv = QuantumVariable(2) + x(qv[0]) + E1 = H.get_measurement(qv) + E1 + # Yields: -1.83858104676077 + + +We calculate the ground state energy with quantum phase estimation. As the results of the phase estimation are modulo :math:`2\pi`, we subtract :math:`2\pi`. + +:: + + U = H.trotterization() + + qpe_res = QPE(qv,U,precision=10,kwargs={"steps":3},iter_spec=True) + + results = qpe_res.get_measurement() + sorted_results= dict(sorted(results.items(), key=lambda item: item[1], reverse=True)) + + phi = list(sorted_results.items())[0][0] + E_qpe = 2*np.pi*(phi-1) # Results are modulo 2*pi, therefore subtract 2*pi + E_qpe + # Yields: -1.8591847149174 \ No newline at end of file diff --git a/documentation/source/reference/Examples/MolecularPotentialEnergyCurve.rst b/documentation/source/reference/Examples/MolecularPotentialEnergyCurve.rst index acc897a3..101ad263 100644 --- a/documentation/source/reference/Examples/MolecularPotentialEnergyCurve.rst +++ b/documentation/source/reference/Examples/MolecularPotentialEnergyCurve.rst @@ -1,18 +1,26 @@ -.. _MolecularPotentialEnergyCurveExample: +.. _MolecularPotentialEnergyCurve: Molecular Potential Energy Curves ================================= .. currentmodule:: qrisp.vqe +Molecular potential energy curves and surfaces are tools in computational quantum chemistry for analysing the molecular geometries and chemical reactions. +The potential energy is defined as + +$$E_{\\text{pot}}(\\{R_A\\})=E_{\\text{elec}}(\\{R_A\\})+E_{\\text{nuc}}(\\{R_A\\})$$ + +where $E_{\text{elec}}(\{R_A\})$ is the electronic energy and $E_{\text{nuc}}(\{R_A\})$ is the nuclear repulsion energy depending on the coordinates of the atoms $\{R_A\}$. + + Example Hydrogen ================ -We caluclate the Potetial Energy Curve for the Hydrogen molecule for varying interatomic distance of the atoms. +We caluclate the potential energy curve for the Hydrogen molecule for varying interatomic distance of the atoms. We implement a function ``problem_data`` that sets up a molecule and -utilizes the `PySCF `_ quantum chemistry library to compute the :meth:`electronic data ` (one- and two-electron integrals, number of orbitals, number of electroins, nuclear repulsion energy, Hartree-Fock energy) +utilizes the `PySCF `_ quantum chemistry library to compute the :meth:`electronic data ` (one- and two-electron integrals, number of orbitals, number of electrons, nuclear repulsion energy, Hartree-Fock energy) for a given molecular geometry. In this example, we vary the **atomic distance** between the two Hydrogen atoms. :: @@ -56,8 +64,8 @@ Finally, we visualize the results. :: import matplotlib.pyplot as plt - plt.scatter(x, y_hf, color='#7d7d7d',marker="o", linestyle='solid',s=10, label='HF energy') - plt.scatter(x, y_qccsd, color='#6929C4',marker="o", linestyle='solid',s=10, label='VQE (QCCSD) energy') + plt.scatter(distances, y_hf, color='#7d7d7d',marker="o", linestyle='solid',s=10, label='HF energy') + plt.scatter(distances, y_qccsd, color='#6929C4',marker="o", linestyle='solid',s=10, label='VQE (QCCSD) energy') plt.xlabel("Atomic distance", fontsize=15, color="#444444") plt.ylabel("Energy", fontsize=15, color="#444444") plt.legend(fontsize=12, labelcolor="#444444") @@ -68,17 +76,86 @@ Finally, we visualize the results. .. figure:: /_static/H2_PEC.png :alt: Hydrogen Potential Energy Curve - :scale: 80% :align: center Example Beryllium hydride ========================= -We caluclate the Potetial Energy Curve for the Beryllium hydride molecule for varying interatomic distance of the atoms. +We caluclate the potential energy curve for the Beryllium hydride molecule for varying interatomic distance of the atoms. +As in the previous example, we set up a function that computes the :meth:`electronic data ` for +the Beryllium hydride molecule for varying **atomic distance** between the Beryllium and Hydrogen atoms. +:: + from pyscf import gto + from qrisp.vqe.problems.electronic_structure import * + from qrisp import QuantumVariable + def problem_data(r): + mol = gto.M( + atom = f'''Be 0 0 0; H 0 0 {r}; H 0 0 {-r}''', + basis = 'sto-3g') + return electronic_data(mol) +We further investigate the problem size: +:: + + data = problem_data(2.0) + print(data['num_orb']) + print(data['num_elec']) + + H = create_electronic_hamiltonian(data).to_pauli_hamiltonian() + print(H.len()) + +In the chosen sto-3g basis, there are 14 molecular orbitals (qubits) and 6 electrons. The problem Hamiltonian has 666 terms. +For `reducing the problem size `_, an active space reduction is applied: We consider 6 **active orbitals** and 2 **active electrons**, that is, + +* 4 electrons occupy the 4 lowest energy molecular orbitals + +* 2 electrons are distributed among the subsequent 4 molecular orbitals (quantum optimization) + +* the remaining (highest energy) 6 orbitals are not occupied + +This lowers the quantum resource requirements to 4 qubits, and a reduced 4 qubit Hamiltonian is considered. + +.. warning:: + + The following code may well take more than 10 minutes to run! + +:: + + distances = np.arange(0.2, 4.0, 0.1) + y_hf = [] + y_qccsd = [] + + for r in distances: + data = problem_data(r) + y_hf.append(data['energy_hf']) + vqe = electronic_structure_problem(data,active_orb=4,active_elec=2) + results = [] + for i in range(5): + res = vqe.run(QuantumVariable(6),depth=2,max_iter=100,optimizer='COBYLA',mes_kwargs={'precision':0.005}) + results.append(res+data['energy_nuc']) + y_qccsd.append(min(results)) + +Finally, we visualize the results. + +:: + + import matplotlib.pyplot as plt + plt.scatter(distances[5:], y_hf[5:], color='#7d7d7d',marker="o", linestyle='solid',s=10, label='HF energy') + plt.scatter(distances[5:], y_qccsd[5:], color='#6929C4',marker="o", linestyle='solid',s=10, label='VQE (QCCSD) energy') + plt.xlabel("Atomic distance", fontsize=15, color="#444444") + plt.ylabel("Energy", fontsize=15, color="#444444") + plt.legend(fontsize=12, labelcolor="#444444") + plt.tick_params(axis='both', labelsize=12) + plt.grid() + plt.show() + + +.. figure:: /_static/BeH2_PEC_4_2.png + :alt: Beryllium Hydrate Potential Energy Curve + :align: center diff --git a/documentation/source/reference/Examples/index.rst b/documentation/source/reference/Examples/index.rst index 7236c68b..a42f6ddc 100644 --- a/documentation/source/reference/Examples/index.rst +++ b/documentation/source/reference/Examples/index.rst @@ -35,9 +35,11 @@ In this section, we provide a glimpse into the diverse range of applications tha - Provides implementations for solving optimization problems using the Quantum Approximate Optimization Algorithm. For a detailed tutorial on how to use QAOA, please refer to the :ref:`in depth QAOA tutorial `. * - :ref:`ShorExample` - Showcases the cryptographic implications of implementating Shor's algorithm and provides insight in how to easily use a custom adder. - * - :ref:`MolecularPotentialEnergyCurveExample` - - Molecular Potetial Energy Curve - + * - :ref:`MolecularPotentialEnergyCurve` + - Calculating molecular potential energy curves with :ref:`VQE`. + * - :ref:`GroundStateEnergyQPE` + - Calculating ground state energies with quantum phase estimation. + .. toctree:: :hidden: @@ -57,3 +59,4 @@ In this section, we provide a glimpse into the diverse range of applications tha QAOA Shor MolecularPotentialEnergyCurve + GroundStateEnergyQPE From 65e80e767423cf0c3714bf1e9441d5e5a330d06e Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Wed, 30 Oct 2024 15:38:53 +0100 Subject: [PATCH 06/13] Restructuring get_measurement --- .../fermionic/fermionic_hamiltonian.py | 57 ++- src/qrisp/operators/hamiltonian.py | 340 +---------------- src/qrisp/operators/hamiltonian_tools.py | 11 +- src/qrisp/operators/measurement.py | 145 -------- .../pauli/bound_pauli_hamiltonian.py | 99 +++++ src/qrisp/operators/pauli/measurement.py | 346 ++++++++++++++++++ .../operators/pauli/pauli_hamiltonian.py | 98 ++++- 7 files changed, 608 insertions(+), 488 deletions(-) delete mode 100644 src/qrisp/operators/measurement.py create mode 100644 src/qrisp/operators/pauli/measurement.py diff --git a/src/qrisp/operators/fermionic/fermionic_hamiltonian.py b/src/qrisp/operators/fermionic/fermionic_hamiltonian.py index 6377ed74..89e203ac 100644 --- a/src/qrisp/operators/fermionic/fermionic_hamiltonian.py +++ b/src/qrisp/operators/fermionic/fermionic_hamiltonian.py @@ -508,9 +508,64 @@ def to_pauli_hamiltonian(self, mapping_type='jordan_wigner', num_qubits=None): return H # - # Partitions + # Measurement # + def get_measurement( + self, + qarg, + precision=0.01, + backend=None, + shots=1000000, + compile=True, + compilation_kwargs={}, + subs_dic={}, + precompiled_qc=None, + _measurement=None # measurement settings + ): + r""" + This method returns the expected value of a Hamiltonian for the state of a quantum argument. + + Parameters + ---------- + qarg : QuantumVariable, QuantumArray or list[QuantumVariable] + The quantum argument to evaluate the Hamiltonian on. + precision: float, optional + The precision with which the expectation of the Hamiltonian is to be evaluated. + The default is 0.01. The number of shots scales quadratically with the inverse precision. + backend : BackendClient, optional + The backend on which to evaluate the quantum circuit. The default can be + specified in the file default_backend.py. + shots : integer, optional + The maximum amount of shots to evaluate the expectation of the Hamiltonian. + The default is 1000000. + compile : bool, optional + Boolean indicating if the .compile method of the underlying QuantumSession + should be called before. The default is True. + compilation_kwargs : dict, optional + Keyword arguments for the compile method. For more details check + :meth:`QuantumSession.compile `. The default + is ``{}``. + subs_dic : dict, optional + A dictionary of Sympy symbols and floats to specify parameters in the case + of a circuit with unspecified, :ref:`abstract parameters`. + The default is {}. + precompiled_qc : QuantumCircuit, optional + A precompiled quantum circuit. + + Raises + ------ + Exception + If the containing QuantumSession is in a quantum environment, it is not + possible to execute measurements. + + Returns + ------- + float + The expected value of the Hamiltonian. + """ + pass + # # Trotterization # diff --git a/src/qrisp/operators/hamiltonian.py b/src/qrisp/operators/hamiltonian.py index aea9ba6a..51207091 100644 --- a/src/qrisp/operators/hamiltonian.py +++ b/src/qrisp/operators/hamiltonian.py @@ -16,9 +16,6 @@ ********************************************************************************/ """ -from qrisp import QuantumVariable, QuantumArray -from qrisp.core.compilation import qompiler -from qrisp.operators.measurement import * from sympy import init_printing # Initialize automatic LaTeX rendering init_printing() @@ -229,15 +226,11 @@ def ground_state_energy(self): """ pass - # - # Evaluate expected value - # - + @abstractmethod def get_measurement( self, qarg, precision=0.01, - method='QWC', backend=None, shots=1000000, compile=True, @@ -245,7 +238,7 @@ def get_measurement( subs_dic={}, circuit_preprocessor=None, precompiled_qc=None, - _measurement=None + _measurement=None # measurement settings ): r""" This method returns the expected value of a Hamiltonian for the state of a quantum argument. @@ -254,13 +247,9 @@ def get_measurement( ---------- qarg : QuantumVariable, QuantumArray or list[QuantumVariable] The quantum argument to evaluate the Hamiltonian on. - method : string, optional - The method for evaluating the expected value of the Hamiltonian. - Available is ``QWC``: Pauli terms are grouped based on qubit-wise commutativity. - The default is ``QWC``. precision: float, optional The precision with which the expectation of the Hamiltonian is to be evaluated. - The default is 0.01. + The default is 0.01. The number of shots scales quadratically with the inverse precision. backend : BackendClient, optional The backend on which to evaluate the quantum circuit. The default can be specified in the file default_backend.py. @@ -278,9 +267,8 @@ def get_measurement( A dictionary of Sympy symbols and floats to specify parameters in the case of a circuit with unspecified, :ref:`abstract parameters`. The default is {}. - circuit_preprocessor : Python function, optional - A function which recieves a QuantumCircuit and returns one, which is applied - after compilation and parameter substitution. The default is None. + precompiled_qc : QuantumCircuit, optional + A precompiled quantum circuit. Raises ------ @@ -301,7 +289,7 @@ def get_measurement( :: from qrisp import QuantumVariable, h - from qrisp.operators import X,Y,Z + from qrisp.operators.pauli import X,Y,Z qv = QuantumVariable(2) h(qv) H = Z(0)*Z(1) @@ -314,7 +302,7 @@ def get_measurement( :: from qrisp import QuantumVariable, QuantumArray, h - from qrisp.operators import X,Y,Z + from qrisp.operators.pauli import X,Y,Z qtype = QuantumVariable(2) q_array = QuantumArray(qtype, shape=(2)) h(q_array) @@ -324,319 +312,7 @@ def get_measurement( #Yields 1.0 """ - - from qrisp import QuantumSession, merge - - if isinstance(qarg,QuantumVariable): - if qarg.is_deleted(): - raise Exception("Tried to get measurement from deleted QuantumVariable") - qs = qarg.qs - - elif isinstance(qarg,QuantumArray): - for qv in qarg.flatten(): - if qv.is_deleted(): - raise Exception( - "Tried to measure QuantumArray containing deleted QuantumVariables" - ) - qs = qarg.qs - elif isinstance(qarg,list): - qs = QuantumSession() - for qv in qarg: - if qv.is_deleted(): - raise Exception( - "Tried to measure QuantumArray containing deleted QuantumVariables" - ) - merge(qs,qv.qs) - - if backend is None: - if qs.backend is None: - from qrisp.default_backend import def_backend - - backend = def_backend - else: - backend = qarg.qs.backend - - if len(qs.env_stack) != 0: - raise Exception("Tried to get measurement within open environment") - - - # Copy circuit in over to prevent modification - if precompiled_qc is None: - if compile: - qc = qompiler( - qs, **compilation_kwargs - ) - else: - qc = qs.copy() - else: - qc = precompiled_qc.copy() - - # Bind parameters - if subs_dic: - qc = qc.bind_parameters(subs_dic) - from qrisp.core.compilation import combine_single_qubit_gates - qc = combine_single_qubit_gates(qc) - - # Execute user specified circuit_preprocessor - if circuit_preprocessor is not None: - qc = circuit_preprocessor(qc) - - qc = qc.transpile() - - from qrisp.misc import get_measurement_from_qc - - if _measurement is None: - pauli_measurement = self.pauli_measurement() - else: - pauli_measurement = _measurement - - meas_circs = pauli_measurement.circuits - meas_qubits = pauli_measurement.qubits - meas_ops = pauli_measurement.operators_int - meas_coeffs = pauli_measurement.coefficients - meas_shots = pauli_measurement.shots - - meas_shots = [round(x/precision**2) for x in meas_shots] - tot_shots = sum(x for x in meas_shots) - if tot_shots>shots: - #meas_shots = [round(x*shots/tot_shots) for x in meas_shots] - raise Warning("Warning: The total amount of shots required: " + str(tot_shots) +" for the target precision: " + str(precision) + " exceeds the allowed maxium amount of shots. Decrease the precision or increase the maxium amount of shots.") - - N = len(meas_circs) - - expectation = 0 - - for k in range(N): - - curr = qc.copy() - curr.append(meas_circs[k].to_gate(), meas_qubits[k]) - - res = get_measurement_from_qc(curr, meas_qubits[k], backend, meas_shots[k]) - - # Groupings - M = len(meas_ops[k]) - for l in range(M): - for outcome,probability in res.items(): - expectation += probability*evaluate_observable(meas_ops[k][l],outcome)*meas_coeffs[k][l] - - return expectation - -# -# NUMBA -# - - def get_measurement_new( - self, - qarg, - precision=0.01, - method='QWC', - backend=None, - shots=1000000, - compile=True, - compilation_kwargs={}, - subs_dic={}, - circuit_preprocessor=None, - precompiled_qc=None, - _measurement=None - ): - r""" - This method returns the expected value of a Hamiltonian for the state of a quantum argument. - - Parameters - ---------- - qarg : QuantumVariable, QuantumArray or list[QuantumVariable] - The quantum argument to evaluate the Hamiltonian on. - method : string, optional - The method for evaluating the expected value of the Hamiltonian. - Available is ``QWC``: Pauli terms are grouped based on qubit-wise commutativity. - The default is ``QWC``. - precision: float, optional - The precision with which the expectation of the Hamiltonian is to be evaluated. - The default is $\epsilon=0.01$. The number of shots scales as $\mathcal O(\epsilon^{-2})$. - backend : BackendClient, optional - The backend on which to evaluate the quantum circuit. The default can be - specified in the file default_backend.py. - shots : integer, optional - The maximum amount of shots to evaluate the expectation of the Hamiltonian. - The default is 1000000. - compile : bool, optional - Boolean indicating if the .compile method of the underlying QuantumSession - should be called before. The default is True. - compilation_kwargs : dict, optional - Keyword arguments for the compile method. For more details check - :meth:`QuantumSession.compile `. The default - is ``{}``. - subs_dic : dict, optional - A dictionary of Sympy symbols and floats to specify parameters in the case - of a circuit with unspecified, :ref:`abstract parameters`. - The default is {}. - circuit_preprocessor : Python function, optional - A function which recieves a QuantumCircuit and returns one, which is applied - after compilation and parameter substitution. The default is None. - - Raises - ------ - Exception - If the containing QuantumSession is in a quantum environment, it is not - possible to execute measurements. - - Returns - ------- - float - The expected value of the Hamiltonian. - - Examples - -------- - - We define a Hamiltonian, and measure its expected value for the state of a :ref:`QuantumVariable`. - - :: - - from qrisp import QuantumVariable, h - from qrisp.operators import X,Y,Z - qv = QuantumVariable(2) - h(qv) - H = Z(0)*Z(1) - res = H.get_measurement(qv) - print(res) - #Yields 0.0 - - We define a Hamiltonian, and measure its expected value for the state of a :ref:`QuantumArray`. - - :: - - from qrisp import QuantumVariable, QuantumArray, h - from qrisp.operators import X,Y,Z - qtype = QuantumVariable(2) - q_array = QuantumArray(qtype, shape=(2)) - h(q_array) - H = Z(0)*Z(1) + X(2)*X(3) - res = H.get_measurement(q_array) - print(res) - #Yields 1.0 - - """ - - from qrisp import QuantumSession, merge - - if isinstance(qarg,QuantumVariable): - if qarg.is_deleted(): - raise Exception("Tried to get measurement from deleted QuantumVariable") - qs = qarg.qs - - elif isinstance(qarg,QuantumArray): - for qv in qarg.flatten(): - if qv.is_deleted(): - raise Exception( - "Tried to measure QuantumArray containing deleted QuantumVariables" - ) - qs = qarg.qs - elif isinstance(qarg,list): - qs = QuantumSession() - for qv in qarg: - if qv.is_deleted(): - raise Exception( - "Tried to measure QuantumArray containing deleted QuantumVariables" - ) - merge(qs,qv.qs) - - if backend is None: - if qs.backend is None: - from qrisp.default_backend import def_backend - - backend = def_backend - else: - backend = qarg.qs.backend - - if len(qs.env_stack) != 0: - raise Exception("Tried to get measurement within open environment") - - - # Copy circuit in over to prevent modification - if precompiled_qc is None: - if compile: - qc = qompiler( - qs, **compilation_kwargs - ) - else: - qc = qs.copy() - else: - qc = precompiled_qc.copy() - - # Bind parameters - if subs_dic: - qc = qc.bind_parameters(subs_dic) - from qrisp.core.compilation import combine_single_qubit_gates - qc = combine_single_qubit_gates(qc) - - # Execute user specified circuit_preprocessor - if circuit_preprocessor is not None: - qc = circuit_preprocessor(qc) - - qc = qc.transpile() - - from qrisp.misc import get_measurement_from_qc - - if _measurement is None: - pauli_measurement = self.pauli_measurement() - else: - pauli_measurement = _measurement - - meas_circs = pauli_measurement.circuits - meas_qubits = pauli_measurement.qubits - meas_ops = pauli_measurement.operators_int - meas_coeffs = pauli_measurement.coefficients - meas_shots = pauli_measurement.shots - - meas_shots = [round(x/precision**2) for x in meas_shots] - tot_shots = sum(x for x in meas_shots) - if tot_shots>shots: - #meas_shots = [round(x*shots/tot_shots) for x in meas_shots] - raise Warning("Warning: The total amount of shots required: " + str(tot_shots) +" for the target precision: " + str(precision) + " exceeds the allowed maxium amount of shots. Decrease the precision or increase the maxium amount of shots.") - - N = len(meas_circs) - - # Perform the measurements - results = [] - for k in range(N): - - curr = qc.copy() - curr.append(meas_circs[k].to_gate(), meas_qubits[k]) - - res = get_measurement_from_qc(curr, meas_qubits[k], backend, meas_shots[k]) - results.append(res) - - # Partition outcomes in uint64 - #outcomes_parts = [partition(result.keys(), len(meas_qubits[k])) for k, result in enumerate(results)] - outcomes_parts = [partition(list(results[k].keys()), len(meas_qubits[k])) if len(meas_qubits[k])>64 else [np.array(list(results[k].keys()), dtype=np.uint64)] for k in range(N)] - - probabilities = [np.array(list(result.values()), dtype=np.float64) for result in results] - - # Partition observables in uint64 - observables_parts = [partition(observable, len(meas_qubits[k])) if len(meas_qubits[k])>64 else [np.array(observable, dtype=np.uint64)] for k, observable in enumerate(meas_ops)] - - coefficients = [np.array(coeffs, dtype=np.float64) for coeffs in meas_coeffs] - - # Evaluate the observable - expectation = 0 - - for k in range(N): - - result = evaluate_observables_parallel(observables_parts[k], outcomes_parts[k], probabilities[k]) - #print(result) - #print(coefficients[k]) - - expectation += result @ coefficients[k] - - - # Groupings - #M = len(meas_ops[k]) - # for outcome,probability in res.items(): - # expectation += probability*evaluate_observable(meas_ops[k][l],outcome)*meas_coeffs[k][l] - - - - return expectation + pass \ No newline at end of file diff --git a/src/qrisp/operators/hamiltonian_tools.py b/src/qrisp/operators/hamiltonian_tools.py index cc5ca708..795c6548 100644 --- a/src/qrisp/operators/hamiltonian_tools.py +++ b/src/qrisp/operators/hamiltonian_tools.py @@ -16,10 +16,6 @@ ********************************************************************************/ """ -from qrisp.operators.hamiltonian import Hamiltonian -from qrisp.operators.pauli.pauli_measurement import PauliMeasurement - - def multi_hamiltonian_measurement( hamiltonians, qarg, @@ -29,7 +25,6 @@ def multi_hamiltonian_measurement( compile=True, compilation_kwargs={}, subs_dic={}, - circuit_preprocessor=None, precompiled_qc=None, _measurements=None ): @@ -62,9 +57,8 @@ def multi_hamiltonian_measurement( A dictionary of Sympy symbols and floats to specify parameters in the case of a circuit with unspecified, :ref:`abstract parameters`. The default is {}. - circuit_preprocessor : Python function, optional - A function which recieves a QuantumCircuit and returns one, which is applied - after compilation and parameter substitution. The default is None. + precompiled_qc : QuantumCircuit, optional + A precompiled quantum circuit. Returns ------- @@ -83,7 +77,6 @@ def multi_hamiltonian_measurement( compile=compile, compilation_kwargs=compilation_kwargs, subs_dic=subs_dic, - circuit_preprocessor=circuit_preprocessor, precompiled_qc=precompiled_qc, _measurement=None if _measurements==None else _measurements[i] )) diff --git a/src/qrisp/operators/measurement.py b/src/qrisp/operators/measurement.py deleted file mode 100644 index 9f4f8297..00000000 --- a/src/qrisp/operators/measurement.py +++ /dev/null @@ -1,145 +0,0 @@ -""" -\******************************************************************************** -* Copyright (c) 2024 the Qrisp authors -* -* This program and the accompanying materials are made available under the -* terms of the Eclipse Public License 2.0 which is available at -* http://www.eclipse.org/legal/epl-2.0. -* -* This Source Code may also be made available under the following Secondary -* Licenses when the conditions for such availability set forth in the Eclipse -* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 -* with the GNU Classpath Exception which is -* available at https://www.gnu.org/software/classpath/license.html. -* -* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 -********************************************************************************/ -""" - -import math -import numpy as np -from numba import njit, prange - - -def partition(values, num_qubits): - """ - Partitions a list of integers into a list of lists of integers with size 64 bit. - - Parameters - ---------- - values : list[int] - A list of integers. - num_qubits : int - The maximal number of bits. - - Returns - ------- - partition : list[np.array] - A list of NumPy numpy.uint64 arrays. - - """ - - - M = math.ceil(num_qubits/64) - N = len(values) - - zeros = [[0]*N for k in range(M)] - partition = [values] - partition.extend(zeros) - - lower_mask = (1<<64) - 1 - - for j in range(1,M): - for i in range(N): - partition[j][i] = partition[j-1][i] & lower_mask - partition[j+1][i] = partition[j-1][i] >> 64 - - if M==1: - return [np.array(partition[0], dtype=np.uint64)] - else: - return [np.array(part, dtype=np.uint64) for part in partition[1:]] - - -def evaluate_observable(observable: int, x: int): - """ - This method evaluates an observable that is a tensor product of Pauli-:math:`Z` operators - with respect to a measurement outcome. - - A Pauli operator of the form :math:`\prod_{i\in I}Z_i`, for some finite set of indices :math:`I\subset \mathbb N`, - is identified with an integer: - We identify the Pauli operator with the binary string that has ones at positions :math:`i\in I` - and zeros otherwise, and then convert this binary string to an integer. - - Parameters - ---------- - - observable : int - The observable represented as integer. - x : int - The measurement outcome represented as integer. - - Returns - ------- - int - The value of the observable with respect to the measurement outcome. - - """ - - if bin(observable & x).count('1') % 2 == 0: - return 1 - else: - return -1 - - -@njit(cache = True) -def evaluate_observable_jitted(observable, x): - """ - - """ - - value = observable & x - count = 0 - while value: - count += value & 1 - value >>= 1 - return 1 if count % 2 == 0 else -1 - - -@njit(parallel = True, cache = True) -def evaluate_observables_parallel(observables_parts, outcome_parts, probabilities): #M=#observables, N=#measurements - """ - - Parameters - ---------- - observables_parts : list[numpy.array[numpy.unit64]] - The observables. - outcome_parts : list[numpy.array[numpy.unit64]] - The measurement outcomes. - probabilities : numpy.array[numpy.float64] - The measurment probabilities. - - Returns - ------- - - """ - - K = len(observables_parts) # Number of observables PARTS - M = len(observables_parts[0]) # Number of observables - N = len(probabilities) # Number of measurement results - - res = np.zeros(M, dtype=np.float64) - - for j in range(M): - - res_array = np.zeros(N, dtype=np.float64) - for i in prange(N): - temp = 1 - for k in range(K): - temp *= evaluate_observable_jitted(observables_parts[k][j], outcome_parts[k][i]) - - res_array[i] += temp - - res[j] = res_array @ probabilities - - return res - diff --git a/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py b/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py index 4f28cd40..5f0cc3c4 100644 --- a/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py +++ b/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py @@ -19,6 +19,7 @@ from qrisp.operators.pauli.helper_functions import * from qrisp.operators.pauli.bound_pauli_term import BoundPauliTerm from qrisp.operators.pauli.pauli_measurement import PauliMeasurement +from qrisp.operators.pauli.measurement import get_measurement import sympy as sp @@ -560,6 +561,104 @@ def commuting_qw_groups(self, show_bases=False): def pauli_measurement(self): return PauliMeasurement(self) + def get_measurement( + self, + qarg, + precision=0.01, + backend=None, + shots=1000000, + compile=True, + compilation_kwargs={}, + subs_dic={}, + precompiled_qc=None, + _measurement=None # measurement settings + ): + r""" + This method returns the expected value of a Hamiltonian for the state of a quantum argument. + + Parameters + ---------- + qarg : QuantumVariable, QuantumArray or list[QuantumVariable] + The quantum argument to evaluate the Hamiltonian on. + precision: float, optional + The precision with which the expectation of the Hamiltonian is to be evaluated. + The default is 0.01. The number of shots scales quadratically with the inverse precision. + backend : BackendClient, optional + The backend on which to evaluate the quantum circuit. The default can be + specified in the file default_backend.py. + shots : integer, optional + The maximum amount of shots to evaluate the expectation of the Hamiltonian. + The default is 1000000. + compile : bool, optional + Boolean indicating if the .compile method of the underlying QuantumSession + should be called before. The default is True. + compilation_kwargs : dict, optional + Keyword arguments for the compile method. For more details check + :meth:`QuantumSession.compile `. The default + is ``{}``. + subs_dic : dict, optional + A dictionary of Sympy symbols and floats to specify parameters in the case + of a circuit with unspecified, :ref:`abstract parameters`. + The default is {}. + circuit_preprocessor : Python function, optional + A function which recieves a QuantumCircuit and returns one, which is applied + after compilation and parameter substitution. The default is None. + precompiled_qc : QuantumCircuit, optional + A precompiled quantum circuit. + + Raises + ------ + Exception + If the containing QuantumSession is in a quantum environment, it is not + possible to execute measurements. + + Returns + ------- + float + The expected value of the Hamiltonian. + + Examples + -------- + + We define a Hamiltonian, and measure its expected value for the state of a :ref:`QuantumVariable`. + + :: + + from qrisp import QuantumVariable, h + from qrisp.operators.pauli import X,Y,Z + qv = QuantumVariable(2) + h(qv) + H = Z(0)*Z(1) + res = H.get_measurement(qv) + print(res) + #Yields 0.0 + + We define a Hamiltonian, and measure its expected value for the state of a :ref:`QuantumArray`. + + :: + + from qrisp import QuantumVariable, QuantumArray, h + from qrisp.operators.pauli import X,Y,Z + qtype = QuantumVariable(2) + q_array = QuantumArray(qtype, shape=(2)) + h(q_array) + H = Z(0)*Z(1) + X(2)*X(3) + res = H.get_measurement(q_array) + print(res) + #Yields 1.0 + + """ + return get_measurement(self, + qarg, + precision=precision, + backend=backend, + shots=shots, + compile=compile, + compilation_kwargs=compilation_kwargs, + subs_dic=subs_dic, + precompiled_qc=precompiled_qc, + _measurement=_measurement) + # # Trotterization # diff --git a/src/qrisp/operators/pauli/measurement.py b/src/qrisp/operators/pauli/measurement.py new file mode 100644 index 00000000..0d45ccbf --- /dev/null +++ b/src/qrisp/operators/pauli/measurement.py @@ -0,0 +1,346 @@ +""" +\******************************************************************************** +* Copyright (c) 2024 the Qrisp authors +* +* This program and the accompanying materials are made available under the +* terms of the Eclipse Public License 2.0 which is available at +* http://www.eclipse.org/legal/epl-2.0. +* +* This Source Code may also be made available under the following Secondary +* Licenses when the conditions for such availability set forth in the Eclipse +* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 +* with the GNU Classpath Exception which is +* available at https://www.gnu.org/software/classpath/license.html. +* +* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 +********************************************************************************/ +""" + +from qrisp import QuantumVariable, QuantumArray +from qrisp.core.compilation import qompiler +import math +import numpy as np +from numba import njit, prange + + +def get_measurement( + hamiltonian, + qarg, + precision=0.01, + backend=None, + shots=1000000, + compile=True, + compilation_kwargs={}, + subs_dic={}, + precompiled_qc=None, + _measurement=None # measurement settings + ): + r""" + This method returns the expected value of a Hamiltonian for the state of a quantum argument. + + Parameters + ---------- + qarg : QuantumVariable, QuantumArray or list[QuantumVariable] + The quantum argument to evaluate the Hamiltonian on. + precision: float, optional + The precision with which the expectation of the Hamiltonian is to be evaluated. + The default is 0.01. The number of shots scales quadratically with the inverse precision. + backend : BackendClient, optional + The backend on which to evaluate the quantum circuit. The default can be + specified in the file default_backend.py. + shots : integer, optional + The maximum amount of shots to evaluate the expectation of the Hamiltonian. + The default is 1000000. + compile : bool, optional + Boolean indicating if the .compile method of the underlying QuantumSession + should be called before. The default is True. + compilation_kwargs : dict, optional + Keyword arguments for the compile method. For more details check + :meth:`QuantumSession.compile `. The default + is ``{}``. + subs_dic : dict, optional + A dictionary of Sympy symbols and floats to specify parameters in the case + of a circuit with unspecified, :ref:`abstract parameters`. + The default is {}. + precompiled_qc : QuantumCircuit, optional + A precompiled quantum circuit. + + Raises + ------ + Exception + If the containing QuantumSession is in a quantum environment, it is not + possible to execute measurements. + + Returns + ------- + float + The expected value of the Hamiltonian. + + """ + + from qrisp import QuantumSession, merge + + if isinstance(qarg,QuantumVariable): + if qarg.is_deleted(): + raise Exception("Tried to get measurement from deleted QuantumVariable") + qs = qarg.qs + + elif isinstance(qarg,QuantumArray): + for qv in qarg.flatten(): + if qv.is_deleted(): + raise Exception( + "Tried to measure QuantumArray containing deleted QuantumVariables" + ) + qs = qarg.qs + elif isinstance(qarg,list): + qs = QuantumSession() + for qv in qarg: + if qv.is_deleted(): + raise Exception( + "Tried to measure QuantumArray containing deleted QuantumVariables" + ) + merge(qs,qv.qs) + + if backend is None: + if qs.backend is None: + from qrisp.default_backend import def_backend + + backend = def_backend + else: + backend = qarg.qs.backend + + if len(qs.env_stack) != 0: + raise Exception("Tried to get measurement within open environment") + + + # Copy circuit in over to prevent modification + if precompiled_qc is None: + if compile: + qc = qompiler( + qs, **compilation_kwargs + ) + else: + qc = qs.copy() + else: + qc = precompiled_qc.copy() + + # Bind parameters + if subs_dic: + qc = qc.bind_parameters(subs_dic) + from qrisp.core.compilation import combine_single_qubit_gates + qc = combine_single_qubit_gates(qc) + + qc = qc.transpile() + + from qrisp.misc import get_measurement_from_qc + + if _measurement is None: + pauli_measurement = hamiltonian.pauli_measurement() + else: + pauli_measurement = _measurement + + meas_circs = pauli_measurement.circuits + meas_qubits = pauli_measurement.qubits + meas_ops = pauli_measurement.operators_int + meas_coeffs = pauli_measurement.coefficients + meas_shots = pauli_measurement.shots + + meas_shots = [round(x/precision**2) for x in meas_shots] + tot_shots = sum(x for x in meas_shots) + if tot_shots>shots: + #meas_shots = [round(x*shots/tot_shots) for x in meas_shots] + raise Warning("Warning: The total amount of shots required: " + str(tot_shots) +" for the target precision: " + str(precision) + " exceeds the allowed maxium amount of shots. Decrease the precision or increase the maxium amount of shots.") + + results = [] + + for index,circuit in enumerate(meas_circs): + + curr = qc.copy() + curr.append(circuit.to_gate(), meas_qubits[index]) + + res = get_measurement_from_qc(curr, meas_qubits[index], backend, meas_shots[index]) + results.append(res) + + expectation = evaluate_expectation(results, meas_ops, meas_coeffs) + #expectation = evaluate_expectation_numba(results, meas_ops, meas_coeffs, meas_qubits) + return expectation + +# +# Evaluate expectation +# + +def evaluate_observable(observable: int, x: int): + """ + This method evaluates an observable that is a tensor product of Pauli-:math:`Z` operators + with respect to a measurement outcome. + + A Pauli operator of the form :math:`\prod_{i\in I}Z_i`, for some finite set of indices :math:`I\subset \mathbb N`, + is identified with an integer: + We identify the Pauli operator with the binary string that has ones at positions :math:`i\in I` + and zeros otherwise, and then convert this binary string to an integer. + + Parameters + ---------- + + observable : int + The observable represented as integer. + x : int + The measurement outcome represented as integer. + + Returns + ------- + int + The value of the observable with respect to the measurement outcome. + + """ + + if bin(observable & x).count('1') % 2 == 0: + return 1 + else: + return -1 + + +def evaluate_expectation(results, operators, coefficients): + """ + Evaluate the expectation. + + """ + + expectation = 0 + + for index1,ops in enumerate(operators): + for index2,op in enumerate(ops): + for outcome,probability in results[index1].items(): + expectation += probability*evaluate_observable(op,outcome)*coefficients[index1][index2] + + return expectation + + +# +# Numba accelearation +# + + +def partition(values, num_qubits): + """ + Partitions a list of integers into a list of lists of integers with size 64 bit. + + Parameters + ---------- + values : list[int] + A list of integers. + num_qubits : int + The maximal number of bits. + + Returns + ------- + partition : list[np.array] + A list of NumPy numpy.uint64 arrays. + + """ + + + M = math.ceil(num_qubits/64) + N = len(values) + + zeros = [[0]*N for k in range(M)] + partition = [values] + partition.extend(zeros) + + lower_mask = (1<<64) - 1 + + for j in range(1,M): + for i in range(N): + partition[j][i] = partition[j-1][i] & lower_mask + partition[j+1][i] = partition[j-1][i] >> 64 + + if M==1: + return [np.array(partition[0], dtype=np.uint64)] + else: + return [np.array(part, dtype=np.uint64) for part in partition[1:]] + + +@njit(cache = True) +def evaluate_observable_jitted(observable, x): + """ + + """ + + value = observable & x + count = 0 + while value: + count += value & 1 + value >>= 1 + return 1 if count % 2 == 0 else -1 + + +@njit(parallel = True, cache = True) +def evaluate_observables_parallel(observables_parts, outcome_parts, probabilities): #M=#observables, N=#measurements + """ + + Parameters + ---------- + observables_parts : list[numpy.array[numpy.unit64]] + The observables. + outcome_parts : list[numpy.array[numpy.unit64]] + The measurement outcomes. + probabilities : numpy.array[numpy.float64] + The measurment probabilities. + + Returns + ------- + + """ + + K = len(observables_parts) # Number of observables PARTS + M = len(observables_parts[0]) # Number of observables + N = len(probabilities) # Number of measurement results + + res = np.zeros(M, dtype=np.float64) + + for j in range(M): + + res_array = np.zeros(N, dtype=np.float64) + for i in prange(N): + temp = 1 + for k in range(K): + temp *= evaluate_observable_jitted(observables_parts[k][j], outcome_parts[k][i]) + + res_array[i] += temp + + res[j] = res_array @ probabilities + + return res + + +def evaluate_expectation_numba(results, operators, coefficients, qubits): + """ + Evaluate the expectation. + + """ + N = len(operators) + + # Partition outcomes in uint64 + #outcomes_parts = [partition(result.keys(), len(meas_qubits[k])) for k, result in enumerate(results)] + outcomes_parts = [partition(list(results[k].keys()), len(qubits[k])) if len(qubits[k])>64 else [np.array(list(results[k].keys()), dtype=np.uint64)] for k in range(N)] + + probabilities = [np.array(list(result.values()), dtype=np.float64) for result in results] + + # Partition observables in uint64 + observables_parts = [partition(observable, len(qubits[k])) if len(qubits[k])>64 else [np.array(observable, dtype=np.uint64)] for k, observable in enumerate(operators)] + + coefficients = [np.array(coeffs, dtype=np.float64) for coeffs in coefficients] + + # Evaluate the observable + expectation = 0 + + + for k in range(N): + + result = evaluate_observables_parallel(observables_parts[k], outcomes_parts[k], probabilities[k]) + #print(result) + #print(coefficients[k]) + + expectation += result @ coefficients[k] + + return expectation \ No newline at end of file diff --git a/src/qrisp/operators/pauli/pauli_hamiltonian.py b/src/qrisp/operators/pauli/pauli_hamiltonian.py index 837d6bae..1eeac821 100644 --- a/src/qrisp/operators/pauli/pauli_hamiltonian.py +++ b/src/qrisp/operators/pauli/pauli_hamiltonian.py @@ -19,6 +19,7 @@ from qrisp.operators.pauli.helper_functions import * from qrisp.operators.pauli.pauli_term import PauliTerm from qrisp.operators.pauli.pauli_measurement import PauliMeasurement +from qrisp.operators.pauli.measurement import get_measurement import sympy as sp @@ -553,11 +554,106 @@ def commuting_qw_groups(self, show_bases=False): return groups # - # Measurement settings + # Measurement settings and measurement # def pauli_measurement(self): return PauliMeasurement(self) + + def get_measurement( + self, + qarg, + precision=0.01, + backend=None, + shots=1000000, + compile=True, + compilation_kwargs={}, + subs_dic={}, + precompiled_qc=None, + _measurement=None # measurement settings + ): + r""" + This method returns the expected value of a Hamiltonian for the state of a quantum argument. + + Parameters + ---------- + qarg : QuantumVariable, QuantumArray or list[QuantumVariable] + The quantum argument to evaluate the Hamiltonian on. + precision: float, optional + The precision with which the expectation of the Hamiltonian is to be evaluated. + The default is 0.01. The number of shots scales quadratically with the inverse precision. + backend : BackendClient, optional + The backend on which to evaluate the quantum circuit. The default can be + specified in the file default_backend.py. + shots : integer, optional + The maximum amount of shots to evaluate the expectation of the Hamiltonian. + The default is 1000000. + compile : bool, optional + Boolean indicating if the .compile method of the underlying QuantumSession + should be called before. The default is True. + compilation_kwargs : dict, optional + Keyword arguments for the compile method. For more details check + :meth:`QuantumSession.compile `. The default + is ``{}``. + subs_dic : dict, optional + A dictionary of Sympy symbols and floats to specify parameters in the case + of a circuit with unspecified, :ref:`abstract parameters`. + The default is {}. + precompiled_qc : QuantumCircuit, optional + A precompiled quantum circuit. + + Raises + ------ + Exception + If the containing QuantumSession is in a quantum environment, it is not + possible to execute measurements. + + Returns + ------- + float + The expected value of the Hamiltonian. + + Examples + -------- + + We define a Hamiltonian, and measure its expected value for the state of a :ref:`QuantumVariable`. + + :: + + from qrisp import QuantumVariable, h + from qrisp.operators.pauli import X,Y,Z + qv = QuantumVariable(2) + h(qv) + H = Z(0)*Z(1) + res = H.get_measurement(qv) + print(res) + #Yields 0.0 + + We define a Hamiltonian, and measure its expected value for the state of a :ref:`QuantumArray`. + + :: + + from qrisp import QuantumVariable, QuantumArray, h + from qrisp.operators.pauli import X,Y,Z + qtype = QuantumVariable(2) + q_array = QuantumArray(qtype, shape=(2)) + h(q_array) + H = Z(0)*Z(1) + X(2)*X(3) + res = H.get_measurement(q_array) + print(res) + #Yields 1.0 + + """ + return get_measurement(self, + qarg, + precision=precision, + backend=backend, + shots=shots, + compile=compile, + compilation_kwargs=compilation_kwargs, + subs_dic=subs_dic, + precompiled_qc=precompiled_qc, + _measurement=_measurement) # # Trotterization From cb6661626a7d24af23751f022e62a608d2756f31 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Wed, 30 Oct 2024 15:46:22 +0100 Subject: [PATCH 07/13] Code cleanup vizualization of PauliHamiltonians --- .../operators/fermionic/visualization.py | 6 +- src/qrisp/operators/pauli/bound_pauli_term.py | 2 +- src/qrisp/operators/pauli/pauli_term.py | 2 +- src/qrisp/operators/pauli/spin.py | 141 ------------------ src/qrisp/operators/pauli/visualization.py | 62 ++++++++ 5 files changed, 66 insertions(+), 147 deletions(-) delete mode 100644 src/qrisp/operators/pauli/spin.py create mode 100644 src/qrisp/operators/pauli/visualization.py diff --git a/src/qrisp/operators/fermionic/visualization.py b/src/qrisp/operators/fermionic/visualization.py index 045fcc1f..0250ec45 100644 --- a/src/qrisp/operators/fermionic/visualization.py +++ b/src/qrisp/operators/fermionic/visualization.py @@ -20,12 +20,10 @@ # ONLY USED FOR LATEX PRINTING # -import sympy as sp -from sympy import Symbol, I -import numpy as np +from sympy import Symbol # -# Pauli symbols (only used for visualization, i.e., LateX printing with SymPy) +# Fermionic symbols (only used for visualization, i.e., LateX printing with SymPy) # class a_(Symbol): diff --git a/src/qrisp/operators/pauli/bound_pauli_term.py b/src/qrisp/operators/pauli/bound_pauli_term.py index 10373d31..e5a56157 100644 --- a/src/qrisp/operators/pauli/bound_pauli_term.py +++ b/src/qrisp/operators/pauli/bound_pauli_term.py @@ -16,7 +16,7 @@ ********************************************************************************/ """ -from qrisp.operators.pauli.spin import X_,Y_,Z_ +from qrisp.operators.pauli.visualization import X_,Y_,Z_ PAULI_TABLE = {("I","I"):("I",1),("I","X"):("X",1),("I","Y"):("Y",1),("I","Z"):("Z",1), ("X","I"):("X",1),("X","X"):("I",1),("X","Y"):("Z",1j),("X","Z"):("Y",-1j), diff --git a/src/qrisp/operators/pauli/pauli_term.py b/src/qrisp/operators/pauli/pauli_term.py index f3108701..e59cb100 100644 --- a/src/qrisp/operators/pauli/pauli_term.py +++ b/src/qrisp/operators/pauli/pauli_term.py @@ -16,7 +16,7 @@ ********************************************************************************/ """ -from qrisp.operators.pauli.spin import X_,Y_,Z_ +from qrisp.operators.pauli.visualization import X_,Y_,Z_ PAULI_TABLE = {("I","I"):("I",1),("I","X"):("X",1),("I","Y"):("Y",1),("I","Z"):("Z",1), ("X","I"):("X",1),("X","X"):("I",1),("X","Y"):("Z",1j),("X","Z"):("Y",-1j), diff --git a/src/qrisp/operators/pauli/spin.py b/src/qrisp/operators/pauli/spin.py deleted file mode 100644 index 95be6cc7..00000000 --- a/src/qrisp/operators/pauli/spin.py +++ /dev/null @@ -1,141 +0,0 @@ -""" -\******************************************************************************** -* Copyright (c) 2024 the Qrisp authors -* -* This program and the accompanying materials are made available under the -* terms of the Eclipse Public License 2.0 which is available at -* http://www.eclipse.org/legal/epl-2.0. -* -* This Source Code may also be made available under the following Secondary -* Licenses when the conditions for such availability set forth in the Eclipse -* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 -* with the GNU Classpath Exception which is -* available at https://www.gnu.org/software/classpath/license.html. -* -* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 -********************************************************************************/ -""" - -# -# ONLY USED FOR LATEX PRINTING -# - -import sympy as sp -from sympy import Symbol, I -import numpy as np - -# -# Helper functions -# - -def delta(i, j): - if i==j: - return 1 - else: - return 0 - -def epsilon(i, j, k): - if (i, j, k) in (("X", "Y", "Z"), ("Y", "Z", "X"), ("Z", "X", "Y")): - return 1 - if (i, j, k) in (("Z", "Y", "X"), ("Y", "X", "Z"), ("X", "Z", "Y")): - return -1 - return 0 - -def mul_helper(P1,P2): - pauli_table = {("X","X"):("I",1),("X","Y"):("Z",1j),("X","Z"):("Y",-1j), - ("Y","X"):("Z",-1j),("Y","Y"):("I",1),("Y","Z"):("X",1j), - ("Z","X"):("Y",1j),("Z","Y"):("X",-1j),("Z","Z"):("I",1)} - - if P1=="I": - return (P2,1) - if P2=="I": - return (P1,1) - return pauli_table[(P1,P2)] -# -# Pauli symbols (only used for visualization, i.e., LateX printing with SymPy) -# - -class X_(Symbol): - - __slots__ = ("axis","index") - - def __new__(cls, index): - obj = Symbol.__new__(cls, "%s%s" %("X",index), commutative=False, hermitian=True) - obj.index = index - return obj - - def _eval_power(b, e): - if e.is_Integer and e.is_positive: - return super().__pow__(int(e) % 2) - - def __mul__(self, other): - if isinstance(other, (X_,Y_,Z_)): - if self.index == other.index: - i = self.axis - j = other.axis - return delta(i, j) \ - + I*epsilon(i, j, "X")*X_(self.index) \ - + I*epsilon(i, j, "Y")*Y_(self.index) \ - + I*epsilon(i, j, "Z")*Z_(self.index) - return super().__mul__(other) - - __rmul__ = __mul__ - -class Y_(Symbol): - - __slots__ = ("axis","index") - - def __new__(cls, index): - obj = Symbol.__new__(cls, "%s%s" %("Y",index), commutative=False, hermitian=True) - obj.index = index - return obj - - def _eval_power(b, e): - if e.is_Integer and e.is_positive: - return super().__pow__(int(e) % 2) - - def __mul__(self, other): - if isinstance(other, (X_,Y_,Z_)): - if self.index == other.index: - i = self.axis - j = other.axis - return delta(i, j) \ - + I*epsilon(i, j, "X")*X_(self.index) \ - + I*epsilon(i, j, "Y")*Y_(self.index) \ - + I*epsilon(i, j, "Z")*Z_(self.index) - return super().__mul__(other) - - __rmul__ = __mul__ - -class Z_(Symbol): - - __slots__ = ("axis","index") - - def __new__(cls, index): - obj = Symbol.__new__(cls, "%s%s" %("Z",index), commutative=False, hermitian=True) - obj.index = index - return obj - - def _eval_power(b, e): - if e.is_Integer and e.is_positive: - return super().__pow__(int(e) % 2) - - def __mul__(self, other): - if isinstance(other, (X_,Y_,Z_)): - if self.index == other.index: - i = self.axis - j = other.axis - return delta(i, j) \ - + I*epsilon(i, j, "X")*X_(self.index) \ - + I*epsilon(i, j, "Y")*Y_(self.index) \ - + I*epsilon(i, j, "Z")*Z_(self.index) - return super().__mul__(other) - - __rmul__ = __mul__ - - - - - - - diff --git a/src/qrisp/operators/pauli/visualization.py b/src/qrisp/operators/pauli/visualization.py new file mode 100644 index 00000000..33a7b711 --- /dev/null +++ b/src/qrisp/operators/pauli/visualization.py @@ -0,0 +1,62 @@ +""" +\******************************************************************************** +* Copyright (c) 2024 the Qrisp authors +* +* This program and the accompanying materials are made available under the +* terms of the Eclipse Public License 2.0 which is available at +* http://www.eclipse.org/legal/epl-2.0. +* +* This Source Code may also be made available under the following Secondary +* Licenses when the conditions for such availability set forth in the Eclipse +* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 +* with the GNU Classpath Exception which is +* available at https://www.gnu.org/software/classpath/license.html. +* +* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 +********************************************************************************/ +""" + +# +# ONLY USED FOR LATEX PRINTING +# + +from sympy import Symbol + +# +# Pauli symbols (only used for visualization, i.e., LateX printing with SymPy) +# + +class X_(Symbol): + + __slots__ = ("axis","index") + + def __new__(cls, index): + obj = Symbol.__new__(cls, "%s%s" %("X",index), commutative=False, hermitian=True) + obj.index = index + return obj + +class Y_(Symbol): + + __slots__ = ("axis","index") + + def __new__(cls, index): + obj = Symbol.__new__(cls, "%s%s" %("Y",index), commutative=False, hermitian=True) + obj.index = index + return obj + +class Z_(Symbol): + + __slots__ = ("axis","index") + + def __new__(cls, index): + obj = Symbol.__new__(cls, "%s%s" %("Z",index), commutative=False, hermitian=True) + obj.index = index + return obj + + + + + + + + From 8b4474cb4372382efc0fef35a45dabd48ea252bc Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Wed, 30 Oct 2024 18:51:43 +0100 Subject: [PATCH 08/13] Restructuring trotterization for PauliHamiltonians --- .../pauli/bound_pauli_hamiltonian.py | 49 ++++++----------- src/qrisp/operators/pauli/bound_pauli_term.py | 24 ++++++++ src/qrisp/operators/pauli/helper_functions.py | 55 ------------------- .../operators/pauli/pauli_hamiltonian.py | 33 ++++------- .../operators/pauli/pauli_measurement.py | 9 ++- src/qrisp/operators/pauli/pauli_term.py | 23 ++++++++ 6 files changed, 83 insertions(+), 110 deletions(-) delete mode 100644 src/qrisp/operators/pauli/helper_functions.py diff --git a/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py b/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py index 5f0cc3c4..45bf6fdf 100644 --- a/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py +++ b/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py @@ -16,10 +16,10 @@ ********************************************************************************/ """ from qrisp.operators.hamiltonian import Hamiltonian -from qrisp.operators.pauli.helper_functions import * from qrisp.operators.pauli.bound_pauli_term import BoundPauliTerm from qrisp.operators.pauli.pauli_measurement import PauliMeasurement from qrisp.operators.pauli.measurement import get_measurement +from qrisp import h, sx, IterationEnvironment, conjugate import sympy as sp @@ -680,7 +680,7 @@ def trotterization(self): This function receives the following arguments: - * qarg : QuantumVariable or list[QuantumVariable] + * qarg : QuantumVariable The quantum argument. * t : float, optional The evolution time $t$. The default is 1. @@ -691,48 +691,31 @@ def trotterization(self): """ - from qrisp import conjugate, rx, ry, rz, cx, h, IterationEnvironment, gphase, QuantumSession, merge + def change_of_basis(pauli_dict): + for qubit, axis in pauli_dict.items(): + if axis=="X": + h(qubit) + if axis=="Y": + sx(qubit) - pauli_measurement = self.pauli_measurement() - bases = pauli_measurement.bases - indices = pauli_measurement.operators_ind # Indices (Qubits) of Z's in PauliTerms (after change of basis) - operators_int = pauli_measurement.operators_int - coeffs = pauli_measurement.coefficients + groups, bases = self.commuting_qw_groups(show_bases=True) def trotter_step(qarg, t, steps): + # qubit to apply gphase for identity term if isinstance(qarg,list): qubit = qarg[0][0] else: qubit = qarg[0] - N = len(bases) - for k in range(N): - basis = bases[k].pauli_dict - - with conjugate(change_of_basis_bound)(basis): - M = len(indices[k]) - - for l in range(M): - if(operators_int[k][l]>0): # Not identity - with conjugate(parity_bound)(indices[k][l]): - rz(-2*coeffs[k][l]*t/steps,indices[k][l][-1]) - else: # Identity - gphase(coeffs[k][l]*t/steps,qubit) + for index,basis in enumerate(bases): + with conjugate(change_of_basis)(basis.pauli_dict): + for term,coeff in groups[index].terms_dict.items(): + term.simulate(coeff*t/steps, qubit) def U(qarg, t=1, steps=1, iter=1): - - if isinstance(qarg,list): - qs = QuantumSession() - for qv in qarg: - merge(qs,qv.qs) - else: - qs = qarg.qs - - with IterationEnvironment(qs, iter): + with IterationEnvironment(qarg.qs, iter): for i in range(steps): trotter_step(qarg, t, steps) - return U - - + return U \ No newline at end of file diff --git a/src/qrisp/operators/pauli/bound_pauli_term.py b/src/qrisp/operators/pauli/bound_pauli_term.py index e5a56157..b3560e7a 100644 --- a/src/qrisp/operators/pauli/bound_pauli_term.py +++ b/src/qrisp/operators/pauli/bound_pauli_term.py @@ -16,6 +16,7 @@ ********************************************************************************/ """ +from qrisp import gphase, rz, cx, conjugate from qrisp.operators.pauli.visualization import X_,Y_,Z_ PAULI_TABLE = {("I","I"):("I",1),("I","X"):("X",1),("I","Y"):("Y",1),("I","Z"):("Z",1), @@ -48,6 +49,29 @@ def __eq__(self, other): def copy(self): return BoundPauliTerm(self.pauli_dict.copy()) + + def is_identity(self): + return len(self.pauli_dict)==0 + + # + # Simulation + # + + # Assume that the operator is diagonal after change of basis + def simulate(self, coeff, qubit): + + def parity(qubits): + n = len(qubits) + for i in range(n-1): + cx(qubits[i],qubits[i+1]) + + if not self.is_identity(): + qubits = list(self.pauli_dict.keys()) + with conjugate(parity)(qubits): + rz(-2*coeff,qubits[-1]) + else: + gphase(coeff,qubit) + # # Printing # diff --git a/src/qrisp/operators/pauli/helper_functions.py b/src/qrisp/operators/pauli/helper_functions.py deleted file mode 100644 index defc86f3..00000000 --- a/src/qrisp/operators/pauli/helper_functions.py +++ /dev/null @@ -1,55 +0,0 @@ -""" -\******************************************************************************** -* Copyright (c) 2024 the Qrisp authors -* -* This program and the accompanying materials are made available under the -* terms of the Eclipse Public License 2.0 which is available at -* http://www.eclipse.org/legal/epl-2.0. -* -* This Source Code may also be made available under the following Secondary -* Licenses when the conditions for such availability set forth in the Eclipse -* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 -* with the GNU Classpath Exception which is -* available at https://www.gnu.org/software/classpath/license.html. -* -* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 -********************************************************************************/ -""" - -def get_integer_from_indices(indices,positions=None): - if positions is not None: - return sum(1 << positions[i] for i in indices) - else: - return sum(1 << i for i in indices) - -# -# Trotterization -# - -from qrisp import conjugate, rx, ry, rz, cx, h, sx, IterationEnvironment, gphase -import numpy as np - -def change_of_basis(qarg, pauli_dict): - for index, axis in pauli_dict.items(): - if axis=="X": - h(qarg[index]) - if axis=="Y": - sx(qarg[index]) - -def parity(qarg, indices): - n = len(indices) - for i in range(n-1): - cx(qarg[indices[i]],qarg[indices[i+1]]) - -def change_of_basis_bound(pauli_dict): - - for qubit, axis in pauli_dict.items(): - if axis=="X": - h(qubit) - if axis=="Y": - sx(qubit) - -def parity_bound(qubits): - n = len(qubits) - for i in range(n-1): - cx(qubits[i],qubits[i+1]) \ No newline at end of file diff --git a/src/qrisp/operators/pauli/pauli_hamiltonian.py b/src/qrisp/operators/pauli/pauli_hamiltonian.py index 1eeac821..7e3a46e5 100644 --- a/src/qrisp/operators/pauli/pauli_hamiltonian.py +++ b/src/qrisp/operators/pauli/pauli_hamiltonian.py @@ -16,10 +16,10 @@ ********************************************************************************/ """ from qrisp.operators.hamiltonian import Hamiltonian -from qrisp.operators.pauli.helper_functions import * from qrisp.operators.pauli.pauli_term import PauliTerm from qrisp.operators.pauli.pauli_measurement import PauliMeasurement from qrisp.operators.pauli.measurement import get_measurement +from qrisp import h, sx, IterationEnvironment, conjugate import sympy as sp @@ -687,29 +687,20 @@ def trotterization(self): """ - from qrisp import conjugate, rx, ry, rz, cx, h, IterationEnvironment, gphase + def change_of_basis(qarg, pauli_dict): + for index, axis in pauli_dict.items(): + if axis=="X": + h(qarg[index]) + if axis=="Y": + sx(qarg[index]) - pauli_measurement = self.pauli_measurement() - bases = pauli_measurement.bases - indices = pauli_measurement.operators_ind # Indices of Z's in PauliTerms (after change of basis) - operators_int = pauli_measurement.operators_int - coefficients = pauli_measurement.coefficients + groups, bases = self.commuting_qw_groups(show_bases=True) def trotter_step(qarg, t, steps): - - N = len(bases) - for k in range(N): - basis = bases[k].pauli_dict - - with conjugate(change_of_basis)(qarg, basis): - M = len(indices[k]) - - for l in range(M): - if operators_int[k][l]>0: # Not identity - with conjugate(parity)(qarg, indices[k][l]): - rz(-2*coefficients[k][l]*t/steps,qarg[indices[k][l][-1]]) - else: # Identity - gphase(coefficients[k][l]*t/steps,qarg[0]) + for index,basis in enumerate(bases): + with conjugate(change_of_basis)(qarg, basis.pauli_dict): + for term,coeff in groups[index].terms_dict.items(): + term.simulate(coeff*t/steps, qarg) def U(qarg, t=1, steps=1, iter=1): with IterationEnvironment(qarg.qs, iter): diff --git a/src/qrisp/operators/pauli/pauli_measurement.py b/src/qrisp/operators/pauli/pauli_measurement.py index b3f29611..3b118d35 100644 --- a/src/qrisp/operators/pauli/pauli_measurement.py +++ b/src/qrisp/operators/pauli/pauli_measurement.py @@ -16,10 +16,17 @@ ********************************************************************************/ """ -from qrisp.operators.pauli.helper_functions import * from qrisp import QuantumVariable, QuantumArray, QuantumCircuit import numpy as np + +def get_integer_from_indices(indices,positions=None): + if positions is not None: + return sum(1 << positions[i] for i in indices) + else: + return sum(1 << i for i in indices) + + class PauliMeasurement: """ diff --git a/src/qrisp/operators/pauli/pauli_term.py b/src/qrisp/operators/pauli/pauli_term.py index e59cb100..8ea98dc0 100644 --- a/src/qrisp/operators/pauli/pauli_term.py +++ b/src/qrisp/operators/pauli/pauli_term.py @@ -16,6 +16,7 @@ ********************************************************************************/ """ +from qrisp import gphase, rz, cx, conjugate from qrisp.operators.pauli.visualization import X_,Y_,Z_ PAULI_TABLE = {("I","I"):("I",1),("I","X"):("X",1),("I","Y"):("Y",1),("I","Z"):("Z",1), @@ -49,6 +50,28 @@ def __eq__(self, other): def copy(self): return PauliTerm(self.pauli_dict.copy()) + def is_identity(self): + return len(self.pauli_dict)==0 + + # + # Simulation + # + + # Assume that the operator is diagonal after change of basis + def simulate(self, coeff, qarg): + + def parity(qarg, indices): + n = len(indices) + for i in range(n-1): + cx(qarg[indices[i]],qarg[indices[i+1]]) + + if not self.is_identity(): + indices = list(self.pauli_dict.keys()) + with conjugate(parity)(qarg, indices): + rz(-2*coeff,qarg[indices[-1]]) + else: + gphase(coeff,qarg[0]) + # # Printing # From 2834454a710f40bb696e7e8a6d3393381fbd62b3 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Wed, 30 Oct 2024 19:20:17 +0100 Subject: [PATCH 09/13] Added lifted decorator to simulate method for PauliTerm --- src/qrisp/operators/pauli/bound_pauli_term.py | 4 +++- src/qrisp/operators/pauli/pauli_term.py | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/qrisp/operators/pauli/bound_pauli_term.py b/src/qrisp/operators/pauli/bound_pauli_term.py index b3560e7a..3d2b14cb 100644 --- a/src/qrisp/operators/pauli/bound_pauli_term.py +++ b/src/qrisp/operators/pauli/bound_pauli_term.py @@ -16,7 +16,7 @@ ********************************************************************************/ """ -from qrisp import gphase, rz, cx, conjugate +from qrisp import gphase, rz, cx, conjugate, lifted from qrisp.operators.pauli.visualization import X_,Y_,Z_ PAULI_TABLE = {("I","I"):("I",1),("I","X"):("X",1),("I","Y"):("Y",1),("I","Z"):("Z",1), @@ -58,6 +58,8 @@ def is_identity(self): # # Assume that the operator is diagonal after change of basis + # Implements exp(i*coeff*\prod_j Z_j) where the product goes over all qubits j in self.pauli_dict + @lifted def simulate(self, coeff, qubit): def parity(qubits): diff --git a/src/qrisp/operators/pauli/pauli_term.py b/src/qrisp/operators/pauli/pauli_term.py index 8ea98dc0..515ef964 100644 --- a/src/qrisp/operators/pauli/pauli_term.py +++ b/src/qrisp/operators/pauli/pauli_term.py @@ -16,7 +16,7 @@ ********************************************************************************/ """ -from qrisp import gphase, rz, cx, conjugate +from qrisp import gphase, rz, cx, conjugate, lifted from qrisp.operators.pauli.visualization import X_,Y_,Z_ PAULI_TABLE = {("I","I"):("I",1),("I","X"):("X",1),("I","Y"):("Y",1),("I","Z"):("Z",1), @@ -58,6 +58,8 @@ def is_identity(self): # # Assume that the operator is diagonal after change of basis + # Implements exp(i*coeff*\prod_j Z_j) where the product goes over all indices j in self.pauli_dict + @lifted def simulate(self, coeff, qarg): def parity(qarg, indices): From 1e0a7ddf2ea57ed39da625effaab2e34fed7135f Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Wed, 30 Oct 2024 20:51:29 +0100 Subject: [PATCH 10/13] Minor update trotterization for PauliHamiltonian --- src/qrisp/operators/pauli/bound_pauli_hamiltonian.py | 8 ++++---- src/qrisp/operators/pauli/pauli_hamiltonian.py | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py b/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py index 45bf6fdf..44a61497 100644 --- a/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py +++ b/src/qrisp/operators/pauli/bound_pauli_hamiltonian.py @@ -19,7 +19,7 @@ from qrisp.operators.pauli.bound_pauli_term import BoundPauliTerm from qrisp.operators.pauli.pauli_measurement import PauliMeasurement from qrisp.operators.pauli.measurement import get_measurement -from qrisp import h, sx, IterationEnvironment, conjugate +from qrisp import h, sx, IterationEnvironment, conjugate, merge import sympy as sp @@ -714,8 +714,8 @@ def trotter_step(qarg, t, steps): term.simulate(coeff*t/steps, qubit) def U(qarg, t=1, steps=1, iter=1): - with IterationEnvironment(qarg.qs, iter): - for i in range(steps): - trotter_step(qarg, t, steps) + merge([qarg]) + with IterationEnvironment(qarg.qs, iter*steps): + trotter_step(qarg, t, steps) return U \ No newline at end of file diff --git a/src/qrisp/operators/pauli/pauli_hamiltonian.py b/src/qrisp/operators/pauli/pauli_hamiltonian.py index 7e3a46e5..da72eb71 100644 --- a/src/qrisp/operators/pauli/pauli_hamiltonian.py +++ b/src/qrisp/operators/pauli/pauli_hamiltonian.py @@ -19,7 +19,7 @@ from qrisp.operators.pauli.pauli_term import PauliTerm from qrisp.operators.pauli.pauli_measurement import PauliMeasurement from qrisp.operators.pauli.measurement import get_measurement -from qrisp import h, sx, IterationEnvironment, conjugate +from qrisp import h, sx, IterationEnvironment, conjugate, merge import sympy as sp @@ -703,8 +703,8 @@ def trotter_step(qarg, t, steps): term.simulate(coeff*t/steps, qarg) def U(qarg, t=1, steps=1, iter=1): - with IterationEnvironment(qarg.qs, iter): - for i in range(steps): - trotter_step(qarg, t, steps) + merge([qarg]) + with IterationEnvironment(qarg.qs, iter*steps): + trotter_step(qarg, t, steps) return U From d412398b45d1f32234bd496cc04a89bb412c3bda Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Wed, 30 Oct 2024 23:34:43 +0100 Subject: [PATCH 11/13] Ising model example --- .../source/_static/Ising_chain_N=6.png | Bin 0 -> 22554 bytes .../source/reference/Examples/IsingModel.rst | 84 ++++++++++++++++++ .../source/reference/Examples/index.rst | 5 +- 3 files changed, 88 insertions(+), 1 deletion(-) create mode 100644 documentation/source/_static/Ising_chain_N=6.png create mode 100644 documentation/source/reference/Examples/IsingModel.rst diff --git a/documentation/source/_static/Ising_chain_N=6.png b/documentation/source/_static/Ising_chain_N=6.png new file mode 100644 index 0000000000000000000000000000000000000000..17b0e56e018f8bc87fddc33398b5ab79e8f59526 GIT binary patch literal 22554 zcmbWf2RPSl-#`AXlm=xqsR)rsQ6w@#Rz}K5N(l*Bp=3uyMj=T^M%hXB%qXNHyX=vy ztcuM1U#IK3pZmF=-+lkazvH-$boqWi^E}_@YrRkJlPYqwYnj$kC=^-+d6`ob%5nt? zWm(?pRro*K&FgmJuOl|c&)A%{xNLJn@0uY+Sk!N_@KpIQ}d|jJ%WPDU#K=7H&=YQ((Zq0#^U)FJB(M zLCvdeX0~tb+O=JyqwBq8E{%@jC~9J>-CXPHPEoB|Ra*U`yV?B_M`P5}r<9VCl5eA{ za3DcRUEQ6J4_kC}O(_Kqk*{`@gnjSs)_b|2siC3wd$irl$HyIiN~TPWx5n9~c?g^B zl6AE@%UqtIz!$Wk0$?Ner$DigN|4I=1Wvt z<1+QaTLwkdJy&FVdwF@u98e%P%TU5W{v}Tp#6kW^zdc7;#vsPT#N_1Z)B2<4Ta`65 zHhmu%QPb9bC!8W_Yn%D5pVt7-!9bttm}`E%=&M(+KD2L&VW3*}@v)$^uI>{J)zw!` ztvYh7RvbHajMF<94=cbg;CofeNFMtbtwabn5X->}8Cr_S?PfcYFT}aH#%*=6F6c3Mx*cP`b zJ|)HD&K>IS-@l(!SEo`{RgK}AZkjX4DOdFgb8OkN#rlRv{IRz7c4eU_)oekruU-l8 zp39+ILA_q4{<-q50Qr~B3$ymwuN{8#Jy(fXd*Q-`oV>hfi_Rip+o1;&B1$q1^6FYz zD=4w25<(rM8V*cNZ7_G^ejVcX4hPf8HrIO>xI4$PdpR3ByRC_ij*d@>JGK7%hx;Ch zJBf(T4zO?Bv?)q8Qas0T=2d{b5(dxni&gwKn?pHPQc*p9`I7F)kt6boiggKQCMJ9R zcZ%Je9B$U{F7>b*$!xalDqbf4>=Z8dle}1ya&@jiW-1O~aFbU+;o7xpQt6jJGBpSZ z3U2!I=g(IDOScn0G&gUFk-JyEL1LcbA*B~@E+pUb<7#ra<9f==Yd?;}?yC_pbs?`h zyiu(E`s^=bo58x{`nO$dYz}Fr>iV4s;oHT=CY6yQ6?O7NG~R+L@)`nrvG z8Qo`HjrSvHeb2W-vxk?L_f4!q&;ipX-lGe@>N!P6ajz3SKKt4D`RT2!tYjrR%KA^! zQK@Ta81K+Aim!RuxLd))!(*T^DMa?Bq$Gv%^yyPc)4*TdcUUe>^t@&2%fB%u^f`!i z#pZ?1#kuO&D=5=qzqiaknVZf!o*S6Z28NM^wETJDax06c!c^DLNS&<2~xz z&HWy$2)W?AhH~bW!J3wAi`30=tg9ahnK9k<@u_`wV#|j}iOzH98^1(J-(d}T&u=$= zG4qGR@6lzHy*xY`E^T2wIfjNCJbl&G)wlNLw6x{l*sy8Srs=+r3mNM>`1bEVT$z^4 zr;r#Q?-~>oB%7O_o^Fya^s~RFB3V0KJ+<%J5jS^t`Yl_iH*DBYIQP){1q}~Zq6yvq z3jskv?rHiJiNX1~xq3gpeZ<5(B{pMtw>K|8zqf+L<+l1>!FTn6n0FZ&8(+M5asQnG zc?m-|=3;1qGqleKas<7Tso0^RshOB%+WYk@$K#7fgoU?$rM`TnpL5{wuB~%Ieak2& z8yS8@ZuR|Ctt>OVrC%dTHtNXLWZ~4cQ4!2F8!V*@T^{$B(wURzu%wcn(!|zq&U5ed!-3QFugy;Ju7S)m#+i!Gsk9U_z!eTEy`T*N*Z%jzVQ`@S#z zzWj|hoK(BSj$F{#e0cZjT>+zKTMY~!M(>GR?!3>?wuznX)n!zoE1gAR>siE;6$4em zg>PZbkn+c*aO&(?ro0$ev>j1}H_b|Qq0lT%h>KR5Rp>J2;JxGr1q zsmRHJyc3JU$%jqOyyx>al63`_%y*&4ck$2f*BGv)@Atn>PE5b~^YYctFPGu%l^?7* z>W+HQbBD!Q^3o-?qS+yy9E|R^9IKogV>`{;^E|XJIBjyB5@p|aHlw)U=FbMKK23}BZl{qkj-QFVBUCoS_yU0wf#({amPRVr4xo{m$|D{))7diCn! zH{OgcB8rdG8S>#ue!Zc=md-Nb-h(lrZaEv7237Bj8P}_X^4!> z!{2geF+S6-kDT4Jac!l6!1wn2EqHZUK<~-6+YUD+EI4l!vDx!nIn1-6L0u!zX6A?C z_M>)m?6SVa8OHS#inJ%qu~VlQUg#F?ZMpO5(@73a&P3h!eb?saf7T>&bd~!tf6lRD z=Tv@t+p4!bVQBLE7e*|IWm6-q{FnhYBQ2YZjg9y4@RZ?|4Eg3O1aX?!QXMPKA1hi% zO-=nF^q`!ExXT~M%tw(9qemQP2jiwEDp>f9>o|l(M22k6QO*?F3u4(#yVj@JXqOuI z#C|O~JWk7hWAt1^2a9$3=diR$mp?}y@}0k%Y1+J;;ygbgM_xD9q2qXluGCyPC}Nd0 zi?82Y{8||>&`~I=Yd^S~SNj#!1?QRFv%d!2A3tWQU*E?fG}SyjIW8^hYNCVB}Jv*xgq;ib1A#c5eK8dQD=Xzn7QmIX*41fa}xadH{AE zg$^;Xu`-jX#`VuvuU#v3`SNoWexENH#*EGjlUn3CjXZ2F;|~c5!N;h8_TrP1Wq^{Z zE4cai%7b~d4(+^d(w--nXEW%Ax^0?qvu#e_e6z#u!QOikvMMV6Lz9?_FU?woD8;CX zU4w(%Eh>A1FZEYPJUVE!{?AN}OLTPf^7cI2vf0Vui8qX@k7CEZoqWLlr$JXu#pC$= zaXFc$6Rj-9YAXsp&i<@;|9%}V;`3`4iP6s2_gY%c)^cBRrB<~x>nPZ!_2Qgbn0|kE z+~VE4cP|L9Hpq zKpaPP>`HZ6NNUOY{{1SceQ~OhO?C!FyA-PYm=9A3;ALt(n{nk<5SJ>hS!(2Zd+*Xf z;NTV&olCQAO1-c@Z`=&W9shW^vn$M^h~Ke-cO&&p#(>;2KUj>#9H!nRY9v|qSMRVH zZc?5+F|D^a*SR=8nm-2Mx;*vuACbpGW@{ZB9KN??uA7@0sUMzhc_MUL#rKSk&N^=y zx8I#E7k(-5JFPmfW%K49%wU-izVqXCC-x_p@$m3Sqi&o#cW#_Ok^=|m{Qdob+60Pb zzAhJHrgXg?)^XBNAKffjZaC388^?N_g;KqtSl4g;YY(nc5qfTA~ zsIj!RE^}MC?oeRBhYw1Xeb)$9#k4;aeBrfoVM)d0@88svS{>V)H*b#r`ju?`DSP*; z&-@Cp`6wj^`Y%L`_ehkIzYkp1BqUgTHt%!G&;NGwB&YHQmh{BLwXN9}&G9;&z4i_cG0I^=@~Wy@AIoVe z+S=OFBiWsFJ9eB;Ez)7oRd^+UwL@g#;;W1}HLrN9x;}+-r+xhVQmg~tK62{{0&5q^qi7&eT0e9{^Hl3m8(|uJPXm4 zx^_+A#f8`Yq}bqon^P}TprjuFc=KfvTuY(byx9}DziT-4wJb{VuAq10UxQRlz8_>~ zUokj1xOD6!t<=P~toTjCY1JtG3V+KRH{?)7u^giCv8X0`0%P(uQ_N?n%) z+oL}}98Bbtp5?6EY|vg24%vA6#8v}FI-}^?^Z{MvIFCSExg4de z+&Nagp5-wM*8SDTG2&mkEIMyLa*gVIo=r)s#cRJ1UHgQR^rNpasZdo)Q3CD<1}>+J zPfi|hesR7rFrfd5!+Oluk{9Rm#}|HgmUVVsIHa8m8skcxO}}A>NNmR35kiNxFh~6pIFoPIo%J|z8jlnBv{oypA=N$T&XP; z{=E3Vz$qRl0Vd+&;(jLD>(4hfHA(S>fS{JeTm3WdSqM$;LqT~#e3feEn z@Xi=-WMexqiE-9b_w0n%-Men;d*4pFrhN3A}6P;J|+P@Rql4$*ed1bdAaH7)UXx%2id0Ns$J%n+Yerb&gL? zTYKB`6)Ti~IKAia9`YM5pw}kMp{>wi3x;{NdAq=k8@{q&PU`BCoI%=VY^y(OhF4L& ztFN!;r~p-;EVb?MO1>k*b_@%g{jUMN(HC36>oz&h^lkYy)cD>cBy`zTuIiy7lkctB z8zMw+2x$vca(M6L;YnBDcWk-qudj|WwBz44?T(gJ+8ca!a=%7``dxflSXA_&3eQOm z4Z}LFR?7Oe=k{G97PxDlM;x{KRuP9m-j6E<4x~L8E)es@RnrG#_CDfB-16e5tc*;S zrvgJGbMsAhHa4BoRx`ZK|E#hY96Hip)ldY{MUYb|uyRk`LnHmZuhtaT;GmjaW;e5L zlEDzm&Z%Hv+0tbHzKU1y^gZxgn+gtZ8h_J5oLCUAV^u|FhnH6+6J!tf-m~OFX>)kj zYDc-@DV=rNlS0d-X`$o3`^DbAo>oSdsrA16)njcRrRD{aAa1_qY`Q#*C)ly93p4)uz2a~Oy6<3nqB zG?R;`Cwd9&gLXl~bZGUR4J(@pfw?mOf{iVZdq$)o1My8O9jV-#o znk|kZ>^v`=W!`=_HkRYV6ES8A`N;UV<%bB-EpqpEHnvRttdWp|7=YtZb8&G|ig{8C z+^&zb=q=1nK00FkT;UQcIn+nCRE}rQpGyHKp;9Lb7#SIHaBy^iqk+xUw`lT(!~oGI z;z2{&>hA4b5$gN?9MS3)&U_7yt&dZU!74L8K7Nwtp1=PpO861$RX{J%W^(|YhVLJC zN8IQr@!4Or8jHd{ovf~}?}DbMS}lf;zL1J?g?v)@@O>Y}H8F87_~9U9AX8rvO|I&R zH^M1C``#F@<6rJu-*P}@#CCEwJ=L)Tl}~8)96h=n-05@vjnuOEO zxLT~+cGfe0n@wk4&rHV7EIM|q1GK#ju?38=`+bxG4Fp14%=SM+ZEL)4_FOg6|ANcH zrZ|-dw}^K4U0TqsWGno4@avb;7#J8_F&iebLO78X`j~6Z%w-<(^gjqFk7=Qis?X?l z;AmQK$yDdbqmfNz|ByC^Y-vi^sLiRDgte0~=ZFfhr)Bid!koN{N_kg_yT%%Z@ve{d z-$IQne|%MFd+pLG9Kjh>S;eby;y#5{>ZJ3usz zL1wwRxjA+d(TV<%kt_uT1@|95e182chyabmpK11%Ow$Tu3e%C-4}1Dmbuh=(q_w|3 zip4hxY=LSdr=~_PB`uw%SGLN)$cX*O5p9kA;^I5Om(|7Iat15%J4{_ISeR-h45_EL zSMBs^_pe_K&gEaHg?1pIC2znAz~!a9`c!`3-^~5;_f5DrC#o-JAF8(s)SSQi%}R1s` ztL$R_5&@L%qa^L&6;il%{kkk*Q($0VM$-|?F4rqFKs!WPVW0HnJtiM38ie{Iu-9jM z5UKkH=KtF&`fq=E`>hJfy*%DncvKileJO``vM!;~nH`t@8XhwuIJDn#i|hamtO)<} z6$@_0LXx2!E_!3V8*77h$ddUW$Z5Y@!L%v41d50R>U7T+qsZgX@^%GTb)UpHP3C_o2Kw(1 z^#sH3id>xM0t%#Q6;62rK`f{2kaZ?S{uUrorSVRL?8VwQ zZ`^nr^HqA%qcZuGz==kNzq$40giw3j^0!#6z(Y0zC!tG$5G+mvd-I+HG9kRaxUKwrBL3tBq$c<%Y%nM0*uLXX`^{nY$hwas;;hXd~WP@FUkXj0?2v$ z&K-Twd@pZr{x^DLO;?f~;}}dzEh=^S^OL17gF4YtgBz#GBOTjSaj}$VkIOWM{E~g z@4f(PPV0~3bs3z(kXg2DnXtnY@7}!` z`zdE3AEz7Duz`d|O5+yjLm%iWoEqyWwE8`Ip7(6VY5?H+&$=g1x?+hT>V?znuiE~T zIFtKA>Psv1Q~-&!6H2 ze*Jp8?|lf0*fx;kQWlr_Td0B3W^)(?ez$K^VPzU@*Cu2CbghGc0#v!d`Z#?;m%znQ zfR|A^IyzRbTbI;aY)(FPLjS&@0`!ZQxWla`H)Q-|u^KftLv{Yp#*7-X(S5(=34;^= z7BQ}hsSl2YObaYdvrT|%=gyt-3JM7~bg-yHejaE$#&#kY^aempo`uOgZxetd|SsP~c!_?8&m#3zTVkzDLf{`DJwYT4D`qO2%;*L$IZ|LLEsg zA;O!wrKU6vgt4|4Eh5Ges6scx7R0(;Z2e}zA2S@s}-v3{r7$M zQnb#Vy&E0P7BQ9b?@*cf;+V@KGf)%<2xxa-pT|xy`-_kcMJ8_zSb1_~l>*aJT(_B) z@k0~%c^yN9MmCkhxaW=c85l_V!ztS7tDti=7zU2+&bmwaQEIlIXgSyZ`Bb#Sm;g6=D@)(!8toqgCXJe+`M^0ADyG zQaK{T95hEJM_RK5s~*6_v>a&>3_p640Vpl`H@`}xxF4u=X<3j4gm-XPK`&<6`!}c8TiZ%5*DM zZ`g3>^=pZ#(e@-w=Rfl|ZN7h|2OJXAz1g?P@A{3=wqvMf-6JEkSaChpZHYeav-wQ6 z`3~nlGe#*R4>tb~&^Nhs_n}%hZytfZdM7ZD!EYU3-lrMGI%v1*-fF;f%a3DQjw!DDbDk|vx$m|W6ao=3&0fMP z6bXt2G}UiP)=q53TLK=2EnPpX>|V{@WAwW~Zl@J>u7YY1u9H_1iac^U)&_ z;o&kt2XXiq$KGwGuOE6^o!#=^Y;9n_W2tei5c19Ls~p{niuTKDp|n(rB6;!R(<}Ly ze?-~W)z!VP_{^kretw>EG-r3*z9j?vl)$@svpguCrEnt-@r4)~7{q7!Kp>MwVa1~5 z4bqZi#anSt%EyuesIyxU|@LJfoc7j4ajN&ua`$<&2i&i@UZG~dP~ ztt9h4?u~1cb+Xs<^YbTiwO+HfR)g*g&XUmLzwO|v*MDYtgsy&eg&n!A*uB*ThLP-@ z^*b_jJxG1Cf7GFHN>OpGkVQxG!>!zlivy=ECzS&)UAlC@q+t(iN*0uf)wHy4A>7=L zi0}jO(z$}<*G0m|VWFX=z3Lc7Z7p|`1@Hw2}6+M;(_%V_= z%lyIuQUB**YQX>AC_X*tznJL_8eM7=xVK0}6BuUm0!7bq6Ta19*o? zs1-3{Ag~SV)g8R*OCv3r#C-vh9!Gia9*m1@;#xZ8E}W7qdWD5z;(w;~r|0H++j7C@ zA!aS7JaL-RUQz+mFVjA;AL5Rf8o@FKLddauw{eO!Ji*#)P;4MD>&Erchf4hy!S;j= zV0{Df{m|3q;*sQJUNWsft0)Z?soY|V^FLRwUtbQ%XY$BKwV05sqs$8WHJ8f5LIMs= z@nV4p%ReC_qfT8I7z^Td$hb~s@@%enYs~PoO!Iaom-&f1;&<3g8i0cWo+OK>9JjwF zj-Nj^0Lj{K$Coh1dhm1-o3TJqT*Je+XWApuA!BBivh=;X*WMU;H){fHev8uz+{nP- zX4?EBL03*zR*H@a2E(yTn@VH{tVY|8L~PF4{YnT<<$rwvw}V6X0+Wc{^ zBP+{x@E{``A6KM3Pz6~7F9My4UyhVDV`2Sc~no}EpA(fyu+{R0Sec6K>w8>9_v{~aV3@-)q+*M9oH3PfLa zM!euq4i~0p@TP)GEe&~uQziV&SQh>bO}8P}T_!7jPn?wL5sQ4>zo&{q;I#$8zG=1> z6|x4~|HdQ%-GdJ}#6(WcSs@)1uP4RUP9H`Wd;R@Bh)St!_}cwi=$>k2kgGfZgBo86{U1WziMO%!;AP({=dpx9 zc!ga4x?7AR+P>EA3NexwAl}B?}`RD`m1_eXELj9p>+`v}n zp0&~{W|iH~yA1SdE4CSwn_=j{66m@o@kbJdC8lL^-pUXZVK=wsST9!sf_n;Yrx$Zu z*TX=+?Ma719VR%zPwJYQnDl_1)6>)OS)F;d>wLECFUM&PTr8$IZIF{;TRM3+24WqS z+zBLJpk(y~Yw4C4MLMm588Za=ej6Bxu;cVg-5tBxMjou|5%(qM8I()qQCOU9ehvTb z4y3AGyu3;o{_1;=U%9g9ReCyVAHCSBEpK+S3GiZ~cDgL?03vCA%*rij*0OQSe%%JV zC-n^w43yw>(#WK?nyvL%j_}nI4o$52iyEnji zZXB09Z(`3(jdiAGR>OI^fGNs|1#ouzUCi6c*k4RRzu*^MT>P^zJrHwm9J8q>A@0d? z3Z6|5)P*;g*%OfSOY7@Ng*zaBhp;83V>mHVPOymI_Pq6(KthA45g6xKF-u8UXZ3RE zD=+OPE}?+GIqAenE=Z)1E9m>Ds1-zyf#baG^}-ahd0Wmcm~vRqqj4}hl#^}{B`=~e z(}(&Vd$=26iV`n+cDoaVMJ9Xd{*{?a4nFpRV;STidLjmCXgteQfrLe4E6i4R8k5s9mfz?Yi#l zERH*W#7Zy7hIJ%$%7llD?Ue>r-OTXi6Xncce8k=+rqlP%*DRQI$<5+G@}7WSjdPDk zl~u&u!7$WZScD}nZ2ir*DKC#hq2eA(1>SqWg9o?4!BMsNDr~qxq8}e58yAf(M^E6&N)e#fjksCIq^#5qv~clzrE(W$QNclE5ao&~RIY z9E;O!j3a2(VN1Hqpi(Law}3$1$%T{F;Kdnct((CsP5NqUYbmY=jj9_$Ir_eR>tR`( zrfzz1KB@leLrkx=9#WHM`9fAuz$%8$nIG4xa@uw+@}bwiV%(PDp%X}`Uc4{1?Ig$) zuv*%UG2P4`PzYo27^Y{3Qzxo;Qm>dv02!hvc0q}x>}(ryK$dVE%V6`(HY{j*OVK%< ztgE+td`_%jL*kF^G|&G170@eAo;nq6RdRgG=9xC<55?HjQILtmFfRU)9I$Z4)g`vPZBLi^OI#N=dkvn{8g zv^HOG*#T}U1B+{@G5Ze6jW?$e7>w*p`N1aVeYCw9d|7mH``)ZRpki7k+&H||O;;f> zHFZCR9p-n`UTbUX6%=F^5~`33xC?c&3)Y0NRS(trokzW(rfb;nEho}`&#`TaBh{I7b|v%i@)ZSp|nF+Rt4xq`X4npEhKty)g;KWW$<1%u zGFvag0EDN##5IT{1;#yqN(XF5MPd+N(&dJo<@~pM|E2Js)p&k#&r9=mkvDJN#6P+E z{^Q5DP>E7uW&u)zG~OLeYB)f-zITYA*r!<95K<+UH$!Wnf~Cc0P`>5 z11SM+Kuf~Btk2mo_45BltiJtYSs_?Y)cO4D(Q`Cw)|9-|D`VQog?sz5v+xB65-dr& z?lIf9Z?_z7+Sh0bh!2VWb+r(&r2#hyaVhvSI}GYw3Ssk8lGc8wJS7?02i7l?MtI;Y z)(2Ns*g{ zH-QZVf(;FL=Ba5mqCfon zY_FsUJ&2Cva;{led9cZoJ;R!jQ?ocyV3d>W9OYBNP3uDz z9osR?q~Tn4xGXyTO%GQM9wIeAZBWos6n^{;5gXS#cc>}x<_SQa@k4OP*FoBxM#|_H zz|#^z2uyY|_Mc?91ZK;}%eTO(sEZS0$(? z0w2fIzpLUsUOvV4VoX-X?=cn|lnLkp8rjWx*GJYKqkD-|5s^M%<8?aE?IX!YMEKlo z9(nQVz~tcY6Z4_9(3P$GoDL`1Y!qvjpqH*fU687{H`E_gX$l1qM0atg*{-3+B)gfu z5PrvLOAanBX~0M%Wds1HDz6zDK4gRLpTaA4BvR{vSG-!?@L(1RXd?P!l3hpCIxL;N zSVjmz1T5jV|8*t9v{?lNF~E?gD&f=FYFO&=X|)dJW{}V)ceJ$jdTh}xUbnyC>DUYJ z_<)JOeCo&)~4b$==OevnP4GarC7M9ZFlateu( zR{XvUTM5cBbK4UyrhE84r<_c~0GLX9j+M0jvewz2<=LoG4UY~DT-P~sW-D1sAb`eS z{z8!7D^8K!3doUehS25%2lWyPIyC+jB#8^$5d(fva2B>*V6b)-zmklBm5pMc{qIq~ zzovUHXdN8c{XjEmsCXzN8<9sm{@Qs#IA4gqk-Uh^HGBRxBZ=GeYf*6EdKkuzC(SlD z!XGg}5rBP*#qIv%$9KnGFYW~5kp$$eO}!w2)nr>N#X!BL5L6;j@rc4i3J)JxMYW8jU4d-=noWl+~OQ_ku&#H$T7r?NOYI{G#u zDZ<8ol`RNVM9bbnA|)W{ zm{=mM$`K-)$eTpKbiBsp4FDdeAKteX27gFu^q2e_{&2fjJ?+Ia`1tG8D_M}@8nMDqzxK6KWEiR~vezU4v7Q3ZA*G=atWuHe zbM1m>TiotT;Z1Kh_Y~Mqy%;hbR?)pY;=})?G?(;P)$Uf2F5H3+AC14I^|mixv@;8f zi2Vn3mLyrQoY*!k04An^sQUeYG*tHX?X!-GD^UMNd0R~P;bLeH6uy9aa`8yCke zB(xbIdR1%cGdYI)p5W%&gluU{%lLr-k+SsqLD5qal}Y^9icqp7t!f*Pg`2$L~3p zXQRN*_@T1$q*>~FL{AT(5dx)C52(V@I=$#Mm3O73a2x}@cUmQfF=FEJ2?-@w^E3>> zt58NjM}3Fhk-BoOz)me?1U*Nj@do)pk`XQ{D&iawr56;hIug9I{7A4@C7VU|+c$46 zqB=d2a1n>jN5{l;mzG(GRaBJu-Me>*_7f8m!gk}VRL4#R({2~0Lzw(6q-!gL!85nz zEi5d+1u?BVgIw#@e`!A19u*Px2m-UgWjwxR#e`v0i;l7G-o28zJQ}4>)JUqn z9vT`-alI!#E3dBuFBS>I08PAzaR$OWkqBHJJy#^hD_YVe*g*NGNc`#& zDRxYImLg^WoAgQgku6-yJKhSmsEomK^5aBzF>HMH-U3A#KSt7R_`m&;q%1=z{Tgg0 zj->lCbT#?$-2Y4Y4PWl+Rk<7}=bFlHwu4KTqz4HmX`dH!7#)en4P7mZy<0g>4zGRx zpSmi&V6l$rkNW$ucfEGr)5mS`$}hzlS^^XCA@|>Yl}+SHykcFFPKO_>BQ{1J^s9g5 zYpWc`MZOB_;2sB+dr)ZEuZnknQsh5|%I`tW2@gC0{nb}Su9EAWsw+?zi74qD#QbpcoucIZ; zl)3g3S2r;!(h9xaZTzy`pIjI{?M()MMvcIIy}{%oQNzEM7l_Ga`D>ByRucI2&j6&Y zc)Mov$PpEO$t>5U-z&YhJe--@Vi~RY%9fl={=tjvey?O&3lc2U%u`WUFic{Wm8>?Q5jxFD(zKy7qW$MWhfL3}n_)06i1`k=exA-4CPyNmDt6h z@GYZgs)*mg!@W0lw3gj&H)hg@(lf)Q+!dk<5r0wvk6F+>annadE=&CICi0L4nDvv_ z8QqQRzPjt?mC82trl}$_BG{5HeUIyA`uXPJr1GTMOXCG^aGIxgY|R= z-UYgOz`t3=#0quIe*XMPK(*b?n`&M=a5g6u2OdVRQw76k8fFzDfW zpwJ40yiFUNN#a7En^9-zapl0Nj{9HDQ3gNEXmJ%x-#bA5rU`3QG(CU5Bc$BMv>Gcw z1kf_#k!UQwbM_|4BBcw29p+)v6)kl5A~HjmC{rtG{{R3b-FNNh=fkw4PC^NdCrwrq zNm$|#^O&~R#5|Z2_$t1S{-0V^wcJf=QOO~7p0j*SQ(v0h>1t~j4_otewMu#7=I@D- z1WV?u=fIs;MGXoR*Tuz5%?#(C>#MyY>S*RptyMQ%=v`4;;BFn^WYJNw?W3sh3cLBa zGo8Qt?5`(W32#n6p;5sQ^Zs{AR8jVBH`djwz6349oxt&60U6QYk#PBQROI@1FSN@r z!5_XWn~%nX5+bTUZ-r*(jM1+iKFvEF&3OB*r;~-ehQ^wo#vcOwBFMs#m$9z!#LzelK%`-MS$Q1CM&x^In zg9K9j%(iXYcBMTpv*}pt(icC(byrTK*_#O&4Dgh+l0meqHth6;VqR8KvJ9DSS15mw z$v@Jshgz3YZfWX$NZTNgYN=%eW5ylY98eH7LP39PbEzV&Xi6p&8mVkf;`Bl#|20@o zO|(l#Cnx0w2i&&F_kR}_MVHh^V(n~8Fv4>%SjT%!iP$m$xnz<3kFoWo0{kB{?caS{ zA4=wOhM(g*Jm8#^n26RtIKS%wPTfaa1OG=}CFsZRsWQ`iA*!HQ5zOx86lnb)+Zb`~ zIt{c?92FB|ZqsEO2)nGEMjE>|ZY*wZ-@)Mh=HguN`##hOS=a+4Z~3o7nq0N?aNTU`2ef`HBD=r2$NCGC~aZ6vk$co?wVkPJ^*Go zVM|uh)dw*E*u6P_0c?ric1RuCX(?I}p~$>NMoebdyt8Nrnx}51rynFK7m}LF%WIJg zXi$z6M%u~eT@V+^^jlh49Rn@kH>}!-WQ$3|NPXPBfg{LmKE1_W7pv+!+??7xyI{94 zYwI`bH+gR9%4GDN4qE?77K3p>a_d-a$eD?G1&!5euDynrguxc>0!Ju4I$%0>`Fop-E49%pvCV~5SMe?T?39$0x2PCIY)5*c7@G4x zdU7FHD(~Z-M=H?=VOG-pJ&nfMadf_KAf3IiPDyh%yUcTyNa21UBIJPdLAQv*2kbKe zsT5|`4j6bej|EMiYraM`mby z{X9MxqhmvA(QM*CfBiWL(~x~So`_-*Iycaez(|?{;mmF)O?G&wcM*k2G5a{1yHVZ< zrFS`n5GoB#HML#npq91!@s4ZwYW~$P0?S-2=clbSzFrlF#gqUsMdPIWH+$)-pW zFE4OXsr9BPfz^#@^@;3w^WZ=G?veIE59ejdn^~xPhOe4d*X*{)7L@#d#A8W&Tj3_cKA*^by`Y@3K!65-Nq1bNv%p1} zmbsQs!sxPzQ4@kn&#KfG{hAa)d6+h{q@-95` z1g8S6VQl*Q3K=RxX^yj57P;G zZ0X*ut~1qfyK4H7rnc%XUFpkwSjA=zjyh5gYZ0BYpiySNEPiT>nLs zlWA>7Ys6ZQo;0&asGBI-h{BZC)TpLBL1ez5urPCH*?p8ZAzZ#YG)=e_ja8zyL$|?t z2zZsmo9Zc@l#6>xw&XPR`5UGvP?0BxsOxoxGX-GQq)alb{HiM_2bFCq} zTYPM2t*}P+1S^TULLF%}*u^Jtazxp55?31ou)p{;vW~%r_@Mf6p4?A zyxOOX#H^k9LOBE_k019KXyLEodb}4||Ev`!Q?n3JWeVYBd66`-)f2<$?NzjlLDG68 zdgDboqZ-Pvmcs&i;k!ePo|GVrGXDGZ;tph3EF{N!z5n$* zgGJHIa}$e6SRTh7?-sV?M)+-*kpq&^4KnC6KT{&($#rxU?FK}$Dbb$PnyIyiVSs+O_VV0^`Z=!!KgMYH->@(mL{qe36v}LRNg(UkNkqA)=WFp%M*oe~`hGxm> zBXYNvi!07PlNR{~Wo24OnGHkjD5OYtbKacK!tL=QJxaC-L6`)CaCow0<4&T~^{ zJM7yIczXR*L&V#}>G1EFL+ZX52c7_`<^DT##`ME}oEu6i8VmYyttcz+`LrXG6aa5e zPY*Vy@P=bk+TLEMJHk|d#Oc9jBgQN9{S2-q^B0^Jy_PAt@SRXp9RF~zri(lgjHcnq zKaYP}zvD<)_K$wM6F0CHCyrdgmR`s`ki8;IuD-9Vlmu$ikcfRSXO}UGL8YdnlQg+B zj)T$r;&y8VHEO~J_~`0kasyTmf=p*;XNOP?KOL?n6$V5nvnXuNk%>(fk$U39kzSU; zJAT|X&x**IhZ?vWVCe=M85j_LD;VjU^-~+2FBrR_Uw0tQ|EJjXK%{-SU|%kA-mGiH?qI!K3Pxib6}E&ZK@`>V|!NeNsp@cp&(vnX&+b zh6s5&PSX zLyxVWW8W?8apq>Z_{T}(3oK#5Z~D16e?G@1ZM$JU$)-4=mFIS4rF`aCCj=KAv`Ll1 zfuLnk+F|3Do}8@Y~wWKTzS+3O^=6 zMQn4m1-tUFBVDFQvT@`?Bnqz)cVn!hH0Eafj}Xtu$jI?pIhHkGqTY1dxAUTV&crc3 z`1Mr;Ww=psWC4D+;DoOiixe(xyzNm`!XyI*Btg z#nd5HzaJ!S4K8p~dlAqfHk3MLhg63L;4le>Ziaxy4K$~RVIb?R5LP2Tr(`!)`1>|m z-quiV)3eL2#=oABJsuKf7h#aCr8ERi!+VVata1wr3nR)+H@cPS07LfVw+-0KPkmJzY^Jwl8+yiNqW=ceaL>Foko*GjirpLM|WWpqO}NT(G!g`Dl4mL55KqZ zm2gVNpJ8Nb5<5^pJg|X=3bxi@*uG;2p+zM0P)g*&yW{Ghs1LYL-v%gJjl3kT-W@K` z-Vb4iRpG=(SMc}GH1s=UVb>qa>(@7AcA8GHo*FVT}rlT-P!5;cUi!Gq*7qOd~Rn?Axb|eJdD1>u*CSUPm^R8MrW0 zUiV3O*rfTD;R%gG?5436J$DReGfj4bg0Dlrw-l<82c}AkL1#9t-$Y7w2Lq}Z*#{A| z`xqFT#{sjJLw-~W329GG9HYC99N4Zsd*lb|p7}H_gv}YH{4ny#a#6%Xf3Vu^A-5Yb zk*;oC+Th~66jlHgU;~@|XXE1I zSA$BEz-gh_;wlcdh4DBlZ6$J_GqQs|z{LRX7l|C!tm^RpM_ z3%R@s6bF7diIifj&S}_3>|BxKHXv-$LXAd_RTx=*`iTErgw04aN)?s$iAhPd$Pf|O zL3ZB4wjS?#IN5s1a}eq+{M}Io7*B^#pQ|f{#Jn&bsVP9`U8oC^*u;kHBtkX{0!}D_ zkU?VLn0$@S@C7=Y=gbiPUWKuCrTN7U3`nbqp5u;aEhQ|NWcfrz#U+&*9ufU3jC4xC zLLmDqp`U3fRC=*G{0Iq%;1Y?W40nM-!TNfKU|!sS^ySN(7!~Ss7&x2o(WP7vnmAm6 z-6>CzF~Z{)CmMRB^X!@$vEShdWgKiPdLoa=7#zgIc&>=FwYHB!RPo{YQMh}vPta*~6c;A&;D+Y>TH#CkNqSv=du%h1^sSE=~Y`8#=ti9S! zIywpkTi+!$nPppaY6hcmZ*yZ~<3LjiGfAGdkh~2x5Us4La+8WWp#Iytl7?7GSEXQd z$WWh#gS&!)^vI#CcLoBKViYG0L+toNP5IuI%Lp8h{PPj$S+w^?gxO#?RdOCRr?YjU z*5-+OXj=cA_H4HKc_TCf+?loY_u^$Brp}H%slpJ-^|7(o4A=(tt6v*Z|uyJeSgOG6p0St}0A!8s<(j2idL{H6eUJyYsxs5C< z!3g6&e@Y{Ej6FP00%3r)F729?o3AKizhu6$gSa}5mrxV)ghbzv)*w9wuAp?JK?~s> z(1+%>b`CZ%1**OnCNdFYOZorrcdGs0&fS0grPm^LzqxAdw Date: Thu, 31 Oct 2024 15:41:18 +0100 Subject: [PATCH 12/13] Update documentation Hamiltonians --- .../source/reference/Hamiltonians/index.rst | 30 ++++++++++++++----- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/documentation/source/reference/Hamiltonians/index.rst b/documentation/source/reference/Hamiltonians/index.rst index 82f2d606..3284bc5b 100644 --- a/documentation/source/reference/Hamiltonians/index.rst +++ b/documentation/source/reference/Hamiltonians/index.rst @@ -7,7 +7,7 @@ Hamiltonians This Hamiltonians submodule of Qrisp provides a unified framework to describe, optimize and simulate quantum Hamiltonians. It provides a collection of different types of Hamiltonians that can be used to model and solve a variety of problems in physics, chemistry, or optimization. -(Up to now, Pauli Hamiltonians have been implemented, but stay tuned for future updates!) +(Up to now, Pauli Hamiltonians and Fermionic Hamiltonians have been implemented, but stay tuned for future updates!) Each type of Hamiltonian comes with comprehensive documentation and brief examples to help you understand its implementation and usage: .. list-table:: @@ -16,12 +16,8 @@ Each type of Hamiltonian comes with comprehensive documentation and brief exampl * - Hamiltonian - USED FOR - * - :ref:`Hamiltonian ` - - unified framework for quantum Hamiltonians * - :ref:`PauliHamiltonian ` - describe Hamiltonians in terms of Pauli operators - * - :ref:`BoundPauliHamiltonian ` - - describe Hamiltonians in terms of Pauli operators, bound to specific QuantumVariables * - :ref:`FermionicHamiltonian ` - describe Hamiltonians in terms of fermionic ladder operators @@ -32,7 +28,27 @@ We encourage you to explore these Hamiltonians, delve into their documentation, .. toctree:: :hidden: - Hamiltonian PauliHamiltonian - BoundPauliHamiltonian FermionicHamiltonian + + +Examples +======== + +.. list-table:: + :widths: 25 25 + :header-rows: 1 + + * - Title + - Description + * - :ref:`VQEHeisenberg` + - Investigating the antiferromagnetic Heisenberg model with :ref:`VQE `. + * - :ref:`VQEElectronicStructure` + - Solving the electronic structure problem with :ref:`VQE `. + * - :ref:`MolecularPotentialEnergyCurve` + - Calculating molecular potential energy curves with :ref:`VQE `. + * - :ref:`GroundStateEnergyQPE` + - Calculating ground state energies with :ref:`quantum phase estimation `. + * - :ref:`IsingModel` + - Hamiltonian simulation of the transverse field Ising model. + From 0a981f1fb47456042e47f0e855e364812bd25c04 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Thu, 31 Oct 2024 16:10:22 +0100 Subject: [PATCH 13/13] Additional tests for Hamiltonians and VQE --- .../test_fermionic_to_pauli.py | 21 +++++++++++++++- tests/test_VQE_electronic_structure.py | 24 +++++++------------ 2 files changed, 29 insertions(+), 16 deletions(-) diff --git a/tests/hamiltonian_tests/test_fermionic_to_pauli.py b/tests/hamiltonian_tests/test_fermionic_to_pauli.py index fad612ed..ec5add51 100644 --- a/tests/hamiltonian_tests/test_fermionic_to_pauli.py +++ b/tests/hamiltonian_tests/test_fermionic_to_pauli.py @@ -18,6 +18,8 @@ from qrisp.operators.fermionic import a, c from qrisp.operators.pauli import X,Y,Z +from qrisp.vqe.problems.electronic_structure import * +from pyscf import gto def test_fermionic_to_pauli(): @@ -38,4 +40,21 @@ def test_fermionic_to_pauli(): + Y(0)*X(1)*X(2)*Y(3) + Y(0)*X(1)*Y(2)*X(3) \ - Y(0)*Y(1)*X(2)*X(3) + Y(0)*Y(1)*Y(2)*Y(3)) - assert str(G1-K)=='0' \ No newline at end of file + assert str(G1-K)=='0' + + +def test_hamiltonian_H2(): + + K = -0.812170607248714 -0.0453026155037992*X(0)*X(1)*Y(2)*Y(3) +0.0453026155037992*X(0)*Y(1)*Y(2)*X(3) +0.0453026155037992*Y(0)*X(1)*X(2)*Y(3) -0.0453026155037992*Y(0)*Y(1)*X(2)*X(3) \ + +0.171412826447769*Z(0) +0.168688981703612*Z(0)*Z(1) +0.120625234833904*Z(0)*Z(2) +0.165927850337703*Z(0)*Z(3) +0.171412826447769*Z(1) \ + +0.165927850337703*Z(1)*Z(2) +0.120625234833904*Z(1)*Z(3) -0.223431536908133*Z(2) +0.174412876122615*Z(2)*Z(3) -0.223431536908133*Z(3) + + mol = gto.M( + atom = '''H 0 0 0; H 0 0 0.74''', + basis = 'sto-3g') + + H = create_electronic_hamiltonian(mol).to_pauli_hamiltonian() + + G = K-H + G.apply_threshold(1e-4) + assert str(G)=='0' \ No newline at end of file diff --git a/tests/test_VQE_electronic_structure.py b/tests/test_VQE_electronic_structure.py index d27bbc54..c4933a96 100644 --- a/tests/test_VQE_electronic_structure.py +++ b/tests/test_VQE_electronic_structure.py @@ -16,22 +16,22 @@ ********************************************************************************/ """ -#from qrisp.vqe.problems.electronic_structure import * +from qrisp.vqe.problems.electronic_structure import * +from pyscf import gto +from qrisp import QuantumVariable +import numpy as np # # H2 molecule # -""" + def test_vqe_electronic_structure_H2(): - - from pyscf import gto - from qrisp import QuantumVariable mol = gto.M( atom = '''H 0 0 0; H 0 0 0.74''', basis = 'sto-3g') - H = create_electronic_hamiltonian(mol) + H = create_electronic_hamiltonian(mol).to_pauli_hamiltonian() assert np.abs(H.ground_state_energy()-(-1.85238817356958)) vqe = electronic_structure_problem(mol) @@ -40,12 +40,10 @@ def test_vqe_electronic_structure_H2(): for i in range(5): res = vqe.run(QuantumVariable(4), depth=1, - max_iter=50, - mes_kwargs={'method':'QWC'}) + max_iter=50) results.append(res) assert np.abs(min(results)-(-1.85238817356958)) < 1e-1 -""" # # BeH2 molecule, active space @@ -53,14 +51,11 @@ def test_vqe_electronic_structure_H2(): """ def test_vqe_electronic_structure_BeH2(): - from pyscf import gto - from qrisp import QuantumVariable - mol = gto.M( atom = f'''Be 0 0 0; H 0 0 3.0; H 0 0 -3.0''', basis = 'sto-3g') - H = create_electronic_hamiltonian(mol,active_orb=6,active_elec=4) + H = create_electronic_hamiltonian(mol,active_orb=6,active_elec=4).to_pauli_hamiltonian() assert np.abs(H.ground_state_energy()-(-16.73195995959339)) # runs for >1 minute @@ -70,8 +65,7 @@ def test_vqe_electronic_structure_BeH2(): for i in range(5): res = vqe.run(QuantumVariable(6), depth=1, - max_iter=50, - mes_kwargs={'method':'QWC'}) + max_iter=50) results.append(res) print(min(results))