Working for D2Q9

lbmweights/exceptions.py

@@ -14,10 +14,17 @@ class ReducibleShellException(Exception):
     """Shell is reducible although it should be irreducible"""
     pass
 
-class NoSolutionError(Exception):
+
+class NoSolutionException(Exception):
     """Equation system cannot be solved"""
     pass
 
-class InfiniteSolutionsError(Exception):
+
+class InfiniteSolutionsException(Exception):
     """Infinite Solutions. May be solveable with Continue.py in https://github.com/BDuenweg/Lattice-Boltzmann-weights"""
+    pass
+
+
+class UninitializedAttributeException(Exception):
+    """Calling a function that is not supposed to be called yet"""
     pass
\ No newline at end of file

lbmweights/lattice.py

 import numpy as np
+import sympy as sp
+from sympy import oo
 from lbmweights.utils.mylog import logger
-from .exceptions import OrderNotImplementedException, ImpossibleVelocityShellException, ReducibleShellException, NoSolutionError, InfiniteSolutionsError
-from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, fill_rhs
+from .exceptions import OrderNotImplementedException, ImpossibleVelocityShellException, ReducibleShellException, NoSolutionException, InfiniteSolutionsException, UninitializedAttributeException
+from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, fill_rhs, approximate_ratio, timer
 
 
 class Lattice:
+
     def __init__(self, dimension=2, order=4, shell_list=[1, 2, 4], seed=None, svd_tolerance=1e-8):
         """
         :param dimension:
@@ -31,6 +34,9 @@ class Lattice:
         self._amount_singular_values = 0
         self._singular_values = []
 
+        self._weight_polynomials = []
+        self._interval = None
+
     def __str__(self):
         string = "D{} - Order: {} - Shells: {}".format(
             self._dimension, self._order, str(self._shell_list)
@@ -106,6 +112,35 @@ class Lattice:
 
         return velocities
 
+    def calculate_valid_interval(self):
+        if not self._weight_polynomials:
+            raise UninitializedAttributeException("Weight Polynomials not give. Calculate them first!")
+
+        interval = sp.Interval(-oo, oo)
+        root = []
+        c_s_sq = sp.Symbol("c_s_sq")
+        for equation in self._weight_polynomials:
+            roots = list(sp.roots(equation, c_s_sq, filter="R").keys())
+            roots.sort()
+            if len(roots) == 0:
+                continue
+            roots = [roots[0] - 1] + roots + [roots[-1] + 1]
+            cur_interval = sp.EmptySet()
+            for i in range(len(roots) - 1):
+                mesh_point = (roots[i + 1] + roots[i])/2
+                if equation.subs(c_s_sq, mesh_point) < 0:
+                    continue
+                # positive weights
+                if i == 0:
+                    cur_interval = sp.Union(cur_interval, sp.Interval(-oo, roots[i]))
+                elif i == (len(roots) - 1):
+                    cur_interval = sp.Union(cur_interval, sp.Interval(roots[i], oo))
+                else:
+                    cur_interval = sp.Union(cur_interval, sp.Interval(roots[i], roots[i+1]))
+            interval = interval.intersect(cur_interval)
+        self._interval = interval
+        return interval
+
     def calculate_weights(self):
         logger.debug("Entering Calculate Weights")
         logger.info(str(self))
@@ -148,7 +183,7 @@ class Lattice:
         for i, shell in enumerate(grand_total_list):
             number_of_velocities = len(shell)
             type = tuple(np.sort(abs(shell[0])))
-            logger.info("Shell number {} with c_i^2 = {} and {} velocities of type {}".format(i+1, (shell[0]**2).sum(), number_of_velocities, type))
+            #logger.info("Shell number {} with c_i^2 = {} and {} velocities of type {}".format(i+1, (shell[0]**2).sum(), number_of_velocities, type))
         # TO HERE: Only Calculate GrandTotalList
 
         # Here: ...
@@ -175,9 +210,11 @@ class Lattice:
 
         # --------------------------------------------------------------------------------------------------
         logger.warning("A Overwritten for test purposes")
+        """
         A = np.array([[2., 4., 8.],
                       [0.50244727, 1.99021093, 8.03915626],
                       [0.47686341, 2.09254637,  7.6298145]])
+                      """
         U, s, V = np.linalg.svd(A)
 
         rows = A.shape[0]
@@ -191,7 +228,7 @@ class Lattice:
             overlap = rows - rank
             if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
                 # NO Solution
-                raise NoSolutionError("Add shells or lower the order")  # change shells ?
+                raise NoSolutionException("Add shells or lower the order")  # change shells ?
             rows = rank
             s.resize(rows)
             rhs.resize((rows, rhs.shape[1]))
@@ -209,9 +246,78 @@ class Lattice:
 
         if cols - rows > 0:
             # TODO Better
-            raise InfiniteSolutionsError("Reduce order, remove shells, or use continue.py in "
+            raise InfiniteSolutionsException("Reduce order, remove shells, or use continue.py in "
                                          "https://github.com/BDuenweg/Lattice-Boltzmann-weights")
 
         # ----------------------------------------------- Unique solution
         solution = np.dot(np.transpose(V), reduced_rhs)
-        print(solution)
\ No newline at end of file
+        solution_columns = self._order // 2
+        assert solution_columns == solution.shape[1]
+
+        coefficients = np.zeros((1 + solution_columns))
+        coefficients[0] = 1
+        for j in range(solution_columns):
+            tmp = 0
+            for i in range(len(grand_total_list)):
+                tmp -= solution[i, j] * len(grand_total_list[i])
+            coefficients[j+1] = tmp
+
+        # -------------------------------------------------- Prettify for sympy
+        weight_polynomials = []
+        c_s_sq = sp.Symbol("c_s_sq")
+        equation = 0
+        for j in range(len(coefficients)):
+            equation += approximate_ratio(coefficients[j]) * c_s_sq ** j
+        weight_polynomials.append(equation)
+        for i in range(len(grand_total_list)):
+            equation = 0
+            for j in range(solution_columns):
+                equation += approximate_ratio(solution[i, j]) * c_s_sq ** (j + 1)
+            weight_polynomials.append(equation)
+        for i in weight_polynomials:
+            print(i)
+
+        # -------------------------------------------------- Calculate Roots
+        """
+        #interval = sp.Interval(-oo, oo)
+        #for equation in weight_polynomials:
+        #    _interval = sp.solveset(equation > 0, c_s_sq, sp.Reals)
+        #    interval = interval.intersect(_interval)
+
+        TOL = 1e-10
+        tmp = []
+        for i in range(len(coefficients)):
+            tmp.append(coefficients[len(coefficients) - 1 - i])
+        tmp_array = np.array(tmp)
+        roots = [root for root in np.roots(tmp_array) if abs(root.imag) < TOL and abs(root.real) > TOL]
+
+        rows = solution.shape[0]
+        cols = solution.shape[1]
+        for j in range(rows):
+            tmp = []
+            for i in range(cols):
+                tmp.append(solution[j, cols - 1 - i])
+            tmp_array = np.array(tmp)
+            roots.append([root for root in np.roots(tmp_array) if abs(root.imag) < TOL and abs(root.real) > TOL])
+        roots.sort()
+        if not len(roots):
+            return []
+        """
+
+        self._weight_polynomials = weight_polynomials
+        self.calculate_valid_interval()
+        logger.info(self._interval)
+
+
+
+
+
+
+
+
+
+
+
+
+
+

lbmweights/utils/utils.py

+from .mylog import logger
 import itertools
 import numpy as np
+import time
 import math
+from fractions import Fraction
+
+times = dict()
+
+
+def timer(func, *args, **kwargs):
+    """
+    Three different methods to time a function.
+
+    Best is used by default, because it's arguably the most important metric of time measurement.
+    """
+
+    def cur_time(func_name, time):
+        return time
+
+    def avg_time(func_name, time):
+        global times
+        if not times.get(func_name):
+            times[func_name] = [1, time]
+        else:
+            avg = times[func_name][1]
+            times[func_name][0] += 1
+            x = times[func_name][0]
+            times[func_name][1] = (x - 1) / x * avg + 1 / x * time
+        return times[func_name][1]
+
+    def best_time(func_name, time):
+        if not times.get(func_name):
+            times[func_name] = time
+        else:
+            if time < times[func_name]:
+                times[func_name] = time
+        return times[func_name]
+
+    def wrapper(*args, **kwargs):
+        start_time = time.time()
+        ret = func(*args, **kwargs)
+        elapsed_time = time.time() - start_time
+
+        logger.info(
+            " - Execution time of {0}: {1}".format(
+                func.__name__, best_time(func.__name__, elapsed_time)
+            )
+        )
+        return ret
+
+    return wrapper
 
 
 def analyse_tensor_dimension(tensor_rank):
@@ -112,6 +161,12 @@ def fill_rhs(max_rank, shells):
     return np.array(rhs)
 
 
+def approximate_ratio(x):
+    if abs(x) < 1e-12:
+        return 0
+    denominator_limit = 1000000 if abs(x) >= 0.1 else int(1. / abs(x)) * 1000000
+    return Fraction(x).limit_denominator(denominator_limit)
+
 
 
 

requirements.txt

 pytest
 pytest-runner
 numpy
-click
\ No newline at end of file
+click
+sympy
\ No newline at end of file