You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
215 lines
17 KiB
215 lines
17 KiB
"""
|
|
unittest for Gram-Schmidt orthogonalization
|
|
|
|
|
|
"""
|
|
|
|
import unittest
|
|
import numpy
|
|
from numpy import array
|
|
np = numpy
|
|
|
|
from wpylib.math.linalg.gram_schmidt import \
|
|
modgs_classic, modgs_fast1, modgs_einsum
|
|
|
|
VERB = 0
|
|
|
|
def Msg(v, msg, *args):
|
|
if VERB >= v:
|
|
print(msg % (args))
|
|
|
|
# Walker 1, alpha sector
|
|
# from file: $GAFQMC/_scrap/test-pyqmc-gafqmc-wlkrs/test01-H2O/rundir1/W_old_run_00000.b0000000005
|
|
# 20150401
|
|
# This one is already orthogonalized
|
|
|
|
AA = \
|
|
array([[ -9.95209228822366981e-02 +4.86671674114909497e-02j, 1.37527009658723297e-02 -1.64322024633423713e-02j, -1.09285124415786812e-02 -4.04641740207571862e-02j, -5.72395447402896350e-03 +5.92599001505337078e-02j, -3.41882963186752311e-03 +3.80981135340452137e-02j],
|
|
[ -7.99330106673113057e-03 -1.20518717234826406e-02j, -8.87492552990528163e-03 +1.72008151178296223e-02j, -6.16159371590404951e-03 +2.09006658805644847e-02j, 2.57068370106492937e-02 -5.06063536781313425e-02j, -7.02797499966956868e-03 -4.23986334680305152e-03j],
|
|
[ -8.60327697198579461e-02 +4.43995763199260401e-02j, -5.43164961056228496e-02 -1.84776851908610533e-02j, 7.80074304097085839e-03 -3.54469899880848552e-02j, -1.76553144832910245e-03 +4.37672272679849594e-02j, 4.62903826016811939e-04 +1.67695032490165104e-02j],
|
|
[ 1.91352887657610937e-02 +3.13442349189221522e-02j, 1.95968816588026444e-02 -4.18043567496832261e-02j, 1.68881629058707690e-02 -1.07635701115742256e-02j, 5.30703802622154858e-02 +8.00220724834391495e-02j, 1.44963777530792044e-02 +1.81593715990712029e-02j],
|
|
[ 1.06904956968698773e-02 -2.12013200511964664e-02j, 1.03749543068047626e-02 +1.53545460523234114e-02j, 5.97814055203535824e-02 +1.27621766600792580e-01j, 7.94812333899607738e-03 -1.92026409769053170e-02j, 5.45038867876157258e-03 -5.02825685665939728e-02j],
|
|
[ 2.48422605809426635e-02 +5.21840592961241231e-02j, -4.72127941479885335e-03 -3.44775219639847849e-02j, 1.94991854284046855e-02 -4.82241195395200078e-02j, -3.53770367411777378e-02 +3.64749093647356165e-02j, 8.90636958127065115e-02 -5.25720694495861035e-02j],
|
|
[ 7.81121506708699984e-03 +1.93009515392476051e-02j, -1.33698552678669976e-02 +3.89448875891383259e-02j, 4.02960279586269121e-02 -3.85277118366297161e-02j, 2.41620662546525339e-01 +1.58986170129235160e-01j, -5.72953229626543652e-03 +4.87180029891744332e-02j],
|
|
[ 1.68192808817195422e-01 -7.51925866264861464e-02j, 9.39468460715941339e-02 +1.06042950621662538e-01j, -1.46229839640162673e-01 -1.26528721151196094e-01j, -1.41880787438075865e-02 +7.81319539429666116e-02j, -1.32877717162429645e-02 -7.21156622637989485e-03j],
|
|
[ -1.83566347139956162e-02 +5.93768751112509936e-03j, 3.58234494797938854e-02 -1.41983055278335130e-01j, 6.84380065166881602e-02 +2.59610136475489292e-02j, -8.53425973710577983e-03 +2.13806696144025338e-02j, 1.56617291707764950e-03 -3.94124827863742977e-02j],
|
|
[ 1.92885776408671791e-03 +4.22283941152761487e-03j, -1.36056970268698114e-03 +8.02821426980370595e-02j, -5.25235074668809095e-04 -6.41254886694271869e-03j, 1.15057562435817336e-02 -3.66991284027769690e-02j, 5.94927628107705798e-04 -1.28246189219457921e-03j],
|
|
[ -2.46284146137509341e-03 -9.58181299132664782e-03j, 2.57055226956389229e-02 -5.71867931529758608e-02j, -3.08989112719658711e-02 +5.53164410626623623e-02j, 6.65266425751886731e-03 -2.70592543440302778e-02j, -6.70201333150303724e-02 -7.16841512197503722e-02j],
|
|
[ 3.69790614682359185e-03 +5.40642032911758767e-03j, 5.97390752568134707e-03 +4.35810083714698471e-03j, -2.41340573086007426e-03 -1.15689738728525018e-02j, 7.52645866701945471e-02 +2.50083947791036425e-02j, 7.09391339921847715e-03 +1.60989723269885053e-02j],
|
|
[ -3.02849303345961848e-03 -1.93031759497581923e-03j, -2.87854508242068215e-02 +5.56890839954279487e-02j, 2.57527789555556073e-02 -1.73670066471564265e-02j, -6.44053868079939688e-03 -3.68364954239200418e-02j, -5.25948922498684565e-02 +1.16130170658301030e-01j],
|
|
[ 3.80721767786398602e-02 -1.48627273926002206e-02j, 4.44815622226979804e-02 -1.80794386092086007e-01j, -6.25426904861647731e-02 -6.02299108605501912e-02j, -2.28222684631696794e-02 -4.37436157161460396e-03j, 1.42789074971693720e-03 +1.67788114866421548e-02j],
|
|
[ 1.95564704058333056e-03 +4.07408151445507166e-03j, -4.77598959544649291e-03 +6.56503893651948084e-02j, 2.59001016226043362e-03 +1.97207772384004783e-02j, 8.67516721745117958e-03 -6.65350695855909668e-02j, 2.85032909316812553e-03 +5.34835212316238894e-03j],
|
|
[ 1.52719721613487819e-01 -6.54545613922360298e-02j, 1.39447115843387669e-02 -4.47694623987207954e-02j, -9.96499866875132367e-02 -7.60225150092632973e-02j, -3.66885688674354130e-03 -4.54640029506012261e-02j, 4.77605181328396005e-03 +9.62701892675258825e-02j],
|
|
[ -1.13875463782783046e-02 -1.24658785761800044e-02j, -3.21550856288821130e-02 +2.65138655099498646e-02j, 2.45651112836127727e-02 -1.35709015294978082e-02j, -1.89717034147060648e-01 +5.32565488020831090e-02j, -3.30091731334601485e-02 -4.41018875091828619e-02j],
|
|
[ 8.49957556161191730e-01 -3.95788224702817271e-01j, 3.00191570945964673e-02 -2.43915619158298266e-04j, 4.18433500608083805e-02 -2.55362907966733606e-03j, -2.15443899642035197e-02 +4.39419955528195863e-02j, -1.90138228123977872e-03 +4.29114093702174307e-02j],
|
|
[ 1.38915979097752815e-02 +2.32683632220901694e-02j, 2.60615584143406194e-03 +4.11362504196044823e-02j, 2.71286187431638397e-02 -3.33158350265078465e-02j, 5.46845668312948008e-01 +1.11124304327869183e-01j, 1.39120897695086451e-02 +7.80247554975158908e-02j],
|
|
[ 1.32162645958668779e-01 -6.85494769394368181e-02j, -2.87063333323404057e-01 -1.27778043884216180e-01j, 4.61932198625061130e-01 +2.92425181002827606e-01j, 2.16668566441828546e-02 -4.99738082504669542e-02j, 1.48390565017132402e-02 -4.34001603159832675e-02j],
|
|
[ 1.74246986372213683e-02 +3.68928267363070883e-02j, -1.37737876798530567e-02 -6.81680986948177570e-03j, 2.84654974587083137e-02 -6.72351600683199213e-02j, -2.74867621016573541e-02 +1.04901216005642894e-01j, 9.56280547122681512e-01 -1.02808878910921442e-01j],
|
|
[ -5.27362415933666090e-02 +1.60582287393329880e-02j, 6.23796461832417992e-01 -1.33455293008663007e-02j, 6.50520697878838927e-01 +1.80640486852616744e-01j, 9.34581267217418632e-03 +1.53807256155453435e-03j, 3.22001822862364546e-03 -4.92179769274162446e-02j],
|
|
[ 8.20254599109718902e-03 +1.34384609529303815e-02j, -1.20484743743531145e-02 +5.68518274116977113e-02j, 3.16427925172714716e-02 -3.95066574613382876e-02j, 6.96598702799640646e-01 +1.10588947910059049e-01j, 7.06059907314447559e-03 +5.81021442176028677e-02j],
|
|
[ -3.49146641949313330e-02 +1.54237222527656131e-02j, -6.16346199195223154e-01 -1.32491257800062351e-01j, 2.96421049204613207e-01 +2.09686320218193417e-01j, 3.41699081552355802e-02 -7.37901017809102638e-02j, 1.18803434309702020e-02 -9.56197919694180311e-03j]])
|
|
|
|
|
|
# Walker 1, alpha sector
|
|
# from file: $GAFQMC/_scrap/test-pyqmc-gafqmc-wlkrs/test01-H2O/rundir2/W_old_run_00000.b0000000005
|
|
# 20150401
|
|
# This one is NOT already orthogonalized
|
|
# AA and BB should result in the same determinant although the detail looks different
|
|
|
|
BB = \
|
|
array([[ -2.59864094736169768e+01 +1.27077292256682153e+01j, -1.93819339245309280e-01 -6.57012599404312780e-01j, 4.95856194837972275e-01 +7.22919527224234026e-01j, -4.85931418409559016e-02 -3.52016600892049336e-01j, 6.55129460126095853e-01 +6.50645351054976007e-01j],
|
|
[ -2.08717110483155643e+00 -3.14692493256434513e+00j, 4.20762460909348754e-02 +1.91445049873742376e-02j, -1.18887860587917504e-01 +1.50342742251630551e-01j, 1.48913741307116904e-01 -1.72199885210744141e-01j, -7.87703149050475726e-02 +8.08730319486135329e-02j],
|
|
[ -2.24644498598059315e+01 +1.15933970193381466e+01j, -4.14774134288470098e-01 -5.84002106654007092e-01j, 5.04727711311872440e-01 +6.28171252793613633e-01j, -4.50752274825731783e-02 -3.35225260056960128e-01j, 5.98125154145126570e-01 +5.04074547042267174e-01j],
|
|
[ 4.99651163656687114e+00 +8.18445106466889527e+00j, -1.21052868850918427e-01 -4.90658249555105674e-02j, 3.13833071832078592e-01 -2.37618053295939435e-01j, -5.43423324575270910e-03 +2.78576676355485553e-01j, 1.98181292398981551e-01 -1.77055917480609615e-01j],
|
|
[ 2.79144918082733184e+00 -5.53598347237528898e+00j, 1.55594626749259973e-01 +1.25090537401377905e-01j, -1.04152913483096313e-02 +3.26565271756842579e-01j, 1.01310804293838916e-01 +3.22527512204784927e-02j, -1.67324899507290931e-01 -1.61978255342079241e-01j],
|
|
[ 6.48668779398844553e+00 +1.36260425807067289e+01j, -3.30481270623765144e-01 +1.45537615170570115e-03j, 5.05862343225446875e-01 -4.24930784435921749e-01j, -3.90649342463602500e-01 +1.25656198930886254e-01j, 5.56917700733675747e-01 -5.04389347686009404e-01j],
|
|
[ 2.03962571227360945e+00 +5.03976867781688309e+00j, -1.62694554314163298e-01 +1.74993753722288453e-01j, 2.98718767553879883e-01 -2.09775323342454412e-01j, 6.61773678253472886e-01 +5.02039054824823028e-01j, 8.89074546998883247e-02 +2.88743298679756538e-02j],
|
|
[ 4.39176715193176364e+01 -1.96339150488831180e+01j, 7.00131212837735051e-01 +1.38374103020257588e+00j, -1.28969000981536430e+00 -1.85914201352741393e+00j, -2.63334573028017119e-02 +1.14792849855204659e+00j, -1.12616096623176176e+00 -9.49352942206428851e-01j],
|
|
[ -4.79319335492863186e+00 +1.55041949493435860e+00j, 9.94200864303528414e-02 -6.10381022658728511e-01j, 2.72530099139980275e-01 +2.40458335441069565e-01j, -9.47117026609448943e-03 -3.05675633129760677e-02j, 1.05684696185104898e-01 -2.55452994832732708e-03j],
|
|
[ 5.03653766688183402e-01 +1.10264686973547921e+00j, -3.01635562346698784e-02 +2.92852809425142524e-01j, 4.29371367396094022e-02 -4.25719095060025920e-02j, 1.03589140406833598e-02 -1.11468053036431688e-01j, 2.35232408329638411e-02 -3.72466186229807594e-02j],
|
|
[ -6.43084939632505570e-01 -2.50195545500392003e+00j, 1.47490676071875326e-01 -2.12245751240811409e-01j, -1.88299209073673784e-01 +2.04574746407425706e-01j, 7.12382282251780852e-02 -7.89775532849481543e-02j, -2.66191618003197628e-01 -1.68176923204563178e-01j],
|
|
[ 9.65578900831461917e-01 +1.41169764497846240e+00j, -1.18880117260648215e-02 +3.46814902011944801e-02j, 3.80415324079927883e-02 -7.57653787499782821e-02j, 2.06892103392281757e-01 +8.56575482596485804e-02j, 4.62629246512387404e-02 +6.00475885255571072e-03j],
|
|
[ -7.90785070880053231e-01 -5.04034950484995137e-01j, -8.91269551008052607e-02 +1.79728803036626306e-01j, 7.31498702189592870e-02 -2.36439904552020880e-02j, -9.01304441913901300e-03 -1.25230248009072148e-01j, -1.64831021587973997e-01 +3.81725786969631387e-01j],
|
|
[ 9.94121785317822848e+00 -3.88088161497335271e+00j, 2.27867824504319999e-01 -4.10510497391186502e-01j, -3.87209518797457186e-01 -5.22728316347241595e-01j, -6.85349834534964425e-02 +1.75535993657036732e-01j, -2.24094898819090366e-01 -1.48334256466603337e-01j],
|
|
[ 5.10648849615284450e-01 +1.06380394591798466e+00j, -4.13639590404624863e-02 +2.41425895974034327e-01j, 4.99749351298805558e-02 +4.06815279160777643e-02j, 2.44927410428529479e-03 -2.05520448837672842e-01j, 2.99401388834034526e-02 -1.54503376679429646e-02j],
|
|
[ 3.98774157796041564e+01 -1.70911702282680267e+01j, 3.67442429849584085e-01 +7.57415697577259017e-01j, -1.05426362723247347e+00 -1.55840662328528157e+00j, -4.29489950024567347e-03 +6.60023511230309889e-01j, -9.43152847569946773e-01 -5.43595898604261163e-01j],
|
|
[ -2.97345959538486149e+00 -3.25502834727855550e+00j, -3.62920265558136640e-02 +3.24360364146693905e-02j, -2.08110858357219639e-02 +7.44017581162866315e-02j, -5.24574479839464147e-01 +1.40374883791183930e-01j, -1.48139452912094488e-01 -1.90471057223798590e-02j],
|
|
[ 2.21936698836037806e+02 -1.03346256988931998e+02j, 2.06031636264815043e+00 +5.10687224765584169e+00j, -4.23530718880099055e+00 -7.30934216538361703e+00j, 1.16664079983079247e-01 +4.69707984636614206e+00j, -5.56233117929216192e+00 -4.47119121207479875e+00j],
|
|
[ 3.62730510400737360e+00 +6.07571952669269244e+00j, -1.31883547178452465e-01 +2.16479989413323776e-01j, 2.87113046792904447e-01 -2.53516269166096042e-01j, 1.59801795501869948e+00 +3.71426414134579386e-01j, 1.54054080016761741e-01 +6.47785418208933761e-02j],
|
|
[ 3.45096542067098824e+01 -1.78992992162899007e+01j, -6.68549361992861302e-01 +3.46089695324148461e-01j, 7.11993272384795484e-01 -1.70285456379856426e-01j, 1.31352010818359177e-01 +5.55606960799409189e-01j, -8.46190822049378788e-01 -8.24996613711939109e-01j],
|
|
[ 4.54985083163883886e+00 +9.63327182308512597e+00j, -2.70419798675604561e-01 +6.21969907623121307e-02j, 4.06228845330372013e-01 -4.03341848532684677e-01j, -2.84449348196536278e-01 +3.37923152699492413e-01j, 3.15417583463721085e+00 -5.60171104618138127e-01j],
|
|
[ -1.37702256817537734e+01 +4.19304499351454574e+00j, 2.12992242034626100e+00 -3.60532203936021645e-01j, 2.24638066039787310e+00 +9.68984023150801232e-01j, 1.07075190287191202e-01 -2.28082116503400323e-01j, 2.40922520207871915e-01 +1.70768453870208420e-01j],
|
|
[ 2.14180810102668762e+00 +3.50898423069580545e+00j, -1.24080442061272236e-01 +2.42878330316865287e-01j, 2.19866401434362174e-01 -2.11281913675311350e-01j, 2.12508004510244408e+00 +3.62657638754537104e-01j, 8.31549289052572055e-02 +7.39399680258190917e-02j],
|
|
[ -9.11674383764443164e+00 +4.02736580871522154e+00j, -2.25086233030583838e+00 -6.77110956632504535e-01j, 1.10545726041556769e+00 +1.03040441466139199e+00j, 8.63289732775991830e-02 -4.44331269872425327e-01j, 3.06459774661179130e-01 +1.60097232824508551e-01j]])
|
|
|
|
|
|
|
|
def check_ortho(V):
|
|
"""Checks orthogonality of the column vectors in matrix V.
|
|
If V is not orthonormal, it will creak.
|
|
"""
|
|
from wpylib.math import ztol
|
|
from numpy.linalg import norm
|
|
prod = np.dot(V.T.conj(), V)
|
|
prod_ref = np.eye(prod.shape[0])
|
|
prod_diff = prod - prod_ref
|
|
NN = norm(prod_diff)
|
|
if NN > 1e-13:
|
|
Msg(1, " Failed check_ortho, norm of orthonormality deviation = %g", NN)
|
|
Msg(1, " prod =\n%s", ztol(prod, 1e-14))
|
|
s = 0
|
|
else:
|
|
Msg(10, " Succeeded check_ortho, norm of orthonormality deviation = %g", NN)
|
|
s = 1
|
|
return NN, ztol(prod, 1e-14)
|
|
|
|
|
|
def ovlp_det(A, B):
|
|
"""If two orthonormal det matrices are identical, then the det
|
|
of the overlap should be one.
|
|
Example:
|
|
ovlp_det(AA, BB3) from above should result in det() that is 1.
|
|
"""
|
|
ovlp = np.dot(A.T.conj(), B)
|
|
Msg(5, "ovlp = \n%s", ovlp)
|
|
mydet = np.linalg.det(ovlp)
|
|
Msg(5, "det = %s", mydet)
|
|
return mydet
|
|
|
|
def Real(Z):
|
|
try:
|
|
return Z.real
|
|
except:
|
|
return Z
|
|
|
|
def Imag(Z):
|
|
try:
|
|
return Z.imag
|
|
except:
|
|
return 0
|
|
|
|
|
|
class ModGS_Tests(unittest.TestCase):
|
|
|
|
def test_modgs_classic_AA(self):
|
|
"""Testing classic modgs on AA.
|
|
Result: Passed.
|
|
"""
|
|
global AA, AA2
|
|
Msg(1, "test_modgs_classic_AA:")
|
|
AA2, AA2n = modgs_classic(AA)
|
|
orth_err, self_prod = check_ortho(AA2)
|
|
self.assertAlmostEqual(orth_err, 0.0, 14)
|
|
det = ovlp_det(AA, AA2)
|
|
self.assertAlmostEqual(Real(det), 1.0, 14)
|
|
self.assertAlmostEqual(Imag(det), 0.0, 14)
|
|
|
|
def test_modgs_classic_BB(self):
|
|
"""Testing classic modgs on BB.
|
|
Result: Passed.
|
|
"""
|
|
global BB, BB2
|
|
Msg(1, "test_modgs_classic_BB:")
|
|
BB2, BB2n = modgs_classic(BB)
|
|
orth_err, self_prod = check_ortho(BB2)
|
|
self.assertAlmostEqual(orth_err, 0.0, 14)
|
|
det = ovlp_det(AA, BB2)
|
|
self.assertAlmostEqual(Real(det), 1.0, 14)
|
|
self.assertAlmostEqual(Imag(det), 0.0, 14)
|
|
|
|
def test_modgs_fast1_AA(self):
|
|
"""Testing fast1 modgs on AA.
|
|
Result: Passed.
|
|
"""
|
|
global AA, AA3
|
|
Msg(1, "test_modgs_fast1_AA:")
|
|
AA3, AA3n = modgs_fast1(AA)
|
|
orth_err, self_prod = check_ortho(AA3)
|
|
self.assertAlmostEqual(orth_err, 0.0, 14)
|
|
det = ovlp_det(AA, AA3)
|
|
self.assertAlmostEqual(Real(det), 1.0, 14)
|
|
self.assertAlmostEqual(Imag(det), 0.0, 14)
|
|
|
|
def test_modgs_fast1_BB(self):
|
|
"""Testing fast1 modgs on BB.
|
|
Result: Passed.
|
|
"""
|
|
global BB, BB3
|
|
Msg(1, "test_modgs_fast1_BB:")
|
|
BB3, BB3n = modgs_fast1(BB)
|
|
orth_err, self_prod = check_ortho(BB3)
|
|
self.assertAlmostEqual(orth_err, 0.0, 14)
|
|
det = ovlp_det(AA, BB3)
|
|
self.assertAlmostEqual(Real(det), 1.0, 14)
|
|
self.assertAlmostEqual(Imag(det), 0.0, 14)
|
|
|
|
def test_modgs_einsum_AA(self):
|
|
"""Testing einsum modgs on AA.
|
|
Result: Passed.
|
|
"""
|
|
global AA, AA4
|
|
Msg(1, "test_modgs_einsum_AA:")
|
|
AA4, AA4n = modgs_einsum(AA)
|
|
orth_err, self_prod = check_ortho(AA4)
|
|
self.assertAlmostEqual(orth_err, 0.0, 14)
|
|
det = ovlp_det(AA, AA4)
|
|
self.assertAlmostEqual(Real(det), 1.0, 14)
|
|
self.assertAlmostEqual(Imag(det), 0.0, 14)
|
|
|
|
def test_modgs_einsum_BB(self):
|
|
"""Testing einsum modgs on BB.
|
|
Result: Passed.
|
|
"""
|
|
global BB, BB4
|
|
Msg(1, "test_modgs_einsum_BB:")
|
|
BB4, BB4n = modgs_einsum(BB)
|
|
orth_err, self_prod = check_ortho(BB4)
|
|
self.assertAlmostEqual(orth_err, 0.0, 14)
|
|
det = ovlp_det(AA, BB4)
|
|
self.assertAlmostEqual(Real(det), 1.0, 14)
|
|
self.assertAlmostEqual(Imag(det), 0.0, 14)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
print("Tests for Gram-Schmidt orthogonalization")
|
|
unittest.main()
|
|
|