final version up to sparse stuff and minor modifications

.gitignore

@@ -5,5 +5,5 @@ odeg/timing/__pycache__
 __pycache__/
 *.pyc
 
-fctExp_results/
+results/
 #asdf
\ No newline at end of file

odeg/potentials.py

 import numpy as np
-from numba import njit
-import odeg.fctexp as fe
-from .utils.resultarchivist import fctExpStored
-from .fctexp import thetaToBeta
+from . import fctexp as fe
 
+default_dict = {
+    'stored': True,
+    }
 
-def grandPotential(r0, rs, pos_charge, p_max, theta_grid, mu_grid):
+def grandPotential(r0, rs, pos_charge, p_max, theta_grid, mu_grid, **kwargs):
 
+    settings = default_dict.copy()
+    settings.update(kwargs)
+    kwargs['part_fct'] = True
     
-    log_part_fct = fctExpStored(rs, r0, pos_charge, p_max, mu_grid,
+    if settings['stored']:
+        log_part_fct = fe.fctExpStored(rs, r0, pos_charge, p_max, mu_grid,
+                     theta_grid, **settings)
+    else:
+        log_part_fct = fe.fctExp(rs, r0, pos_charge, p_max, mu_grid,
                      theta_grid)
-    
+        
+        
     gr_pot = - log_part_fct / np.outer(np.ones( np.size(mu_grid) ), theta_grid) 
     
     if np.size(theta_grid) > 1 and np.size(mu_grid) > 1:
@@ -21,27 +29,46 @@ def grandPotential(r0, rs, pos_charge, p_max, theta_grid, mu_grid):
         return float(gr_pot)
         
 
-def freeEnergy(r0, rs, pos_charge, p_max, theta_grid):
+def freeEnergy(r0, rs, pos_charge, p_max, theta_grid, **kwargs):
     
-    mu_min = -0.5
-    mu_max = 0.1
-    mu_grid = np.linspace(mu_min, mu_max, 1000)
+    settings = default_dict.copy()
+    settings.update(kwargs)
+    kwargs['part_fct'] = True
+    kwargs['mu'] = True
     
-    log_part_fct, number_exp = fctExpStored(rs, r0, pos_charge, p_max,
-                                         mu_grid, theta_grid, fe.numberFct, part_fct = True)
+    if settings['stored']:
+        mu, log_part_fct = fe.eqFctExpStored(rs, r0, pos_charge, p_max, theta_grid, **kwargs)
+    else:
+        mu, log_part_fct = fe.eqFctExp(rs, r0, pos_charge, p_max, theta_grid, **kwargs)
     
-    gr_pot = - log_part_fct / np.outer(np.ones( np.size(mu_grid) ), thetaToBeta(theta_grid, rs))
-    free_energy = np.zeros(np.size(theta_grid))
+    free_energy = - log_part_fct/fe.thetaToBeta(theta_grid, rs) + mu * pos_charge
     
-    for b_idx in range(np.size(theta_grid)):
-        
-        m_idx = np.argmin(np.abs( number_exp[:,b_idx] - pos_charge))
-        result_mu[ = mu_grid[m_idx]
-        free_energy[b_idx] = gr_pot[m_idx, b_idx] - mu_grid[m_idx] * pos_charge
+    if np.size(theta_grid) > 1:
+        return free_energy
+    else:
+        return float(free_energy)
     
-    return result_mu
+    
+def entropy(r0, rs, pos_charge, p_max, theta_grid, **kwargs):
+    
+    settings = default_dict.copy()
+    settings.update(kwargs)
+    kwargs['part_fct'] = True
+    kwargs['mu'] = True
+    
+    if settings['stored']:
+        mu, log_part_fct, e_exp = fe.eqFctExpStored(rs, r0, pos_charge, p_max, theta_grid,
+                                                    fe.energyFct, **kwargs)
+    else:
+        mu, log_part_fct, e_exp = fe.eqFctExp(rs, r0, pos_charge, p_max, theta_grid,
+                                              fe.energyFct, **kwargs)
+    
+    entr = log_part_fct + fe.thetaToBeta(theta_grid, rs) * (e_exp - mu * pos_charge)
     
     if np.size(theta_grid) > 1:
-        return free_energy
+        return entr
     else:
-        return float(free_energy)
\ No newline at end of file
+        return float(entr)
+    
+    
+    
\ No newline at end of file

odeg/utils/fctexpcore.py

 import numpy as np
-from numba import njit, int64
+from numba import njit
+import scipy.sparse as sp
 
-from . import bitoperations as bo
 from . import hamiltonian as ham
-from . import neighbourhood as nbh
-import scipy.sparse as sp
+from . import iterators as it
 
 
-debug_mode = False
+fctexpblock_debug = False
 
 # (parallel) worker function that receives assignments from a queue (-> queue)
 # calls core algorithm and writes the results in the shared arrays part_fct and fct_exp
 # p_max, pos_charge: parameters of the hamiltonian
 # v_tilde: fourier transform of the potential, calculated in advance
 # kin_coeff: coefficient of the kinetic contribution
-def fctExpWorker(queue, p_max, pos_charge, spin_boundary_bits, v_tilde, kin_coeff, part_fct,
-                 fct_exp, exponent_shift, shift_init, shift_thr, mu_grid, beta_grid,
-                 fcts, fct_count, wlock):
+def fctExpWorker(mp_queue, p_max, pos_charge, spin_boundary_bits, v_tilde, kin_coeff,
+                 mp_part_fct, mp_fct_exp, mp_exponent_shift, shift_init, shift_thr,
+                 mu_grid, beta_grid, mp_lock, *args):
     
     
     local_part_fct = np.zeros((mu_grid.shape[0], beta_grid.size))
-    local_fct_exp = np.zeros((fct_count, mu_grid.shape[0], beta_grid.size))
+    local_fct_exp = np.zeros((len(args), mu_grid.shape[0], beta_grid.size))
     local_exponent_shift = np.ones((mu_grid.shape[0], beta_grid.size)) * shift_init
     
     while True:
         
-        spin_dist, symmetry_factor, spin_block_size = queue.get()
+        init_fstate, spin_dist, p_total, symmetry, dia_shift = mp_queue.get()
        
-        
         if spin_dist is None:
             break
         
-        spin_dist = np.array(spin_dist)
-        fctExpBlock(spin_dist, spin_block_size, p_max, pos_charge, spin_boundary_bits,
-                    v_tilde, kin_coeff, local_part_fct, local_fct_exp, local_exponent_shift,
-                    shift_thr, mu_grid, beta_grid, fct_count, symmetry_factor, *fcts)
-        
+        spin_dist = np.array(spin_dist)        
+        fctExpBlock(init_fstate, spin_dist, p_total, dia_shift, spin_boundary_bits,
+                    v_tilde, kin_coeff, p_max, local_part_fct, local_fct_exp, 
+                    local_exponent_shift, shift_thr, mu_grid, beta_grid, symmetry, *args)
 
-    with wlock:
+
+    with mp_lock:
         for m_idx in range(mu_grid.shape[0]):
             for b_idx in range(beta_grid.size):
                 
-                global_shift = exponent_shift[m_idx + b_idx * mu_grid.shape[0]]
+                global_shift = mp_exponent_shift[m_idx + b_idx * mu_grid.shape[0]]
                 local_shift = local_exponent_shift[m_idx, b_idx]
                 
                 if global_shift > local_shift:
                     local_part_fct[m_idx, b_idx] *= np.exp(local_shift - global_shift)
                     local_fct_exp[:,m_idx, b_idx] *= np.exp(local_shift - global_shift)
                 else:
-                    part_fct[m_idx + b_idx * mu_grid.shape[0]] *= np.exp(global_shift - local_shift)
+                    mp_part_fct[m_idx + b_idx * mu_grid.shape[0]] *= np.exp(global_shift - local_shift)
                     
-                    for f_idx in range(fct_count):
-                        fct_exp[m_idx + b_idx * mu_grid.shape[0] + f_idx * mu_grid.shape[0] * beta_grid.size
+                    for f_idx in range(len(args)):
+                        mp_fct_exp[m_idx + b_idx * mu_grid.shape[0] + f_idx * mu_grid.shape[0] * beta_grid.size
                                 ] *= np.exp(global_shift - local_shift)
                         
-                    exponent_shift[m_idx + b_idx * mu_grid.shape[0]] = local_shift
+                    mp_exponent_shift[m_idx + b_idx * mu_grid.shape[0]] = local_shift
                     
-                part_fct[m_idx + b_idx * mu_grid.shape[0]] += local_part_fct[m_idx, b_idx]
+                mp_part_fct[m_idx + b_idx * mu_grid.shape[0]] += local_part_fct[m_idx, b_idx]
                 
-                for f_idx in range(fct_count):
-                    fct_exp[m_idx + b_idx * mu_grid.shape[0] + f_idx * mu_grid.shape[0] * beta_grid.size
+                for f_idx in range(len(args)):
+                    mp_fct_exp[m_idx + b_idx * mu_grid.shape[0] + f_idx * mu_grid.shape[0] * beta_grid.size
                             ] += local_fct_exp[f_idx, m_idx, b_idx]
 
 
@@ -73,14 +71,15 @@ def fctExpWorker(queue, p_max, pos_charge, spin_boundary_bits, v_tilde, kin_coef
 # kin_coeff: coefficient of the kinetic contribution
 def fctExpBlock(init_fstate, spin_dist, p_total, dia_shift, spin_boundary_bits,
                 v_tilde, kin_coeff, p_max, part_fct, fct_exp, exponent_shift,
-                shift_thr, mu_grid, beta_grid, fct_count, symmetry, *fcts):
+                shift_thr, mu_grid, beta_grid, symmetry, *fcts):
     
     
     row, col, data, block_size = cooIngredients(init_fstate, spin_boundary_bits, v_tilde, kin_coeff, p_max)
     coo_matrix = sp.coo_array((data, (row, col)), shape=(block_size, block_size))
     
     
-    if debug_mode:
+    ############ DEBUG SECTION START #######################
+    if fctexpblock_debug:
         csr_matrix = coo_matrix.tocsr()
         
         dia_max = csr_matrix[0,0]
@@ -116,13 +115,13 @@ def fctExpBlock(init_fstate, spin_dist, p_total, dia_shift, spin_boundary_bits,
         
         return
     
+    ######### DEBUG SECTION END ################
+    
     matrix = coo_matrix.toarray()
     evs = np.linalg.eigvalsh(matrix) + dia_shift
     occ = spin_dist.sum()
            
     updatePartFct(evs, occ, part_fct, fct_exp, exponent_shift, shift_thr, mu_grid, beta_grid, symmetry)
-    
-
     for f_idx, fct in enumerate(fcts):
         updateFctExp(evs, occ, spin_dist, p_total, fct, f_idx, fct_exp, exponent_shift, shift_thr,
                          mu_grid, beta_grid, symmetry)
@@ -147,23 +146,6 @@ def updatePartFct(evs, occ, part_fct, fct_exp, exponent_shift, shift_thr, mu_gri
                 
             part_fct[m_idx, b_idx] += symmetry * np.sum(np.exp(exponents - shift))
 
-            
-@njit
-def updatePartFctV2(evs, occ, part_fct, fct_exp, exponent_shift, shift_thr, mu_grid, beta_grid, symmetry):
-    
-    for ev in evs:
-        exponents = - beta_grid * (ev - mu_grid * occ)     
-        shift_update = np.maximum(0, exponents - exponent_shift)
-            
-            
-        if np.max(shift_update) > shift_thr:
-            
-            part_fct *= np.exp(-shift_update)
-            fct_exp[:] *= np.exp(-shift_update)
-            
-            exponent_shift += shift_update
-            
-        part_fct += symmetry * np.exp(exponents - exponent_shift)
         
 @njit
 def updateFctExp(evs, occ, spin_dist, p_total, fct, f_idx, fct_exp, exponent_shift, shift_thr,
@@ -174,15 +156,6 @@ def updateFctExp(evs, occ, spin_dist, p_total, fct, f_idx, fct_exp, exponent_shi
             exponents = - beta_grid[b_idx] * (evs - mu_grid[m_idx, b_idx] * occ)
             fct_exp[f_idx, m_idx, b_idx] += symmetry * np.sum(np.exp(exponents - exponent_shift[m_idx, b_idx]) *
                                                     fct(occ, spin_dist, p_total, evs) )
-            
-            
-@njit
-def updateFctExpV2(evs, occ, spin_dist, p_total, fct, f_idx, fct_exp, exponent_shift, shift_thr,
-                 mu_grid, beta_grid, symmetry):
-    
-    for ev in evs:
-        exponents = - beta_grid * (ev - mu_grid * occ)
-        fct_exp[f_idx] += symmetry * np.exp(exponents - exponent_shift) * fct(occ, spin_dist, p_total, evs)
     
 
 
@@ -209,13 +182,8 @@ def cooIngredients(init_fstate, spin_boundary_bits, v_tilde, kin_coeff, p_max):
         fstate = stack_mimic.pop()
         ind = fstate_to_ind[fstate]
         
-        dia = ham.exchangeInteraction(fstate, spin_boundary_bits, v_tilde)
-        bit = 1
-        for spin_idx in range(bo.bitCount(spin_boundary_bits)):
-            for p in range(-p_max, p_max + 1):
-                if bit & fstate:
-                    dia += p*p*kin_coeff
-                bit *= 2
+        dia = ham.exchangeInteraction(fstate, spin_boundary_bits, v_tilde, p_max, kin_coeff)
+        
 
 
         if entry_count == max_size:            
@@ -230,7 +198,7 @@ def cooIngredients(init_fstate, spin_boundary_bits, v_tilde, kin_coeff, p_max):
         entry_count += 1
         
     
-        for neighbour in nbh.neighbourIterator(fstate, spin_boundary_bits):
+        for neighbour in it.neighbourGen(fstate, spin_boundary_bits):
             
             if not(neighbour in fstate_to_ind):
                 n_ind = len(fstate_to_ind)
@@ -265,32 +233,4 @@ def cooIngredients(init_fstate, spin_boundary_bits, v_tilde, kin_coeff, p_max):
                 
                 entry_count += 2
         
-        
-    return row, col, data, fstate_count
-
-# builds the full block and returns all of its eigenvalues
-@njit
-def blockEVs(fstates, v_tilde, spin_boundary_bits, spin_dist, pos_charge, kin_coeff, p_max):
-    # build the matrix
-    n = spin_dist.sum()
-    matrix = np.eye(fstates.size) * (v_tilde[0] * (n*(n-1)//2 - n * pos_charge))
-    
-    for i in range(fstates.size):
-        matrix[i,i] += ham.exchangeInteraction(fstates[i], spin_boundary_bits, v_tilde)
-        
-        bit = 1
-        for spin_idx in range(spin_dist.size):
-            for p in range(-p_max, p_max + 1):
-                if bit & fstates[i]:
-                    matrix[i,i] += p*p*kin_coeff
-                bit *= 2
-    
-    for i in range(fstates.size):
-        for j in range(i):
-            dyn_int = ham.dynamicInteraction(fstates[i], fstates[j], spin_boundary_bits, v_tilde)
-            # dyn_int = j
-            matrix[i,j] = dyn_int
-            matrix[j,i] = dyn_int
-    
-    # calculate eigenvalues
-    return np.linalg.eigvalsh(matrix)
\ No newline at end of file
+    return row, col, data, fstate_count
\ No newline at end of file

odeg/utils/hamiltonian.py

@@ -2,9 +2,18 @@ from numba import njit, int64
 import numpy as np
 from . import bitoperations as bo
 
+
+# calculates interactions that leave the state unchanged (q = 0)
+# plus electron proton interaction
+@njit
+def diaShift(n, vt0, pos_charge):
+    return vt0*(n*(n-1)//2 - n*pos_charge)
+
+
 # calculates interactions that leave the state unchanged (q != 0)
+# plus the kinetic contribution
 @njit
-def exchangeInteraction(fstate, spin_boundary_bits, v_tilde):
+def exchangeInteraction(fstate, spin_boundary_bits, v_tilde, p_max, kin_coeff):
     
     res = 0
     
@@ -17,6 +26,13 @@ def exchangeInteraction(fstate, spin_boundary_bits, v_tilde):
             if state2 & fstate:
                 res -= v_tilde[q]
                 
+    bit = 1
+    for spin_idx in range(bo.bitCount(spin_boundary_bits)):
+        for p in range(-p_max, p_max + 1):
+            if bit & fstate:
+                res += p*p*kin_coeff
+            bit *= 2
+            
     return res
 
 # calculates off-diagonal matrix elements
@@ -38,12 +54,10 @@ def dynamicInteraction(fstate1, fstate2, spin_boundary_bits, v_tilde):
         ones_count += bo.bitCount( (involved_bits[1]-2*involved_bits[0]) & fstate1 )
         sign = 1 - 2 * (ones_count % 2)
         
-        mom_diff1 = bo.bitCount(involved_bits[1] - involved_bits[0])
-        
-        
-        
+        p_diff1 = bo.bitCount(involved_bits[1] - involved_bits[0])
+
         if (involved_bits[2] - involved_bits[1]) & spin_boundary_bits:
-            return sign * v_tilde[mom_diff1]
+            return sign * v_tilde[p_diff1]
         else:
-            mom_diff2 = bo.bitCount(involved_bits[2] - involved_bits[0])
-            return sign * (v_tilde[mom_diff1] - v_tilde[mom_diff2])
\ No newline at end of file
+            p_diff2 = bo.bitCount(involved_bits[2] - involved_bits[0])
+            return sign * (v_tilde[p_diff1] - v_tilde[p_diff2])
\ No newline at end of file

odeg/utils/iterators.py

@@ -101,4 +101,41 @@ def pTotalGen(spin_dist, spin_boundary_bits, p_state_count, h_is_p_symmetric = T
             left_bits -= bit
             right_bits += (bit >> 1)
         else:
-            left_bits -= (bit >> 1)
\ No newline at end of file
+            left_bits -= (bit >> 1)
+            
+            
+@njit
+def neighbourGen(fstate, spin_boundary_bits):
+    
+    
+    for occ_bit_1 in bo.bitIteration(fstate):
+        for occ_bit_2 in bo.bitIteration(fstate):
+            
+            if occ_bit_2 >= occ_bit_1:
+                break
+            
+            vac_bit_1 = occ_bit_1
+            vac_bit_2 = occ_bit_2 >> 1
+            
+            while vac_bit_2 and not( (vac_bit_1 + vac_bit_2) & spin_boundary_bits):
+                        
+                vac_bit_1 <<= 1
+                
+                if not( (vac_bit_1 + vac_bit_2) & fstate ):
+                    
+                    yield (occ_bit_1 + occ_bit_2 + vac_bit_1 + vac_bit_2) ^ fstate
+
+                vac_bit_2 >>= 1    
+                
+                    
+            vac_bit_1 = occ_bit_1 >> 1
+            vac_bit_2 = occ_bit_2 << 1
+            
+            while vac_bit_1 > vac_bit_2 and not( (vac_bit_1 + (vac_bit_2>>1)) & spin_boundary_bits):
+                    
+                if not( (vac_bit_1 + vac_bit_2) & fstate ):
+                    
+                    yield (occ_bit_1 + occ_bit_2 + vac_bit_1 + vac_bit_2) ^ fstate
+                    
+                vac_bit_1 >>= 1
+                vac_bit_2 <<= 1
\ No newline at end of file

odeg/utils/neighbourhood.py

-import numpy as np
-from numba import njit, int64
-from . import bitoperations as bo
-
-
-@njit
-def connectedComponent(fstate, spin_boundary_bits):
-    
-    block_fstates = np.hstack((np.array([fstate]),neighbours(fstate, spin_boundary_bits)))
-    count = block_fstates.size
-    size = count
-
-    ind = 1
-    while ind < count:
-        
-        # print(ind,count,size)
-        local_neighbours = neighbours(block_fstates[ind], spin_boundary_bits)
-        
-        for candidate in local_neighbours:
-            if not( candidate in block_fstates ):
-                if count == size:
-                    block_fstates = np.hstack((block_fstates, np.zeros(size, dtype = int64)))
-                    size *= 2
-                
-                block_fstates[count] = candidate
-                count += 1
-        
-        ind += 1
-            
-    return block_fstates[:count]
-
-
-
-@njit
-def neighbours(fstate, spin_boundary_bits):
-    
-    max_size = 100
-    neigh = np.zeros(max_size, dtype = int64)
-    
-    # common_neighbours = np.zeros(bo.bitLength(spin_boundary_bits), dtype = int)
-    count = 0
-    for occ_bit_1 in bo.bitIteration(fstate):
-        for occ_bit_2 in bo.bitIteration(fstate):
-            
-            if occ_bit_2 >= occ_bit_1:
-                break
-            
-            vac_bit_1 = occ_bit_1
-            vac_bit_2 = occ_bit_2 >> 1
-            
-            while vac_bit_2 and not( (vac_bit_1 + vac_bit_2) & spin_boundary_bits):
-                        
-                vac_bit_1 <<= 1
-                
-                if not( (vac_bit_1 + vac_bit_2) & fstate ):
-                    
-                    neigh[count] = (occ_bit_1 + occ_bit_2 + vac_bit_1 + vac_bit_2) ^ fstate
-                    count += 1
-                    
-                    if count == max_size:
-                        neigh = np.hstack((neigh, np.zeros(max_size, dtype = int64)))
-                        max_size *= 2
-
-                vac_bit_2 >>= 1    
-                
-                    
-            vac_bit_1 = occ_bit_1 >> 1
-            vac_bit_2 = occ_bit_2 << 1
-            
-            while vac_bit_1 > vac_bit_2 and not( (vac_bit_1 + (vac_bit_2>>1)) & spin_boundary_bits):
-                    
-                if not( (vac_bit_1 + vac_bit_2) & fstate ):
-                    
-                    neigh[count] = (occ_bit_1 + occ_bit_2 + vac_bit_1 + vac_bit_2) ^ fstate
-                    count += 1
-                    
-                    if count == max_size:
-                        neigh = np.hstack((neigh, np.zeros(max_size, dtype = int64)))
-                        max_size *= 2
-                    
-                vac_bit_1 >>= 1
-                vac_bit_2 <<= 1
-                
-    return neigh[:count]
-
-
-@njit
-def neighbourIterator(fstate, spin_boundary_bits):
-    
-    
-    for occ_bit_1 in bo.bitIteration(fstate):
-        for occ_bit_2 in bo.bitIteration(fstate):
-            
-            if occ_bit_2 >= occ_bit_1:
-                break
-            
-            vac_bit_1 = occ_bit_1
-            vac_bit_2 = occ_bit_2 >> 1
-            
-            while vac_bit_2 and not( (vac_bit_1 + vac_bit_2) & spin_boundary_bits):
-                        
-                vac_bit_1 <<= 1
-                
-                if not( (vac_bit_1 + vac_bit_2) & fstate ):
-                    
-                    yield (occ_bit_1 + occ_bit_2 + vac_bit_1 + vac_bit_2) ^ fstate
-
-                vac_bit_2 >>= 1    
-                
-                    
-            vac_bit_1 = occ_bit_1 >> 1
-            vac_bit_2 = occ_bit_2 << 1
-            
-            while vac_bit_1 > vac_bit_2 and not( (vac_bit_1 + (vac_bit_2>>1)) & spin_boundary_bits):
-                    
-                if not( (vac_bit_1 + vac_bit_2) & fstate ):
-                    
-                    yield (occ_bit_1 + occ_bit_2 + vac_bit_1 + vac_bit_2) ^ fstate
-                    
-                vac_bit_1 >>= 1
-                vac_bit_2 <<= 1
-
-
-@njit
-def commonNeighbours(fstate1, fstate2, spin_boundary_bits):
-    
-    if fstate1 == fstate2:
-        return neighbours(fstate1, spin_boundary_bits)
-    
-    
-    assert(bo.binDigitSum(fstate1^fstate2) == 4)
-    
-    common_neighbours = np.zeros(3 * bo.bitLength(spin_boundary_bits), dtype = int64)
-    # common_neighbours = np.zeros(bo.bitLength(spin_boundary_bits), dtype = int)
-    count = 0
-    
-    differences = fstate1 ^ fstate2
-    bystanders = fstate1 & (~differences)
-    
-    bits = np.zeros(4, dtype = int64)
-    # bits = np.zeros(4, dtype = int)
-    for ind in range(4):
-        bits[ind] = differences & (~differences + 1)
-        differences ^= bits[ind]
-    
-    differences = fstate1 ^ fstate2
-    dist1 = bo.binDigitSum(bits[3] - bits[2])
-    dist2 = bo.binDigitSum(bits[3] - bits[1])
-    
-    
-    
-    # search for seesaw types corresponding to dist1
-    
-    for left_border in bo.bitIteration(spin_boundary_bits):
-        
-        left_bit = left_border
-        right_bit = left_bit >> dist1
-        
-        while right_bit and not( (left_bit - right_bit) & spin_boundary_bits):
-            
-            if ((left_bit & bystanders) >> dist1) ^ (right_bit & bystanders) and not(
-                    (left_bit + right_bit) & differences):
-                
-                busy_bits = bits[1]+bits[3] if left_bit&bystanders else bits[0]+bits[2]
-                common_neighbours[count] = (busy_bits + bystanders) ^ (left_bit + right_bit)
-                count += 1
-                
-            left_bit >>= 1
-            right_bit >>= 1
-            
-            
-    
-    # search for seesaw types corresponding to dist2
-    
-    if (bits[2]-bits[1]) & spin_boundary_bits:
-        
-        left_bit = bits[3] << 1
-        right_bit = bits[1] << 1
-        while not( (right_bit + left_bit) & (2*spin_boundary_bits)):
-            
-            if ((left_bit & bystanders) >> dist2) ^ (right_bit & bystanders) and not(
-                    (left_bit + right_bit) & differences):
-                
-                busy_bits = bits[3]+bits[2] if left_bit&fstate1 else bits[1] + bits[0]
-                common_neighbours[count] = (busy_bits + bystanders) ^ (left_bit + right_bit)
-                count += 1
-                
-            left_bit <<= 1
-            right_bit <<= 1
-            
-        left_bit = bits[3] >> 1
-        right_bit = bits[1] >> 1
-        while right_bit and not( (left_bit + right_bit) & spin_boundary_bits ):
-            
-            if ((left_bit & bystanders) >> dist2) ^ (right_bit & bystanders) and not(
-                    (left_bit + right_bit) & differences):
-                
-                busy_bits = bits[3]+bits[2] if left_bit&fstate1 else bits[1] + bits[0]
-                common_neighbours[count] = (busy_bits + bystanders) ^ (left_bit + right_bit)
-                count += 1
-                
-            left_bit >>= 1
-            right_bit >>= 1
-        
-    else:
-        for left_border in bo.bitIteration(spin_boundary_bits):
-            
-            left_bit = left_border
-            right_bit = left_bit >> dist2
-            
-            while right_bit and not( (left_bit - right_bit) & spin_boundary_bits):
-                
-                if ((left_bit & bystanders) >> dist2) ^ (right_bit & bystanders) and not(
-                        (left_bit + right_bit) & differences):
-                    
-                    busy_bits = bits[2]+bits[3] if left_bit & bystanders else bits[0]+bits[1]
-                    common_neighbours[count] = busy_bits + (bystanders ^ (left_bit + right_bit))
-                    count += 1
-                    
-                left_bit >>= 1
-                right_bit >>= 1
-    
-    
-    # search for stopover types
-    left_bit = bits[3]
-    right_bit = bits[0] >> 1
-    
-    while right_bit and not( (left_bit + right_bit) & spin_boundary_bits):
-                
-        left_bit <<= 1
-        
-        if left_bit & differences:
-            right_bit >>= 1
-            continue
-        
-        if bo.binDigitSum( (left_bit + right_bit) & bystanders ) == 2:
-            
-            common_neighbours[count] = bits.sum() + bystanders ^ (left_bit + right_bit)
-            # print("c1 found", common_neighbours[count])
-            count += 1
-            
-        elif (left_bit + right_bit) & bystanders == 0:
-            
-            common_neighbours[count] = bystanders ^ (left_bit + right_bit)
-            # print("c2 found", common_neighbours[count])
-            count += 1
-
-        right_bit >>= 1    
-        
-            
-    left_bit = bits[3] >> 1
-    right_bit = bits[0] << 1
-    
-    while right_bit < left_bit and not( (left_bit + (right_bit>>1)) & spin_boundary_bits):
-        
-        if left_bit & differences:
-            left_bit >>= 1
-            right_bit <<= 1
-            continue
-            
-        if bo.binDigitSum( (left_bit + right_bit) & bystanders ) == 2:
-            
-            common_neighbours[count] = bits.sum() + bystanders ^ (left_bit + right_bit)
-            # print("c3 found", common_neighbours[count], "   left_bit, right_bit = ", left_bit, right_bit)
-            count += 1
-            
-        elif (left_bit + right_bit) & bystanders == 0:
-            
-            common_neighbours[count] = bystanders ^ (left_bit + right_bit)
-            # print("c4 found", common_neighbours[count], "   left_bit, right_bit = ", left_bit, right_bit)
-            count += 1
-            
-        left_bit >>= 1
-        right_bit <<= 1
-    
-    return common_neighbours[:count]
-
-
-
-@njit
-def randomNeighbour(fstate, spin_boundary_bits):
-    
-    max_size = bo.bitLength(spin_boundary_bits)
-    neigh = np.zeros(max_size, dtype = int64)
-    count = 0
-    for occ_bit_1 in bo.bitIteration(fstate):
-        for occ_bit_2 in bo.bitIteration(fstate):
-            
-            if occ_bit_2 >= occ_bit_1:
-                break
-            
-            vac_bit_1 = occ_bit_1
-            vac_bit_2 = occ_bit_2 >> 1
-            
-            while vac_bit_2 and not( (vac_bit_1 + vac_bit_2) & spin_boundary_bits):
-                        
-                vac_bit_1 <<= 1
-                
-                if not( (vac_bit_1 + vac_bit_2) & fstate ):
-                    
-                    neigh[count] = (occ_bit_1 + occ_bit_2 + vac_bit_1 + vac_bit_2) ^ fstate
-                    count += 1
-                    
-                    if count == max_size:
-                        neigh = np.hstack((neigh, np.zeros(max_size, dtype = int64)))
-                        max_size *= 2
-                        
-                vac_bit_2 >>= 1    
-                
-                    
-            vac_bit_1 = occ_bit_1 >> 1
-            vac_bit_2 = occ_bit_2 << 1
-            
-            while vac_bit_1 > vac_bit_2 and not( (vac_bit_1 + (vac_bit_2>>1)) & spin_boundary_bits):
-                    
-                if not( (vac_bit_1 + vac_bit_2) & fstate ):
-                    
-                    neigh[count] = (occ_bit_1 + occ_bit_2 + vac_bit_1 + vac_bit_2) ^ fstate
-                    count += 1
-                    
-                    if count == max_size:
-                        neigh = np.hstack((neigh, np.zeros(max_size, dtype = int64)))
-                        max_size *= 2
-                        
-                vac_bit_1 >>= 1
-                vac_bit_2 <<= 1
-                
-                
-    if count == 0:
-        return -1
-    else:
-        return neigh[np.random.randint(count)]
-
-
-@njit
-def randomCommonNeighbour(fstate1, fstate2, spin_boundary_bits):
-    
-    if fstate1 == fstate2:
-        return randomNeighbour(fstate1, spin_boundary_bits)
-    
-    assert(bo.binDigitSum(fstate1^fstate2) == 4)
-    
-    neigh = np.zeros(3 * bo.bitLength(spin_boundary_bits), dtype = int64)
-    count = 0
-    
-    differences = fstate1 ^ fstate2
-    bystanders = fstate1 & (~differences)
-    
-    bits = np.zeros(4, dtype = int64)
-    # bits = np.zeros(4, dtype = int)
-    for ind in range(4):
-        bits[ind] = differences & (~differences + 1)
-        differences ^= bits[ind]
-    
-    differences = fstate1 ^ fstate2
-    dist1 = bo.binDigitSum(bits[3] - bits[2])
-    dist2 = bo.binDigitSum(bits[3] - bits[1])
-    
-    
-    
-    # search for seesaw types corresponding to dist1
-    
-    for left_border in bo.bitIteration(spin_boundary_bits):
-        
-        left_bit = left_border
-        right_bit = left_bit >> dist1
-        
-        while right_bit and not( (left_bit - right_bit) & spin_boundary_bits):
-            
-            if ((left_bit & bystanders) >> dist1) ^ (right_bit & bystanders) and not(
-                    (left_bit + right_bit) & differences):
-                
-                busy_bits = bits[1]+bits[3] if left_bit&bystanders else bits[0]+bits[2]
-                neigh[count] = (busy_bits + bystanders) ^ (left_bit + right_bit)
-                count += 1
-                
-            left_bit >>= 1
-            right_bit >>= 1
-            
-            
-    
-    # search for seesaw types corresponding to dist2
-    
-    if (bits[2]-bits[1]) & spin_boundary_bits:
-        
-        left_bit = bits[3] << 1
-        right_bit = bits[1] << 1
-        while not( (right_bit + left_bit) & (2*spin_boundary_bits)):
-            
-            if ((left_bit & bystanders) >> dist2) ^ (right_bit & bystanders) and not(
-                    (left_bit + right_bit) & differences):
-                
-                busy_bits = bits[3]+bits[2] if left_bit&fstate1 else bits[1] + bits[0]
-                neigh[count] = (busy_bits + bystanders) ^ (left_bit + right_bit)
-                count += 1
-                
-            left_bit <<= 1
-            right_bit <<= 1
-            
-        left_bit = bits[3] >> 1
-        right_bit = bits[1] >> 1
-        while right_bit and not( (left_bit + right_bit) & spin_boundary_bits ):
-            
-            if ((left_bit & bystanders) >> dist2) ^ (right_bit & bystanders) and not(
-                    (left_bit + right_bit) & differences):
-                
-                busy_bits = bits[3]+bits[2] if left_bit&fstate1 else bits[1] + bits[0]
-                neigh[count] = (busy_bits + bystanders) ^ (left_bit + right_bit)
-                count += 1
-                
-            left_bit >>= 1
-            right_bit >>= 1
-        
-    else:
-        for left_border in bo.bitIteration(spin_boundary_bits):
-            
-            left_bit = left_border
-            right_bit = left_bit >> dist2
-            
-            while right_bit and not( (left_bit - right_bit) & spin_boundary_bits):
-                
-                if ((left_bit & bystanders) >> dist2) ^ (right_bit & bystanders) and not(
-                        (left_bit + right_bit) & differences):
-                    
-                    busy_bits = bits[2]+bits[3] if left_bit & bystanders else bits[0]+bits[1]
-                    neigh[count] = busy_bits + (bystanders ^ (left_bit + right_bit))
-                    count += 1
-                    
-                left_bit >>= 1
-                right_bit >>= 1
-    
-    
-    # search for stopover types
-    left_bit = bits[3]
-    right_bit = bits[0] >> 1
-    
-    while right_bit and not( (left_bit + right_bit) & spin_boundary_bits):
-                
-        left_bit <<= 1
-        
-        if left_bit & differences:
-            right_bit >>= 1
-            continue
-        
-        if bo.binDigitSum( (left_bit + right_bit) & bystanders ) == 2:
-            
-            neigh[count] = bits.sum() + bystanders ^ (left_bit + right_bit)
-            count += 1
-            
-        elif (left_bit + right_bit) & bystanders == 0:
-            
-            neigh[count] = bystanders ^ (left_bit + right_bit)
-            count += 1
-
-        right_bit >>= 1    
-        
-            
-    left_bit = bits[3] >> 1
-    right_bit = bits[0] << 1
-    
-    while right_bit < left_bit and not( (left_bit + (right_bit>>1)) & spin_boundary_bits):
-        
-        if left_bit & differences:
-            left_bit >>= 1
-            right_bit <<= 1
-            continue
-            
-        if bo.binDigitSum( (left_bit + right_bit) & bystanders ) == 2:
-            
-            neigh[count] = bits.sum() + bystanders ^ (left_bit + right_bit)
-            count += 1
-            
-        elif (left_bit + right_bit) & bystanders == 0:
-            
-            neigh[count] = bystanders ^ (left_bit + right_bit)
-            count += 1
-            
-        left_bit >>= 1
-        right_bit <<= 1
-    
-    
-    
-    if count == 0:
-        return -1
-    else:
-        return neigh[np.random.randint(count)]
\ No newline at end of file

odeg/utils/resultencoder.py

+import numpy as np
+
+
+def isLinspace(array, tol = 1e-8):
+    if len(np.shape(array)) != 1 or array.size == 1:
+        return False
+    
+    if tol < np.linalg.norm(array - np.linspace(array[0], array[-1], array.size)):
+        return False
+    else:
+        return True
+       
+        
+
+def arrayEncoder(name, array, svc_max):
+    
+    if isLinspace(array):
+        return "_" + name + "LS_%g_%g_%d" % (array[0], array[-1], array.size)
+    
+    elif array.size <= svc_max:
+        res = "_" + name + "SV"
+        for value in array:
+            res += "_%g" % value
+        return res
+    
+    else:
+        return "toomanynonlinspacevalues"
\ No newline at end of file

odeg/utils/tobenamed.py

-from numba import njit
-import numpy as np
-from math import comb
-
-@njit
-def isLinspace(array):
-    if len(np.shape(array)) != 1:
-        return False
-    
-    if 1e-4 < np.linalg.norm(array - np.linspace(array[0], array[-1], array.size)):
-        return False
-    else:
-        return True
-    
-    
-@njit
-def diaShift(n, vt0, pos_charge):
-    return vt0*(n*(n-1)//2 - n*pos_charge)    
-
-
-@njit
-def numberExpApp(n_max, vt0, pos_charge, beta_grid, mu_grid):
-    
-    n_exp = np.zeros((mu_grid.size, beta_grid.size))
-    part_fct = np.zeros((mu_grid.size, beta_grid.size))
-    shifts = np.zeros((mu_grid.size, beta_grid.size))
-    bin_coeff = 1
-    
-    for n in range(n_max+1):
-        
-        exponents = np.outer(diaShift(n, vt0, pos_charge) - mu_grid * n, -beta_grid)
-        
-        for m_idx in range(mu_grid.size):
-            for b_idx in range(beta_grid.size):
-                
-                em = exponents[m_idx, b_idx]
-                if em > shifts[m_idx, b_idx]:
-                    n_exp[m_idx, b_idx] *= np.exp(shifts[m_idx, b_idx] - em)
-                    part_fct[m_idx, b_idx] *= np.exp(shifts[m_idx, b_idx] - em)
-                    shifts[m_idx, b_idx] = em
-        
-        
-        n_exp += n * bin_coeff * np.exp(exponents-shifts)
-        part_fct += bin_coeff * np.exp(exponents-shifts)
-        bin_coeff *= (n_max - n)
-        bin_coeff //= (n + 1)
-    
-    return n_exp/part_fct
-
-
-# @njit
-# def eqMuGuess(n_max, vt0, pos_charge, beta_grid, match_tol, mu_res):
-    
-#     res = np.ones(beta_grid.size, mu_res + 1)
-    
-#     mu_grid = np.array([])
-#     n_exp_target = np.array([[pos_charge - match_tol], [pos_charge + match_tol]])
-    
-    
-#     n_exp = numberExpApp(n_max, vt0, pos_charge, beta_grid, mu)
-    
-    
-#     while np.max(res) > mu_res:
-        
-        
-        
-        
-        
-#         n_exp = numberExpApp(n_max, vt0, pos_charge, beta_grid, mu)
-        
-        
-        
-#     return mu
\ No newline at end of file

test_fctexp.py

 import numpy as np
 import matplotlib.pyplot as plt
 from time import time
-from math import comb
 
 plt.rcParams['figure.dpi'] = 300
 plt.rcParams['savefig.dpi'] = 300
@@ -9,53 +8,52 @@ plt.rcParams['savefig.dpi'] = 300
 
 # import odeg.old.fctexp_05 as fe
 import odeg.fctexp as fe
-from odeg.utils.tobenamed import numberExpApp
 
 
 if __name__ == '__main__':
     
-    r_0 = 0.5
-    r_s = 4
-    p_max = 3
-    pos_charge = 7
-    # pos_charge = 7
+    r_0 = 0.6
+    r_s = 3
+    p_max = 4
+    pos_charge = 8
+    
     
-    # method = "matchtol"
     method = "min"
     
     # mu_grid = np.array([2])
-    mu_grid = np.linspace(-0.3, 0.1, 300)
-    theta_grid1 = np.linspace(0.01, 2, 100)
-    theta_grid = np.array([0.05])
-    pm = 1
+    mu_grid = np.linspace(-3, 3, 300)
+    theta_grid1 = np.linspace(0.5, 4, 50)
+    theta_grid = np.array([0.01])
+    mode = 1
     
     length = 2 * r_s * pos_charge
     vt0 = 2/length*(2 + np.log(length/(2*r_0)))
     
     start_time = time()
-    e1 = fe.fctExp(r_s, r_0, pos_charge, p_max, mu_grid, theta_grid, fe.numberFct,
-                    parallel = False)
+    # e1 = fe.fctExp(r_s, r_0, pos_charge, p_max, mu_grid, theta_grid, fe.numberFct,
+    #                 parallel = False)
     
-    if pm == 0:
+    if mode == 0:
         e1 = fe.fctExpStored(r_s, r_0, pos_charge, p_max, mu_grid, theta_grid,
-                                fe.numberFct, parallel = False)
-    elif pm == 1:
-        mu, pf, n_exp, e_exp = fe.eqFctExp(r_s, r_0, pos_charge, p_max, theta_grid1,
+                                fe.numberFct, parallel = False, process_count = 4)
+    elif mode == 1:
+        mu, pf, n_exp, e_exp, n2_exp = fe.eqFctExpStored(r_s, r_0, pos_charge, p_max, theta_grid1,
                                 fe.numberFct,
                                 fe.energyFct,
+                                fe.number2Fct,
                                 mu = True,
                                 part_fct = True,
                                 method = method,
                                 tol = 1e-2,
                                 match_tol = 1e-8,
-                                zoom = (100,100,100))
+                                zoom = (300,300,300,300))
                            
                            
     print('time used: ' + str(time()-start_time))
     
     # print(ex.numberToGenMu(evs, beta, pos_charge_count))
     
-    if pm == 0:
+    if mode == 0:
         plt.rcParams['text.usetex'] = True
         lbl  = label = [(r'$\Theta$ = ' + str(t)) for t in theta_grid] if theta_grid.size > 1 else r'$\Theta$ = ' + str(theta_grid[0]) 
         fig, ax = plt.subplots()
@@ -70,13 +68,13 @@ if __name__ == '__main__':
         # app = np.outer(mu_grid, beta_grid)
         # app = nmax/(1+np.exp(-app))
         start = time()
-        app = numberExpApp(nmax, vt0, pos_charge, beta_grid, mu_grid)
+        app = fe.numberExpApp(nmax, vt0, pos_charge, beta_grid, mu_grid)
         print("appTime:", time()-start)
         ax. plot(mu_grid, app, c = 'red')
     
-    elif pm == 1:
+    elif mode == 1:
         
-        # free_energy = - pf/fe.thetaToBeta(theta_grid1, r_s) + mu * pos_charge
+        free_energy = - pf/fe.thetaToBeta(theta_grid1, r_s) + mu * pos_charge
         
         fig, ax = plt.subplots()
         ax.plot(theta_grid1, np.transpose(mu), label = method)

test_potentials.py

@@ -36,7 +36,7 @@ if __name__ == '__main__':
     print('time used: ' + str(time()-start_time))
     
     
-    np.save('results/testrun_fe', fr_en_dpp)
+    # np.save('results/testrun_fe', fr_en_dpp)
     
     
     # fig, ax = plt.subplots()