Debuging swsym.bonds

SpinW · Apr 30, 2018 · 160ecab · 160ecab
1 parent 16cdd09
commit 160ecab
Show file tree

Hide file tree

Showing 8 changed files with 105 additions and 38 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -6,7 +6,7 @@ set(MAKE_IDIVIDUAL_TESTS "FALSE")
 
 set(CMAKE_CXX_STANDARD 11)
 #Make highly optimised
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
+#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
 
 # I've had problems with gcc < 5.3 and OpenMP does not work correctly.
 if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
@@ -60,7 +60,7 @@ add_library(
         swsym/src/symOperator.cpp swsym/include/symOperator.h
         src/atom.cpp include/atom.h
         src/templateFuncs.tpp include/templateFuncs.h # Template functions
-)
+        swsym/src/oporder.cpp swsym/include/oporder.h)
 
 # Link with the required libraries.
 target_link_libraries(SpinW

diff --git a/src/templateFuncs.tpp b/src/templateFuncs.tpp
@@ -158,7 +158,10 @@ template <typename T> arma::Mat<typename T::elem_type> cmodM(T &inMat, double to
 
 template <typename T> arma::Cube<typename T::elem_type> cmodC(T &inCube, double toll){
     arma::Cube<typename T::elem_type> retCube = armaModC(inCube,1);
-    retCube.elem(arma::find(retCube > (arma::ones(retCube.n_rows,retCube.n_cols, retCube.n_slices) - toll))) -= 1;
+    arma::umat indices = arma::find(retCube > (1 - toll));
+    if (indices.n_elem > 0) {
+        retCube.elem(arma::find(retCube > (1 - toll))) -= 1;
+    }
     return retCube;
 }
 
@@ -214,13 +217,12 @@ template <typename T> arma::Cube<typename T::elem_type> repCube(T &inCube,int i,
     return retCube;
 }
 
-template <typename T> std::tuple<arma::uvec, arma::uvec> isnewUC(T &A, T &B, double toll){
+template <typename T, typename TT> std::tuple<arma::urowvec, arma::uvec> isnewUC(T &A, TT &B, double toll){
 
-//    arma::uword nA[2] = {A.n_rows, A.n_cols};
-//    arma::uword nB[2] = {B.n_rows, B.n_cols};
+    std::cout << A << B << std::endl;
 
-    arma::Cube<typename T::elem_type> AC(A.n_rows, A.n_cols,1); AC.slice(0) = A;
-    arma::Cube<typename T::elem_type> BC(B.n_rows, B.n_cols,1); BC.slice(0) = B;
+    arma::Cube<typename T::elem_type> AC(A.memptr(), A.n_rows, A.n_cols, 1);
+    arma::Cube<typename T::elem_type> BC(B.memptr(), B.n_rows, B.n_cols, 1);
 
     int perm1[] = {2, 3, 1};
     int perm2[] = {3, 2, 1};
@@ -231,18 +233,20 @@ template <typename T> std::tuple<arma::uvec, arma::uvec> isnewUC(T &A, T &B, dou
     AC = repCube(AC, 1, (int)B.n_cols, 1);
     BC = repCube(BC, (int)A.n_cols, 1, 1);
 
-    arma::Cube<typename T::elem_type> temp = cmodC(arma::abs(AC - BC),toll);
+    arma::cube temp = cmodC(arma::abs(AC - BC),toll);
     temp %= temp; // Remember element wise multiplication....
-    arma::Mat<typename T::elem_type> preFind = arma::sum(temp,2);
+    arma::mat preFind = arma::sum(temp, 2);
 
-    arma::umat postFind = arma::find(preFind > toll);
-    arma::uvec isNew = arma::all(postFind, 0);
+    preFind.elem(arma::find(preFind > toll)).ones();
+    preFind.elem(arma::find(preFind < toll)).zeros();
+    arma::urowvec isNew = arma::all(preFind > 0, 0);
 
-    arma::uvec idx = arma::linspace<arma::uvec>(0, B.n_cols - 1, B.n_cols);
+    //TODO FIX this linspace....
+    arma::urowvec idx = arma::linspace<arma::urowvec>(0, B.n_cols - 1, B.n_cols);
+    !preFind(arma::span(),idx(isNew))*
 
-    std::cout << postFind << std::endl << isNew << std::endl;
-    //    symIdx = max(bsxfun(@times,~notequal(:,idx(~isnew)),(1:size(notequal,1))'),[],1);
 
+    std::cout << preFind << std::endl << isNew << std::endl << idx << std::endl;
 
     return std::make_tuple(isNew, idx);
 };

diff --git a/swsym/include/oporder.h b/swsym/include/oporder.h
@@ -0,0 +1,8 @@
+//
+// Created by ward_s on 4/25/18.
+//
+
+#ifndef SPINW_OPORDER_H
+#define SPINW_OPORDER_H
+
+#endif //SPINW_OPORDER_H
diff --git a/swsym/include/swsym.h b/swsym/include/swsym.h
@@ -39,6 +39,7 @@ class swsym{
 
     int searchSym(std::string searchString);
     arma::cube symOperator(int symNumber);
+    int oporder(arma::mat symOP);
 };
 
 #endif //SPINW_SWSYM_H
diff --git a/swsym/src/bond.cpp b/swsym/src/bond.cpp
@@ -16,11 +16,11 @@ std::tuple<arma::mat, arma::umat> swsym::bond(arma::mat r,
     double tolDist = 1e-5;
 
 
-    arma::colvec r1 = r(arma::span(),bond(3));
-    arma::colvec r2 = r(arma::span(),bond(4));
+    arma::colvec r1 = r.col(bond(3));
+    arma::colvec r2 = r.col(bond(4));
     arma::colvec dl = bond(arma::span(0,2));
 
-    arma::cube thisSym = symOp(symN);
+    arma::cube thisSym = symOperator(symN);
     arma::cube redThisSym = thisSym(arma::span(), arma::span(0, 2), arma::span());
     arma::cube addThisSym = thisSym(arma::span(), arma::span(3, 3), arma::span());
 
@@ -56,16 +56,16 @@ std::tuple<arma::mat, arma::umat> swsym::bond(arma::mat r,
         r2New = arma::join_rows(r2New,r2TempC.slice(0));
         dlNew = arma::join_rows(dlNew,dlTempC.slice(0));
     }
-    std::tuple<arma::uvec, arma::uvec> atom1 = isnewUC(r,r1New,tolDist);
-    arma::uvec atom1Select = std::get<1>(atom1);
+    std::tuple<arma::urowvec, arma::uvec> thisAtom1 = isnewUC(r,r1New,tolDist);
+    arma::uvec atom1Select = std::get<1>(thisAtom1);
 
-    if (arma::any(std::get<0>(atom1))){
+    if (arma::any(std::get<0>(thisAtom1))){
         throw std::runtime_error("The generated positions for atom1 are wrong!");
     }
 
-    std::tuple<arma::uvec, arma::uvec> atom2 = isnewUC(r,r2New,tolDist);
-    arma::uvec atom2Select = std::get<1>(atom2);
-    if (arma::any(std::get<0>(atom2))){
+    std::tuple<arma::urowvec, arma::uvec> thisAtom2 = isnewUC(r,r2New,tolDist);
+    arma::uvec atom2Select = std::get<1>(thisAtom2);
+    if (arma::any(std::get<0>(thisAtom2))){
         throw std::runtime_error("The generated positions for atom2 are wrong!");
     }
 

diff --git a/swsym/src/oporder.cpp b/swsym/src/oporder.cpp
@@ -0,0 +1,22 @@
+//
+// Created by ward_s on 4/25/18.
+//
+
+#include "../include/swsym.h"
+#include "../../include/templateFuncs.h"
+
+int swsym::oporder(arma::mat symOP){
+    arma::mat RN = symOP.cols(0, 2);
+    arma::colvec TN = arma::round(symOP.col(3) * 12);
+
+    int N = 0;
+
+    while ((arma::norm(RN - arma::eye(3, 3)) > 10 * std::numeric_limits<double>::epsilon() ||
+            arma::norm(TN)
+            ) && (N < 10)){
+        RN *= symOP.cols(0, 2);
+        TN = armaModM(arma::round(symOP.cols(0, 2)*TN + symOP.col(3)), 12);
+        N++;
+    }
+    return N;
+}
diff --git a/swsym/src/position.cpp b/swsym/src/position.cpp
@@ -23,7 +23,6 @@ std::tuple<arma::mat, arma::urowvec> swsym::position(arma::mat &r0, int symN, do
         for (arma::uword s = 0; s < thisSym.n_slices; s++) {
             arma::mat mTemp = armaModM(mmat(redThisSym.slice(s), r0(arma::span(), i), order) + addThisSym.slice(s), 1);
             arma::cube rTemp(mTemp.memptr(), mTemp.n_rows, mTemp.n_cols, 1);
-//            rTemp.slice(0) = mTemp;
             rTemp = permute(rTemp, perm);
             if (rTemp.n_elem > 3) {
                 rTemp = armaModC(rTemp, 1);

diff --git a/swsym/src/symOperator.cpp b/swsym/src/symOperator.cpp
@@ -6,27 +6,60 @@
 #include "../../include/templateFuncs.h"
 
 arma::cube swsym::symOperator(int symNumber){
-    arma::cube genOP = getSym(symNumber - 1);
-    arma::uword nGen  = genOP.n_slices;
-    arma::cube symOp;
+
+
+    arma::cube genOP = getSym(symNumber);
+
+    genOP.each_slice([](arma::mat&X){
+        X.col(3) = X.col(3)*12;
+        X.col(3) = arma::round(X.col(3));
+    });
+
+    arma::cube symOp(3,4,1);
     symOp.slice(0)= arma::join_rows(arma::eye(3,3), arma::zeros(3,1));
 
-    genOP.each_slice([symOp](const arma::mat& X){
-        arma::mat R0 = X(arma::span(),arma::span(0,2));
-        arma::mat T0 = X(arma::span(),3);
+    for (arma::uword s = 0; s < genOP.n_slices; s++){
+        arma::mat X = genOP.slice(s);
+        arma::mat R0 = X.cols(0,2);
+        arma::colvec T0 = X.col(3);
 
         arma::mat R = arma::eye(3,3);
-        arma::mat T = arma::zeros(3,1);
+        arma::colvec T = arma::zeros(3,1);
 
-        //TODO Change P from 1 to the proper value
-        int P = 1;
+        int P = oporder(arma::join_rows(R0, T0));
 
         for (int i = 0; i < P; i++){
-            R %= R0;
-            T %= R0;
-            T += T0;
+            R = R0*R;
+            T = R0*T + T0;
 
+            arma::uword nSym = symOp.n_slices;
+            for (arma::uword j = 0; j < nSym; j++){
+                arma::mat temp = symOp.slice(j);
+                arma::mat RS = R * temp.cols(0, 2);
+                arma::mat TS = armaModM(arma::round(R*temp.col(3) + T), 12);
+
+                arma::cube newCube = symOp;
+                newCube.each_slice([RS, TS](arma::mat& Y) {
+                    Y(0, 0) = arma::accu(arma::abs(Y.cols(0,2) - RS));
+                    Y(0, 1) = arma::any(arma::vectorise(Y.col(3) - TS) > 0);
+                });
+
+                arma::colvec idxR = newCube.subcube(0, 0, 0, 0, 0, newCube.n_slices-1);
+                arma::colvec idxT = newCube.subcube(0, 1, 0, 0, 1, newCube.n_slices-1);
+
+                bool add = true;
+                for (int a = 0; a < idxR.n_elem; a++){
+                    if ((idxR(a) == 0) && (idxT(a) == 0)){
+                        add = false;
+                        break;
+                    }
+                }
+
+                if (add){
+                    symOp = arma::join_slices(symOp,arma::join_rows(RS,TS/12));
+                }
+            }
         }
-    });
+    }
     return  symOp;
 }