Test of Renku

script for computing 1D spin-1/2 correlation and magnetization

not optimal, it is a script (it does the job) ;-)

README.md

 # 1D XY spin-half model
-A Renku project.
+
+The longitudinal (SzSz) correlation is the same than the
+density-density correlation in a free fermionic chain (up to a shift
+in the reciprocal space)
+
+The transverse (SxSx) correlation is however computed using a pfaffian
+operation over an antisymetric matrix due to vick theorem
+
+All results can be found in Cruz and Goncalvez publication (maped on
+hole)
\ No newline at end of file

pfapack/README.org

+https://arxiv.org/pdf/1102.3440.pdf
+
+code @
+https://michaelwimmer.org/downloads.html
+

pfapack/__init__.py

+import pfaffian

pfapack/pfaffian.py

+"""A package for computing Pfaffians"""
+
+from __future__ import division
+import numpy as np
+import scipy.linalg as la
+import scipy.sparse as sp
+import math, cmath
+
+def householder_real(x):
+    """(v, tau, alpha) = householder_real(x)
+
+    Compute a Householder transformation such that
+    (1-tau v v^T) x = alpha e_1
+    where x and v a real vectors, tau is 0 or 2, and
+    alpha a real number (e_1 is the first unit vector)
+    """
+
+    assert x.shape[0]>0
+
+    sigma=np.dot(x[1:],x[1:])
+
+    if sigma==0:
+        return (np.zeros(x.shape[0]), 0, x[0])
+    else:
+        norm_x=math.sqrt(x[0]**2+sigma)
+
+        v=x.copy();
+
+        #depending on whether x[0] is positive or negatvie
+        #choose the sign
+        if x[0]<=0:
+            v[0]-=norm_x
+            alpha=+norm_x
+        else:
+            v[0]+=norm_x
+            alpha=-norm_x
+
+        v/=np.linalg.norm(v)
+
+        return (v, 2, alpha)
+
+def householder_complex(x):
+    """(v, tau, alpha) = householder_real(x)
+
+    Compute a Householder transformation such that
+    (1-tau v v^T) x = alpha e_1
+    where x and v a complex vectors, tau is 0 or 2, and
+    alpha a complex number (e_1 is the first unit vector)
+    """
+    assert x.shape[0]>0
+
+    sigma=np.dot(np.conj(x[1:]), x[1:])
+
+    if sigma==0:
+        return (np.zeros(x.shape[0]), 0, x[0])
+    else:
+        norm_x=cmath.sqrt(x[0].conjugate()*x[0]+sigma)
+
+        v=x.copy();
+
+        phase=cmath.exp(1j*math.atan2(x[0].imag, x[0].real))
+
+        v[0]+=phase*norm_x
+
+        v/=np.linalg.norm(v)
+
+    return (v, 2, -phase*norm_x)
+
+def skew_tridiagonalize(A, overwrite_a=False, calc_q=True):
+    """ T, Q = skew_tridiagonalize(A, overwrite_a, calc_q=True)
+
+    or
+
+    T = skew_tridiagonalize(A, overwrite_a, calc_q=False)
+
+    Bring a real or complex skew-symmetric matrix (A=-A^T) into
+    tridiagonal form T (with zero diagonal) with a orthogonal
+    (real case) or unitary (complex case) matrix U such that
+    A = Q T Q^T
+    (Note that Q^T and *not* Q^dagger also in the complex case)
+
+    A is overwritten if overwrite_a=True (default: False), and
+    Q only calculated if calc_q=True (default: True)
+    """
+
+    #Check if matrix is square
+    assert A.shape[0] == A.shape[1] > 0
+    #Check if it's skew-symmetric
+    assert abs((A+A.T).max())<1e-14
+
+    n = A.shape[0]
+    A = np.asarray(A)  #the slice views work only properly for arrays
+
+    #Check if we have a complex data type
+    if np.issubdtype(A.dtype, np.complexfloating):
+        householder = householder_complex
+    elif not np.issubdtype(A.dtype, np.number):
+        raise TypeError("pfaffian() can only work on numeric input")
+    else:
+        householder = householder_real
+
+    if not overwrite_a:
+        A = A.copy()
+
+    if calc_q:
+        Q = np.eye(A.shape[0], dtype=A.dtype)
+
+    for i in xrange(A.shape[0]-2):
+        #Find a Householder vector to eliminate the i-th column
+        v, tau, alpha = householder(A[i+1:,i])
+        A[i+1, i] = alpha
+        A[i, i+1] = -alpha
+        A[i+2:, i] = 0
+        A[i, i+2:] = 0
+
+        #Update the matrix block A(i+1:N,i+1:N)
+        w = tau*np.dot(A[i+1:, i+1:], v.conj());
+        A[i+1:,i+1:]+=np.outer(v,w)-np.outer(w,v)
+
+        if calc_q:
+            #Accumulate the individual Householder reflections
+            #Accumulate them in the form P_1*P_2*..., which is
+            # (..*P_2*P_1)^dagger
+            y = tau*np.dot(Q[:, i+1:], v)
+            Q[:, i+1:]-=np.outer(y,v.conj())
+
+    if calc_q:
+        return (np.asmatrix(A), np.asmatrix(Q))
+    else:
+        return np.asmatrix(A)
+
+def skew_LTL(A, overwrite_a=False, calc_L=True, calc_P=True):
+    """ T, L, P = skew_LTL(A, overwrite_a, calc_q=True)
+
+    Bring a real or complex skew-symmetric matrix (A=-A^T) into
+    tridiagonal form T (with zero diagonal) with a lower unit
+    triangular matrix L such that
+    P A P^T= L T L^T
+
+    A is overwritten if overwrite_a=True (default: False),
+    L and P only calculated if calc_L=True or calc_P=True,
+    respectively (default: True).
+    """
+
+    #Check if matrix is square
+    assert A.shape[0] == A.shape[1] > 0
+    #Check if it's skew-symmetric
+    assert abs((A+A.T).max())<1e-14
+
+    n = A.shape[0]
+    A = np.asarray(A)  #the slice views work only properly for arrays
+
+    if not overwrite_a:
+        A = A.copy()
+
+    if calc_L:
+        L = np.eye(n, dtype=A.dtype)
+
+    if calc_P:
+        Pv = np.arange(n)
+
+    for k in xrange(n-2):
+        #First, find the largest entry in A[k+1:,k] and
+        #permute it to A[k+1,k]
+        kp = k+1+np.abs(A[k+1:,k]).argmax()
+
+        #Check if we need to pivot
+        if kp != k+1:
+            #interchange rows k+1 and kp
+            temp = A[k+1,k:].copy()
+            A[k+1,k:] = A[kp,k:]
+            A[kp,k:] = temp
+
+            #Then interchange columns k+1 and kp
+            temp = A[k:,k+1].copy()
+            A[k:,k+1] = A[k:,kp]
+            A[k:,kp] = temp
+
+            if calc_L:
+                #permute L accordingly
+                temp = L[k+1,1:k+1].copy()
+                L[k+1,1:k+1] = L[kp,1:k+1]
+                L[kp,1:k+1] = temp
+
+            if calc_P:
+                #accumulate the permutation matrix
+                temp = Pv[k+1]
+                Pv[k+1] = Pv[kp]
+                Pv[kp] = temp
+
+        #Now form the Gauss vector
+        if A[k+1,k] != 0.0:
+            tau = A[k+2:,k].copy()
+            tau /= A[k+1,k]
+
+            #clear eliminated row and column
+            A[k+2:,k] = 0.0
+            A[k,k+2:] = 0.0
+
+            #Update the matrix block A(k+2:,k+2)
+            A[k+2:,k+2:] += np.outer(tau, A[k+2:,k+1])
+            A[k+2:,k+2:] -= np.outer(A[k+2:,k+1], tau)
+
+            if calc_L:
+                L[k+2:,k+1] = tau
+
+    if calc_P:
+        #form the permutation matrix as a sparse matrix
+        P = sp.csr_matrix( (np.ones(n), (np.arange(n), Pv)) )
+
+    if calc_L:
+        if calc_P:
+            return (np.asmatrix(A), np.asmatrix(L), P)
+        else:
+            return (np.asmatrix(A), np.asmatrix(L))
+    else:
+        if calc_P:
+            return (np.asmatrix(A), P)
+        else:
+            return np.asmatrix(A)
+
+def pfaffian(A, overwrite_a=False, method='P'):
+    """ pfaffian(A, overwrite_a=False, method='P')
+
+    Compute the Pfaffian of a real or complex skew-symmetric
+    matrix A (A=-A^T). If overwrite_a=True, the matrix A
+    is overwritten in the process. This function uses
+    either the Parlett-Reid algorithm (method='P', default),
+    or the Householder tridiagonalization (method='H')
+    """
+    #Check if matrix is square
+    assert A.shape[0] == A.shape[1] > 0
+    #Check if it's skew-symmetric
+    assert abs((A+A.T).max())<1e-14
+    #Check that the method variable is appropriately set
+    assert method == 'P' or method == 'H'
+
+    # Make sure the matrix is using floating point algebra
+    if not np.issubdtype(A.dtype, np.inexact):
+        A = 1.0 * A
+
+    if method == 'P':
+        return pfaffian_LTL(A, overwrite_a)
+    else:
+        return pfaffian_householder(A, overwrite_a)
+
+def pfaffian_LTL(A, overwrite_a=False):
+    """ pfaffian_LTL(A, overwrite_a=False)
+
+    Compute the Pfaffian of a real or complex skew-symmetric
+    matrix A (A=-A^T). If overwrite_a=True, the matrix A
+    is overwritten in the process. This function uses
+    the Parlett-Reid algorithm.
+    """
+    #Check if matrix is square
+    assert A.shape[0] == A.shape[1] > 0
+    #Check if it's skew-symmetric
+    assert abs((A+A.T).max())<1e-14
+
+    n = A.shape[0]
+    A = np.asarray(A)  #the slice views work only properly for arrays
+
+    #Quick return if possible
+    if n%2==1:
+        return 0
+
+    if not overwrite_a:
+        A = A.copy()
+
+    pfaffian_val = 1.0
+
+    for k in xrange(0, n-1, 2):
+        #First, find the largest entry in A[k+1:,k] and
+        #permute it to A[k+1,k]
+        kp = k+1+np.abs(A[k+1:,k]).argmax()
+
+        #Check if we need to pivot
+        if kp != k+1:
+            #interchange rows k+1 and kp
+            temp = A[k+1,k:].copy()
+            A[k+1,k:] = A[kp,k:]
+            A[kp,k:] = temp
+
+            #Then interchange columns k+1 and kp
+            temp = A[k:,k+1].copy()
+            A[k:,k+1] = A[k:,kp]
+            A[k:,kp] = temp
+
+            #every interchange corresponds to a "-" in det(P)
+            pfaffian_val *= -1
+
+        #Now form the Gauss vector
+        if A[k+1,k] != 0.0:
+            tau = A[k,k+2:].copy()
+            tau /= A[k,k+1]
+
+            pfaffian_val *= A[k,k+1]
+
+            if k+2<n:
+                #Update the matrix block A(k+2:,k+2)
+                A[k+2:,k+2:] += np.outer(tau, A[k+2:,k+1])
+                A[k+2:,k+2:] -= np.outer(A[k+2:,k+1], tau)
+        else:
+            #if we encounter a zero on the super/subdiagonal, the
+            #Pfaffian is 0
+            return 0.0
+
+    return pfaffian_val
+
+
+def pfaffian_householder(A, overwrite_a=False):
+    """ pfaffian(A, overwrite_a=False)
+
+    Compute the Pfaffian of a real or complex skew-symmetric
+    matrix A (A=-A^T). If overwrite_a=True, the matrix A
+    is overwritten in the process. This function uses the
+    Householder tridiagonalization.
+
+    Note that the function pfaffian_schur() can also be used in the
+    real case. That function does not make use of the skew-symmetry
+    and is only slightly slower than pfaffian_householder().
+    """
+
+    #Check if matrix is square
+    assert A.shape[0] == A.shape[1] > 0
+    #Check if it's skew-symmetric
+    assert abs((A+A.T).max())<1e-14
+
+    n = A.shape[0]
+
+    #Quick return if possible
+    if n%2==1:
+        return 0
+
+    #Check if we have a complex data type
+    if np.issubdtype(A.dtype, np.complexfloating):
+        householder=householder_complex
+    elif not np.issubdtype(A.dtype, np.number):
+        raise TypeError("pfaffian() can only work on numeric input")
+    else:
+        householder=householder_real
+
+    A = np.asarray(A)  #the slice views work only properly for arrays
+
+    if not overwrite_a:
+        A = A.copy()
+
+    pfaffian_val = 1.
+
+    for i in xrange(A.shape[0]-2):
+        #Find a Householder vector to eliminate the i-th column
+        v, tau, alpha = householder(A[i+1:,i])
+        A[i+1, i] = alpha
+        A[i, i+1] = -alpha
+        A[i+2:, i] = 0
+        A[i, i+2:] = 0
+
+        #Update the matrix block A(i+1:N,i+1:N)
+        w = tau*np.dot(A[i+1:, i+1:], v.conj());
+        A[i+1:,i+1:]+=np.outer(v,w)-np.outer(w,v)
+
+        if tau!=0:
+            pfaffian_val *= 1-tau
+        if i%2==0:
+            pfaffian_val *= -alpha
+
+    pfaffian_val *= A[n-2,n-1]
+
+    return pfaffian_val
+
+def pfaffian_schur(A, overwrite_a=False):
+    """Calculate Pfaffian of a real antisymmetric matrix using
+    the Schur decomposition. (Hessenberg would in principle be faster,
+    but scipy-0.8 messed up the performance for scipy.linalg.hessenberg()).
+
+    This function does not make use of the skew-symmetry of the matrix A,
+    but uses a LAPACK routine that is coded in FORTRAN and hence faster
+    than python. As a consequence, pfaffian_schur is only slightly slower
+    than pfaffian().
+    """
+
+    assert np.issubdtype(A.dtype, np.number) and not np.issubdtype(A.dtype, np.complexfloating)
+
+    assert A.shape[0] == A.shape[1] > 0
+
+    assert abs(A + A.T).max() < 1e-14
+
+    #Quick return if possible
+    if A.shape[0]%2 == 1:
+        return 0
+
+    (t, z) = la.schur(A, output='real', overwrite_a=overwrite_a)
+    l = np.diag(t, 1)
+    return np.prod(l[::2]) * la.det(z)

requirements.txt

+numpy
+scipy
+pfapack

sxsx.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy
+import pfapack.pfaffian as pf
+
+##################
+# spin-1/2 chain #
+##################
+
+# chain length
+N = 40
+
+# <Sx(xi,ti)Sx(xj,0)>
+xi = numpy.array([20, 21])
+ti = numpy.array([-0.2, 0., 0.2, 0.4])
+xj = numpy.array([20])
+# note that the code is not threaded thus slow. Not a problem, because
+# it is a script and we want fast solution rather than efficient
+# coding at this point
+
+# ----------------------------------------------------------
+
+Jxy = 1. # energy scale
+hz = 0.8 # external magnetic field
+T = 0.002 # temperature
+B = 1./T # 1/T
+
+dk = numpy.pi/(N+1.);
+k = dk*numpy.linspace(1,N,N); # numpy.linspace(dk,N*dk,N);
+
+# ----------------------------------------------------------
+
+
+# dispersion
+E = lambda kn: Jxy*numpy.cos(kn) - hz;
+
+AitAj = lambda i,t,j: (2./(N+1.))*( numpy.sum( \
+        ( numpy.sin(k*j) * numpy.sin(k*i) ) * \
+        ( +(1.+0.j) * numpy.cos(E(k)*t) - (0.+1.j) * numpy.sin(E(k)*t) * numpy.tanh(B*E(k)/2.) ) \
+                  )     )
+AitBj = lambda i,t,j: (2./(N+1.))*( numpy.sum( \
+        ( numpy.sin(k*j) * numpy.sin(k*i) ) * \
+        ( -(0.+1.j) * numpy.sin(E(k)*t) + (1.+0.j) * numpy.cos(E(k)*t) * numpy.tanh(B*E(k)/2.) ) \
+                  )     )
+BitBj = lambda i,t,j: -AitAj(i,t,j)
+BitAj = lambda i,t,j: -AitBj(i,t,j)
+
+
+#-----------------------------------------------------------
+
+sxsxitj = numpy.zeros((len(xi),len(ti),len(xj)),dtype='complex128')
+
+for tindex,time in enumerate(ti):
+    for iindex,i in enumerate(xi.astype(int)):
+        for jindex,j in enumerate(xj.astype(int)):
+            # print("i,j",i,j)
+            
+            dim = 2*(i+j-1)
+            PfafA = numpy.zeros((dim,dim),dtype='complex128')
+
+            ### 2i-1 x 2i-2
+            for matiindex in xrange(0,i):
+                if matiindex != i-1:
+                    for matjindex in xrange(matiindex,i-1):
+                        # print('@[{0},1+{1}] : <A_{2}(0)B_{3}>'.format(2*matiindex,2*matjindex,matiindex+1,matjindex+1))
+                        PfafA[2*matiindex,1+2*matjindex] = AitBj(matiindex+1,0,matjindex+1)
+                for matjindex in xrange(matiindex+1,i):
+                    # print('@[{0},1+{1}] : <A_{2}(0)A_{3}>'.format(2*matiindex,2*matjindex-1,matiindex+1,matjindex+1))
+                    PfafA[2*matiindex,1+2*matjindex-1] = AitAj(matiindex+1,0,matjindex+1)            
+
+            for matiindex in xrange(1,i):
+                for matjindex in xrange(matiindex,i):
+                    # print('@[{0},1+{1}] : <B_{2}(0)A_{3}>'.format(2*matiindex-1,2*matjindex-1,matiindex,matjindex+1))
+                    PfafA[2*matiindex-1,1+2*matjindex-1] = BitAj(matiindex,0,matjindex+1)
+                if matiindex != i-1:
+                    for matjindex in xrange(matiindex,i-1):
+                        # print('@[{0},1+{1}] : <B_{2}(0)B_{3}>'.format(2*matiindex-1,2*matjindex,matiindex,matjindex+1))
+                        PfafA[2*matiindex-1,1+2*matjindex] = BitBj(matiindex,0,matjindex+1)
+                        
+            ### 2i-1 x 2j-1
+            for matiindex in xrange(0,i):
+                for matjindex in xrange(0,j):
+                    # print('@[{0},1+2*i-2+{1}={2}] : <A_{3}(t)A_{4}>'.format(2*matiindex,2*matjindex,1+2*i-2+2*matjindex,matiindex+1,matjindex+1))
+                    PfafA[2*matiindex,1+2*i-2+2*matjindex] = AitAj(matiindex+1,time,matjindex+1)
+                if matiindex != i:
+                    for matjindex in xrange(1,j):
+                        # print('@[{0},1+2*i-2+{1}={2}] : <A_{3}(t)B_{4}>'.format(2*matiindex,2*matjindex-1,1+2*i-2+2*matjindex-1,matiindex+1,matjindex))
+                        PfafA[2*matiindex,1+2*i-2+2*matjindex-1] = AitBj(matiindex+1,time,matjindex)
+            for matiindex in xrange(1,i):
+                for matjindex in xrange(0,j):
+                    # print('@[{0},1+2*i-2+{1}={2}] : <B_{3}(t)A_{4}>'.format(2*matiindex-1,2*matjindex,1+2*i-2+2*matjindex,matiindex,matjindex+1))
+                    PfafA[2*matiindex-1,1+2*i-2+2*matjindex] = BitAj(matiindex,time,matjindex+1)
+                if matiindex != i:
+                    for matjindex in xrange(1,j):
+                        # print('@[{0},1+2*i-2+{1}={2}] : <B_{3}(t)B_{4}>'.format(2*matiindex-1,2*matjindex-1,1+2*i-2+2*matjindex-1,matiindex,matjindex))
+                        PfafA[2*matiindex-1,1+2*i-2+2*matjindex-1] = BitBj(matiindex,time,matjindex)
+
+
+            ### 2j-2 x 2j-1
+            for matiindex in xrange(0,j-1):
+                for matjindex in xrange(matiindex,j):
+                    if matjindex!=j-1:
+                        # print('@[2*i-1+{0}={1},1+2*i-1+{2}={3}] : <A_{4}(0)B_{5}>'.format(2*matiindex,2*i-1+2*matiindex,2*matjindex,1+2*i-1+2*matjindex,matiindex+1,matjindex+1))
+                        PfafA[2*i-1+2*matiindex,1+2*i-1+2*matjindex] = AitBj(matiindex+1,0,matjindex+1)
+                for matjindex in xrange(matiindex+1,j):
+                    # print('@[2*i-1+{0}={1},1+2*i-1+{2}={3}] : <A_{4}(0)A_{5}>'.format(2*matiindex,2*i-1+2*matiindex,2*matjindex-1,1+2*i-1+2*matjindex-1,matiindex+1,matjindex+1))
+                    PfafA[2*i-1+2*matiindex,1+2*i-1+2*matjindex-1] = AitAj(matiindex+1,0,matjindex+1)
+            for matiindex in xrange(1,j):
+                for matjindex in xrange(matiindex,j):
+                    # print('@[2*i-1+{0}={1},1+2*i-1+{2}={3}] : <B_{4}(0)A_{5}>'.format(2*matiindex-1,2*i-1+2*matiindex-1,2*matjindex-1,1+2*i-1+2*matjindex-1,matiindex,matjindex+1))
+                    PfafA[2*i-1+2*matiindex-1,1+2*i-1+2*matjindex-1] = BitAj(matiindex,0,matjindex+1)
+                for matjindex in xrange(matiindex,j):
+                    if matjindex!=j-1:
+                        # print('@[2*i-1+{0}={1},1+2*i-1+{2}={3}] : <B_{4}(0)B_{5}>'.format(2*matiindex-1,2*i-1+2*matiindex-1,2*matjindex,1+2*i-1+2*matjindex,matiindex,matjindex+1))
+                        PfafA[2*i-1+2*matiindex-1,1+2*i-1+2*matjindex] = BitBj(matiindex,0,matjindex+1)
+
+            # sxsxitj[iindex,tindex,jindex] = numpy.sqrt(numpy.linalg.det(PfafA-numpy.transpose(PfafA)))/4. # sign problem ?!
+            sxsxitj[iindex,tindex,jindex] = pf.pfaffian(PfafA-numpy.transpose(PfafA))/4.
+
+
+
+# save to filename
+filename = 'data/sxsx_N{N}_sh-xy_Jxy{Jxy}_hz{hz}_B{B}_szsz_xi{ximin}s{ximax}_ti{timin}s{timax}_xj{xjmin}s{xjmax}'.format(N=N,Jxy=Jxy,hz=hz,B=B,ximin=numpy.min(xi),ximax=numpy.max(xi),timin=numpy.min(ti),timax=numpy.max(ti),xjmin=numpy.min(xj),xjmax=numpy.max(xj))
+numpy.savez(filename,N=N,Jxy=Jxy,hz=hz,B=B,xi=xi,ti=ti,xj=xj,sxsxitj=sxsxitj)

szsz.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy
+
+##################
+# spin-1/2 chain #
+##################
+
+# chain length
+N = 40
+
+# <Sz(xi,ti)Sz(xj,0)>
+xi = numpy.array([20,21])
+ti = numpy.array([-0.2, 0., 0.2, 0.4])
+xj = numpy.array([20])
+# note that the code is not threaded thus slow. Not a problem, because
+# it is a script and we want fast solution rather than efficient
+# coding at this point
+
+# ----------------------------------------------------------
+
+Jxy = 1. # energy scale
+hz = 0.8 # external magnetic field
+T = 0.002 # temperature
+B = 1./T # 1/T
+
+dk = numpy.pi/(N+1.);
+k = dk*numpy.linspace(1,N,N); # numpy.linspace(dk,N*dk,N);
+
+# ----------------------------------------------------------
+
+
+# dispersion
+E = lambda kn: Jxy*numpy.cos(kn) - hz;
+
+
+AitAj = lambda i,t,j: (2./(N+1.))*( numpy.sum( \
+        ( numpy.sin(k*j) * numpy.sin(k*i) ) * \
+        ( +(1.+0.j) * numpy.cos(E(k)*t) - (0.+1.j) * numpy.sin(E(k)*t) * numpy.tanh(B*E(k)/2.) ) \
+                  )     )
+AitBj = lambda i,t,j: (2./(N+1.))*( numpy.sum( \
+        ( numpy.sin(k*j) * numpy.sin(k*i) ) * \
+        ( -(0.+1.j) * numpy.sin(E(k)*t) + (1.+0.j) * numpy.cos(E(k)*t) * numpy.tanh(B*E(k)/2.) ) \
+                  )     )
+BitBj = lambda i,t,j: -AitAj(i,t,j)
+BitAj = lambda i,t,j: -AitBj(i,t,j)
+
+
+# ----------------------------------------------------------