Commit cb526af7 authored by MischaD's avatar MischaD
Browse files

Implemented SVD

parent 3559138d
......@@ -12,4 +12,12 @@ class ImpossibleVelocityShellException(Exception):
class ReducibleShellException(Exception):
"""Shell is reducible although it should be irreducible"""
pass
class NoSolutionError(Exception):
"""Equation system cannot be solved"""
pass
class InfiniteSolutionsError(Exception):
"""Infinite Solutions. May be solveable with Continue.py in https://github.com/BDuenweg/Lattice-Boltzmann-weights"""
pass
\ No newline at end of file
import numpy as np
import random
from lbmweights.utils.mylog import logger
from .exceptions import OrderNotImplementedException, ImpossibleVelocityShellException, ReducibleShellException
from .utils.utils import analyse_tensor_dimension, get_group, get_subshells
from .exceptions import OrderNotImplementedException, ImpossibleVelocityShellException, ReducibleShellException, NoSolutionError, InfiniteSolutionsError
from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, fill_rhs
class Lattice:
def __init__(self, dimension=2, order=4, shell_list=[1, 2, 4], seed=None):
def __init__(self, dimension=2, order=4, shell_list=[1, 2, 4], seed=None, svd_tolerance=1e-8):
"""
:param dimension:
:param order:
:param shell_list:
:param seed:
:param tolerance: tolerance to set the singular values to zero after svd
"""
self._dimensions = dimension
self._dimension = dimension
self._order = order
self._shell_list = shell_list
self._seed = seed
self._svd_tolerance = svd_tolerance
np.random.seed(self._seed)
# LES
self._tensor_space_dimension = 0
self._possible_tensors = []
self._tensor_dimensions = []
self._velocities = 0
self._shells = 0
# SVD
self._amount_singular_values = 0
self._singular_values = []
def __str__(self):
string = "D{} - Order: {} - Shells: {}".format(
self._dimensions, self._order, str(self._shell_list)
self._dimension, self._order, str(self._shell_list)
)
if self._seed:
string += " - Seed: {}".format(self._seed)
......@@ -111,15 +119,16 @@ class Lattice:
for k in np.arange(2, self._order + 2, 2):
possible_tensors = analyse_tensor_dimension(k)
self._possible_tensors.append(possible_tensors)
self._tensor_dimensions = [len(x) for x in self._possible_tensors]
################################################
# FROM HERE ...
velocities_amount = 1
grand_total_list = []
total_subshells = []
group = get_group(self._dimensions) # TODO test output
group = get_group(self._dimension) # TODO test output
for shell in range(len(self._shell_list)):
velocities = self.velocities_for_shell(dimension=self._dimensions, shell=self._shell_list[shell])
velocities = self.velocities_for_shell(dimension=self._dimension, shell=self._shell_list[shell])
if len(velocities) == 0:
raise ImpossibleVelocityShellException("No velocity sets found for velocity shell {}.".format(shell))
......@@ -140,6 +149,69 @@ class Lattice:
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))
random.seed(self._seed)
# TO HERE: Only Calculate GrandTotalList
# Here: ...
# SpacialDimension = self._dimension
# MaxTensorRank = self._order
# ListOfTensordimensions = [len(x) for x in self._possible_tensors]
# GrandTotalList = grand_total_list == list shells which are list of arrays with the actual values
# Arguments = <class 'dict'>: {'d': 2, 'm': 4, 'c': [1, 2, 4], 's': 44, 'y': 'test': 'quiet': 'write_latex':}
lhs = []
for i, rank in enumerate(np.arange(2, self._order + 2, 2)):
tensor_space_dimension = self._tensor_dimensions[i]
for j in range(tensor_space_dimension):
vector = 2 * np.random.rand(self._dimension) - 1
vector = vector / np.linalg.norm(vector)
rows = []
for shell_velocities in grand_total_list:
shell_sum = lattice_sum(vector, shell_velocities, rank)
rows.append(shell_sum)
lhs.append(rows)
A = np.array(lhs)
b = fill_rhs(self._order, self._tensor_dimensions)
# --------------------------------------------------------------------------------------------------
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]
cols = A.shape[1]
self._singular_values = len(s)
s = np.array([sv if sv > self._svd_tolerance else 0 for sv in s])
rank = np.count_nonzero(s)
rhs = np.dot(np.transpose(U), b)
if rank < rows:
overlap = rows - rank
if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
# NO Solution
raise NoSolutionError("Add shells or lower the order") # change shells ?
rows = rank
s.resize(rows)
rhs.resize((rows, rhs.shape[1]))
self._singular_values = len(s)
reduced_rhs = np.zeros((rows, rhs.shape[1]))
for i, SingularValue in enumerate(s):
for j in range(rhs.shape[1]):
reduced_rhs[i, j] = rhs[i, j] / SingularValue
x = np.copy(reduced_rhs)
reduced_rhs = np.zeros((rows, rhs.shape[1]))
for i, singular_value in enumerate(s):
reduced_rhs[i, :] = rhs[i, :] / singular_value
if cols - rows > 0:
# TODO Better
raise InfiniteSolutionsError("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
......@@ -62,9 +62,11 @@ def get_group(dimension):
def contains(arr, list_of_arr):
return any((arr == x).all() for x in list_of_arr)
def contains_in_sublist(arr, list_of_lists):
return any(contains(arr, list_of_arr) for list_of_arr in list_of_lists)
def get_subshells(shell, group):
subshells = []
for velocity in shell:
......@@ -78,3 +80,38 @@ def get_subshells(shell, group):
subshells.append(subshell)
return subshells
def double_factorial(n):
assert n >= 0
if n == 0 or n == 1:
return 1
return n * double_factorial(n-2)
def lattice_sum(vector, velocities, rank):
ret = 0
for velocity in velocities:
scalar_product = 0
for dim in range(len(vector)):
scalar_product += vector[dim] * velocity[dim]
ret += scalar_product ** rank
return ret / double_factorial(rank - 1)
def fill_rhs(max_rank, shells):
rhs = []
cols = max_rank // 2
for k, rank in enumerate(np.arange(2, max_rank + 2, 2)):
for j in range(shells[k]):
rows = []
for i in range(cols):
c_s_power = 2 * i + 2
rows.append(1 if c_s_power == rank else 0)
rhs.append(rows)
return np.array(rhs)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment