from qiskit_nature.second_q.operators import BosonicOp, FermionicOp

# Create a creation operator
bosonic_op = BosonicOp({'+_0': 1}, num_modes=1)
# Create an annihilation operator
bosonic_op = BosonicOp({'-_0': 1}, num_modes=1)

bosonic_op1 = BosonicOp({'+_0 -_0 -_1 +_1': 2}, num_modes=2)
bosonic_op2 = BosonicOp({'+_2 -_2': 1.5, '+_0': 3}, num_modes=3)
# Sum
print("Sum")
print(bosonic_op1 + bosonic_op2)
# Multiplication
print("\nMultiplication")
print(0.5 * bosonic_op1 @ bosonic_op2)

bosonic_op = BosonicOp({'+_0 -_0 -_1 +_0': 1}, num_modes=2)
# Normal order (reorder all operators such that all creation op are in front). Due to the commutation relations, when inverting -_0 +_0 we get 1 + +_0 -_0, which explains the two terms
print("Normal Ordering")
print(bosonic_op.normal_order())
# Index order (reorder all operators such that the lowest indices come first)
print("\nIndex Ordering")
print(bosonic_op.index_order())

# Simplify
print("Simplify")
# Simplify can be used to reduce numerical noise:
print(BosonicOp({'+_0 -_0': 1}) + BosonicOp({'+_0 -_0': 1e-10}).simplify(atol=1e-9))

print("\nSimplify with two creation operators")
print("\nBosons:")
print(BosonicOp({'+_0 -_0 -_1 +_0 +_0': 1}, num_modes=2).simplify())
print("\nFermions:")
print(FermionicOp({'+_0 -_0 -_1 +_0 +_0': 1}, num_spin_orbitals=2).simplify())

FermionicOp({'+_0 -_0 -_1 +_0': 1}, num_spin_orbitals=2).simplify()

from qiskit_nature.second_q.mappers import BosonicLinearMapper
from qiskit_nature import settings
settings.use_pauli_sum_op = False

mapper_occ1 = BosonicLinearMapper(max_occupation=1)
boson1 = mapper_occ1.map(BosonicOp({'+_0 -_0': 1}))
print(boson1)

print("Max occupation = 1")
boson_p1 = mapper_occ1.map(BosonicOp({'+_0': 1}))
boson_m1 = mapper_occ1.map(BosonicOp({'-_0': 1}))
print(boson_p1)

print("\nMax occupation = 2")
mapper_occ2 = BosonicLinearMapper(max_occupation=2)
boson_p2 = mapper_occ2.map(BosonicOp({'+_0': 1}))
boson_m2 = mapper_occ2.map(BosonicOp({'-_0': 1}))
print(boson_p2)

print("\nMax occupation = 3")
mapper_occ3 = BosonicLinearMapper(max_occupation=3)
boson_p3 = mapper_occ3.map(BosonicOp({'+_0': 1}))
boson_m3 = mapper_occ3.map(BosonicOp({'-_0': 1}))
print(boson_p3)

