entropy_new.py 5.57 KB
Newer Older
mgassner's avatar
mgassner committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import numpy as np
import sys
import math
from numpy.core.arrayprint import format_float_positional

from numpy.core.fromnumeric import trace

from joblib import delayed, Parallel
import numpy as np
import scipy.special
import pickle

import numba
from numba import uint64, float32, float64
import sys
import os

import lattice_symmetries
from lattice_symmetries import *

@numba.jit("uint64(uint64)", nogil=True, nopython=True)
def _hamming_weight(n: int):
    """The Hamming weight counts the non zero-values in a bit string. 
    In the case of spin-chains this is the amount of up-spins.
    Example: 10 -> 1010 -> hamming weight = 2
    See https://stackoverflow.com/a/9830282
    """
    n = (n & 0x5555555555555555) + ((n & 0xAAAAAAAAAAAAAAAA) >> 1)
    n = (n & 0x3333333333333333) + ((n & 0xCCCCCCCCCCCCCCCC) >> 2)
    n = (n & 0x0F0F0F0F0F0F0F0F) + ((n & 0xF0F0F0F0F0F0F0F0) >> 4)
    n = (n & 0x00FF00FF00FF00FF) + ((n & 0xFF00FF00FF00FF00) >> 8)
    n = (n & 0x0000FFFF0000FFFF) + ((n & 0xFFFF0000FFFF0000) >> 16)
    n = (n & 0x00000000FFFFFFFF) + ((n & 0xFFFFFFFF00000000) >> 32)
    return n

mgassner's avatar
mgassner committed
36
def sector_dm(trace_states, sub_states, amplitudes, spin_inversion):
mgassner's avatar
mgassner committed
37
38
39
    distinct_sub_states = list(set(sub_states))
    distinct_trace_states = list(set(trace_states))
    n = len(distinct_sub_states)
mgassner's avatar
mgassner committed
40
41
    # print('Trace states:', trace_states)
    # print('Sub_states:', sub_states)
mgassner's avatar
mgassner committed
42
43
44
45
46
47
48
49
50
    matrix = np.zeros((n,n))
    for state in distinct_trace_states:
        
        single_trace_result = [(amplitudes[i], distinct_sub_states.index(sub_states[i])) 
                                for i, x in enumerate(trace_states) if x == state]
        
        psi = np.zeros(n)
        for x in single_trace_result:  
            psi[x[1]] += x[0]
mgassner's avatar
mgassner committed
51
52
53
54
55
        matrix_add = np.outer(psi,psi)
        if spin_inversion is not None:
            matrix += 0.5 * (matrix_add + matrix_add[::-1,::-1])
        else:
            matrix += matrix_add
mgassner's avatar
mgassner committed
56
57
    return matrix

mgassner's avatar
mgassner committed
58
def reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes, spin_inversion=None, first_trace_spin=None):
mgassner's avatar
mgassner committed
59
    """
mgassner's avatar
mgassner committed
60
61
62
63
64
    :sub_dim: amount of spins in the subsystem (uint)
    :dim: amount of spins of total system (uint)
    :hamming: amount of upspins of total system(uint or None)
    :amplitudes: amplitudes of input state (1d-array of floats)
    :states_of_amplitudes: state which belongs to amplitude (1d-array of uints)
mgassner's avatar
mgassner committed
65
    :first_spin: first spin that's traced out; start counting at 0 (rightmost spin)!
mgassner's avatar
mgassner committed
66
    :return: List of sector density matrices
mgassner's avatar
mgassner committed
67
    """
mgassner's avatar
mgassner committed
68
69
    if spin_inversion is not None: 
        dim -= 1
mgassner's avatar
mgassner committed
70
71
72
73
74
75
76
77
78
79
    if first_trace_spin is None:
        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
    
mgassner's avatar
mgassner committed
80
    hamming_of_sub_states = np.array([_hamming_weight(n) for n in sub_states])
mgassner's avatar
mgassner committed
81
82
    # print(sub_states, trace_states, hamming_of_sub_states, amplitudes)

mgassner's avatar
mgassner committed
83
    rho = []
mgassner's avatar
mgassner committed
84
    if hamming is None:
mgassner's avatar
mgassner committed
85
        rho.append(sector_dm(trace_states, sub_states, amplitudes, spin_inversion))
mgassner's avatar
mgassner committed
86
    else:
mgassner's avatar
mgassner committed
87
        for sub_hamming in range(max(0, hamming - (dim - sub_dim)), min(hamming, sub_dim) + 1):
mgassner's avatar
mgassner committed
88
            # print('Sub hamming is', sub_hamming)
mgassner's avatar
mgassner committed
89
90
91
            indices = np.where(hamming_of_sub_states == sub_hamming)[0].tolist()
            sector_sub_states = [sub_states[i] for i in indices]
            sector_trace_states = [trace_states[i] for i in indices]
mgassner's avatar
mgassner committed
92
            # print(sector_trace_states)
mgassner's avatar
mgassner committed
93
94
            sector_amplitudes = [amplitudes[i] for i in indices]
            rho.append(sector_dm(sector_trace_states, sector_sub_states, sector_amplitudes, spin_inversion))
mgassner's avatar
mgassner committed
95
96
97
98
99
    return rho

    
def lambda_log_lambda(x):
    """
mgassner's avatar
mgassner committed
100
    Von Neumann Entropy (log with basis 2)
mgassner's avatar
mgassner committed
101
    :x: array of eigenvalues of a sector density matrix
mgassner's avatar
mgassner committed
102
    :return: negative Bipartite Entanglement entropy of a sector density_matrix
mgassner's avatar
mgassner committed
103
    """
mgassner's avatar
mgassner committed
104
    y = np.where(x > 1e-7, np.log(x) / np.log(2), 1e-17)
mgassner's avatar
mgassner committed
105
106
107
108
109
110
    y *= x
    return y

def compute_entropy(rho):
    """
    :rho: list of sector density matrices
mgassner's avatar
mgassner committed
111
    :return: Bipartite Entanglement Entropy (von Neumann)
mgassner's avatar
mgassner committed
112
113
114
115
116
117
118
    """
    entropy = 0
    for i in range(len(rho)):
            spectrum, _ = np.linalg.eigh(rho[i])
            entropy -= lambda_log_lambda(spectrum).sum()
    return entropy

mgassner's avatar
mgassner committed
119
120
121
def compute_MI(sub_dim, dim, hamming, amplitudes, states_of_amplitudes, spin_inversion=None, first_trace_spin=None):
    # 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
mgassner's avatar
mgassner committed
122
    first_trace_spin_left = first_trace_spin - sub_dim
mgassner's avatar
mgassner committed
123
124
125
126
127
128
129
130
131
132
133
134
    #S_AB:
    rhos_AB = reduced_dm(sub_dim*2, dim, hamming, amplitudes, states_of_amplitudes, 
                            spin_inversion=spin_inversion, first_trace_spin=first_trace_spin)
    S_AB = compute_entropy(rhos_AB)
    #S_B (right)
    rhos_B = reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes, 
                            spin_inversion=spin_inversion, first_trace_spin=first_trace_spin)
    S_B = compute_entropy(rhos_B)
    #S_A (left)
    rhos_A = reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes, 
                            spin_inversion=spin_inversion, first_trace_spin=first_trace_spin_left)
    S_A = compute_entropy(rhos_A)
mgassner's avatar
mgassner committed
135
    return [S_A, S_B, S_AB] #S_A + S_B -S_AB