Migrated sim algorithm.

diff --git a/.gitmodules b/.gitmodules
new file mode 100644 (file)
index 0000000..6386119
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "latlib"]
+       path = latlib
+       url = treefish.org:~/public_git/phys/latlib.git

diff --git a/latlib b/latlib
new file mode 160000 (submodule)
index 0000000..c08c4b7
--- /dev/null
+++ b/latlib
@@ -0,0 +1 @@
+Subproject commit c08c4b78cd61f072d566aefaa872c5348ab129a0

diff --git a/u1casc-ordinary/.gitignore b/u1casc-ordinary/.gitignore
new file mode 100644 (file)
index 0000000..332a25a
--- /dev/null
+++ b/u1casc-ordinary/.gitignore
@@ -0,0 +1,6 @@
+CMakeFiles
+*~
+Makefile
+cmake_install.cmake
+CMakeCache.txt
+u1casc-ordinary

diff --git a/u1casc-ordinary/CMakeLists.txt b/u1casc-ordinary/CMakeLists.txt
new file mode 100644 (file)
index 0000000..40225bd
--- /dev/null
+++ b/u1casc-ordinary/CMakeLists.txt
@@ -0,0 +1,20 @@
+cmake_minimum_required(VERSION 2.8)
+
+project(u1casc)
+
+# #PROFILING
+# set(CMAKE_CXX_FLAGS -pg)
+# set(CMAKE_EXE_LINKER_FLAGS -pg)
+# set(CMAKE_SHARED_LINKER_FLAGS -pg)
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "./")
+
+FIND_PACKAGE(GSL REQUIRED)
+include_directories(${GSL_INCLUDE_DIRS})
+
+SET(CMAKE_BUILD_TYPE Release)
+
+add_subdirectory(latlib)
+
+add_executable(u1casc-ordinary u1casc-ordinary.cpp)
+target_link_libraries(u1casc-ordinary o815 ${GSL_LIBRARIES} lat_neigh)

diff --git a/u1casc-ordinary/FindGSL.cmake b/u1casc-ordinary/FindGSL.cmake
new file mode 100644 (file)
index 0000000..486be73
--- /dev/null
+++ b/u1casc-ordinary/FindGSL.cmake
@@ -0,0 +1,135 @@
+# Try to find gnu scientific library GSL
+# See
+# http://www.gnu.org/software/gsl/  and
+# http://gnuwin32.sourceforge.net/packages/gsl.htm
+#
+# Based on a script of Felix Woelk and Jan Woetzel
+# (www.mip.informatik.uni-kiel.de)
+#
+# It defines the following variables:
+#  GSL_FOUND - system has GSL lib
+#  GSL_INCLUDE_DIRS - where to find headers
+#  GSL_LIBRARIES - full path to the libraries
+#  GSL_LIBRARY_DIRS, the directory where the PLplot library is found.
+ 
+#  CMAKE_GSL_CXX_FLAGS  = Unix compiler flags for GSL, essentially "`gsl-config --cxxflags`"
+#  GSL_LINK_DIRECTORIES = link directories, useful for rpath on Unix
+#  GSL_EXE_LINKER_FLAGS = rpath on Unix
+ 
+set( GSL_FOUND OFF )
+set( GSL_CBLAS_FOUND OFF )
+ 
+# Windows, but not for Cygwin and MSys where gsl-config is available
+if( WIN32 AND NOT CYGWIN AND NOT MSYS )
+  # look for headers
+  find_path( GSL_INCLUDE_DIR
+    NAMES gsl/gsl_cdf.h gsl/gsl_randist.h
+    )
+  if( GSL_INCLUDE_DIR )
+    # look for gsl library
+    find_library( GSL_LIBRARY
+      NAMES gsl
+    ) 
+    if( GSL_LIBRARY )
+      set( GSL_INCLUDE_DIRS ${GSL_INCLUDE_DIR} )
+      get_filename_component( GSL_LIBRARY_DIRS ${GSL_LIBRARY} PATH )
+      set( GSL_FOUND ON )
+    endif( GSL_LIBRARY )
+ 
+    # look for gsl cblas library
+    find_library( GSL_CBLAS_LIBRARY
+        NAMES gslcblas
+      )
+    if( GSL_CBLAS_LIBRARY )
+      set( GSL_CBLAS_FOUND ON )
+    endif( GSL_CBLAS_LIBRARY )
+       
+    set( GSL_LIBRARIES ${GSL_LIBRARY} ${GSL_CBLAS_LIBRARY} )
+  endif( GSL_INCLUDE_DIR )
+   
+  mark_as_advanced(
+    GSL_INCLUDE_DIR
+    GSL_LIBRARY
+    GSL_CBLAS_LIBRARY
+  )
+else( WIN32 AND NOT CYGWIN AND NOT MSYS )
+  if( UNIX OR MSYS )
+    find_program( GSL_CONFIG_EXECUTABLE gsl-config
+      /usr/bin/
+      /usr/local/bin
+    )
+     
+    if( GSL_CONFIG_EXECUTABLE )
+      set( GSL_FOUND ON )
+       
+      # run the gsl-config program to get cxxflags
+      execute_process(
+        COMMAND sh "${GSL_CONFIG_EXECUTABLE}" --cflags
+        OUTPUT_VARIABLE GSL_CFLAGS
+        RESULT_VARIABLE RET
+        ERROR_QUIET
+        )
+      if( RET EQUAL 0 )
+        string( STRIP "${GSL_CFLAGS}" GSL_CFLAGS )
+        separate_arguments( GSL_CFLAGS )
+ 
+        # parse definitions from cflags; drop -D* from CFLAGS
+        string( REGEX MATCHALL "-D[^;]+"
+          GSL_DEFINITIONS  "${GSL_CFLAGS}" )
+        string( REGEX REPLACE "-D[^;]+;" ""
+          GSL_CFLAGS "${GSL_CFLAGS}" )
+ 
+        # parse include dirs from cflags; drop -I prefix
+        string( REGEX MATCHALL "-I[^;]+"
+          GSL_INCLUDE_DIRS "${GSL_CFLAGS}" )
+        string( REPLACE "-I" ""
+          GSL_INCLUDE_DIRS "${GSL_INCLUDE_DIRS}")
+        string( REGEX REPLACE "-I[^;]+;" ""
+          GSL_CFLAGS "${GSL_CFLAGS}")
+ 
+        message("GSL_DEFINITIONS=${GSL_DEFINITIONS}")
+        message("GSL_INCLUDE_DIRS=${GSL_INCLUDE_DIRS}")
+        message("GSL_CFLAGS=${GSL_CFLAGS}")
+      else( RET EQUAL 0 )
+        set( GSL_FOUND FALSE )
+      endif( RET EQUAL 0 )
+ 
+      # run the gsl-config program to get the libs
+      execute_process(
+        COMMAND sh "${GSL_CONFIG_EXECUTABLE}" --libs
+        OUTPUT_VARIABLE GSL_LIBRARIES
+        RESULT_VARIABLE RET
+        ERROR_QUIET
+        )
+      if( RET EQUAL 0 )
+        string(STRIP "${GSL_LIBRARIES}" GSL_LIBRARIES )
+        separate_arguments( GSL_LIBRARIES )
+ 
+        # extract linkdirs (-L) for rpath (i.e., LINK_DIRECTORIES)
+        string( REGEX MATCHALL "-L[^;]+"
+          GSL_LIBRARY_DIRS "${GSL_LIBRARIES}" )
+        string( REPLACE "-L" ""
+          GSL_LIBRARY_DIRS "${GSL_LIBRARY_DIRS}" )
+      else( RET EQUAL 0 )
+        set( GSL_FOUND FALSE )
+      endif( RET EQUAL 0 )
+       
+      MARK_AS_ADVANCED(
+        GSL_CFLAGS
+      )
+      message( STATUS "Using GSL from ${GSL_PREFIX}" )
+    else( GSL_CONFIG_EXECUTABLE )
+      message( STATUS "FindGSL: gsl-config not found.")
+    endif( GSL_CONFIG_EXECUTABLE )
+  endif( UNIX OR MSYS )
+endif( WIN32 AND NOT CYGWIN AND NOT MSYS )
+ 
+if( GSL_FOUND )
+  if( NOT GSL_FIND_QUIETLY )
+    message( STATUS "FindGSL: Found both GSL headers and library" )
+  endif( NOT GSL_FIND_QUIETLY )
+else( GSL_FOUND )
+  if( GSL_FIND_REQUIRED )
+    message( FATAL_ERROR "FindGSL: Could not find GSL headers or library" )
+  endif( GSL_FIND_REQUIRED )
+endif( GSL_FOUND )

diff --git a/u1casc-ordinary/latlib b/u1casc-ordinary/latlib
new file mode 120000 (symlink)
index 0000000..7076779
--- /dev/null
+++ b/u1casc-ordinary/latlib
@@ -0,0 +1 @@
+../latlib
\ No newline at end of file

diff --git a/u1casc-ordinary/sim.hpp b/u1casc-ordinary/sim.hpp
new file mode 100644 (file)
index 0000000..af989c9
--- /dev/null
+++ b/u1casc-ordinary/sim.hpp
@@ -0,0 +1,141 @@
+#ifndef SIM_HPP
+#define SIM_HPP
+
+#include <gsl/gsl_rng.h>
+#include <complex>
+#include <math.h>
+
+#include "latlib/neigh.h"
+
+#define EPSILONU 1
+#define EPSILONPHI 0.5
+
+class sim : public o815::sim {
+public:
+  sim(o815 *_O815);
+  unsigned int lsize4;
+  neigh *nb;
+  complex<double> *U, *phi;
+  double kappa[2], lambda[2], beta;
+
+private:
+  void _makeSweep();
+  void _newParas();
+  gsl_rng* rangsl;
+  double rhoPhi(const int& iphi, const int& x0, const complex<double>& candPhi);
+  double rhoU(const int& x0, const int& nu0, const complex<double>& candU);
+  int updatePhi(const int& iphi, const int& x0);
+  int updateU(const int& x0, const int& nu0);
+};
+
+sim::sim(o815 *_O815) : o815::sim( _O815, 
+                                  sizeof(complex<double>)*
+                                  _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1]*(2+4) ) {
+
+  lsize4 = _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1];
+
+  nb = new neigh(4, _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[1]);
+
+  phi = (complex<double>*)confMem;
+  U = (complex<double>*)(confMem + sizeof(complex<double>)*lsize4*2);
+
+  rangsl = gsl_rng_alloc(gsl_rng_ranlxs0);
+  gsl_rng_set(rangsl, time(NULL));
+}
+
+void sim::_makeSweep() {  
+  for( int ix=0; ix<lsize4; ix++ ) {
+    for( int inu=0; inu<4; inu++) updateU(ix, inu);
+    for( int iphi=0; iphi<2; iphi++) updatePhi(iphi, ix);
+  }
+}
+
+void sim::_newParas() {
+  kappa[0] = (*O815->paraQ)["kappaone"];
+  kappa[1] = (*O815->paraQ)["kappatwo"];
+  lambda[0] = (*O815->paraQ)["lambdaone"];
+  lambda[1] = (*O815->paraQ)["lambdatwo"];
+  beta = (*O815->paraQ)["beta"];
+
+  for(int ix=0; ix<lsize4; ix++) {
+    for(int i=0; i<2; i++) phi[ i*lsize4 + ix ] = 0;
+    for(int nu=0; nu<4; nu++) U[ ix*4 + nu ] = 1;
+  }
+}
+
+int sim::updateU(const int& x0, const int& nu0)
+{
+  complex<double> candU = U[x0*4+nu0] * polar(1.0, 2*EPSILONU*( 0.5 - gsl_rng_uniform(rangsl) ));
+
+  if ( gsl_rng_uniform(rangsl) <= rhoU(x0, nu0, candU) ) {
+    U[x0*4 + nu0] = candU;
+    return 1;
+  }
+
+  return 0;
+}
+
+int sim::updatePhi(const int& iphi, const int& x0)
+{
+  complex<double> candPhi = phi[ iphi*lsize4 + x0 ] + 
+    complex<double>(2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ), 
+                   2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ));
+
+  if ( gsl_rng_uniform(rangsl) <= rhoPhi(iphi, x0, candPhi) ) {
+    phi[ iphi*lsize4 + x0 ] = candPhi;
+    return 1;
+  }
+
+  return 0;
+}
+
+double sim::rhoPhi(const int& iphi, const int& x0, const complex<double>& candPhi)
+{
+  double deltaS=0;
+
+  for( int mu=0; mu<4; mu++) {
+    if( iphi == 0 ) {
+      deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * U[ x0*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
+      deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * conj(U[ (*nb)[x0*8+mu+4]*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
+      deltaS -= 2 * real( conj(candPhi) * U[ x0*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
+      deltaS -= 2 * real( conj(candPhi) * conj(U[ (*nb)[x0*8+mu+4]*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
+    }
+    else if( iphi == 1 ) {
+      deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * conj(U[ x0*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
+      deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * U[ (*nb)[x0*8+mu+4]*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
+      deltaS -= 2 * real( conj(candPhi) * conj(U[ x0*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
+      deltaS -= 2 * real( conj(candPhi) * U[ (*nb)[x0*8+mu+4]*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
+    }
+  }
+
+  deltaS -= kappa[iphi] * norm(phi[ iphi*lsize4 + x0 ]);
+  deltaS += kappa[iphi] * norm(candPhi);
+
+  deltaS -= lambda[iphi] * pow(norm(phi[ iphi*lsize4 + x0 ]),2);
+  deltaS += lambda[iphi] * pow(norm(candPhi),2);
+
+  return exp(-deltaS);
+}
+
+double sim::rhoU(const int& x0, const int& nu0, const complex<double>& candU)
+{
+  double deltaS=0;
+
+  for( int nu=0; nu<4; nu++ ) {
+    if( nu == nu0 ) continue;
+    deltaS += beta * real( U[x0*4+nu0] * U[ (*nb)[x0*8+nu0]*4 + nu ] * conj(U[ (*nb)[x0*8+nu]*4 + nu0 ]) * conj(U[ x0*4 + nu ]) );
+    deltaS += beta * real( U[ (*nb)[x0*8+nu+4]*4 + nu0 ] * U[ (*nb)[ (*nb)[x0*8+nu+4]*8+nu0 ]*4 + nu ] * conj(U[ x0*4 + nu0 ]) * conj(U[ (*nb)[x0*8+nu+4]*4 + nu ]) );
+    deltaS -= beta * real( candU * U[ (*nb)[x0*8+nu0]*4 + nu ] * conj(U[ (*nb)[x0*8+nu]*4 + nu0 ]) * conj(U[ x0*4 + nu ]) );
+    deltaS -= beta * real( U[ (*nb)[x0*8+nu+4]*4 + nu0 ] * U[ (*nb)[ (*nb)[x0*8+nu+4]*8+nu0 ]*4 + nu ] * conj(candU) * conj(U[ (*nb)[x0*8+nu+4]*4 + nu ]) );
+  }
+
+  deltaS += 2 * real( conj(phi[ 0*lsize4 + x0 ]) * U[ x0*4 + nu0 ] * phi[ 0*lsize4 + (*nb)[x0*8+nu0] ]  );
+  deltaS -= 2 * real( conj(phi[ 0*lsize4 + x0 ]) * candU * phi[ 0*lsize4 + (*nb)[x0*8+nu0] ]  );
+  
+  deltaS += 2 * real( conj(phi[ 1*lsize4 + x0 ]) * conj(U[ x0*4 + nu0 ]) * phi[ 1*lsize4 + (*nb)[x0*8+nu0] ]  );
+  deltaS -= 2 * real( conj(phi[ 1*lsize4 + x0 ]) * conj(candU) * phi[ 1*lsize4 + (*nb)[x0*8+nu0] ]  );
+
+  return exp(-deltaS);
+}
+
+#endif

diff --git a/u1casc-ordinary/u1casc-ordinary.cpp b/u1casc-ordinary/u1casc-ordinary.cpp
new file mode 100644 (file)
index 0000000..f6f05fd
--- /dev/null
+++ b/u1casc-ordinary/u1casc-ordinary.cpp
@@ -0,0 +1,77 @@
+#include <complex>
+#include <cstring>
+
+#include "latlib/o815/o815.h"
+
+#include "sim.hpp"
+
+o815 *O815;
+sim *Sim;
+
+//const complex<double> _i_ = complex<double>(0.0,1.0);
+
+//#include "obs_phi2.hpp"
+
+o815::comoption specOps[] = {
+  { "kappaone", required_argument, NULL, 'k', "set inverse mass kappa_1", "min:max:inc" },
+  { "kappatwo", required_argument, NULL, 'm', "set inverse mass kappa_2", "min:max:inc" },
+  { "beta", required_argument, NULL, 'b', "set inverse gauge coupling beta", "min:max:inc" },
+  { "lambdaone", required_argument, NULL, 'p', "set quartic coupling lambda_1", "min:max:inc" },
+  { "lambdatwo", required_argument, NULL, 'q', "set quartic coupling lambda_2", "min:max:inc" },
+  { "", 0, NULL, 0, "", "" }
+};
+
+void parseSpecOps() 
+{
+  for (int isopt = 0; isopt < O815->parsedSpecOps.size(); isopt++)
+    switch(O815->parsedSpecOps[isopt].first) {
+    case 'm':
+      O815->paraQ->addRange("mass", O815->parsedSpecOps[isopt].second);
+      break;
+    }
+}
+
+void helpHeader() 
+{ 
+  cout << "Usage: ./u1casc-ordinary [OPTIONS] [obs1] ... [obsN]" << endl << endl;
+}
+
+void parseLonelyArgs()
+{
+  for (vector<char*>::iterator lonit = O815->lonelyArgs.begin(); lonit != O815->lonelyArgs.end(); ++lonit) {
+    /*
+    if ( strcmp(*lonit, "phi2") == 0 ) {
+      *O815->out->log << "MASTER: registered observable: phi2" << endl << flush;
+      O815->observables.push_back(new obs_phi2(O815));
+    }
+    else if ( strcmp(*lonit, "phi2_hist") == 0 ) {
+      *O815->out->log << "MASTER: registered observable: phi2_hist" << endl << flush;
+      O815->observables.push_back(new obs_phi2_hist(O815));
+    }
+    else if ( strcmp(*lonit, "phip2_hist") == 0 ) {
+      *O815->out->log << "MASTER: registered observable: phip2_hist" << endl << flush;
+      O815->observables.push_back(new obs_phip2_hist(O815));
+    }
+    */
+  }
+}
+
+int main (int argc, char *argv[])
+{
+  O815 = new o815(argc, argv, "u1casc-ordinary", specOps, &helpHeader);
+
+  O815->addPara("beta", 1);
+  O815->addPara("lambdaone", 1);
+  O815->addPara("lambdatwo", 1);
+  O815->addPara("kappaone", 1);
+  O815->addPara("kappatwo", 1);
+
+  parseSpecOps();
+  O815->postParaInit();
+  O815->Sim = new sim(O815);
+  parseLonelyArgs();
+  O815->mainLoop();
+
+  delete O815;
+  return 0;
+}