Commit 3fdbebc5 authored by mgassner's avatar mgassner
Browse files

added computation of MI

parent 799e42bf
...@@ -2,8 +2,6 @@ import lattice_symmetries as ls ...@@ -2,8 +2,6 @@ import lattice_symmetries as ls
import numpy as np import numpy as np
import yaml import yaml
from entropy_new import idxs_rearrangement
# Spin operators # Spin operators
sigma_x = np.array([[0., 1.], [1., 0.]]) sigma_x = np.array([[0., 1.], [1., 0.]])
sigma_y = np.array([[0., -1.0j], [1.0j, 0.]]) sigma_y = np.array([[0., -1.0j], [1.0j, 0.]])
...@@ -84,13 +82,13 @@ class spin_model: ...@@ -84,13 +82,13 @@ class spin_model:
def create_hamiltonian(self): def create_hamiltonian(self):
self.hamiltonian = ls.Operator(self.basis, [ls.Interaction(self.__matrix, self.edges)]) self.hamiltonian = ls.Operator(self.basis, [ls.Interaction(self.__matrix, self.edges)])
"""
def update_basis_and_hamiltonian(self, rearrangement): def update_basis_and_hamiltonian(self, rearrangement):
# TODO: find out if I need to update edges as well!!! # TODO: find out if I need to update edges as well!!!
new_states = idxs_rearrangement(self.basis.states, rearrangement) new_states = idxs_rearrangement(self.basis.states, rearrangement)
self.basis.build(representatives=new_states) self.basis.build(representatives=new_states)
self.hamiltonian = ls.Operator(self.basis, [ls.Interaction(self.__matrix, self.edges)]) self.hamiltonian = ls.Operator(self.basis, [ls.Interaction(self.__matrix, self.edges)])
"""
def compute_ew_and_ev(self, print_gs_energy=True): def compute_ew_and_ev(self, print_gs_energy=True):
self.eigenvalues, self.eigenstates = ls.diagonalize(self.hamiltonian, k=1) self.eigenvalues, self.eigenstates = ls.diagonalize(self.hamiltonian, k=1)
if print_gs_energy: if print_gs_energy:
......
...@@ -37,7 +37,8 @@ def sector_dm(trace_states, sub_states, amplitudes, spin_inversion): ...@@ -37,7 +37,8 @@ def sector_dm(trace_states, sub_states, amplitudes, spin_inversion):
distinct_sub_states = list(set(sub_states)) distinct_sub_states = list(set(sub_states))
distinct_trace_states = list(set(trace_states)) distinct_trace_states = list(set(trace_states))
n = len(distinct_sub_states) n = len(distinct_sub_states)
print(trace_states, sub_states) # print('Trace states:', trace_states)
# print('Sub_states:', sub_states)
matrix = np.zeros((n,n)) matrix = np.zeros((n,n))
for state in distinct_trace_states: for state in distinct_trace_states:
...@@ -54,31 +55,41 @@ def sector_dm(trace_states, sub_states, amplitudes, spin_inversion): ...@@ -54,31 +55,41 @@ def sector_dm(trace_states, sub_states, amplitudes, spin_inversion):
matrix += matrix_add matrix += matrix_add
return matrix return matrix
def reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes, spin_inversion=False): def reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes, spin_inversion=None, first_trace_spin=None):
""" """
:sub_dim: amount of spins in the subsystem (uint) :sub_dim: amount of spins in the subsystem (uint)
:dim: amount of spins of total system (uint) :dim: amount of spins of total system (uint)
:hamming: amount of upspins of total system(uint or None) :hamming: amount of upspins of total system(uint or None)
:amplitudes: amplitudes of input state (1d-array of floats) :amplitudes: amplitudes of input state (1d-array of floats)
:states_of_amplitudes: state which belongs to amplitude (1d-array of uints) :states_of_amplitudes: state which belongs to amplitude (1d-array of uints)
:first_spin: first spin that's traced out; start counting at 0 (rightmost spin)!
:return: List of sector density matrices :return: List of sector density matrices
""" """
if spin_inversion is not None: if spin_inversion is not None:
dim -= 1 dim -= 1
sub_states = states_of_amplitudes >> (dim - sub_dim) if first_trace_spin is None:
trace_states = states_of_amplitudes & (2**(dim-sub_dim)-1) sub_states = states_of_amplitudes >> (dim - sub_dim)
trace_states = states_of_amplitudes & (2**(dim-sub_dim) - 1)
else:
right_sub_states = states_of_amplitudes & (2**first_trace_spin - 1)
left_sub_states = states_of_amplitudes >> (dim - (sub_dim - first_trace_spin))
sub_states = right_sub_states + (left_sub_states << first_trace_spin)
trace_states = states_of_amplitudes & (2**(dim - sub_dim + first_trace_spin) - 2**first_trace_spin )
trace_states >>= first_trace_spin
hamming_of_sub_states = np.array([_hamming_weight(n) for n in sub_states]) hamming_of_sub_states = np.array([_hamming_weight(n) for n in sub_states])
#print(sub_states, trace_states, hamming_of_sub_states, amplitudes) # print(sub_states, trace_states, hamming_of_sub_states, amplitudes)
rho = [] rho = []
if hamming is None: if hamming is None:
rho.append(sector_dm(trace_states, sub_states, amplitudes, spin_inversion)) rho.append(sector_dm(trace_states, sub_states, amplitudes, spin_inversion))
else: else:
for sub_hamming in range(max(0, hamming - (dim - sub_dim)), min(hamming, sub_dim) + 1): for sub_hamming in range(max(0, hamming - (dim - sub_dim)), min(hamming, sub_dim) + 1):
print('Sub hamming is', sub_hamming) # print('Sub hamming is', sub_hamming)
indices = np.where(hamming_of_sub_states == sub_hamming)[0].tolist() indices = np.where(hamming_of_sub_states == sub_hamming)[0].tolist()
sector_sub_states = [sub_states[i] for i in indices] sector_sub_states = [sub_states[i] for i in indices]
sector_trace_states = [trace_states[i] for i in indices] sector_trace_states = [trace_states[i] for i in indices]
#print(sector_trace_states) # print(sector_trace_states)
sector_amplitudes = [amplitudes[i] for i in indices] sector_amplitudes = [amplitudes[i] for i in indices]
rho.append(sector_dm(sector_trace_states, sector_sub_states, sector_amplitudes, spin_inversion)) rho.append(sector_dm(sector_trace_states, sector_sub_states, sector_amplitudes, spin_inversion))
return rho return rho
...@@ -105,27 +116,23 @@ def compute_entropy(rho): ...@@ -105,27 +116,23 @@ def compute_entropy(rho):
entropy -= lambda_log_lambda(spectrum).sum() entropy -= lambda_log_lambda(spectrum).sum()
return entropy return entropy
def compute_MI(sub_dim, dim, hamming, amplitudes, states_of_amplitudes, spin_inversion=None, first_trace_spin=None):
def idxs_rearrangement(old_states, rearrangement): # so far left and right sub-dimension is the same and equal to sub_dim
""" assert (first_trace_spin + sub_dim) <= dim and first_trace_spin >= sub_dim
Rearranges the old states according to the spin number list rearrangement first_trace_spin_left = first_trace_spin + (dim - 2*sub_dim)
:old_states: index of old states #S_AB:
:rearrangement: list of new spin positions print('--------------A and B---------------')
:return: index of rearranged spins rhos_AB = reduced_dm(sub_dim*2, dim, hamming, amplitudes, states_of_amplitudes,
""" spin_inversion=spin_inversion, first_trace_spin=first_trace_spin)
new_states = [] S_AB = compute_entropy(rhos_AB)
batch_size = 1000000 #S_B (right)
number_spins = len(rearrangement) print('-----------------B------------------')
rhos_B = reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes,
# Rearranges the columns in the old spin matrix and computes the index of the rearranged spins spin_inversion=spin_inversion, first_trace_spin=first_trace_spin)
for i in range(len(old_states) // batch_size + 1): S_B = compute_entropy(rhos_B)
print(i * 1. / (len(old_states) // batch_size + 1), "not zero if len(old_states) is large", flush=True) #S_A (left)
new_states.append( print('-----------------A------------------')
spin_to_index( rhos_A = reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes,
index_to_spin( spin_inversion=spin_inversion, first_trace_spin=first_trace_spin_left)
old_states[i * batch_size:np.min([(i + 1) * batch_size, len(old_states)])], number_spins S_A = compute_entropy(rhos_A)
)[:, rearrangement].astype(bool), number_spins return [S_A, S_B, S_AB]
) \ No newline at end of file
)
return np.concatenate(new_states)
\ No newline at end of file
...@@ -2,13 +2,17 @@ from heisenberg import test_xxz, test_area_law ...@@ -2,13 +2,17 @@ from heisenberg import test_xxz, test_area_law
from tfim import test_tfim, test_area_law_tfim from tfim import test_tfim, test_area_law_tfim
from package_fails import test_symmetries_xxx from package_fails import test_symmetries_xxx
from ls_example_from_github import github_test from ls_example_from_github import github_test
from mutual_information import test_mutual_information
if __name__ == "__main__": if __name__ == "__main__":
#test_xxz() #test_xxz()
#test_area_law() #test_area_law()
#test_area_law_tfim() #test_area_law_tfim()
#test_tfim() #test_tfim()
test_symmetries_xxx() #test_symmetries_xxx('tfim')
test_mutual_information()
"""print('------------------------------------------------') """print('------------------------------------------------')
print('---With symmetries and spin inversion---') print('---With symmetries and spin inversion---')
github_test() github_test()
......
import numpy as np
import yaml
import time
import matplotlib.pyplot as plt
from create_model import spin_model
from entropy_new import compute_MI, compute_entropy, reduced_dm
def test_mutual_information():
model_name = 'xxz'
number_spins = 10
periodic = False
spin_inversion = None
param = 1.0
if param < -1.0:
hamming_weight = None
else:
hamming_weight = number_spins // 2
model = spin_model(model_name=model_name, number_spins=number_spins, periodic=periodic,
param=param, hamming_weight=hamming_weight ,spin_inversion=spin_inversion)
model.compute_ew_and_ev()
print('EIGENSTATE at h/J = ', param, 'is: ', model.eigenstates[:,0])
sub_dim = 4
first_trace_spin = 4
number_spins = model.number_spins()
basis_states = model.basis.states
gs = model.eigenstates[:,0]
print('Number Spins: ', model.basis.number_spins)
print('States', model.basis.states)
print('Sub Dimension is ', sub_dim)
entropies = compute_MI(sub_dim, number_spins, hamming_weight, gs, basis_states, spin_inversion=None, first_trace_spin=first_trace_spin)
MI = entropies[0] + entropies[1] - entropies[2]
print(MI)
print(entropies)
print('--------------------Test-----------------------')
rhos4 = reduced_dm(4, number_spins, hamming_weight, gs, basis_states)
rhos8 = reduced_dm(4, number_spins, hamming_weight, gs, basis_states)
print(compute_entropy(rhos4), compute_entropy(rhos8))
\ No newline at end of file
...@@ -6,24 +6,23 @@ import matplotlib.pyplot as plt ...@@ -6,24 +6,23 @@ import matplotlib.pyplot as plt
from create_model import spin_model from create_model import spin_model
from entropy_new import _hamming_weight, reduced_dm, compute_entropy from entropy_new import _hamming_weight, reduced_dm, compute_entropy
def test_symmetries_xxx(): def test_symmetries_xxx(model_name='xxz'):
model_name = 'xxz'
number_spins = [8,10,12] number_spins = [8,10,12]
periodic = True periodic = True
param = 1.0 param = 1.0
print('Test if gs-energy stays the same if we use symmetries for the ls package:') print('Test if gs-energy stays the same if we use symmetries for the ls package:')
for length in number_spins: for length in number_spins:
hamming_weight = length // 2 hamming_weight = length // 2 if model_name=='xxz' else None
print('\n' , '------- Number spins is: ' + str(length) + '-----------' ) print('\n' , '------- Number spins is: ' + str(length) + '-----------' )
print('---With symmetries and spin inversion---') print('---With symmetries and spin inversion---')
model = spin_model(model_name=model_name, number_spins=length, periodic=periodic, model = spin_model(model_name=model_name, number_spins=length, periodic=periodic,
param=param, hamming_weight=hamming_weight, use_symmetries=True, spin_inversion=-1) param=param, hamming_weight=hamming_weight, use_symmetries=True, spin_inversion=-1 if model_name=='xxz' else 1)
model.compute_ew_and_ev() model.compute_ew_and_ev()
#print('EIGENSTATE at Delta/J = ', param, 'is: ', model.eigenstates[:,0]) #print('EIGENSTATE at Delta/J = ', param, 'is: ', model.eigenstates[:,0])
del model del model
print('---Without symmetries, but spin inversion---') print('---Without symmetries, but spin inversion---')
model = spin_model(model_name=model_name, number_spins=length, periodic=periodic, model = spin_model(model_name=model_name, number_spins=length, periodic=periodic,
param=param, hamming_weight=hamming_weight, spin_inversion=-1) param=param, hamming_weight=hamming_weight, spin_inversion=-1 if model_name=='xxz' else 1)
model.compute_ew_and_ev() model.compute_ew_and_ev()
#print('EIGENSTATE at Delta/J = ', param, 'is: ', model.eigenstates[:,0]) #print('EIGENSTATE at Delta/J = ', param, 'is: ', model.eigenstates[:,0])
del model del model
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment