diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..27c0fda --- /dev/null +++ b/.travis.yml @@ -0,0 +1,28 @@ +notifications: + email: false + +matrix: + include: + - sudo: required + services: + - docker + env: DOCKER_IMAGE=quay.io/pypa/manylinux1_x86_64 + PLAT=manylinux1_x86_64 + - sudo: required + services: + - docker + env: DOCKER_IMAGE=quay.io/pypa/manylinux1_i686 + PRE_CMD=linux32 + PLAT=manylinux1_i686 + - sudo: required + services: + - docker + env: DOCKER_IMAGE=quay.io/pypa/manylinux2010_x86_64 + PLAT=manylinux2010_x86_64 + +install: + - docker pull $DOCKER_IMAGE + +script: + - docker run --rm -e PLAT=$PLAT -v `pwd`:/io $DOCKER_IMAGE $PRE_CMD /io/travis/build-wheels.sh + - ls wheelhouse/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt deleted file mode 100644 index c9ebe92..0000000 --- a/CMakeLists.txt +++ /dev/null @@ -1,55 +0,0 @@ -cmake_minimum_required(VERSION 3.4...3.18) -include(CMakePrintHelpers) -project(pyket LANGUAGES CXX) -set(CMAKE_BUILD_TYPE RelWithDebInfo) -#set(CMAKE_BUILD_TYPE Debug) -set(DEFAULT_NUM_MODES 2) -set(DEFAULT_NUM_BITS 8) - -set(rapidjson_DIR home/ethan/libraries/builds/rapidjson/include) -set(Eigen3_DIR /home/ethan/libraries/builds/eigen_3_install) -FIND_PACKAGE( Boost REQUIRED ) -find_package(Eigen3 3.3 NO_MODULE REQUIRED) - - - -if(NOT DEFINED NUM_MODES) - message(STATUS "Using Default Number of modes") - set(NUM_MODES ${DEFAULT_NUM_MODES}) -endif() -if(NOT DEFINED NUM_BITS) - message(STATUS "Using Default Number of bits per mode") - set(NUM_BITS ${DEFAULT_NUM_BITS}) -endif() -#-------DEFINE MODULE NAME------------ -set(pyket pyket${NUM_MODES}_${NUM_BITS}) - -add_compile_definitions(NUM_MODES=${NUM_MODES} NUM_BITS=${NUM_BITS}) -cmake_print_variables(NUM_MODES) -cmake_print_variables(NUM_BITS) - -# specify the C++ standard -set(CMAKE_CXX_STANDARD 17) -set(CMAKE_CXX_STANDARD_REQUIRED True) -add_library(HamiltonianFuncs STATIC "") - -add_subdirectory(src) - - -add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/extern/pybind11) -pybind11_add_module(${pyket} SHARED ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cpp) - - -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/src/pyket_config.hpp.in ${CMAKE_CURRENT_SOURCE_DIR}/src/pyket_config.hpp @ONLY) -include_directories(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src) -target_include_directories(HamiltonianFuncs SYSTEM PUBLIC ${rapidjson_DIR}/include) -target_link_libraries (HamiltonianFuncs Eigen3::Eigen) -target_compile_options(HamiltonianFuncs PRIVATE -fPIC) -target_link_libraries(${pyket} PRIVATE HamiltonianFuncs) -find_package(OpenMP) -if(OpenMP_CXX_FOUND) - target_link_libraries(${pyket} PUBLIC OpenMP::OpenMP_CXX) -endif() - -install(TARGETS ${pyket} DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/pyket) - diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..ba7c18a --- /dev/null +++ b/environment.yml @@ -0,0 +1,18 @@ +name: adaptive-spin-evolution +channels: + - conda-forge + - defaults +dependencies: + - python>=3.6 + - pybind11=2.10 + - pip>=21.3 #For editable installationhttps://pip.pypa.io/en/latest/reference/build-system/pyproject-toml/?highlight=project%20toml#editable-installation + - numpy>=1.16 + - scipy>=1.0 + - matplotlib>=1.2.1 + - pytest>=5.3 + - pytest-rerunfailures>=5.3 + - pymeep=1.24 + - pymeep-extras=1.24 + - pyscf=2.0.1 + - pip: + - pyscf[geomopt]==2.0.1 \ No newline at end of file diff --git a/extern/eigen b/extern/eigen index 4a56359..3147391 160000 --- a/extern/eigen +++ b/extern/eigen @@ -1 +1 @@ -Subproject commit 4a56359406eca3b041b06641656aad06315d4461 +Subproject commit 3147391d946bb4b6c68edd901f2add6ac1f31f8c diff --git a/extern/pybind11 b/extern/pybind11 index f1abf5d..aa304c9 160000 --- a/extern/pybind11 +++ b/extern/pybind11 @@ -1 +1 @@ -Subproject commit f1abf5d9159b805674197f6bc443592e631c9130 +Subproject commit aa304c9c7d725ffb9d10af08a3b34cb372307020 diff --git a/extern/rapidjson b/extern/rapidjson index 232389d..f54b0e4 160000 --- a/extern/rapidjson +++ b/extern/rapidjson @@ -1 +1 @@ -Subproject commit 232389d4f1012dddec4ef84861face2d2ba85709 +Subproject commit f54b0e47a08782a6131cc3d60f94d038fa6e0a51 diff --git a/extern/test b/extern/test index 5ac090e..d2895eb 160000 --- a/extern/test +++ b/extern/test @@ -1 +1 @@ -Subproject commit 5ac090e82cc44f50ef34a990524879a9d196a1c2 +Subproject commit d2895ebfdfdf16074c58c9801d53e190c4654fcb diff --git a/pyket/compare.py b/pyket/compare.py deleted file mode 100644 index 9813bca..0000000 --- a/pyket/compare.py +++ /dev/null @@ -1,107 +0,0 @@ -import pyket -import qutip -import matplotlib.pyplot as plt - -import numpy as np -fock_N1 = 40 -fock_N2 = 40 -#basic enery parameters -w0 = 3 -v0 = 4 -g0 = .2 -w1 = 4 -v1 = 5 -g1 = .07 - -#step size -dt = .01 - - - -p = pyket.Params() -pinit_state = pyket.StateVector() - -ket = pyket.StateKet() -ket.set_spin(True) -ket.set_amp(complex('1')) -pinit_state.push_back(ket) -ket = pyket.StateKet() -ket.set_mode(0,1) -ket.set_amp(complex(0)) -pinit_state.push_back(ket) -ket = pyket.StateKet() - - -d = { - "run_id":1420, - "output_directory":"~/code/scratch", - "initial_state":pinit_state, - "mode_energies":[v0,v1], - "mode_couplings":[g0,g1], - "up_energy":.5*w0, - "down_energy":-.5*w0, - "energy_cutoff":-1, - "dt":.01 - -} - -p = pyket.Params() -p.load_dict(d) -h = pyket.H(p) -print(pinit_state) -print(qinit_state) -dt = .01 -N = min(fock_N1,fock_N2)-1 -up_list_qu = np.empty(N) -up_list_py = np.empty(N) - -t = np.empty(N) -for i in range(min(fock_N1,fock_N2)-1): - h.run_grow() - -for i in range(min(fock_N1,fock_N2)-1): - print('--------------ITERATION ' +str(i) + '--------------------') - h.blas_evolve(complex(0,-dt)) - qinit_state = qinit_state +complex(0,-dt)*H*qinit_state - up,down = h.get_blas_spin_pop() - up_list_py[i] = up - up_list_qu[i] = abs(qinit_state.unit().ptrace(0)[0][0][0]) - t[i] = dt*i - -h.reset() -qinit_state = qutip.tensor(spin,vac0,vac1).unit() -for i in range(min(fock_N1,fock_N2)-1): - h.run_grow() - -for i in range(min(fock_N1,fock_N2)-1): - print('--------------ITERATION ' +str(i) + '--------------------') - h.blas_evolve(complex(0,-dt)) - qinit_state = qinit_state +complex(0,-dt)*H*qinit_state - up,down = h.get_blas_spin_pop() - up_list_py[i] = up - up_list_qu[i] = abs(qinit_state.unit().ptrace(0)[0][0][0]) - t[i] = dt*i - - -h.reset_with_state(pinit_state) -qinit_state = qutip.tensor(spin,vac0,vac1).unit() -for i in range(min(fock_N1,fock_N2)-1): - h.run_grow() - -for i in range(min(fock_N1,fock_N2)-1): - print('--------------ITERATION ' +str(i) + '--------------------') - h.blas_evolve(complex(0,-dt)) - qinit_state = qinit_state +complex(0,-dt)*H*qinit_state - up,down = h.get_blas_spin_pop() - up_list_py[i] = up - up_list_qu[i] = abs(qinit_state.unit().ptrace(0)[0][0][0]) - t[i] = dt*i -print(h.get_size()) - - -fig, ax = plt.subplots(2,1,figsize=(9,4)) -ax[0].plot(t,up_list_py,label="pyket_blas_evolve") -ax[0].legend() -ax[1].plot(t,up_list_qu,label="qutip") -ax[1].legend() -plt.show() diff --git a/pyket/just_compare.py b/pyket/just_compare.py deleted file mode 100644 index dd15afc..0000000 --- a/pyket/just_compare.py +++ /dev/null @@ -1,50 +0,0 @@ -import pyket -fock_N1 = 3 -fock_N2 = 3 - -#basic enery parameters -w0 = 3 -v0 = 4 -g0 = .2 -w1 = 4 -v1 = 5 -g1 = .07 -p = pyket.Params() -pinit_state = pyket.StateVector() - -ket = pyket.StateKet() -ket.set_spin(True) -ket.set_amp(complex('1')) -pinit_state.push_back(ket) -ket = pyket.StateKet() -ket.set_mode(0,1) -ket.set_amp(complex(0)) -pinit_state.push_back(ket) -ket = pyket.StateKet() - - -d = { - "run_id":1420, - "output_directory":"~/code/scratch", - "initial_state":pinit_state, - "mode_energies":[v0,v1], - "mode_couplings":[g0,g1], - "up_energy":.5*w0, - "down_energy":-.5*w0, - "energy_cutoff":-1, - "dt":.01 - -} - -p = pyket.Params() -p.load_dict(d) -h = pyket.H(p) - - - -for i in range(min(fock_N1,fock_N2)-1): - print('--------------ITERATION ' +str(i) + '--------------------') - h.run_grow() - h.run_step(complex(1,0)) - print(h.get_state()) - diff --git a/pyket/one_compare.py b/pyket/one_compare.py deleted file mode 100644 index 77afae1..0000000 --- a/pyket/one_compare.py +++ /dev/null @@ -1,85 +0,0 @@ -import pyket -import qutip - - -import numpy as np -fock_N1 = 2 - -#basic enery parameters -w0 = 3 -v0 = 4 -g0 = .2 - -#step size -dt = .01 - -#--------qutip code---------- -#basic operators - -s_eye = qutip.operators.identity(2) -n_eye1 = qutip.operators.identity(fock_N1) - - -oz = qutip.tensor(qutip.operators.sigmaz(),n_eye1) -oplus = qutip.tensor(qutip.operators.sigmap(),n_eye1) -ox = qutip.tensor(qutip.operators.sigmax(),n_eye1) -a0 = qutip.tensor(s_eye,qutip.operators.destroy(fock_N1)) - - - -#make initial state -vac0 = qutip.states.basis(fock_N1) -spin = qutip.states.basis(2) -qinit_state = qutip.tensor(spin,vac0).unit() - - - -#hamiltonian -H = 0.5*w0*oz + v0*a0.dag()*a0 + g0*ox*(a0.dag()+a0) - - - -p = pyket.Params() -pinit_state = pyket.StateVector() - -ket = pyket.StateKet() -ket.set_spin(True) -ket.set_amp(complex('1')) -pinit_state.push_back(ket) -ket = pyket.StateKet() -ket.set_mode(0,1) -ket.set_amp(complex(0)) -pinit_state.push_back(ket) -ket = pyket.StateKet() - - -d = { - "run_id":1420, - "output_directory":"~/code/scratch", - "initial_state":pinit_state, - "mode_energies":[v0], - "mode_couplings":[g0], - "up_energy":.5*w0, - "down_energy":-.5*w0, - "energy_cutoff":0, - "dt":.01 - -} - -p = pyket.Params() -p.load_dict(d) -h = pyket.H(p) -print(pinit_state) -print(qinit_state) -dt = .1 -for i in range(fock_N1-1): - print('--------------ITERATION ' +str(i) + '--------------------') - - print(h.get_matrix()) - print(H) - h.run_step(complex(1,0)) - qinit_state = qinit_state +complex(1,0)*H*qinit_state - h.run_grow() - print(h.get_state()) - print(qinit_state) - print(qinit_state.unit()) diff --git a/pyket/plot_2_modes.py b/pyket/plot_2_modes.py new file mode 100644 index 0000000..7787e0d --- /dev/null +++ b/pyket/plot_2_modes.py @@ -0,0 +1,80 @@ +import matplotlib.pyplot as plt +import numpy as np + +import pyket2_8 as pyket + +# basic enery parameters +w0 = 1.0 +v0 = 1.0 +v1 = 0.0 +g0 = 0.5 +g1 = 0.0 + +# setup run +dt = 0.01 +p = pyket.Params() + +ket1 = pyket.StateKet() +ket1.set_spin(True) +ket1.set_amp(complex(1, 0)) +ket1.set_mode(0, 0) +ket2 = pyket.StateKet() +ket2.set_spin(True) +ket2.set_mode(0, 1) +ket2.set_amp(complex(0, 0)) +ket3 = pyket.StateKet() +ket3.set_spin(False) +ket3.set_amp(complex(0,0)) +ket3.set_mode(0,1) +ket4 = pyket.StateKet() +ket4.set_spin(False) +ket4.set_mode(0,0) +ket4.set_amp(complex(0,0)) + +pinit_state = pyket.StateVector() +pinit_state.push_back(ket1) +pinit_state.push_back(ket2) +pinit_state.push_back(ket3) +pinit_state.push_back(ket4) + +d = { + "run_id": 1420, + "output_directory": "~/code/scratch", + "initial_state": pinit_state, + "mode_energies": [v0, v1], + "mode_couplings": [g0, g1], + "up_energy": 0.5 * w0, + "down_energy": -0.5 * w0, + "energy_cutoff": 10, + "dt": dt, +} +N = 10000 # number of timesteps +p = pyket.Params() +p.load_dict(d) +h = pyket.H(p) + +up_list = np.zeros(N) +down_list = np.zeros(N) +num_states = np.zeros(N) +t = np.zeros(N) +print(h.get_state()) +for i in range(N): + h.grow() + h.run_step(complex(0, -dt)) + up, down = h.get_spin_pop() + num_states[i] = h.get_size() + t[i] = dt * i + up_list[i], down_list[i] = up, down +fig, ax = plt.subplots( + 2, 1, figsize=(9, 4) +) # Create a figure containing a single axes. +ax[0].plot(t, up_list, label="up") +ax[0].set_title("up and down vs time") +ax[0].set_ylabel("Probability") +ax[0].plot(t, down_list, label="down") +ax[0].legend() +ax[1].plot(t, num_states, label="number of states") +ax[1].set_ylabel("Number of States") +ax[0].set_xlabel("Time (arbitrary)") +fig.suptitle(f"w0={w0},v0={v0},v1={v1},g0={g0},g1={g1}") +plt.savefig("plot_2_modes." + f"w0={w0}.v0={v0}.v1={v1}.g0={g0}.g1={g1}.png") diff --git a/pyket/recursion_tests.py b/pyket/recursion_tests.py deleted file mode 100644 index bf0e908..0000000 --- a/pyket/recursion_tests.py +++ /dev/null @@ -1,540 +0,0 @@ -import numpy as np -import argparse -import matplotlib.pyplot as plt -from enum import Enum,unique - - -def make_tuples_hermetian(tuple_list): - rows = [x[0] for x in tuple_list] - cols = [x[1] for x in tuple_list] - vals = [x[2] for x in tuple_list] - for r,c,v in tuple_list: - if r !=c: - rows.append(c) - cols.append(r) - vals.append(np.conjugate(v)) - return(np.array(rows,int),np.array(cols,int),np.array(vals,complex)) - -class pTests(Enum): - p100g = 1 - p2q = 2 - p100q = 3 - p1000 = 4 -if __name__ == '__main__': - parser = argparse.ArgumentParser(description='Test different cases for recursion (repeating behaviour)') - parser.add_argument('tests', metavar='arg', type=lambda val:pTests[val], nargs='+', - help='tests to choose from') - parser.add_argument('--bits', dest='bits',type=int, - help='number of bits to use', default=8) - - - args = parser.parse_args() - - - for bt in args.tests: - if(bt==pTests.p100g and args.bits==8): - - import pyket100_8 as pyket - w0 = 2 - v0 = 2 - g0 = .2 - tf = .1 - dt = .01 - k1 = pyket.StateKet() - k2 = pyket.StateKet() - k1.set_mode(99,1) - k1.set_spin(False) - k2.set_spin(True) - k1.set_amp(complex(.001,0)) - k2.set_amp(complex(1,0)) - init_state = pyket.StateVector() - init_state.push_back(k1) - init_state.push_back(k2) - d = { - "run_id":1420, - "output_directory":"~/code/scratch", - "initial_state":init_state, - "mode_energies":np.append(np.linspace(0,10,99), v0), - "mode_couplings": np.append(np.linspace(.0001,.001,99), g0), - "up_energy":-.5*w0, - "down_energy":.5*w0, - "energy_cutoff":.00001, - "dt":dt - } - N = int(tf/dt) - p = pyket.Params() - p.load_dict(d) - h = pyket.H(p) - up =[0 for x in range(N)] - num =[0 for x in range(N)] - t = [0 for x in range(N)] - for i in range(N): - h.grow() - h.run_step(complex(0,-dt)) - up[i],down = h.get_spin_pop() - num[i] = h.get_size() - t[i] = dt*i - fig,(ax1, ax2) = plt.subplots(2, 1) - # make a little extra space between the subplots - fig.subplots_adjust(hspace=0.5) - - ax1.plot(t, up) - ax1.set_xlim(0, t[-1]) - ax1.set_xlabel('Time') - ax1.set_ylabel('Spin Population') - ax1.grid(True) - - - ax2.plot(t, num) - ax2.set_xlim(0, t[-1]) - ax2.set_xlabel('Time') - ax2.set_ylabel('State_Size') - ax2.grid(True) - plt.show() - - if(bt==pTests.p2q and args.bits==8): - - import pyket2_8 as pyket - import qutip - fock_N1 = 40 - fock_N2 = 40 - - w0 = 3 - v0 = 4 - g0 = .2 - w1 = 7 - v1 = 6 - g1 = .007 - pinit_state = pyket.StateVector() - - ket = pyket.StateKet() - ket.set_spin(True) - ket.set_amp(complex('1')) - pinit_state.push_back(ket) - ket = pyket.StateKet() - ket.set_mode(0,1) - ket.set_amp(complex(0)) - pinit_state.push_back(ket) - ket = pyket.StateKet() - - d = { - "run_id":1420, - "output_directory":"~/code/scratch", - "initial_state":pinit_state, - "mode_energies":[v0,v1], - "mode_couplings":[g0,g1], - "up_energy":.5*w0, - "down_energy":-.5*w0, - "energy_cutoff":-1, - "dt":.01 - } - - p = pyket.Params() - p.load_dict(d) - h = pyket.H(p) - - #--------qutip code---------- - #basic operators - - s_eye = qutip.operators.identity(2) - n_eye1 = qutip.operators.identity(fock_N1) - n_eye2 = qutip.operators.identity(fock_N2) - - oz = qutip.tensor(qutip.operators.sigmaz(),n_eye1,n_eye2) - oplus = qutip.tensor(qutip.operators.sigmap(),n_eye1,n_eye2) - ox = qutip.tensor(qutip.operators.sigmax(),n_eye1,n_eye2) - a0 = qutip.tensor(s_eye,qutip.operators.destroy(fock_N1),n_eye2) - a1 = qutip.tensor(s_eye,n_eye1,qutip.operators.destroy(fock_N2)) - - - #make initial state - vac0 = qutip.states.basis(fock_N1) - vac1 = qutip.states.basis(fock_N2) - spin = qutip.states.basis(2) - qinit_state = qutip.tensor(spin,vac0,vac1).unit() - - - #hamiltonian - H = 0.5*w0*oz + v0*a0.dag()*a0 + v1*a1.dag()*a1 + g0*ox*(a0.dag()+a0) + g1*ox*(a1.dag()+a1) - - - tf1 = .003 - tf2 = 12 - dt = .0005 - NPTS = 100 - - N1 = int(tf1/dt) - N2 = int(tf2/dt) - def get_qtip_pop( qstate): - return(abs(qstate.unit().ptrace(0)[0][0][0]),abs(qstate.unit().ptrace(0)[0][0][1])) - - up1 =[0 for x in range(N1)] - up_qutip1 = [0 for x in range(N1)] - num1 =[0 for x in range(N1)] - num_qutip1 = [fock_N1*fock_N2*2 for x in range(N1)] - t1 = [0 for x in range(N1)] - for i in range(N1): - print(i) - h.grow() - h.run_step(complex(0,-dt)) - qinit_state = qinit_state +complex(0,-dt)*H*qinit_state - up1[i],d = h.get_spin_pop() - up_qutip1[i],d = get_qtip_pop(qinit_state) - num1[i] = h.get_size() - - t1[i] = dt*i - - up2 =[0 for x in range(N2)] - down =[0 for x in range(N2)] - up_qutip2 = [0 for x in range(N2)] - num2 =[0 for x in range(N2)] - up_vac_cav_qutip2= [fock_N1*fock_N2*2 for x in range(N2)] - - t2 = [0 for x in range(N2)] - h.set_zero_except_init() - from scipy.sparse import csc_matrix - from scipy.sparse.linalg import expm, expm_multiply - from scipy.linalg import norm - from scipy.sparse import diags - v = h.get_tuples() - n_size = h.get_size() - row1 = [x[0] for x in v] - col1 = [x[1] for x in v] - data1 = [x[2] for x in v] - for r,c,dat, in v: - if r != c: - row1.append(c) - col1.append(r) - data1.append(dat) - - istate = np.zeros((n_size,1)) - istate[1] = 1 - A=complex(0,-1)*csc_matrix((data1, (row1, col1)), shape=(n_size,n_size)) - timeSci = np.linspace(0,tf2,NPTS,endpoint=True) - res =expm_multiply(A,istate, start=0, stop=tf2, num=NPTS, endpoint=True) - spin_idxs = h.get_spin_idxs() - - v_up = np.zeros(n_size) - for i in range(n_size): - v_up[i] = spin_idxs[i][0] - - upListSci = np.zeros(len(res)) - up_vac_cav = np.zeros(len(res)) - for i in range(len(res)): - r = diags(v_up,format="csc").dot(res[i]) - val = norm(r)**2 - upListSci[i] = val - up_vac_cav[i] = norm(res[i][1])**2 - - print(upListSci) - qinit_state = qutip.tensor(spin,vac0,vac1).unit() - for i in range(N2): - print(i) - qinit_state = (qinit_state +complex(0,-dt)*H*qinit_state ).unit() - up_qutip2[i],d = get_qtip_pop(qinit_state) - up_vac_cav_qutip2[i] = abs(qinit_state.overlap(qutip.tensor(spin,vac0,vac1).unit()))**2 - t2[i] = dt*i - - up2 = upListSci - fig, ax = plt.subplots(2,2 ,figsize=(12,5) ) - # make a little extra space between the subplots - fig.subplots_adjust(hspace=0.5) - - ax[0][0].plot(t1, up1, label="my code") - ax[0][0].plot(t1, up_qutip1, label = "qutip") - ax[0][0].set_xlim(0, t1[-1]) - ax[0][0].set_xlabel('Time') - ax[0][0].set_ylabel('Spin Population during Growth') - ax[0][0].grid(True) - ax[0][0].legend() - - - ax[1][0].plot(t1, num1,label='my code') - ax[1][0].plot(t1, num_qutip1,label='qutip') - ax[1][0].set_xlim(0, t1[-1]) - ax[1][0].set_xlabel('Time') - ax[1][0].set_ylabel('State_Size') - ax[1][0].grid(True) - ax[1][0].legend() - - - ax[0][1].plot(timeSci, up2, label="my code population") - ax[0][1].plot(t2, up_qutip2, label = "qutip") - ax[0][1].set_xlabel('Time') - ax[0][1].set_ylabel('Spin Population during Growth') - ax[0][1].grid(True) - ax[0][1].legend() - - - ax[1][1].plot(timeSci, up_vac_cav,label='my code') - ax[1][1].plot(t2, up_vac_cav_qutip2,label='qutip') - ax[1][1].set_xlabel('Time') - ax[1][1].set_ylabel('Spin Up Vacuum Cavity Probability') - ax[1][1].grid(True) - ax[1][1].legend() - fig.savefig('/home/ethan/code/run_output/outp2q.png',dpi=300) - - - - - if(bt==pTests.p100q and args.bits==8): - from scipy.sparse import csc_matrix - from scipy.sparse.linalg import expm, expm_multiply - from scipy.linalg import norm - from scipy.sparse import diags - import pyket100_8 as pyket - from numpy.random import default_rng - rng = default_rng() - w0 = 3 - v0 = 3 - g0 = 2 - - m = np.linspace(0,5,100) - g = np.linspace(0,.01,100) - g[0] = g0; - m[0] = v0; - tf1 = .001 - tf2 = 1000 - dt = .001 - N2 = 10000 - N1 = int(tf1/dt) - - #grow initial space - psi_init = pyket.StateVector() - k = pyket.StateKet() - k.set_mode(0,1) - k.set_spin(False) - k.set_amp(complex(0,0)) - psi_init.push_back(k) - k = pyket.StateKet() - k.set_spin(True) - k.set_amp(complex(1,0)) - psi_init.push_back(k) - d = { - "initial_state":psi_init, - "mode_energies":m, - "mode_couplings":g, - "up_energy":.5*w0, - "down_energy":-.5*w0, - "energy_cutoff":.00001, - } - p = pyket.Params() - p.load_dict(d) - h = pyket.H(p) - - - t1 = np.zeros(N1) - t2 = np.zeros(N2) - pop1 = np.zeros((N1,2)) - num1 = np.zeros(N1) - - for i in range(N1): - print(i) - h.grow() - h.run_step(complex(0,-dt)) - pop1[i][0],pop1[i][1] = h.get_spin_pop() - num1[i] = h.get_size() - t1[i] = i*dt - h.set_zero_except_init() - rows,cols,data = make_tuples_hermetian(h.get_tuples()) - data = data*complex(0,-1) - h_size = h.get_size() - psi_init2 = np.zeros((h.get_size(),1),complex) - for ket in psi_init: - idx = h.get_ket_idx(ket) - assert(idx > -1) - psi_init2[idx] = ket.get_amp() - - - pop2 = np.zeros((N2,2)) - up_vacs = np.zeros(N2) - H_Sparse = csc_matrix((data,(rows,cols)),shape=(h_size,h_size)) - results = expm_multiply(H_Sparse,psi_init2, start=0,stop=tf2,num=N2,endpoint=True) - spin_idxs = h.get_spin_idxs() - up_mask = np.array([x[0] for x in spin_idxs]) - down_mask = np.array([x[1] for x in spin_idxs]) - ket = pyket.StateKet() - ket.set_spin(True) - up_vac_idx = h.get_ket_idx(ket) - if(up_vac_idx < 0): - print("Error!") - print(h.get_state()) - t2 = np.linspace(0,tf2,N2) - for i in range(len(results)): - r_up = diags(up_mask,format="csc").dot(results[i]) - r_down = diags(down_mask,format="csc").dot(results[i]) - pop2[i] = (norm(r_up)**2,norm(r_down)**2) - up_vacs[i] = abs(results[i][up_vac_idx])**2 - - - - - import matplotlib.ticker as ticker - fig, ax = plt.subplots(2,2 ,figsize=(12,5)) - fig.subplots_adjust(hspace=.5,wspace=.5) - - ax[0][0].plot(t1,pop1[:,0],marker="*",label="Spin Up") - ax[0][0].plot(t1,pop1[:,1],marker="1",label="Spin Down") - ax[0][0].set_xlabel('Time') - ax[0][0].set_ylabel('Spin Population During Growth') - ax[0][0].grid(True) - ax[0][0].set_xlim(left=0) - ax[0][0].xaxis.set_major_locator(ticker.LinearLocator(5)) - ax[0][0].xaxis.set_minor_locator(ticker.LinearLocator(20)) - ax[0][0].legend() - - ax[1][0].plot(t1,num1,marker="*") - ax[1][0].set_xlabel('Time') - ax[1][0].set_ylabel('State Space Size') - ax[1][0].set_xlim(left=0) - ax[1][0].xaxis.set_major_locator(ticker.LinearLocator(5)) - ax[1][0].xaxis.set_minor_locator(ticker.LinearLocator(20)) - ax[1][0].grid(True) - - ax[0][1].plot(t2,pop2[:,0],label='Spin Up') - ax[0][1].plot(t2,pop2[:,1],label='Spin Down') - ax[0][1].set_xlabel('Time') - ax[0][1].set_ylabel('Spin Population During Sparse Exp. Mat.') - ax[0][1].grid(True) - ax[0][1].legend() - - ax[1][1].plot(t2,up_vacs) - ax[1][1].set_ylabel('Initial Spin Up State Probability') - ax[1][1].set_xlabel('Time') - ax[1][1].grid(True) - plt.show() - - if(bt==pTests.p1000 and args.bits==8): - from scipy.sparse import csc_matrix - from scipy.sparse.linalg import expm, expm_multiply - from scipy.linalg import norm - from scipy.sparse import diags - import pyket1000_8 as pyket - from numpy.random import default_rng - - rng = default_rng() - w0 = 3 - v0 = 3 - g0 = .3 - - m = np.linspace(2,4,pyket.num_modes()) - g = np.linspace(0,.01,pyket.num_modes()) - g[0] = g0; - m[0] = v0; - tf1 = .11 - tf2 = 70 - dt = .001 - N2 = 1000 - N1 = int(tf1/dt) - - #grow initial space - psi_init = pyket.StateVector() - k = pyket.StateKet() - k.set_mode(0,1) - k.set_spin(False) - k.set_amp(complex(0,0)) - psi_init.push_back(k) - k = pyket.StateKet() - k.set_spin(True) - k.set_amp(complex(1,0)) - psi_init.push_back(k) - d = { - "initial_state":psi_init, - "mode_energies":m, - "mode_couplings":g, - "up_energy":.5*w0, - "down_energy":-.5*w0, - "energy_cutoff":.00001, - } - p = pyket.Params() - p.load_dict(d) - h = pyket.H(p) - - - t1 = np.zeros(N1) - t2 = np.zeros(N2) - pop1 = np.zeros((N1,2)) - num1 = np.zeros(N1) - - for i in range(N1): - print(i) - h.grow() - h.run_step(complex(0,-dt)) - pop1[i][0],pop1[i][1] = h.get_spin_pop() - num1[i] = h.get_size() - t1[i] = i*dt - h.set_zero_except_init() - rows,cols,data = make_tuples_hermetian(h.get_tuples()) - data = data*complex(0,-1) - h_size = h.get_size() - psi_init2 = np.zeros((h.get_size(),1),complex) - for ket in psi_init: - idx = h.get_ket_idx(ket) - assert(idx > -1) - psi_init2[idx] = ket.get_amp() - - - pop2 = np.zeros((N2,2)) - up_vacs = np.zeros((N2,2)) - H_Sparse = csc_matrix((data,(rows,cols)),shape=(h_size,h_size)) - results = expm_multiply(H_Sparse,psi_init2, start=0,stop=tf2,num=N2,endpoint=True) - spin_idxs = h.get_spin_idxs() - up_mask = np.array([x[0] for x in spin_idxs]) - down_mask = np.array([x[1] for x in spin_idxs]) - ket = pyket.StateKet() - ket.set_spin(True) - up_vac_idx = h.get_ket_idx(ket) - ket.set_spin(False) - ket.set_mode(0,1) - down_vac_idx = h.get_ket_idx(ket) - if(up_vac_idx < 0): - print("Error!") - print(h.get_state()) - t2 = np.linspace(0,tf2,N2) - for i in range(len(results)): - r_up = diags(up_mask,format="csc").dot(results[i]) - r_down = diags(down_mask,format="csc").dot(results[i]) - pop2[i] = (norm(r_up)**2,norm(r_down)**2) - up_vacs[i] = abs(results[i][up_vac_idx])**2,abs(results[i][down_vac_idx])**2 - - - - - import matplotlib.ticker as ticker - fig, ax = plt.subplots(2,2 ,figsize=(12,5)) - fig.subplots_adjust(hspace=.5,wspace=.5) - - ax[0][0].plot(t1,pop1[:,0],marker="*",label="Spin Up") - ax[0][0].plot(t1,pop1[:,1],marker="1",label="Spin Down") - ax[0][0].set_xlabel('Time') - ax[0][0].set_ylabel('Spin Population During Growth') - ax[0][0].grid(True) - ax[0][0].set_xlim(left=0) - ax[0][0].xaxis.set_major_locator(ticker.LinearLocator(5)) - ax[0][0].xaxis.set_minor_locator(ticker.LinearLocator(20)) - ax[0][0].legend() - - ax[1][0].plot(t1,num1,marker="*") - ax[1][0].set_xlabel('Time') - ax[1][0].set_ylabel('State Space Size') - ax[1][0].set_xlim(left=0) - ax[1][0].xaxis.set_major_locator(ticker.LinearLocator(5)) - ax[1][0].xaxis.set_minor_locator(ticker.LinearLocator(20)) - ax[1][0].grid(True) - - ax[0][1].plot(t2,pop2[:,0],label='Spin Up') - ax[0][1].plot(t2,pop2[:,1],label='Spin Down') - ax[0][1].set_xlabel('Time') - ax[0][1].set_ylabel('Spin Population During Sparse Exp. Mat.') - ax[0][1].grid(True) - ax[0][1].legend() - - ax[1][1].plot(t2,up_vacs) - ax[1][1].set_ylabel('Initial Spin Up State Probability') - ax[1][1].set_xlabel('Time') - ax[1][1].grid(True) - plt.show() - - - #calculate emission rate - diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..6f7026c --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,19 @@ +[build-system] +requires = ["setuptools>=42", "wheel", "pybind11~=2.6.1"] +build-backend = "setuptools.build_meta" + +[project] +name = "pyket" +description = "Adaptive code for pyket" +dynamic = ["version"] +dependencies = [ + "test", + "pytest" +] + + +[tool.setuptools] +zip-safe=false + + + diff --git a/scripts/example_benzen_in_cavity.ipynb b/scripts/example_benzen_in_cavity.ipynb new file mode 100644 index 0000000..07a0094 --- /dev/null +++ b/scripts/example_benzen_in_cavity.ipynb @@ -0,0 +1,94 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Adaptive Calculation of Benzene in a Nano-Tip Cavity\n", + "\n", + "The purpose of this example is to show how an adaptive code can be used to find the *'active space'* of a light-matter system. The adaptive code builds *'marries*' two methods/calculations:\n", + "\n", + "1. A quasi-normal mode (QNM) calculation of a nano-cavity using **MEEP**. MEEP is an open-source ectrodynamics python/C++ package. [link](https://meep.readthedocs.io/en/latest/)\n", + "2. An excited state calculation of a molecule using **PySCF**. The excited state calculation yeilds the electronic energy levels and electric transition dipole moments. PySCF is an open-source electronic structure software for python. [link](https://pyscf.org/)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Quasi-Normal Mode Calculation\n", + "\n", + "We hope to calculate the quasi-normal modes (QNM) of a nano-tip cavity. A quasi-normal mode is a 'resonant' solution of Maxwell's wave equation for an open cavity. This boils down to an eigenvalue problem where solutions(QNMs) are pairs of:\n", + "- `Eigen spatial electric fields` - A vector field $\\vec E_m(\\vec r)$ for QNM $m$ \n", + "- `Eigen complex frequencies` - A complex frequency: $\\omega + i\\Gamma$\n", + " - $\\omega$ represents oscillation frequency\n", + " - $\\Gamma$ represents decay rate \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Excited State Calculation of Molecule\n", + "\n", + "Given some molecule (n this case benzene) we hope to calculate \n", + "`Excited states` - A list of $N_e$ Excited States, which are used to find\n", + "- `Excited state energies` - A list of $N_e$ excited state energies \n", + "- `Electric Dipole Transition Moments` - A list of transition moments such that the dipole moment $\\vec \\mu$ is $$ \\vec \\mu = \\sum_{i j} d_{i,j}\\langle i | \\vec r | j \\rangle \\; \\bigl [ | i \\rangle | \\langle j |\\bigr]$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Use adaptive code to 'marry' the excited state and quasi-normal mode calculations\n", + "\n", + "In this section we\n", + "- Prepare parameters\n", + " - Place the molecule (benzene) in the cavity at $\\vec R_{mol}$\n", + " - Calculate the couplings $g_{i j m}$ which is given by the dipole interaction $$H_{int} = -\\vec \\mu \\cdot \\vec E = - \\sum_{i j} d_{i,j}\\langle i | \\vec r | j \\rangle \\; \\bigl [ | i \\rangle | \\langle j |\\bigr] \\cdot \\sum_{m\\in \\text{QNMs}} \\vec E_m(\\vec R_{mol})$$ \n", + " - For each QNM m, represent loss by distributing couplings $g_{i j m}$ on a discrete Lorentzian of width $\\Gamma_m$ and area $g_{i j m}$\n", + "- Find 'active space' with adaptive code\n", + " - Choose some initial excited state $\\psi$\n", + " - Find the transitions $\\mathcal T_q (\\psi)$ where $|\\mathcal T_q(\\psi)| > \\epsilon$ \n", + " - Expand active space so $\\psi$ can transition into the $\\mathcal T_q$'s\n", + " - Update $\\psi$ to the next time step $\\psi(t) \\rightarrow \\psi(t+\\Delta t)$ \n", + " - repeat until you have hit an iteration limit\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.10.5 ('.venv': venv)", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.10.5" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "f1a5ab9681ce3ca299cb303913c7832bf9396572c6e1d32e3bab8d7a0afaa5f4" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..b7e9559 --- /dev/null +++ b/setup.py @@ -0,0 +1,44 @@ +#taken heavily from: +# https://github.com/pybind/python_example +# https://pybind11.readthedocs.io/en/stable/compiling.html#building-with-setuptools +from pathlib import Path +import sys + +from pybind11 import get_cmake_dir +# Available at setup time due to pyproject.toml +from pybind11.setup_helpers import Pybind11Extension, build_ext +from setuptools import setup + +__version__ = "0.0.1" +NUM_MODES = 8000 +NUM_BITS = 8 +include_dirs = Path('./ext/boost/') +# The main interface is through Pybind11Extension. +# * You can add cxx_std=11/14/17, and then build_ext can be removed. +# * You can set include_pybind11=false to add the include directory yourself, +# say from a submodule. +# +# Note: +# Sort input source files if you glob sources to ensure bit-for-bit +# reproducible builds (https://github.com/pybind/python_example/pull/53) + +ext_modules = [ + Pybind11Extension( + name="pyket", + sources=sorted(Path('.').glob('**/*.cpp')), + # Example: passing in the version to the compiled code + define_macros = [ + ('VERSION_INFO', __version__), + ('NUM_MODES', NUM_MODES), + ('NUM_BITS',NUM_BITS), + ('MODULE_NAME',__name__) + ], + ), + include_dirs = [ + Path('./extern/rapidjson/include'), + Path('./extern/boost'), + Path('./extern/eigen'), + ] +] + +setup(ext_modules=ext_modules,cmdclass={"build_ext": build_ext}) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt deleted file mode 100644 index 37623a3..0000000 --- a/src/CMakeLists.txt +++ /dev/null @@ -1 +0,0 @@ -add_subdirectory(${CMAKE_CURRENT_LIST_DIR}/hamiltonian) diff --git a/src/hamiltonian.cpp b/src/hamiltonian.cpp new file mode 100644 index 0000000..e9878ea --- /dev/null +++ b/src/hamiltonian.cpp @@ -0,0 +1,45 @@ +#include "pyket_config.hpp" +#include +#include +#include +#include "hamiltonian.hpp" +#include "input_tools.hpp" + +PYBIND11_MAKE_OPAQUE(State_Ket); +PYBIND11_MAKE_OPAQUE(std::vector>); +PYBIND11_MAKE_OPAQUE(hamiltonian); + +namespace py = pybind11; + + +using state_ket = State_Ket; +using state_vector = std::vector; + +void init_hamiltonian(py::module_ &m){ + + py::class_(m,"H") + .def(py::init([](const std::string &json_arg){ + return std::unique_ptr(new hamiltonian(json_arg)); + })) + .def(py::init([](param_vals &p){ + return std::unique_ptr(new hamiltonian(p)); + })) + .def("get_ket_idx",&hamiltonian::get_ket_idx,"returns -1 if it cannot get idx.") + .def("store_matrix",&hamiltonian::store_matrix,py::call_guard()) + .def("store_vector",&hamiltonian::store_vector,py::call_guard()) + .def("evolve_state",&hamiltonian::evolve_state,py::call_guard()) + .def("evolve_step",&hamiltonian::evolve_step,py::call_guard()) + .def("get_tuples",&hamiltonian::get_tuples) + .def("get_spin_idxs",&hamiltonian::get_spin_idxs) + .def("reset",&hamiltonian::reset) + .def("reset_with_state",&hamiltonian::reset_with_state) + .def("run_step", &hamiltonian::run_step) + .def("grow",&hamiltonian::run_grow) + .def("get_state", &hamiltonian::get_state) + .def("get_size", &hamiltonian::get_psi_size) + .def("get_spin_pop",&hamiltonian::get_spin_pop) + .def("get_emitter_cavity_prob",&hamiltonian::get_emitter_cavity_prob) + .def("get_mode_pop",&hamiltonian::get_modeLbl_quanta_pop) + .def("set_zero_except_init",&hamiltonian::set_zero_except_init) + .def("set_epsilon",&hamiltonian::set_epsilon); +} \ No newline at end of file diff --git a/src/hamiltonian/CMakeLists.txt b/src/hamiltonian/CMakeLists.txt deleted file mode 100644 index dee4ecf..0000000 --- a/src/hamiltonian/CMakeLists.txt +++ /dev/null @@ -1,17 +0,0 @@ -target_sources(HamiltonianFuncs PRIVATE - append_connections.cpp - merge_states.cpp - output.cpp - runtime.cpp - core.cpp - grow_configuration_space.cpp - par_runtime.cpp - setup.cpp - io_tools/input_tools.hpp - io_tools/input_tools.cpp - PUBLIC - configuration.hpp - ) - - -target_include_directories(HamiltonianFuncs PUBLIC ${CMAKE_CURRENT_LIST_DIR} ${CMAKE_CURRENT_LIST_DIR}/io_tools ) diff --git a/src/hamiltonian/io_tools/.#output_tools.hpp b/src/hamiltonian/io_tools/.#output_tools.hpp deleted file mode 120000 index e1a3616..0000000 --- a/src/hamiltonian/io_tools/.#output_tools.hpp +++ /dev/null @@ -1 +0,0 @@ -ethan@ethan-x570.22715:1609009941 \ No newline at end of file diff --git a/src/hamiltonian/io_tools/CMakeLists.txt b/src/hamiltonian/io_tools/CMakeLists.txt deleted file mode 100644 index f7f5471..0000000 --- a/src/hamiltonian/io_tools/CMakeLists.txt +++ /dev/null @@ -1,11 +0,0 @@ -add_library(IO_Tools - input_tools.cpp - output_tools.cpp - ) -target_include_directories(IO_Tools - PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} - ) - -# install rules -install(TARGETS IO_Tools DESTINATION lib) -install(FILES input_tools.hpp output_tools.hpp DESTINATION include) diff --git a/src/hamiltonian/par_runtime.cpp b/src/hamiltonian/par_runtime.cpp deleted file mode 100644 index 72c4929..0000000 --- a/src/hamiltonian/par_runtime.cpp +++ /dev/null @@ -1,18 +0,0 @@ -#include "hamiltonian.hpp" - -#include -#include -#include - - - - - - -void hamiltonian::par_test_one(){ - -} - -void hamiltonian::par_test_two(){ - -} diff --git a/src/hamiltonian/par_runtime.cpp.old b/src/hamiltonian/par_runtime.cpp.old deleted file mode 100644 index 6159640..0000000 --- a/src/hamiltonian/par_runtime.cpp.old +++ /dev/null @@ -1,114 +0,0 @@ -#include "hamiltonian.hpp" -#include "configuration_char.hpp" - - - -#include -#include -#include -/*#include -#include -#include -#include -#include -#include -#include -#include */ -#include - -#define VEXCL_BACKEND_COMPUTE -#include -#include -#include - -#include -#include - -namespace compute = boost::compute; - - - -struct basic_config{ - - bool spin; - uint8_t rep[63]; - uint idx; - - basic_config():spin(false),rep{1},idx(0){} - bool operator<(const basic_config c)const{ - - if(c.spin > spin){ - return true; - }else if( spin == c.spin){ - for(int i = 0; i < 63; i++){ - if(c.rep[i] > rep[i]){ - return true; - } - return false; - } - } - - return false; - } - -}; - - BOOST_COMPUTE_ADAPT_STRUCT(basic_config,basic_config, (spin,rep,idx)) - -void hamiltonian::par_test_one(){ - - - - BOOST_COMPUTE_FUNCTION(bool, cmp_config, (basic_config b, basic_config a), - { - - if(a.spin > b.spin){ - return true; - }else if(b.spin == a.spin){ - for(int i = 0; i < 63; i++){ - if(a.rep[i] > b.rep[i]){ - return true; - } - return false; - } - }else{ - return false; - } - }); - - - - - compute::device device = compute::system::default_device(); - compute::context context(device); - compute::command_queue queue(context,device); - - - compute::vector state_vec(context); - for(int i = 0; i < 100000; i++){ - basic_config k; - k.rep[i%64] = i; - state_vec.push_back(k,queue); - } - - compute::sort(state_vec.begin(),state_vec.end(),cmp_config,queue); - - - - -} - -void hamiltonian::par_test_two(){ - using namespace boost::numeric::ublas; - symmetric_matrix ml (3, 3); - for (unsigned i = 0; i < ml.size1 (); ++ i) - for (unsigned j = 0; j <= i; ++ j) - ml (i, j) = 3 * i + j; - std::cout << ml << std::endl; - symmetric_matrix mu (3, 3); - for (unsigned i = 0; i < mu.size1 (); ++ i) - for (unsigned j = i; j < mu.size2 (); ++ j) - mu (i, j) = 3 * i + j; - std::cout << mu << std::endl; - -} diff --git a/src/hamiltonian/py_generate_config_class.py b/src/hamiltonian/py_generate_config_class.py deleted file mode 100644 index 8dc20ab..0000000 --- a/src/hamiltonian/py_generate_config_class.py +++ /dev/null @@ -1,6 +0,0 @@ - -num_modes = 10 -result = "struct config_manual" -for i in range(num_modes): - result - diff --git a/src/hamiltonian/py_str_gen.py b/src/hamiltonian/py_str_gen.py deleted file mode 100644 index 91969c0..0000000 --- a/src/hamiltonian/py_str_gen.py +++ /dev/null @@ -1,17 +0,0 @@ -from string import Template -if __name__=="__main__": - out = "" - s = Template('%|${mode}$$-d|:%|${level}$$-d|%|${space}t|') - for i in range(10): - d = dict(mode=str(2*i+1),level=str(2*i+2),space=str(8*(i+1))) - out = out + s.substitute(d) - out = out + "\\n" - print(out) - - out = "" - - s = Template('fmtr % mode_level_array[${idx}].first;\nfmtr % mode_level_array[${idx}].second;\n') - for i in range(10): - d = dict(idx ="i*modes_per_line +"+ str(i)) - out = out + s.substitute(d) - print(out) diff --git a/src/hamiltonian/runtime.cpp b/src/hamiltonian/runtime.cpp index 7734909..394eef4 100644 --- a/src/hamiltonian/runtime.cpp +++ b/src/hamiltonian/runtime.cpp @@ -37,48 +37,6 @@ void hamiltonian::run_grow(){ } -void hamiltonian::get_decay_diags(int delta_e){ - - vector> decay_diags(psi_lbl.size()); - - //for each vector in psi, get the connected states - for(int i=0; i < int(psi_amp.size()); i++){ - double E0 = std::real(matrix_diag(psi_amp[i])); - int idx = psi_amp[i].idx; - for(int j = 0; j < NUM_MODES; j++){ - ket_pair kp = get_connected_states(psi_amp[i],j); - int level = psi_amp[i].get_mode(j); - if( (kp.raised.i != state_ket::null_i) && (magnitude* g[level+1][j]> params.energy_cutoff) ){ - //Check to see if energy is within delta_e - double E_k = std::real(matrix_diag(kp.raised )); - double V_ik = g[level+1][j]; - if(std::abs(E_k-E0) <= delta_e){ - if(binary_search_state(kp.raised, psi_lbl) < 0){ - //state is not in current space, so we add it as decay term - decay_diags[idx] += complex(0,-(V_ik*V_ik)/delta_e); - } - } - } - if( (kp.lowered.i != state_ket::null_i) && (magnitude*g[level][j] > params.energy_cutoff)){ - //Check to see if energy is within delta_e - double E_k = std::real(matrix_diag(kp.lowered )); - double V_ik = g[level][j]; - if(std::abs(E_k-E0) <= delta_e){ - //state is within energy threshold. - if(binary_search_state(kp.lowered,psi_lbl) < 0){ - //State is not in current space, add it as complex energy - //Given by fermi's golden rule - decay_diags[idx] += complex(0,(V_ik*V_ik)/delta_e); - } - } - } - - } - } - - return(decay_diags); -} - void hamiltonian::set_zero_except_init(){ using namespace std; diff --git a/src/main.cpp b/src/main.cpp index c89d4fe..922b524 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -15,9 +15,8 @@ PYBIND11_MAKE_OPAQUE(std::vector); PYBIND11_MAKE_OPAQUE(std::vector); -PYBIND11_MAKE_OPAQUE(State_Ket); -PYBIND11_MAKE_OPAQUE(std::vector>); -PYBIND11_MAKE_OPAQUE(hamiltonian); + + using namespace std; @@ -27,7 +26,23 @@ using state_ket = State_Ket; using state_vector = std::vector; +/* include definitions for module initialization +defined in main_*****.cpp ` +*/ +void init_hamiltonian(py::module_ &); +void init_state_ket(py::module_ &); +void init_state_vector(py::module_ &); +void init_params(py::module_ &); + + PYBIND11_MODULE(PYKET,m){ + + /* initialize separately defined module components */ + init_state_ket(m); + init_state_vector(m); + init_hamiltonian(m); + init_params(m); + m.doc() = "adaptive spin evolution ported to python"; m.def("num_modes",[](){return int( NUM_MODES);},"Number of allowed modes in a state"); m.def("num_bits", [](){return int(NUM_BITS);},"Number of bits that can store a level"); @@ -41,81 +56,7 @@ PYBIND11_MODULE(PYKET,m){ [](auto &it1,const auto &it2){return norm(it1.amp)>norm(it2.amp);}); }); - py::class_(m,"StateVector") - .def(py::init<>()) - .def("__iter__",[](state_vector &v){ - return py::make_iterator(v.begin(),v.end()); - },py::keep_alive<0,1>()) - .def("__len__",[](const state_vector &v){return v.size();}) - .def("clear", &state_vector::clear) - .def("push_back", [](state_vector &v, const state_ket &k)mutable{v.push_back(k);} ) - .def("__repr__",[](state_vector &v){ - string result = "[\n "; - for(const auto &d: v){ - result = result + "{" + d.str() + "}\n"; - } - result += "]"; - return(py::str(result)); - }) - .def("__getitem__",[](const state_vector &v,uint i) - {return v[i];}); - - py::class_(m,"StateKet") - .def(py::init<>()) - .def("get_mode",&state_ket::get_mode) - .def("set_mode",&state_ket::set_mode) - .def("increment_mode",&state_ket::increment_mode) - .def("decrement_mode",&state_ket::decrement_mode) - .def("get_spin",[](const state_ket &k){return k.spin;}) - .def("set_spin",[](state_ket &k,bool spin)mutable{ k.spin = spin;}) - .def("__repr__",&state_ket::str) - .def("__eq__",&state_ket::operator==) - .def("__lt__",&state_ket::operator<) - .def("set_amp",[](state_ket &k,int i)mutable{k.idx = i;}) - .def("get_idx",[](const state_ket &k){return k.idx;}) - .def("get_amp",[](const state_ket &k){return k.amp;}) - .def("set_amp",[](state_ket &k,std::complex a)mutable{k.amp = a;}); - py::class_(m,"H") - .def(py::init([](const std::string &json_arg){ - return std::unique_ptr(new hamiltonian(json_arg)); - })) - .def(py::init([](param_vals &p){ - return std::unique_ptr(new hamiltonian(p)); - })) - .def("get_ket_idx",&hamiltonian::get_ket_idx,"returns -1 if it cannot get idx.") - .def("store_matrix",&hamiltonian::store_matrix,py::call_guard()) - .def("store_vector",&hamiltonian::store_vector,py::call_guard()) - .def("evolve_state",&hamiltonian::evolve_state,py::call_guard()) - .def("evolve_step",&hamiltonian::evolve_step,py::call_guard()) - .def("get_tuples",&hamiltonian::get_tuples) - .def("get_spin_idxs",&hamiltonian::get_spin_idxs) - .def("reset",&hamiltonian::reset) - .def("reset_with_state",&hamiltonian::reset_with_state) - .def("run_step", &hamiltonian::run_step) - .def("grow",&hamiltonian::run_grow) - .def("get_state", &hamiltonian::get_state) - .def("get_size", &hamiltonian::get_psi_size) - .def("get_spin_pop",&hamiltonian::get_spin_pop) - .def("get_emitter_cavity_prob",&hamiltonian::get_emitter_cavity_prob) - .def("get_mode_pop",&hamiltonian::get_modeLbl_quanta_pop) - .def("set_zero_except_init",&hamiltonian::set_zero_except_init) - .def("set_epsilon",&hamiltonian::set_epsilon) - .def("get_matrix",[](hamiltonian &h){ - // Redirect cout. - streambuf* oldCoutStreamBuf = cout.rdbuf(); - ostringstream strCout; - cout.rdbuf( strCout.rdbuf() ); - - // This goes to the string stream. - cout << "(TODO) not implimented" ; - - // Restore old cout. - cout.rdbuf( oldCoutStreamBuf ); - - - return strCout.str(); - }); - + py::class_>(m, "IntVector") .def(py::init<>()) .def("clear", &std::vector::clear) @@ -145,68 +86,4 @@ PYBIND11_MODULE(PYKET,m){ }); - py::class_(m, "Params") - .def(py::init<>()) - .def("load_dict",[](param_vals &p,py::dict dict){ - for(auto item:dict){ - string key = string(py::str(item.first)); - if(key.compare(string{"run_id"})==0){ - p.run_id = item.second.cast(); - }else if(key.compare(string{"output_directory"})==0){ - p.output_directory = item.second.cast(); - }else if(key.compare(string{"energy_cutoff"})==0){ - p.energy_cutoff = item.second.cast(); - }else if(key.compare(string{"energy_spectral_density"})==0){ - p.energy_spectral_density = item.second.cast(); - }else if(key.compare(string{"mode_couplings"})==0){ - p.mode_couplings.clear(); - for(auto val:item.second){ - p.mode_couplings.push_back(val.cast()); - } - assert(p.mode_couplings.size() == NUM_MODES); - }else if(key.compare(string{"mode_energies"})==0){ - p.mode_energies.clear(); - - for(auto val:item.second){ - p.mode_energies.push_back(val.cast()); - } - assert(p.mode_energies.size() == NUM_MODES); - }else if(key.compare(string{"up_energy"})==0){ - p.up_energy = item.second.cast(); - }else if(key.compare(string{"down_energy"})==0){ - p.down_energy = item.second.cast(); - }else if(key.compare(string{"t0"})==0){ - p.t0 = item.second.cast(); - }else if(key.compare(string{"tf"})==0){ - p.tf = item.second.cast(); - }else if(key.compare(string{"dt"})==0){ - p.dt = item.second.cast(); - }else if(key.compare(string{"N"})==0){ - p.N = item.second.cast(); - }else if(key.compare(string{"initial_state"})==0){ - state_vector *v = item.second.cast(); - p.initial_state.resize(v->size()); - copy(v->begin(),v->end(),p.initial_state.begin()); - }else{ - py::print(py::str("Key not able to be assigned: " + key)); - } - - } - - }) - .def_readwrite("run_id",¶m_vals::run_id) - .def_readwrite("output_directory",¶m_vals::output_directory) - .def_readwrite("energy_cutoff",¶m_vals::energy_cutoff) - .def_readwrite("energy_spectral_density",¶m_vals::energy_spectral_density) - .def_readwrite("mode_couplings",¶m_vals::mode_couplings) - .def_readwrite("mode_energies",¶m_vals::mode_energies) - .def_readwrite("up_energy",¶m_vals::up_energy) - .def_readwrite("down_energy",¶m_vals::down_energy) - .def_readwrite("t0",¶m_vals::t0) - .def_readwrite("tf",¶m_vals::tf) - .def_readwrite("dt",¶m_vals::dt) - .def_readwrite("N",¶m_vals::N) - .def_readwrite("initial_state",¶m_vals::initial_state); - - } diff --git a/src/main_hamiltonian.cpp b/src/main_hamiltonian.cpp new file mode 100644 index 0000000..e9878ea --- /dev/null +++ b/src/main_hamiltonian.cpp @@ -0,0 +1,45 @@ +#include "pyket_config.hpp" +#include +#include +#include +#include "hamiltonian.hpp" +#include "input_tools.hpp" + +PYBIND11_MAKE_OPAQUE(State_Ket); +PYBIND11_MAKE_OPAQUE(std::vector>); +PYBIND11_MAKE_OPAQUE(hamiltonian); + +namespace py = pybind11; + + +using state_ket = State_Ket; +using state_vector = std::vector; + +void init_hamiltonian(py::module_ &m){ + + py::class_(m,"H") + .def(py::init([](const std::string &json_arg){ + return std::unique_ptr(new hamiltonian(json_arg)); + })) + .def(py::init([](param_vals &p){ + return std::unique_ptr(new hamiltonian(p)); + })) + .def("get_ket_idx",&hamiltonian::get_ket_idx,"returns -1 if it cannot get idx.") + .def("store_matrix",&hamiltonian::store_matrix,py::call_guard()) + .def("store_vector",&hamiltonian::store_vector,py::call_guard()) + .def("evolve_state",&hamiltonian::evolve_state,py::call_guard()) + .def("evolve_step",&hamiltonian::evolve_step,py::call_guard()) + .def("get_tuples",&hamiltonian::get_tuples) + .def("get_spin_idxs",&hamiltonian::get_spin_idxs) + .def("reset",&hamiltonian::reset) + .def("reset_with_state",&hamiltonian::reset_with_state) + .def("run_step", &hamiltonian::run_step) + .def("grow",&hamiltonian::run_grow) + .def("get_state", &hamiltonian::get_state) + .def("get_size", &hamiltonian::get_psi_size) + .def("get_spin_pop",&hamiltonian::get_spin_pop) + .def("get_emitter_cavity_prob",&hamiltonian::get_emitter_cavity_prob) + .def("get_mode_pop",&hamiltonian::get_modeLbl_quanta_pop) + .def("set_zero_except_init",&hamiltonian::set_zero_except_init) + .def("set_epsilon",&hamiltonian::set_epsilon); +} \ No newline at end of file diff --git a/src/main_params.cpp b/src/main_params.cpp new file mode 100644 index 0000000..146b9c4 --- /dev/null +++ b/src/main_params.cpp @@ -0,0 +1,89 @@ +#include "pyket_config.hpp" +#include +#include +#include +#include +#include +#include"configuration.hpp" +#include "hamiltonian.hpp" +#include "input_tools.hpp" + + +PYBIND11_MAKE_OPAQUE(State_Ket); +PYBIND11_MAKE_OPAQUE(std::vector>); +PYBIND11_MAKE_OPAQUE(hamiltonian); + +namespace py = pybind11; +using namespace std; + +using state_ket = State_Ket; +using state_vector = std::vector; + +void init_params(py::module_ &m){ + + + py::class_(m, "Params") + .def(py::init<>()) + .def("load_dict",[](param_vals &p,py::dict dict){ + for(auto item:dict){ + string key = string(py::str(item.first)); + if(key.compare(string{"run_id"})==0){ + p.run_id = item.second.cast(); + }else if(key.compare(string{"output_directory"})==0){ + p.output_directory = item.second.cast(); + }else if(key.compare(string{"energy_cutoff"})==0){ + p.energy_cutoff = item.second.cast(); + }else if(key.compare(string{"energy_spectral_density"})==0){ + p.energy_spectral_density = item.second.cast(); + }else if(key.compare(string{"mode_couplings"})==0){ + p.mode_couplings.clear(); + for(auto val:item.second){ + p.mode_couplings.push_back(val.cast()); + } + assert(p.mode_couplings.size() == NUM_MODES); + }else if(key.compare(string{"mode_energies"})==0){ + p.mode_energies.clear(); + + for(auto val:item.second){ + p.mode_energies.push_back(val.cast()); + } + assert(p.mode_energies.size() == NUM_MODES); + }else if(key.compare(string{"up_energy"})==0){ + p.up_energy = item.second.cast(); + }else if(key.compare(string{"down_energy"})==0){ + p.down_energy = item.second.cast(); + }else if(key.compare(string{"t0"})==0){ + p.t0 = item.second.cast(); + }else if(key.compare(string{"tf"})==0){ + p.tf = item.second.cast(); + }else if(key.compare(string{"dt"})==0){ + p.dt = item.second.cast(); + }else if(key.compare(string{"N"})==0){ + p.N = item.second.cast(); + }else if(key.compare(string{"initial_state"})==0){ + state_vector *v = item.second.cast(); + p.initial_state.resize(v->size()); + copy(v->begin(),v->end(),p.initial_state.begin()); + }else{ + py::print(py::str("Key not able to be assigned: " + key)); + } + + } + + }) + .def_readwrite("run_id",¶m_vals::run_id) + .def_readwrite("output_directory",¶m_vals::output_directory) + .def_readwrite("energy_cutoff",¶m_vals::energy_cutoff) + .def_readwrite("energy_spectral_density",¶m_vals::energy_spectral_density) + .def_readwrite("mode_couplings",¶m_vals::mode_couplings) + .def_readwrite("mode_energies",¶m_vals::mode_energies) + .def_readwrite("up_energy",¶m_vals::up_energy) + .def_readwrite("down_energy",¶m_vals::down_energy) + .def_readwrite("t0",¶m_vals::t0) + .def_readwrite("tf",¶m_vals::tf) + .def_readwrite("dt",¶m_vals::dt) + .def_readwrite("N",¶m_vals::N) + .def_readwrite("initial_state",¶m_vals::initial_state); + + +} \ No newline at end of file diff --git a/src/main_state_ket.cpp b/src/main_state_ket.cpp new file mode 100644 index 0000000..19550da --- /dev/null +++ b/src/main_state_ket.cpp @@ -0,0 +1,37 @@ +#include "pyket_config.hpp" +#include +#include +#include +#include"configuration.hpp" + +PYBIND11_MAKE_OPAQUE(State_Ket); +PYBIND11_MAKE_OPAQUE(std::vector>); + + +namespace py = pybind11; + +using state_ket = State_Ket; + +namespace py = pybind11; + +void init_state_ket(py::module_ &m){ + + py::class_(m,"StateKet") + .def(py::init<>()) + .def("get_mode",&state_ket::get_mode) + .def("set_mode",&state_ket::set_mode) + .def("increment_mode",&state_ket::increment_mode) + .def("decrement_mode",&state_ket::decrement_mode) + .def("get_spin",[](const state_ket &k){return k.spin;}) + .def("set_spin",[](state_ket &k,bool spin)mutable{ k.spin = spin;}) + .def("__repr__",&state_ket::str) + .def("__eq__",&state_ket::operator==) + .def("__lt__",&state_ket::operator<) + .def("set_amp",[](state_ket &k,int i)mutable{k.idx = i;}) + .def("get_idx",[](const state_ket &k){return k.idx;}) + .def("get_amp",[](const state_ket &k){return k.amp;}) + .def("set_amp",[](state_ket &k,std::complex a)mutable{k.amp = a;}); + + + +} \ No newline at end of file diff --git a/src/main_state_vector.cpp b/src/main_state_vector.cpp new file mode 100644 index 0000000..bc8044d --- /dev/null +++ b/src/main_state_vector.cpp @@ -0,0 +1,37 @@ +#include "pyket_config.hpp" +#include +#include +#include +#include +#include +#include"configuration.hpp" + +PYBIND11_MAKE_OPAQUE(State_Ket); +PYBIND11_MAKE_OPAQUE(std::vector>); + +namespace py = pybind11; + +using state_ket = State_Ket; +using state_vector = std::vector; + +void init_state_vector(py::module_ &m){ + py::class_(m,"StateVector") + .def(py::init<>()) + .def("__iter__",[](state_vector &v){ + return py::make_iterator(v.begin(),v.end()); + },py::keep_alive<0,1>()) + .def("__len__",[](const state_vector &v){return v.size();}) + .def("clear", &state_vector::clear) + .def("push_back", [](state_vector &v, const state_ket &k)mutable{v.push_back(k);} ) + .def("__repr__",[](state_vector &v){ + std::string result = "[\n "; + for(const auto &d: v){ + result = result + "{" + d.str() + "}\n"; + } + result += "]"; + return(py::str(result)); + }) + .def("__getitem__",[](const state_vector &v,uint i) + {return v[i];}); +} + \ No newline at end of file diff --git a/src/params.cpp b/src/params.cpp new file mode 100644 index 0000000..146b9c4 --- /dev/null +++ b/src/params.cpp @@ -0,0 +1,89 @@ +#include "pyket_config.hpp" +#include +#include +#include +#include +#include +#include"configuration.hpp" +#include "hamiltonian.hpp" +#include "input_tools.hpp" + + +PYBIND11_MAKE_OPAQUE(State_Ket); +PYBIND11_MAKE_OPAQUE(std::vector>); +PYBIND11_MAKE_OPAQUE(hamiltonian); + +namespace py = pybind11; +using namespace std; + +using state_ket = State_Ket; +using state_vector = std::vector; + +void init_params(py::module_ &m){ + + + py::class_(m, "Params") + .def(py::init<>()) + .def("load_dict",[](param_vals &p,py::dict dict){ + for(auto item:dict){ + string key = string(py::str(item.first)); + if(key.compare(string{"run_id"})==0){ + p.run_id = item.second.cast(); + }else if(key.compare(string{"output_directory"})==0){ + p.output_directory = item.second.cast(); + }else if(key.compare(string{"energy_cutoff"})==0){ + p.energy_cutoff = item.second.cast(); + }else if(key.compare(string{"energy_spectral_density"})==0){ + p.energy_spectral_density = item.second.cast(); + }else if(key.compare(string{"mode_couplings"})==0){ + p.mode_couplings.clear(); + for(auto val:item.second){ + p.mode_couplings.push_back(val.cast()); + } + assert(p.mode_couplings.size() == NUM_MODES); + }else if(key.compare(string{"mode_energies"})==0){ + p.mode_energies.clear(); + + for(auto val:item.second){ + p.mode_energies.push_back(val.cast()); + } + assert(p.mode_energies.size() == NUM_MODES); + }else if(key.compare(string{"up_energy"})==0){ + p.up_energy = item.second.cast(); + }else if(key.compare(string{"down_energy"})==0){ + p.down_energy = item.second.cast(); + }else if(key.compare(string{"t0"})==0){ + p.t0 = item.second.cast(); + }else if(key.compare(string{"tf"})==0){ + p.tf = item.second.cast(); + }else if(key.compare(string{"dt"})==0){ + p.dt = item.second.cast(); + }else if(key.compare(string{"N"})==0){ + p.N = item.second.cast(); + }else if(key.compare(string{"initial_state"})==0){ + state_vector *v = item.second.cast(); + p.initial_state.resize(v->size()); + copy(v->begin(),v->end(),p.initial_state.begin()); + }else{ + py::print(py::str("Key not able to be assigned: " + key)); + } + + } + + }) + .def_readwrite("run_id",¶m_vals::run_id) + .def_readwrite("output_directory",¶m_vals::output_directory) + .def_readwrite("energy_cutoff",¶m_vals::energy_cutoff) + .def_readwrite("energy_spectral_density",¶m_vals::energy_spectral_density) + .def_readwrite("mode_couplings",¶m_vals::mode_couplings) + .def_readwrite("mode_energies",¶m_vals::mode_energies) + .def_readwrite("up_energy",¶m_vals::up_energy) + .def_readwrite("down_energy",¶m_vals::down_energy) + .def_readwrite("t0",¶m_vals::t0) + .def_readwrite("tf",¶m_vals::tf) + .def_readwrite("dt",¶m_vals::dt) + .def_readwrite("N",¶m_vals::N) + .def_readwrite("initial_state",¶m_vals::initial_state); + + +} \ No newline at end of file diff --git a/src/pyket_config.hpp b/src/pyket_config.hpp deleted file mode 100644 index fe810ff..0000000 --- a/src/pyket_config.hpp +++ /dev/null @@ -1 +0,0 @@ -#define PYKET pyket1_8 diff --git a/src/pyket_config.hpp.in b/src/pyket_config.hpp.in deleted file mode 100644 index a996709..0000000 --- a/src/pyket_config.hpp.in +++ /dev/null @@ -1 +0,0 @@ -#define PYKET @pyket@ diff --git a/src/state_ket.cpp b/src/state_ket.cpp new file mode 100644 index 0000000..19550da --- /dev/null +++ b/src/state_ket.cpp @@ -0,0 +1,37 @@ +#include "pyket_config.hpp" +#include +#include +#include +#include"configuration.hpp" + +PYBIND11_MAKE_OPAQUE(State_Ket); +PYBIND11_MAKE_OPAQUE(std::vector>); + + +namespace py = pybind11; + +using state_ket = State_Ket; + +namespace py = pybind11; + +void init_state_ket(py::module_ &m){ + + py::class_(m,"StateKet") + .def(py::init<>()) + .def("get_mode",&state_ket::get_mode) + .def("set_mode",&state_ket::set_mode) + .def("increment_mode",&state_ket::increment_mode) + .def("decrement_mode",&state_ket::decrement_mode) + .def("get_spin",[](const state_ket &k){return k.spin;}) + .def("set_spin",[](state_ket &k,bool spin)mutable{ k.spin = spin;}) + .def("__repr__",&state_ket::str) + .def("__eq__",&state_ket::operator==) + .def("__lt__",&state_ket::operator<) + .def("set_amp",[](state_ket &k,int i)mutable{k.idx = i;}) + .def("get_idx",[](const state_ket &k){return k.idx;}) + .def("get_amp",[](const state_ket &k){return k.amp;}) + .def("set_amp",[](state_ket &k,std::complex a)mutable{k.amp = a;}); + + + +} \ No newline at end of file diff --git a/src/state_vector.cpp b/src/state_vector.cpp new file mode 100644 index 0000000..bc8044d --- /dev/null +++ b/src/state_vector.cpp @@ -0,0 +1,37 @@ +#include "pyket_config.hpp" +#include +#include +#include +#include +#include +#include"configuration.hpp" + +PYBIND11_MAKE_OPAQUE(State_Ket); +PYBIND11_MAKE_OPAQUE(std::vector>); + +namespace py = pybind11; + +using state_ket = State_Ket; +using state_vector = std::vector; + +void init_state_vector(py::module_ &m){ + py::class_(m,"StateVector") + .def(py::init<>()) + .def("__iter__",[](state_vector &v){ + return py::make_iterator(v.begin(),v.end()); + },py::keep_alive<0,1>()) + .def("__len__",[](const state_vector &v){return v.size();}) + .def("clear", &state_vector::clear) + .def("push_back", [](state_vector &v, const state_ket &k)mutable{v.push_back(k);} ) + .def("__repr__",[](state_vector &v){ + std::string result = "[\n "; + for(const auto &d: v){ + result = result + "{" + d.str() + "}\n"; + } + result += "]"; + return(py::str(result)); + }) + .def("__getitem__",[](const state_vector &v,uint i) + {return v[i];}); +} + \ No newline at end of file diff --git a/travis/build_wheels.sh b/travis/build_wheels.sh new file mode 100644 index 0000000..e6788e2 --- /dev/null +++ b/travis/build_wheels.sh @@ -0,0 +1,38 @@ +#!/bin/bash +set -e -u -x + +function repair_wheel { + wheel="$1" + if ! auditwheel show "$wheel"; then + echo "Skipping non-platform wheel $wheel" + else + auditwheel repair "$wheel" --plat "$PLAT" -w /io/wheelhouse/ + fi +} + + +# Install headers required for the extension +curl -sSL https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz | tar -xvz +curl -sSL https://github.com/Tencent/rapidjson/archive/refs/tags/v1.1.0.tar.gz | tar -xvz +curl -sSL https://boostorg.jfrog.io/artifactory/main/release/1.80.0/source/boost_1_80_0.tar.gz | tar -xvz +mkdir include_dir +cp -r boost_*/boost include_dir +cp -r eigen-*/Eigen include_dir +cp -r rapidjson-*/include/rapidjson include_dir + +# Compile wheels +for PYBIN in /opt/python/*/bin; do + "${PYBIN}/pip" install -r /io/dev-requirements.txt + "${PYBIN}/pip" wheel /io/ --no-deps -w wheelhouse/ +done + +# Bundle external shared libraries into the wheels +for whl in wheelhouse/*.whl; do + repair_wheel "$whl" +done + +# Install packages and test +for PYBIN in /opt/python/*/bin/; do + "${PYBIN}/pip" install pyket -I ./include_dir --no-index --f /io/wheelhouse + (cd "$HOME"; "${PYBIN}/nosetests" pyket) +done \ No newline at end of file diff --git a/unit_tests/main_test.cpp b/unit_tests/main_test.cpp deleted file mode 100644 index a14d778..0000000 --- a/unit_tests/main_test.cpp +++ /dev/null @@ -1,282 +0,0 @@ -#define BOOST_TEST_MODULE dataset_example -#include -#include -#include -#include <../configuration.hpp> -#include <../hamiltonian.hpp> -//input utilities -#include <../input_tools.hpp> -#include -#include -#include -#include - -namespace data = boost::unit_test::data; -namespace utf = boost::unit_test; - - -const int modes = 500; -const int bits = 4; -typedef partial_config p_conf; -BOOST_TEST_DECORATOR(* utf::label("partial_state_comparison")) -BOOST_DATA_TEST_CASE(partial_config1 /*suite name*/, - data::xrange(500)*data::xrange(500)^data::random(1,15), - m1,m2,level) -{ - - p_conf c1; - p_conf c2; - c1.set_mode(m1,level); - c2.set_mode(m2,level); - BOOST_TEST((c1 < c2)==(m1 > m2)); -} - -typedef State_Ket, modes,bits> Ket; -BOOST_TEST_DECORATOR(* utf::label("ket_comparison")) -BOOST_DATA_TEST_CASE(State_Ket_1 /*suite name*/, - data::xrange(500)*data::xrange(500)^data::random(1,15)^data::random(0,1)^data::random(0,1), - m1,m2,level,spin1,spin2) -{ - - Ket k1; - Ket k2; - k1.spin = spin1; - k1.set_mode(m1,level); - - k2.spin = spin2; - k2.set_mode(m2,level); - //(k1 < k2)==(m1>m2)&&(spin1 < spin2)); - bool kcheck = k1m2)&&(spin1 < spin2))); -} - - - - -struct F { -public: - - const int num_modes=100; - const int bits = 4; - typedef vector::iterator vIt; - typedef State_Ket,100,4> State_Ket_T; - typedef vector State_Vector; - typedef hamiltonian hamiltonian_T; - - - F(){ - BOOST_TEST_MESSAGE( "setup fixture" ); - //read in params - int argc = 1; - char argv[] = {'p','r','\0'}; - char * argptr[] = {argv}; - bool exit = get_params(params,argc, argptr); - BOOST_TEST_REQUIRE(!exit); - - Spin_Params spin_params; - spin_params.energy = 1; - - double energy_cutoff = .1; - - State_Vector initial_state; - State_Ket_T ket; - ket.amp=1; - initial_state.push_back(ket); - } - - ~F() - { BOOST_TEST_MESSAGE( "teardown fixture" ); } - void setup(){ - - } - void teardown(){ - - } - - param_vals params; - hamiltonian_T h(); - - -}; - -BOOST_TEST_DECORATOR(*utf::tolerance(0.00001) *utf::label("operators")) -BOOST_DATA_TEST_CASE_F(F, vector_normalize, - data::xrange(300)^data::random()^data::random()^data::random(1,500), - idx,re,im,mul) -{ - - State_Vector v(300); - for( int i = 1; i < 300; i++){ - v[i].amp += 1; - v[i].amp *=mul; - } - BOOST_TEST_MESSAGE( v[idx]); - v[idx].amp =0 ; - - F::hamiltonian_T::normalize_state(v); - double norm = F::hamiltonian_T::state_magnitude(v); - - BOOST_TEST( norm == 1 ); - -} - - - - -struct FSMALL { -public: - - const int num_modes=2; - const int bits = 4; - typedef vector::iterator vIt; - typedef State_Ket,2,4> State_Ket_T; - typedef vector State_Vector_T; - typedef hamiltonian hamiltonian_T; - - param_vals params; - State_Vector_T delta; - State_Ket_T ket; - hamiltonian_T h; - - FSMALL(){ - BOOST_TEST_MESSAGE( "setup fixture" ); - //read in params - int argc = 1; - char argv[] = {'p','r','\0'}; - char * argptr[] = {argv}; - bool exit = get_params(params,argc, argptr); - BOOST_TEST_REQUIRE(!exit); - - Spin_Params spin_params; - spin_params.energy = 1; - - double energy_cutoff = .00000000001; - - State_Vector_T initial_state(1); - initial_state[0].amp=1; - - h.setup(params.mode_energies.begin(),params.mode_couplings.begin(), - spin_params,initial_state,energy_cutoff); - - - h.h_dipole(delta,ket,-1,1); - h.simple_run(1); - - - } - - ~FSMALL() - { BOOST_TEST_MESSAGE( "teardown fixture" ); } - void setup(){ - - } - void teardown(){ - - } - - - -}; - - -BOOST_TEST_DECORATOR(*utf::tolerance(0.00001) *utf::label("operators")) -BOOST_DATA_TEST_CASE_F(FSMALL,single_mode, - data::xrange(1), - r) -{ - - State_Vector_T my_v = delta; - - BOOST_TEST(my_v[r].get_mode(0)==1); - -} - - -const int mode_one_expected[] = {0,1,0}; -const int mode_two_expected[] = {0,0,1}; - -BOOST_TEST_DECORATOR(*utf::tolerance(0.00001) *utf::label("auxillary")) -BOOST_DATA_TEST_CASE_F(FSMALL, single_run, - data::xrange(3)^data::make(mode_one_expected)^data::make(mode_two_expected), - idx,m1,m2) -{ - - State_Vector_T v = h.get_psi_lbl(); - BOOST_TEST_MESSAGE(v[idx]); - BOOST_REQUIRE(v.size() == 3); - BOOST_TEST(v[idx].get_mode(0) == m1); - BOOST_TEST(v[idx].get_mode(1) == m2); - -} - - - -struct FBIG { -public: - - const int num_modes=100; - const int bits = 4; - typedef vector::iterator vIt; - typedef State_Ket,2,4> State_Ket_T; - typedef vector State_Vector_T; - typedef hamiltonian hamiltonian_T; - - param_vals params; - State_Vector_T delta; - State_Ket_T ket; - hamiltonian_T h; - - FBIG(){ - BOOST_TEST_MESSAGE( "setup fixture" ); - //read in params - int argc = 1; - char argv[] = {'p','r','\0'}; - char * argptr[] = {argv}; - bool exit = get_params(params,argc, argptr); - BOOST_TEST_REQUIRE(!exit); - - Spin_Params spin_params; - spin_params.energy = params.atom_levels[1]; - - double energy_cutoff = .00000000001; - - State_Vector_T initial_state(1); - initial_state[0].amp=1; - - h.setup(params.mode_energies.begin(),params.mode_couplings.begin(), - spin_params,initial_state,energy_cutoff); - - - h.h_dipole(delta,ket,-1,1); - h.simple_run(1/20); - - - } - - ~FBIG() - { BOOST_TEST_MESSAGE( "teardown fixture" ); } - void setup(){ - - } - void teardown(){ - - } - - - -}; - - -BOOST_TEST_DECORATOR(*utf::tolerance(0.00001) *utf::label("leakage")) -BOOST_DATA_TEST_CASE_F(FBIG, single_run_big, - data::xrange(3), - idx) -{ - - State_Vector_T v = h.get_psi_lbl(); - BOOST_TEST_MESSAGE(v[0]); - BOOST_REQUIRE(v.size() > 1); - -} - diff --git a/unit_tests/makefile b/unit_tests/makefile deleted file mode 100644 index 389060d..0000000 --- a/unit_tests/makefile +++ /dev/null @@ -1,15 +0,0 @@ -objects = input_tools.o main_test.o -flags = -g -std=c++17 -#-DBOOST_TEST_TOOLS_UNDER_DEBUGGER -boost_test: ${objects} - g++ ${flags} -o boost_test ${objects} -lboost_unit_test_framework -lboost_program_options - -main_test.o: ../configuration.hpp ../input_tools.hpp ../hamiltonian.hpp main_test.cpp - g++ ${flags} -c main_test.cpp - -input_tools.o: ../input_tools.hpp ../input_tools.cpp - g++ ${flags} -c ../input_tools.cpp -clean: - rm *.o boost_test - -