From 7d6a435bd0015bbed9c512e7517876842a9ae646 Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Fri, 17 Apr 2020 14:28:39 +0200 Subject: [PATCH 1/8] first stab --- src/main.cpp | 90 +++++++++++++++++++++++++++++----------------------- 1 file changed, 51 insertions(+), 39 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index b6f385f..316b483 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -159,6 +159,52 @@ void fillsims(const std::vector& args, std::vector& sims, Callable c) emplace_from_tuple(sims, c(p)); } +template +void loop(const Lattice& latt, const std::vector& params, Callable& filter, int nsweeps) +{ + // model + std::vector sims; + + // simulation vector + sims.reserve(params.size());//MARQOV has issues with copying -> reuires reserve in vector + + fillsims( params, sims, filter); + +//init + #pragma omp parallel for + for(std::size_t i = 0; i < sims.size(); ++i) //for OMP + { + auto myid = i; + auto& marqov = sims[i]; + + // number of cluster updates and metropolis sweeps + const int ncluster = 0; + + // number of EMCS during relaxation and measurement + const int nrlx = 5000; + // perform simulation + marqov.init_hot(); + marqov.wrmploop(nrlx, ncluster, nsweeps, myid); + } + +std::cout<<"warmup done"< id(nreplicas); @@ -295,17 +341,8 @@ int main() // lattice RegularLattice latt(L, dim); - - // model - std::vector >> sims; - - // simulation vector - sims.reserve(parameters.size());//MARQOV has issues with copying -> reuires reserve in vector - - - fillsims( parameters, - sims, - [&latt, &outdir, L]( decltype(parameters[0]) p) + + auto xxzfilter = [&latt, &outdir, L]( decltype(parameters[0]) p) { // write a filter to determine output file path and name std::string str_beta = "beta"+std::to_string(std::get<0>(p)); @@ -316,33 +353,8 @@ int main() std::string outsubdir = outdir+"/"+std::to_string(L)+"/"; return std::tuple_cat(std::forward_as_tuple(latt), std::make_tuple(outsubdir+outname), p); - } - ); - - - - - // ------------- execute ------------- - - #pragma omp parallel for - for(std::size_t i = 0; i < sims.size(); ++i) //for OMP - { - auto myid = i; - auto& marqov = sims[i]; + }; - // number of cluster updates and metropolis sweeps - const int ncluster = 0; - const int nsweeps = 2*L; - - // number of EMCS during relaxation and measurement - const int nrlx = 5000; - const int nmsr = 15000; - - - // perform simulation - marqov.init_hot(); - marqov.wrmploop(nrlx, ncluster, nsweeps, myid); - marqov.gameloop(nmsr, ncluster, nsweeps, myid); - } + loop >>(latt, parameters, xxzfilter, 2*L); } } -- GitLab From 3c3486d11b219b812e035f4f47b6a3eea562aede Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Fri, 17 Apr 2020 14:35:58 +0200 Subject: [PATCH 2/8] tidy up inclusions and move algebra helpers to marqov header --- src/main.cpp | 64 +++------------------------------------------------- src/marqov.h | 51 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 61 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 316b483..a6615da 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -11,6 +11,7 @@ #include "helpers.h" #include "cartprod.h" #include "registry.h" +#include "marqov.h" //helper, delete later @@ -65,59 +66,6 @@ class Grid std::vector getnbr(int i); }; -// ------- elementary state vector calculus - -template -StateVector operator + (StateVector lhs, StateVector rhs) -{ - StateVector res(lhs); - for(int i = 0; i < std::tuple_size::value; ++i) - res[i] += rhs[i]; - return res; -} - -template -StateVector operator - (StateVector lhs, StateVector rhs) -{ - StateVector res(lhs); - for(int i = 0; i < std::tuple_size::value; ++i) - res[i] -= rhs[i]; - return res; -} - -inline double dot(const double& a, const double& b) -{ - return a*b; -} - -template -inline typename VecType::value_type dot(const VecType& a, const VecType& b) -{ - typedef typename VecType::value_type FPType; - return std::inner_product(begin(a), end(a), begin(b), 0.0); -} - - -template -inline void reflect(StateVector& vec, const StateVector mirror) -{ - const int SymD = std::tuple_size::value; - - const double dotp = dot(vec,mirror); - - for (int i=0; i -inline void normalize(Container& a) -{ - typename Container::value_type tmp_abs=std::sqrt(dot(a, a)); - - for (int i = 0; i < a.size(); ++i) a[i] /= tmp_abs; -} - - - template inline void coutsv(StateVector& vec) { @@ -163,7 +111,7 @@ template void loop(const Lattice& latt, const std::vector& params, Callable& filter, int nsweeps) { // model - std::vector sims; + std::vector > sims; // simulation vector sims.reserve(params.size());//MARQOV has issues with copying -> reuires reserve in vector @@ -217,8 +165,6 @@ const int myid = 0; // remove once a parallelization is available #include "XXZAntiferro.h" #include "XXZAntiferroSingleAniso.h" -#include "marqov.h" - int main() { // read config files @@ -333,10 +279,6 @@ int main() for (auto& param_tuple : parameters) std::swap(std::get<0>(param_tuple), std::get<1>(param_tuple)); - - - - // ----------- set up simulations ------------ // lattice @@ -355,6 +297,6 @@ int main() return std::tuple_cat(std::forward_as_tuple(latt), std::make_tuple(outsubdir+outname), p); }; - loop >>(latt, parameters, xxzfilter, 2*L); + loop >(latt, parameters, xxzfilter, 2*L); } } diff --git a/src/marqov.h b/src/marqov.h index 2fa3dc8..41983f7 100644 --- a/src/marqov.h +++ b/src/marqov.h @@ -29,6 +29,57 @@ auto _call(Function f, Object& obj, Tuple t) { return _call(f, obj, t, std::make_index_sequence{}); } +// ------- elementary state vector calculus + +template +StateVector operator + (StateVector lhs, StateVector rhs) +{ + StateVector res(lhs); + for(int i = 0; i < std::tuple_size::value; ++i) + res[i] += rhs[i]; + return res; +} + +template +StateVector operator - (StateVector lhs, StateVector rhs) +{ + StateVector res(lhs); + for(int i = 0; i < std::tuple_size::value; ++i) + res[i] -= rhs[i]; + return res; +} + +inline double dot(const double& a, const double& b) +{ + return a*b; +} + +template +inline typename VecType::value_type dot(const VecType& a, const VecType& b) +{ + typedef typename VecType::value_type FPType; + return std::inner_product(begin(a), end(a), begin(b), 0.0); +} + + +template +inline void reflect(StateVector& vec, const StateVector mirror) +{ + const int SymD = std::tuple_size::value; + + const double dotp = dot(vec,mirror); + + for (int i=0; i +inline void normalize(Container& a) +{ + typename Container::value_type tmp_abs=std::sqrt(dot(a, a)); + + for (int i = 0; i < a.size(); ++i) a[i] /= tmp_abs; +} + // --------------------------- MARQOV CLASS ------------------------------- template -- GitLab From f7b64e31300dbe227f7aacbaf6999a533fffdfa4 Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Fri, 17 Apr 2020 14:36:20 +0200 Subject: [PATCH 3/8] remove unused grid dummy --- src/main.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index a6615da..612a749 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -61,11 +61,6 @@ struct GetTrait >//helper trait to break up a string at } }; -class Grid -{ - std::vector getnbr(int i); -}; - template inline void coutsv(StateVector& vec) { -- GitLab From d675c408000aa4b373c84721b2287523f08ee3b7 Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Fri, 17 Apr 2020 15:25:10 +0200 Subject: [PATCH 4/8] use registryfunctions for retrieving keys --- src/main.cpp | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 612a749..1247079 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -160,12 +160,35 @@ const int myid = 0; // remove once a parallelization is available #include "XXZAntiferro.h" #include "XXZAntiferroSingleAniso.h" +void selectsim(RegistryDB& reg) +{ + const std::string ham = reg.Get("mc", "General", "Hamiltonian" ); + if (ham == "Ising") + { + } + else if (ham == "Heisenberg") + { + } + else if (ham == "Phi4") + { + } + else if (ham == "BlumeCapel") + { + } + else if(ham == "XXZAntiferro") + { + } + else if(ham == "XXZAntiferroSingleAniso") + { + } +} + int main() { // read config files RegistryDB registry("../src/config"); - + //FIXME: NEVER DELETE USER DATA // remove old output and prepare new one std::string outdir = registry.Get("mc", "IO", "outdir" ); std::string logdir = registry.Get("mc", "IO", "logdir" ); @@ -211,8 +234,13 @@ int main() double lvfinal = registry.Get("mc", "General", "lvfinal" ); int lvsteps = registry.Get("mc", "General", "lvsteps" ); - auto parnames = registry.Get>("mc", "Hamiltonian", "names"); - + auto parnames = registry.GetBlock("mc", "Hamiltonian").GetKeys(); + + for(auto t: out) + std::cout<> par; for (int i=0; i Date: Fri, 17 Apr 2020 15:26:09 +0200 Subject: [PATCH 5/8] remove debugging things --- src/main.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 1247079..888c7f5 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -235,12 +235,6 @@ int main() int lvsteps = registry.Get("mc", "General", "lvsteps" ); auto parnames = registry.GetBlock("mc", "Hamiltonian").GetKeys(); - - for(auto t: out) - std::cout<> par; for (int i=0; i Date: Fri, 17 Apr 2020 21:57:48 +0200 Subject: [PATCH 6/8] unify and actually call ALL Hamiltonians --- src/BlumeCapel.h | 2 +- src/config/mc | 4 +- src/main.cpp | 181 +++++++++++++++++++++++++---------------------- 3 files changed, 98 insertions(+), 89 deletions(-) diff --git a/src/BlumeCapel.h b/src/BlumeCapel.h index 71c23c9..0d20cc5 100644 --- a/src/BlumeCapel.h +++ b/src/BlumeCapel.h @@ -150,7 +150,7 @@ class BlumeCapel template inline void wolff_flip(StateVector& sv, const A a=0) { - sv[0] *= -static_cast(1.0); + sv[0] *= -static_cast(1.0); } }; diff --git a/src/config/mc b/src/config/mc index c18e1a3..0ffef70 100644 --- a/src/config/mc +++ b/src/config/mc @@ -8,10 +8,10 @@ lvstart = 0.91 lvfinal = 0.93 lvsteps = 10 loopstyle = lin +Hamiltonian = XXZAntiferroSingleAniso - -[Hamiltonian] // model parameters other than temperature +[XXZAntiferroSingleAniso] // model parameters other than temperature names = extfield, aniso, singleaniso diff --git a/src/main.cpp b/src/main.cpp index 888c7f5..7ce8680 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include "rndwrapper.h" #include "regular_lattice.h" #include "vectorhelpers.h" @@ -102,17 +103,17 @@ void fillsims(const std::vector& args, std::vector& sims, Callable c) emplace_from_tuple(sims, c(p)); } -template -void loop(const Lattice& latt, const std::vector& params, Callable& filter, int nsweeps) +template +void loop(const std::vector& params, Callable filter, int nsweeps) { // model std::vector > sims; // simulation vector - sims.reserve(params.size());//MARQOV has issues with copying -> reuires reserve in vector + sims.reserve(params.size());//MARQOV has issues with copying -> requires reserve in vector fillsims( params, sims, filter); - + //init #pragma omp parallel for for(std::size_t i = 0; i < sims.size(); ++i) //for OMP @@ -160,72 +161,81 @@ const int myid = 0; // remove once a parallelization is available #include "XXZAntiferro.h" #include "XXZAntiferroSingleAniso.h" -void selectsim(RegistryDB& reg) +template +void RegularLatticeloop(RegistryDB& reg, const std::string outdir, std::string logdir, const std::vector& parameters, Callable filter) { - const std::string ham = reg.Get("mc", "General", "Hamiltonian" ); + using namespace std::placeholders; + const auto dim = reg.Get("mc", "General", "dim" ); + const auto nL = reg.Get >("mc", "General", "nL" ); + + // ---------------- main loop ----------------- + + cout << endl << "The dimension is " << dim << endl; + + // lattice size loop + for (int j=0; j(parameters, f, 2*L); + } +} + +void selectsim(RegistryDB& registry, std::string outdir, std::string logdir) +{ + const std::string ham = registry.Get("mc", "General", "Hamiltonian" ); + const int nreplicas = registry.Get("mc", "General", "nreplicas" ); + auto defaultfilter = [](RegularLattice& latt, std::string outdir, int L, auto p) + { + // write a filter to determine output file path and name + std::string str_beta = "beta"+std::to_string(std::get<0>(p)); + std::string outname = str_beta+".h5"; + std::string outsubdir = outdir+"/"+std::to_string(L)+"/"; + return std::tuple_cat(std::forward_as_tuple(latt), std::make_tuple(outsubdir+outname), p); + }; if (ham == "Ising") { + auto betas = registry.Get >("mc", ham, "betas"); + std::vector myj = {1.0}; + auto parameters = cart_prod(betas, myj); + RegularLatticeloop >(registry, outdir, logdir, parameters, defaultfilter); } else if (ham == "Heisenberg") { + auto betas = registry.Get >("mc", ham, "betas"); + std::vector myj = {1.0}; + auto parameters = cart_prod(betas, myj); + RegularLatticeloop >(registry, outdir, logdir, parameters, defaultfilter); } else if (ham == "Phi4") { + auto betas = registry.Get >("mc", ham, "betas"); + std::vector myj = {1.0}; + auto parameters = cart_prod(betas, myj, myj); + RegularLatticeloop >(registry, outdir, logdir, parameters, defaultfilter); } else if (ham == "BlumeCapel") { + auto betas = registry.Get >("mc", ham, "betas"); + std::vector myj = {1.0}; + auto parameters = cart_prod(betas, myj); + RegularLatticeloop >(registry, outdir, logdir, parameters, defaultfilter); } else if(ham == "XXZAntiferro") { + auto betas = registry.Get >("mc", ham, "betas"); + std::vector myj = {1.0}; + auto parameters = cart_prod(betas, myj, myj,myj); + RegularLatticeloop >(registry, outdir, logdir, parameters, defaultfilter); } else if(ham == "XXZAntiferroSingleAniso") { - } -} - -int main() -{ - // read config files - RegistryDB registry("../src/config"); - - //FIXME: NEVER DELETE USER DATA - // remove old output and prepare new one - std::string outdir = registry.Get("mc", "IO", "outdir" ); - std::string logdir = registry.Get("mc", "IO", "logdir" ); - - std::string command; - command = "rm -r " + outdir; - system(command.c_str()); - command = "rm -r " + logdir; - system(command.c_str()); - - makeDir(outdir); - makeDir(logdir); - - - // ------------------ live view ----------------------- - /* - - const int L = 30; - RegularLattice lattice(L, 3); - std::string outfile = outdir+"temp.h5"; - Marqov > marqov(lattice, outfile, 1/0.36); - marqov.init_hot(); - const int ncluster = 0; - const int nsweeps = L/2; - marqov.wrmploop(50, ncluster, nsweeps); - marqov.gameloop_liveview(); - - */ - // ---------------------------------------------------- - - - - // --------- unpack configuration file --------- - - const auto dim = registry.Get("mc", "General", "dim" ); - const auto nL = registry.Get >("mc", "General", "nL" ); - const int nreplicas = registry.Get("mc", "General", "nreplicas" ); + // --------- unpack configuration file --------- auto lvname = registry.Get("mc", "General", "loopvar" ); auto loopstyle = registry.Get("mc", "General", "loopstyle" ); @@ -234,7 +244,7 @@ int main() double lvfinal = registry.Get("mc", "General", "lvfinal" ); int lvsteps = registry.Get("mc", "General", "lvsteps" ); - auto parnames = registry.GetBlock("mc", "Hamiltonian").GetKeys(); + auto parnames = registry.GetBlock("mc", ham).GetKeys(); std::vector> par; for (int i=0; i parval = {0}; - parval[0] = registry.Get("mc", "Hamiltonian", parnames[i]); + parval[0] = registry.Get("mc", ham, parnames[i]); par.push_back(parval); } } // create range for loop variable std::vector loopvar = create_range(lvstart, lvfinal, lvsteps); - - // write values in external fields in logfile - std::string logfile = registry.Get("mc", "IO", "logfile" ); - std::ofstream os(logdir+"/"+logfile); - os << std::setprecision(7); - for (int i=0; i(param_tuple), std::get<1>(param_tuple)); - // ----------- set up simulations ------------ - // lattice - RegularLattice latt(L, dim); - - auto xxzfilter = [&latt, &outdir, L]( decltype(parameters[0]) p) + // write values in external fields in logfile + std::string logfile = registry.Get("mc", "IO", "logfile" ); + std::ofstream os(logdir+"/"+logfile); + os << std::setprecision(7); + for (int i=0; i(p)); @@ -313,7 +301,28 @@ int main() return std::tuple_cat(std::forward_as_tuple(latt), std::make_tuple(outsubdir+outname), p); }; + RegularLatticeloop >(registry, outdir, logdir, parameters, xxzfilter); + } +} - loop >(latt, parameters, xxzfilter, 2*L); - } +int main() +{ + // read config files + RegistryDB registry("../src/config"); + + // remove old output and prepare new one + std::string outdir = registry.Get("mc", "IO", "outdir" ); + std::string logdir = registry.Get("mc", "IO", "logdir" ); + + //FIXME: NEVER DELETE USER DATA + std::string command; + command = "rm -r " + outdir; + system(command.c_str()); + command = "rm -r " + logdir; + system(command.c_str()); + + makeDir(outdir); + makeDir(logdir); + +selectsim(registry, outdir, logdir); } -- GitLab From 45b9d7db32ef3b31c1b86927d95f2120f0f56877 Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Sat, 18 Apr 2020 00:41:49 +0200 Subject: [PATCH 7/8] Add Neighbours class --- src/main.cpp | 19 +++++++++++++++++-- src/marqov.h | 34 ++++++++++++++++++++++++---------- 2 files changed, 41 insertions(+), 12 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 7ce8680..6c19fb2 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -13,6 +13,7 @@ #include "cartprod.h" #include "registry.h" #include "marqov.h" +#include "neighbourclass.h" //helper, delete later @@ -164,7 +165,6 @@ const int myid = 0; // remove once a parallelization is available template void RegularLatticeloop(RegistryDB& reg, const std::string outdir, std::string logdir, const std::vector& parameters, Callable filter) { - using namespace std::placeholders; const auto dim = reg.Get("mc", "General", "dim" ); const auto nL = reg.Get >("mc", "General", "nL" ); @@ -190,7 +190,7 @@ void selectsim(RegistryDB& registry, std::string outdir, std::string logdir) { const std::string ham = registry.Get("mc", "General", "Hamiltonian" ); const int nreplicas = registry.Get("mc", "General", "nreplicas" ); - auto defaultfilter = [](RegularLattice& latt, std::string outdir, int L, auto p) + auto defaultfilter = [](auto& latt, std::string outdir, int L, auto p) { // write a filter to determine output file path and name std::string str_beta = "beta"+std::to_string(std::get<0>(p)); @@ -303,6 +303,21 @@ void selectsim(RegistryDB& registry, std::string outdir, std::string logdir) }; RegularLatticeloop >(registry, outdir, logdir, parameters, xxzfilter); } + else if(ham == "IregularIsing") + { + std::vector > dummy; + Neighbours nbrs(dummy); + auto betas = registry.Get >("mc", ham, "betas"); + std::vector myj = {1.0}; + auto parameters = cart_prod(betas, myj); + // extract lattice size and prepare directories + int L = nbrs.size(); + cout << endl << "L = " << L << endl << endl; + makeDir(outdir+"/"+std::to_string(L)); + // lattice + auto f = [&defaultfilter, &nbrs, &outdir, L](auto p){return defaultfilter(nbrs, outdir, L, p);};//partially apply filter + loop, Neighbours >(parameters, f, 2*L); + } } int main() diff --git a/src/marqov.h b/src/marqov.h index 41983f7..39372f7 100644 --- a/src/marqov.h +++ b/src/marqov.h @@ -169,20 +169,36 @@ struct ObsTupleToObsCacheTuple Marqov(Marqov&& other) = default; Marqov& operator=(Marqov&& other) = default; - template + //For reference: this template is just to cool to forget.... +// template +// inline typename std::enable_if_t +// marqov_measure(std::tuple& t, Args... args) {} +// +// template +// inline typename std::enable_if_t +// marqov_measure(std::tuple& t, Args... args) +// { +// auto retval = _call(&std::tuple_element >::type::template measure, +// std::get(t), +// std::make_tuple(std::forward(args)...) ); +// marqov_measure(t, args...); +// std::get(obscache)< inline typename std::enable_if_t - marqov_measure(std::tuple& t, Args... args) {} + marqov_measure(std::tuple& t, S& s, G&& grid) {} - template + template inline typename std::enable_if_t - marqov_measure(std::tuple& t, Args... args) + marqov_measure(std::tuple& t, S& s, G&& grid) { auto retval = _call(&std::tuple_element >::type::template measure, + std::tuple >::type::template measure, std::get(t), - std::make_tuple(args...) ); - marqov_measure(t, args...); - std::get(obscache)<(t, s, grid); + std::get(obscache)< metro;//C++11 - // obs now handled differently - // number of EMCS static constexpr int nstep = 250; }; -- GitLab From 9c3065b4505be65608bb55c1c9cfee75ff139d01 Mon Sep 17 00:00:00 2001 From: Florian Goth Date: Sat, 18 Apr 2020 00:54:01 +0200 Subject: [PATCH 8/8] forgot neighbourclass --- src/neighbourclass.h | 312 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 312 insertions(+) create mode 100644 src/neighbourclass.h diff --git a/src/neighbourclass.h b/src/neighbourclass.h new file mode 100644 index 0000000..d5dfb53 --- /dev/null +++ b/src/neighbourclass.h @@ -0,0 +1,312 @@ +/* contact + Copyright (C) 2018 - 2020 Manuel Schrauth, Jefferson S.E. Portela, Florian Goth + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +*/ + +#ifndef NEIGHBOURCLASS_H +#define NEIGHBOURCLASS_H +#include +#include +#include + +/** This function prepares the input data gridx, gridy and n + * so that n is stored in a format usable by the neighbour class. + * This task can not be done by the neighbour class since it only encodes + * the neighbourhood relations but not the coordinates. Nevertheless the coordinates need to + * be changed, too. + * Again: After exection of this function gridx, gridy, and n have changed data! + * low amount of neighbours can be found at the beginning + * large amounts at the end of n. + * For a particular point the neighbours are sorted in ascending order. + * @param gridx The array of the x coordinats + * @param gridy the array of the y components + * @param n the vector containing the neighbour relations that needs to be transformed. + * */ +void prepare_for_neighbours_class(std::vector& gridx, std::vector& gridy, + std::vector< std::vector >& n) +{ + std::vector mymap; + for(const std::vector& v : n) + { + mymap.push_back(v.size()); + } + auto ret = std::minmax_element(mymap.begin(), mymap.end()); + int minneighbours = *ret.first; + int classes = 1 + *ret.second - minneighbours; + std::vector* arr = new std::vector[classes]; + for(uint i = 0; i < n.size(); ++i) + { + arr[n[i].size() - minneighbours].push_back(i); + } + //arr is now an array where we can find the indices of neighbours with a particular number of neighbours + std::vector newarr; + for(uint i = 0; i < classes; ++i) + for(uint j = 0; j < arr[i].size(); ++j) + newarr.push_back(arr[i][j]); + //newarr is the mapping which new index corresponds to which old index + std::vector invnewarr(newarr.size()); + for(int i = 0; i < newarr.size(); ++i) + invnewarr[newarr[i]] = i; + //invnewarr holds the info which old index corresponds to which new index + // now we should have all data to map the indices to the new format. + delete [] arr; + std::vector gridxnew(n.size()); + std::vector gridynew(n.size()); + std::vector< std::vector > nnew(n.size()); + for (uint i = 0; i < n.size(); ++i) + { + gridxnew[i] = gridx[newarr[i]]; + gridynew[i] = gridy[newarr[i]]; + const std::vector& oldpoint = n[newarr[i]]; + nnew[i].resize(oldpoint.size()); + for(uint j = 0; j < oldpoint.size(); ++j) + nnew[i][j] = invnewarr[oldpoint[j]]; + //impose some more structure by sorting the neighbour indices. + std::sort(nnew[i].begin(), nnew[i].end()); + } + gridx = gridxnew; + gridy = gridynew; + n = nnew; + //The following commented out block can be used to dump the resulting structure + // std::ofstream os("test.dat"); + // for(uint i = 0; i < n.size(); ++i) + // { + // os<<2+n[i].size()<<" "< +class Neighbours +{ +public: + /** This function returns the position of a particular point in the neighbour array. + * can be found + * @param c The class of the point + * @param pos The index of the point + * */ + inline std::size_t mypos(int c, const IntType& pos) const + { + int res = bookkeeping[c].arr_idx_beg + bookkeeping[c].nbr_class*(pos - bookkeeping[c].nbr_idx_beg); + return res; + } + /** This function queries to which neighour class an index belongs. + * Hence it can be used to query the number of neighbours. + * @param pos The position of the index that we ask + * @return The number of neighbours of point pos. + * */ + inline int get_class(const IntType& pos) const + { + int c = 0; + while (((c) < nr_of_starts) && (pos >= bookkeeping[c].nbr_idx_beg )) ++c; + return c-1; + } +// inline IntType& operator[](const IntType& pos) const +// { +// return data[pos*nbsize]; +// } + + inline auto getnbrs(int fam, int pos) const + { + //Determine neighbour class (i.e. has it three or four or five neighbours) + int c = get_class(pos); + //determine the offset of the relevant data into the array + std::size_t seg = mypos(c, pos); + return std::vector(data + seg, data + seg + c); + } + + /** This function gives a particular neighbour of a site at a particular position. + * @param pos The index of the position. + * @param idx The index of the neighbour + * @return The index that the neighbour has in this array. + * */ + inline IntType& get_neighbour(const IntType& pos, const IntType& idx) const + { + //Determine neighbour class (i.e. has it three or four or five neighbours) + int c = get_class(pos); + //determine the offset of the relevant data into the array + std::size_t seg = mypos(c, pos); + return data[seg + idx]; + } + + inline IntType getRandomNeighbour(std::size_t idx, uint16_t number_of_nn, const IntType& pos, int& rndnumber) const +// inline IntType getRandomNeighbour(std::size_t idx, uint16_t number_of_nn, const IntType& pos, int rndnumber) const + { +// /* + switch (number_of_nn) + { + case 1: + break; + case 2: + idx += (rndnumber)%2; + break; + case 3: + idx += (rndnumber)%3; + break; + case 4: + idx += (rndnumber)%4; + break; + default: + idx += rndnumber%number_of_nn; + break; + } +// */ +// idx += rndnumber; + return data[idx]; + } + + inline Neighbours(const std::vector >& n) : Neighbours(&n) {} + inline Neighbours(const std::vector >* n) + { + /*A helper structure to store some properties that we need for the construction of the data array*/ + struct Props { + uint occurences; ///< Used to store how often a particular amount of neighbours occured + uint min; ///< used to store the minimum neighbour that occured for a particular class + uint max; ///< used to store the maximum neighbour that occured for a particular class + uint lowest_idx; ///< used to find out where the first occurence of a particular class was + Props () : occurences(0), min(-1), max(0), lowest_idx(-1) {} + }; + std::map mymap; + for (uint i = 0; i < n->size(); ++i ) + { + uint myneighbours_count = (*n)[i].size(); + Props& prop = mymap[myneighbours_count];//if not found the [] operator calls the default copy constructor. + ++prop.occurences; + //we know that the neighbours are sorted in ascending order, therefore we easily find min and max + const std::vector& point = (*n)[i]; + if (point.size() > 0) + { + if (point[0] < prop.min) prop.min = point[0]; + if (point.back() > prop.max) prop.max = point.back(); + } + prop.lowest_idx = std::min(i, prop.lowest_idx); + } + //calculate required amount of memory elements + std::size_t memsize = 0; + for (auto it = mymap.cbegin(); it != mymap.cend(); ++it) + memsize += it->first*it->second.occurences; + //layout: first we store all with one neighbour, then all with two and so on... + nr_of_starts = mymap.size(); + int ret = posix_memalign((void**)&bookkeeping, 64, nr_of_starts*sizeof(typename Neighbours::Info));//request cache line aligned memory + if(ret != 0) throw("Can't allocate memory required for bookkeeping!!"); + // fill our array with the number of starts + auto it = mymap.cbegin(); + bookkeeping[0].nbr_class = it->first; + bookkeeping[0].arr_idx_beg = 0; + bookkeeping[0].nbr_idx_beg = it->second.lowest_idx; + for (uint i = 1; i < nr_of_starts; ++i) + { + auto prev = it; + ++it; + bookkeeping[i].nbr_class = it->first; + bookkeeping[i].arr_idx_beg = prev->second.occurences*prev->first + bookkeeping[i-1].arr_idx_beg; + bookkeeping[i].nbr_idx_beg = it->second.lowest_idx; + } + ret = posix_memalign((void**)&data, 4096, memsize*sizeof(IntType));//request page-aligned memory + if(ret != 0) throw("Can't allocate memory for neighbour array!!"); + if (memsize < 1E6) + { +// std::cout<<"Allocated "<size(); ++i) + { + int nbrs = (*n)[i].size(); + for(int j = 0; j < nbrs; ++j) + data[cur_pos++] = static_cast((*n)[i][j]); + } + nr_of_points = n->size(); + } + inline ~Neighbours() { + if (data != NULL) free(data); + if (bookkeeping != NULL) free(bookkeeping); + } + std::size_t size() const {return nr_of_points;} + + uint16_t neighbours_in_class(uint c) const {return bookkeeping[c].nbr_class;} + + /*This friend declaration enables a function from readwrite.h to write to this data structure directly. This is needed for a more memory efficient I/O.*/ + friend void import_geometry_directly(const int N, std::vector& gridx, std::vector& gridy, Neighbours& outn, const std::string path ); + + inline Neighbours() : data(NULL), bookkeeping(NULL), nr_of_starts(0), nr_of_points(0) {} + + /** + * A move assignment operator to avoid temporary copies. + * @param n the temporary object. + * @return An object with the state of n. + */ + inline Neighbours& operator=(Neighbours&& n) + { + if (this != &n) + { + //tidy up our own data + free(data); + free(bookkeeping); + //copy over state from other temporary object + data = n.data; + bookkeeping = n.bookkeeping; + nr_of_starts = n.nr_of_starts; + nr_of_points = n.nr_of_points; + //invalidate the OTHER class + n.nr_of_starts = 0; + n.nr_of_points = 0; + n.bookkeeping = NULL; + n.data = NULL; + } + return *this; + } + /** + * A move constructor to avoid temporary copies. + * @param n the temporary object. + * @return An object with the state of n. + */ + inline Neighbours(Neighbours&& n) + { + //copy over state from other temporary object + data = n.data; + bookkeeping = n.bookkeeping; + nr_of_starts = n.nr_of_starts; + nr_of_points = n.nr_of_points; + //invalidate the OTHER class + n.nr_of_starts = 0; + n.nr_of_points = 0; + n.bookkeeping = NULL; + n.data = NULL; + } +private: + struct Info { + std::size_t nbr_class;///< How many neighbours do we have? + std::size_t arr_idx_beg;///< where does this class start in the neighbour array? + std::size_t nbr_idx_beg;///< At which neighbour index does this start? + }; + IntType* data; ///< The data array start stores the actual neighbours + Info* bookkeeping;///< properties of the data array. The array is sorted in ascending order. + int16_t nr_of_starts; ///< The length of the bookkeeping array + std::size_t nr_of_points; ///< with this we can easily access the total amount of neighbours we have +}; +#endif -- GitLab