hydra log

.gitignore

@@ -7,4 +7,5 @@ jobs/*
 *.prof
 *.h5
 __pycache__
-dist
+dist
+outputs

Makefile

@@ -0,0 +1,2 @@
+clear:
+	rm -rf output/*

conf/gfpeps/default.yaml

@@ -1,25 +1,26 @@
 params:
-  Nv: 2 # Number of virtual fermions on each bond
-  seed: 42 # Random seed
+  Nv: 1 # Number of virtual fermions on each bond
+  seed: 44 # Random seed
 
 lattice:
   Lx: 101 # System size in x-direction
   Ly: 101 # System size in y-direction
 
 hamiltonian:
   ht: 1.0 # Hoping amplitude
-  DeltaX: 1.0 # d-wave paring amplitude in x-direction
-  DeltaY: -1.0 # d-wave paring amplitude in y-direction
+  DeltaX: 0.0 # d-wave paring amplitude in x-direction
+  DeltaY: 0.0 # d-wave paring amplitude in y-direction
   delta: 0.0 # desity of holes
   Mu: 0.0 # chemical potential
   solve_mu_from_delta: false # whether to solve chemical potential from hole density
 
 file:
-  LoadFile: "./data/default.h5" # Load initial T from this file, if it exists
+  # LoadFile: "./data/default.h5" # Load initial T from this file, if it exists
+  LoadFile: ""
   WriteFile: "./data/default.h5" # Write information about final Gaussian fPEPS to this file
 
 optimizer:
   MaxIter: 100
   gtol: 1E-7 # gtol for optimizer
 
-backend: gpu
+backend: cpu

gfpeps_app.py

@@ -0,0 +1,15 @@
+from gfpeps import gaussian_fpeps
+
+# hydra config
+import hydra
+from omegaconf import DictConfig
+
+import logging
+log = logging.getLogger(__name__)
+
+@hydra.main(version_base=None, config_path="conf/gfpeps", config_name="default")
+def main_app(cfg: DictConfig) -> None:
+    return gaussian_fpeps(cfg)
+
+if __name__ == '__main__':
+    main_app()

src/gfpeps/GaussianLinearMap.py

@@ -1,6 +1,6 @@
 import jax.numpy as jnp
 from jax.scipy.linalg import inv
-from ABD import getABD
+from .ABD import getABD
 
 def GaussianLinearMap(Glocal,Gin):
     r""" Glocal: Local project, Gin, Gamma in"""

src/gfpeps/__init__.py

@@ -2,4 +2,5 @@
 #
 # SPDX-License-Identifier: MIT
 
-from main import gaussian_fpeps
+__all__ = ['gaussian_fpeps']
+from .gaussian_fpeps import gaussian_fpeps

src/gfpeps/deltatomu.py

@@ -2,7 +2,7 @@
 import numpy as np
 import jax.numpy as jnp
 import jax
-from jax.config import config
+from jax import config
 config.update("jax_enable_x64", True)
 from scipy.optimize import root_scalar
 

src/gfpeps/exact.py

@@ -1,7 +1,7 @@
 import jax.numpy as jnp
 from jax.scipy.linalg import eigh
 import jax
-from Gin import BatchK
+from .Gin import BatchK
 
 def exact(k,ht,D1X,D1Y,Mu):
     D = D1X*jnp.cos(k[0])+D1Y*jnp.cos(k[1])

src/gfpeps/gaussian_fpeps.py

@@ -0,0 +1,87 @@
+import jax
+jax.config.update("jax_enable_x64", True)
+import numpy as np
+import jax.numpy as jnp
+from jax import jit
+from .loss import optimize_runtime_loss
+from .loadwrite import initialT,savelog,savelog_trivial
+from .exact import eg
+from .deltatomu import solve_mu
+import logging
+
+# import Manopt
+import pymanopt
+from pymanopt.manifolds import Stiefel
+from pymanopt import Problem
+from pymanopt.optimizers import TrustRegions,ConjugateGradient
+
+def gaussian_fpeps(cfg):
+    # unpack cfg
+    np.random.seed(cfg.params.seed)
+    Nv = cfg.params.Nv
+    Lx, Ly = cfg.lattice.Lx, cfg.lattice.Ly
+    LoadKey, WriteKey  = cfg.file.LoadFile, cfg.file.WriteFile
+
+    cfgh = cfg.hamiltonian
+    ht = cfgh.ht
+    DeltaX, DeltaY = cfgh.DeltaX, cfgh.DeltaY
+    delta, Mu = cfgh.delta, cfgh.Mu
+
+    Tsize = 8*Nv+4
+    T = initialT(LoadKey,Tsize)
+    U,S,V = np.linalg.svd(T)
+    T = U @ V
+
+    if cfgh.solve_mu_from_delta:
+        logging.info("Overwrite Origin Mu")
+        Mu = solve_mu(DeltaX,delta)
+
+    lossT = jit(optimize_runtime_loss(Nv=Nv,Lx=Lx,Ly=Ly,
+    hoping=ht,DeltaX=DeltaX,DeltaY=DeltaY,Mu=Mu), backend=cfg.backend)
+
+    def egrad(x): return np.array(jax.grad(lossT)(jnp.array(x)))
+
+    @jax.jit
+    def hvp_loss(primals, tangents): return jax.jvp(jax.grad(lossT), primals, tangents)[1]
+
+    def ehessa(x,v): return np.array(hvp_loss((jnp.array(x),), (jnp.array(v),)))
+
+    Eg = eg(Lx,Ly,ht,DeltaX,DeltaY,Mu) # Will use solved Mu to calculate Eg
+    logging.info("Eg = {}\n".format(Eg))
+    # Optimizer
+
+    manifold = Stiefel(Tsize, Tsize)
+    @pymanopt.function.numpy(manifold)
+    def cost(x):
+        return lossT(x)
+
+    @pymanopt.function.numpy(manifold)
+    def euclidean_gradient(x):
+        return egrad(x)
+
+    @pymanopt.function.numpy(manifold)
+    def euclidean_hessian(x,y):
+        return ehessa(x,y)
+
+
+    problem = Problem(manifold=manifold, 
+                      cost=cost,
+                      euclidean_gradient=euclidean_gradient,
+                      euclidean_hessian=euclidean_hessian)
+    solver = ConjugateGradient(log_verbosity=1, max_iterations=cfg.optimizer.MaxIter)
+
+    result = solver.run(problem, initial_point=T)
+    log_cost = np.array(result.log["iterations"]["cost"])
+    log_gnorm = np.array(result.log["iterations"]["gradient_norm"])
+
+    logging.info("Iterations \t Cost \t Gradient Norm")
+    for iter in range(len(log_cost)):
+        logging.info(f"{iter} \t {log_cost[iter]} \t {log_gnorm[iter]}")
+
+    logging.info(f"Optimization done!, final cost: {result.cost}, gnorm: {result.gradient_norm }")
+
+    Xopt = result.point
+    args = {"Mu":Mu,"DeltaX":DeltaX,"DeltaY":DeltaY,"delta":delta,
+            "ht":ht,"Lx":Lx,"Ly":Ly,"Nv":Nv,"seed":cfg.params.seed}
+    savelog_trivial(WriteKey,Xopt,lossT(Xopt),Eg,args)
+    return Xopt

src/gfpeps/loadwrite.py

@@ -1,24 +1,25 @@
 import numpy as np
 import jax.numpy as jnp
 import h5py,os
+import logging
 
 def initialT(loadfile,Tsize):
     if loadfile != None and os.path.isfile(loadfile):
-        print('Try to initialize T from',loadfile)
+        logging.info(f'Try to initialize T from {loadfile}')
         with h5py.File(loadfile, 'r') as f:
             if "/transformer/T" in f.keys():
                 T = f["/transformer/T"][:]
                 f.close()
                 return jnp.reshape(T,(Tsize,Tsize))
             else:
-                print("Load Failed! No /transformer/T in",loadfile," switch to random initialize!")
+                logging.debug(f"Load Failed! No /transformer/T in {loadfile} switch to random initialize!")
                 f.close()
                 return jnp.array(np.random.rand(Tsize,Tsize))
     return jnp.array(np.random.rand(Tsize,Tsize))
 
 def savelog(writefile,results):
     if writefile != None:
-        print('Save T to',writefile)
+        logging.info(f'Save T to {writefile}')
         with h5py.File(writefile, 'w') as f:
             f["/transformer/T"] = results.x
             f["/energy/EABD"] = results.fun
@@ -30,7 +31,7 @@ def savelog(writefile,results):
 
 def savelog_trivial(writefile,x,fun,Eg,args):
     if writefile != None:
-        print('Save T to',writefile)
+        logging.info(f'Save T to {writefile}')
         with h5py.File(writefile, 'w') as f:
             f["/transformer/T"] = x
             f["/energy/EABD"] = fun

src/gfpeps/loss.py

@@ -1,8 +1,8 @@
 from jax import vmap
 import jax.numpy as jnp
-from ABD import getGammaProjector
-from Gin import BatchGammaIn,BatchK
-from GaussianLinearMap import GaussianLinearMap
+from .ABD import getGammaProjector
+from .Gin import BatchGammaIn,BatchK
+from .GaussianLinearMap import GaussianLinearMap
 
 def energy_function(hoping=1.0,DeltaX=0.0,DeltaY=0.0,Mu=0.0,Lx=100,Ly=100):
     batch_k = BatchK(Lx,Ly)