<reponame>Jakl3/Heat-Diffusion
"""
Author: <NAME>
"""

import numpy as np
import heapq
import sys
from tqdm import tqdm
from collections import defaultdict
from mesh_reader import MeshReader
import itertools


class Node:
    """
    A vertex/Node on the grid

    :Attributes:
    - Status: "FAR" - unvisited Node, "NAB" - narrow band/considered Node, "FRZ" - frozen Node
    - Value: distance at the vertex, computed by the T function
    """

    def __init__(self, status, value):
        self.status = status
        self.value = value

    def set_status_far(self):
        self.status = "FAR"

    def set_status_nb(self):
        self.status = "NAB"

    def set_status_frozen(self):
        self.status = "FRZ"

    def get_value(self):
        return self.value

    def set_value(self, value):
        self.value = value

    def is_far(self):
        return self.status == "FAR"

    def is_nb(self):
        return self.status == "NAB"

    def is_frozen(self):
        return self.status == "FRZ"

    def __str__(self):
        return '(%.1f %s)' % (self.value, self.status)

    def __repr__(self):
        return '(%.1f %s)' % (self.value, self.status)


class FMMDistance:
    """
    Creates the distance map from an initial condition using the FMM,
    computing distances via the Eikonal equation
    """

    def __init__(self, vertex_adj, initial_point):
        """
        Initialize the FMM solver
        :param vertex_adj: the vertex adjacency list generated by the mesh reader
        :param initial_point: the index of the initial point/vertex, selected by the user
        """
        self.mat = defaultdict()
        for key in vertex_adj:
            self.mat[key] = Node("FAR", np.inf)

        self.adj_list = vertex_adj

        self.F = 1  # there is a uniform propagation speed across the entire mesh
        self.list_nb = []

        self.mat[initial_point].set_value(0)
        self.mat[initial_point].set_status_frozen()
        # print(self.adj_list[initial_point])
        for vn in self.non_frozen_neighbors(initial_point):
            """
            Since the first set of neighbors only has one source vertex, we calculate its distance
            using the euclidean distance.
            """
            a = self.adj_list[initial_point].coord
            b = self.adj_list[vn].coord
            dist = np.linalg.norm(a - b)
            # print(dist, type(dist))
            self.mat[vn].set_value(dist)
            heapq.heappush(self.list_nb, (dist, vn))
            self.mat[vn].set_status_nb()

    def non_frozen_neighbors(self, point_idx):
        """
        Returns list of non frozen neighboring Nodes
        :param point_idx: index of the point
        :return result: list of non frozen neighboring Node indices
        """
        result = []
        children = self.adj_list[point_idx].children_vertices
        for child in children:
            if not self.mat[child].is_frozen():
                result.append(child)
        return result

    def set_distance(self, v):
        """
        Calculate/recalculate the distance T at the Node v(i, j)
        :param v: the Node to be updated; this Node should either be narrow band or far
        """
        # Find neighboring vertices that are frozen
        frozen = []
        children = self.adj_list[v].children_vertices
        for child in children:
            if self.mat[child].is_frozen():
                frozen.append(child)
        # Break the loop if there are less than two frozen neighboring vertices
        if len(frozen) < 2:
            return

        # Loop through all possible pairs of frozen; in other words, calculate T
        # value for each triangle that contains the desired vertex
        t_list = []
        for i in range(len(frozen)):
            for j in range(len(frozen)):
                if i == j:
                    continue

                ''' Pull out index of the frozen vertices '''
                choice_i = frozen[i]
                choice_j = frozen[j]

                # Make sure that the condition T(b) > T(a) is met
                if self.mat[choice_i].get_value() > self.mat[choice_j].get_value():
                    vert_b = choice_i
                    vert_a = choice_j
                else:
                    vert_a = choice_i
                    vert_b = choice_j

                # print('>', vert_a, vert_b)
                # print('>>', self.mat[vert_a].get_value(), self.mat[vert_b].get_value())

                va_adj = self.adj_list[vert_a]
                vb_adj = self.adj_list[vert_b]
                vc_adj = self.adj_list[v]

                va = self.mat[vert_a]
                vb = self.mat[vert_b]
                vc = self.mat[v]

                # Get the coordinates info from adjacency list
                vec_a = va_adj.coord
                vec_b = vb_adj.coord
                vec_c = vc_adj.coord

                # Get the length of the edges
                len_b = np.linalg.norm(vec_a - vec_c)
                len_a = np.linalg.norm(vec_b - vec_c)

                ''' Calculate the angles between the vectors '''
                # Normalize the vectors into unit vectors
                def normalize(vec):
                    norm = np.linalg.norm(vec)
                    if norm == 0:
                        return vec
                    return vec / norm

                vec_a = normalize(vec_a)
                vec_b = normalize(vec_b)

                # Precompute the angles for later
                theta = np.arccos(np.dot(vec_a, vec_b))
                sin_theta = np.sin(theta)
                cos_theta = np.dot(vec_a, vec_b)

                ''' Get the u value, which will satisfy the condition u >= 0 '''
                u = vb.get_value() - va.get_value()

                ''' Create and solve the quadratic function '''
                # Equation will be in the form of at^2 + bt + c = 0
                a = len_a**2 + len_b**2 - 2 * len_a * len_b * cos_theta
                b = 2 * len_b * u * (len_a * cos_theta - len_b)
                c = len_b**2 * (u**2 - self.F**2 * len_a**2 * sin_theta**2)

                # Discriminant in quadratic function
                discriminant = b**2 - 4 * a * c

                # There are several cases based on the result of the discriminant
                if discriminant < 0:
                    t = min(vc.get_value(), len_b + va.get_value(), len_a + vb.get_value())
                    t_list.append(t)
                else:
                    t1 = (-b + np.sqrt(discriminant)) / (2 * a)
                    t2 = (-b - np.sqrt(discriminant)) / (2 * a)

                    if t1 != 0:
                        pending1 = len_b*(t1 - u)/t1
                        if cos_theta != 0 and t1 > 0 and t1 > u and len_a * cos_theta < pending1 < len_a / cos_theta:
                            t = min(vc.get_value(), t1+va.get_value())
                            t_list.append(t)
                        else:
                            t = min(vc.get_value(), len_b + va.get_value(), len_a + vb.get_value())
                            t_list.append(t)

                    if t2 != 0:
                        pending2 = len_b*(t2 - u)/t2
                        if cos_theta != 0 and t2 > 0 and t2 > u and len_a * cos_theta < pending2 < len_a / cos_theta:
                            t = min(vc.get_value(), t2+va.get_value())
                            t_list.append(t)
                        else:
                            t = min(vc.get_value(), len_b + va.get_value(), len_a + vb.get_value())
                            t_list.append(t)

            ''' Update the narrow band heap '''
            t_fin = min(t_list)
            self.mat[v].set_value(t_fin)
            # If the current node is marked as far, add it to the heap
            if self.mat[v].is_far():
                heapq.heappush(self.list_nb, (t_fin, v))
                self.mat[v].set_status_nb()
            # If the current node is marked as narrow band, update its position in the heap
            else:
                # Find previous index
                idx = -1
                for nb_item in range(len(self.list_nb)):
                    val, vert = self.list_nb[nb_item]
                    if vert == v:
                        idx = nb_item
                        break
                self.list_nb[idx] = (t_fin, v)
                heapq.heapify(self.list_nb)

    def iterate(self):
        """
        Runs one iteration of the loop.
        """
        u = heapq.heappop(self.list_nb)
        u_idx = u[1]

        self.mat[u_idx].set_status_frozen()
        for vn in self.non_frozen_neighbors(u_idx):
            self.set_distance(vn)

    def test(self, n):
        for i in range(n):
            self.iterate()

    def calculate_distance(self):
        """
        Runs iterate until all distances are calculated.
        """
        iteration = 0
        dist_map = []

        def generator():
            while self.list_nb:
                yield

        for _ in tqdm(generator(), desc="Calculating Distance"):
            iteration += 1
            self.iterate()
            dist_map.append({key: self.mat[key].get_value() for key in self.mat})

        print(f"Calculations finished with {iteration} iterations.")
        return dist_map


if __name__ == "__main__":
    file_path = "/home/jack/Code/GitHub/Heat-Diffusion/VTK-Mesh-Render/data/"
    file_name = input("Give filename: ")
    # file_name = 'InternalNodes_2.vtk'

    reader = MeshReader(file_path=file_path, file_name=file_name)
    reader.setup()

    fmm = FMMDistance(vertex_adj=reader.vertex_adj, initial_point=0)
    fmm.calculate_distance()

    print(fmm.mat)
"""
Author: <NAME>
"""

import numpy as np
import heapq
import sys
from tqdm import tqdm
from collections import defaultdict
from mesh_reader import MeshReader
import itertools


class Node:
    """
    A vertex/Node on the grid

    :Attributes:
    - Status: "FAR" - unvisited Node, "NAB" - narrow band/considered Node, "FRZ" - frozen Node
    - Value: distance at the vertex, computed by the T function
    """

    def __init__(self, status, value):
        self.status = status
        self.value = value

    def set_status_far(self):
        self.status = "FAR"

    def set_status_nb(self):
        self.status = "NAB"

    def set_status_frozen(self):
        self.status = "FRZ"

    def get_value(self):
        return self.value

    def set_value(self, value):
        self.value = value

    def is_far(self):
        return self.status == "FAR"

    def is_nb(self):
        return self.status == "NAB"

    def is_frozen(self):
        return self.status == "FRZ"

    def __str__(self):
        return '(%.1f %s)' % (self.value, self.status)

    def __repr__(self):
        return '(%.1f %s)' % (self.value, self.status)


class FMMDistance:
    """
    Creates the distance map from an initial condition using the FMM,
    computing distances via the Eikonal equation
    """

    def __init__(self, vertex_adj, initial_point):
        """
        Initialize the FMM solver
        :param vertex_adj: the vertex adjacency list generated by the mesh reader
        :param initial_point: the index of the initial point/vertex, selected by the user
        """
        self.mat = defaultdict()
        for key in vertex_adj:
            self.mat[key] = Node("FAR", np.inf)

        self.adj_list = vertex_adj

        self.F = 1  # there is a uniform propagation speed across the entire mesh
        self.list_nb = []

        self.mat[initial_point].set_value(0)
        self.mat[initial_point].set_status_frozen()
        # print(self.adj_list[initial_point])
        for vn in self.non_frozen_neighbors(initial_point):
            """
            Since the first set of neighbors only has one source vertex, we calculate its distance
            using the euclidean distance.
            """
            a = self.adj_list[initial_point].coord
            b = self.adj_list[vn].coord
            dist = np.linalg.norm(a - b)
            # print(dist, type(dist))
            self.mat[vn].set_value(dist)
            heapq.heappush(self.list_nb, (dist, vn))
            self.mat[vn].set_status_nb()

    def non_frozen_neighbors(self, point_idx):
        """
        Returns list of non frozen neighboring Nodes
        :param point_idx: index of the point
        :return result: list of non frozen neighboring Node indices
        """
        result = []
        children = self.adj_list[point_idx].children_vertices
        for child in children:
            if not self.mat[child].is_frozen():
                result.append(child)
        return result

    def set_distance(self, v):
        """
        Calculate/recalculate the distance T at the Node v(i, j)
        :param v: the Node to be updated; this Node should either be narrow band or far
        """
        # Find neighboring vertices that are frozen
        frozen = []
        children = self.adj_list[v].children_vertices
        for child in children:
            if self.mat[child].is_frozen():
                frozen.append(child)
        # Break the loop if there are less than two frozen neighboring vertices
        if len(frozen) < 2:
            return

        # Loop through all possible pairs of frozen; in other words, calculate T
        # value for each triangle that contains the desired vertex
        t_list = []
        for i in range(len(frozen)):
            for j in range(len(frozen)):
                if i == j:
                    continue

                ''' Pull out index of the frozen vertices '''
                choice_i = frozen[i]
                choice_j = frozen[j]

                # Make sure that the condition T(b) > T(a) is met
                if self.mat[choice_i].get_value() > self.mat[choice_j].get_value():
                    vert_b = choice_i
                    vert_a = choice_j
                else:
                    vert_a = choice_i
                    vert_b = choice_j

                # print('>', vert_a, vert_b)
                # print('>>', self.mat[vert_a].get_value(), self.mat[vert_b].get_value())

                va_adj = self.adj_list[vert_a]
                vb_adj = self.adj_list[vert_b]
                vc_adj = self.adj_list[v]

                va = self.mat[vert_a]
                vb = self.mat[vert_b]
                vc = self.mat[v]

                # Get the coordinates info from adjacency list
                vec_a = va_adj.coord
                vec_b = vb_adj.coord
                vec_c = vc_adj.coord

                # Get the length of the edges
                len_b = np.linalg.norm(vec_a - vec_c)
                len_a = np.linalg.norm(vec_b - vec_c)

                ''' Calculate the angles between the vectors '''
                # Normalize the vectors into unit vectors
                def normalize(vec):
                    norm = np.linalg.norm(vec)
                    if norm == 0:
                        return vec
                    return vec / norm

                vec_a = normalize(vec_a)
                vec_b = normalize(vec_b)

                # Precompute the angles for later
                theta = np.arccos(np.dot(vec_a, vec_b))
                sin_theta = np.sin(theta)
                cos_theta = np.dot(vec_a, vec_b)

                ''' Get the u value, which will satisfy the condition u >= 0 '''
                u = vb.get_value() - va.get_value()

                ''' Create and solve the quadratic function '''
                # Equation will be in the form of at^2 + bt + c = 0
                a = len_a**2 + len_b**2 - 2 * len_a * len_b * cos_theta
                b = 2 * len_b * u * (len_a * cos_theta - len_b)
                c = len_b**2 * (u**2 - self.F**2 * len_a**2 * sin_theta**2)

                # Discriminant in quadratic function
                discriminant = b**2 - 4 * a * c

                # There are several cases based on the result of the discriminant
                if discriminant < 0:
                    t = min(vc.get_value(), len_b + va.get_value(), len_a + vb.get_value())
                    t_list.append(t)
                else:
                    t1 = (-b + np.sqrt(discriminant)) / (2 * a)
                    t2 = (-b - np.sqrt(discriminant)) / (2 * a)

                    if t1 != 0:
                        pending1 = len_b*(t1 - u)/t1
                        if cos_theta != 0 and t1 > 0 and t1 > u and len_a * cos_theta < pending1 < len_a / cos_theta:
                            t = min(vc.get_value(), t1+va.get_value())
                            t_list.append(t)
                        else:
                            t = min(vc.get_value(), len_b + va.get_value(), len_a + vb.get_value())
                            t_list.append(t)

                    if t2 != 0:
                        pending2 = len_b*(t2 - u)/t2
                        if cos_theta != 0 and t2 > 0 and t2 > u and len_a * cos_theta < pending2 < len_a / cos_theta:
                            t = min(vc.get_value(), t2+va.get_value())
                            t_list.append(t)
                        else:
                            t = min(vc.get_value(), len_b + va.get_value(), len_a + vb.get_value())
                            t_list.append(t)

            ''' Update the narrow band heap '''
            t_fin = min(t_list)
            self.mat[v].set_value(t_fin)
            # If the current node is marked as far, add it to the heap
            if self.mat[v].is_far():
                heapq.heappush(self.list_nb, (t_fin, v))
                self.mat[v].set_status_nb()
            # If the current node is marked as narrow band, update its position in the heap
            else:
                # Find previous index
                idx = -1
                for nb_item in range(len(self.list_nb)):
                    val, vert = self.list_nb[nb_item]
                    if vert == v:
                        idx = nb_item
                        break
                self.list_nb[idx] = (t_fin, v)
                heapq.heapify(self.list_nb)

    def iterate(self):
        """
        Runs one iteration of the loop.
        """
        u = heapq.heappop(self.list_nb)
        u_idx = u[1]

        self.mat[u_idx].set_status_frozen()
        for vn in self.non_frozen_neighbors(u_idx):
            self.set_distance(vn)

    def test(self, n):
        for i in range(n):
            self.iterate()

    def calculate_distance(self):
        """
        Runs iterate until all distances are calculated.
        """
        iteration = 0
        dist_map = []

        def generator():
            while self.list_nb:
                yield

        for _ in tqdm(generator(), desc="Calculating Distance"):
            iteration += 1
            self.iterate()
            dist_map.append({key: self.mat[key].get_value() for key in self.mat})

        print(f"Calculations finished with {iteration} iterations.")
        return dist_map


if __name__ == "__main__":
    file_path = "/home/jack/Code/GitHub/Heat-Diffusion/VTK-Mesh-Render/data/"
    file_name = input("Give filename: ")
    # file_name = 'InternalNodes_2.vtk'

    reader = MeshReader(file_path=file_path, file_name=file_name)
    reader.setup()

    fmm = FMMDistance(vertex_adj=reader.vertex_adj, initial_point=0)
    fmm.calculate_distance()

    print(fmm.mat)
Author: <NAME> A vertex/Node on the grid      :Attributes:     - Status: "FAR" - unvisited Node, "NAB" - narrow band/considered Node, "FRZ" - frozen Node     - Value: distance at the vertex, computed by the T function Creates the distance map from an initial condition using the FMM,     computing distances via the Eikonal equation Initialize the FMM solver         :param vertex_adj: the vertex adjacency list generated by the mesh reader         :param initial_point: the index of the initial point/vertex, selected by the user # there is a uniform propagation speed across the entire mesh # print(self.adj_list[initial_point]) Since the first set of neighbors only has one source vertex, we calculate its distance             using the euclidean distance. # print(dist, type(dist)) Returns list of non frozen neighboring Nodes         :param point_idx: index of the point         :return result: list of non frozen neighboring Node indices Calculate/recalculate the distance T at the Node v(i, j)         :param v: the Node to be updated; this Node should either be narrow band or far # Find neighboring vertices that are frozen # Break the loop if there are less than two frozen neighboring vertices # Loop through all possible pairs of frozen; in other words, calculate T # value for each triangle that contains the desired vertex Pull out index of the frozen vertices # Make sure that the condition T(b) > T(a) is met # print('>', vert_a, vert_b) # print('>>', self.mat[vert_a].get_value(), self.mat[vert_b].get_value()) # Get the coordinates info from adjacency list # Get the length of the edges Calculate the angles between the vectors # Normalize the vectors into unit vectors # Precompute the angles for later Get the u value, which will satisfy the condition u >= 0 Create and solve the quadratic function # Equation will be in the form of at^2 + bt + c = 0 # Discriminant in quadratic function # There are several cases based on the result of the discriminant Update the narrow band heap # If the current node is marked as far, add it to the heap # If the current node is marked as narrow band, update its position in the heap # Find previous index Runs one iteration of the loop. Runs iterate until all distances are calculated. # file_name = 'InternalNodes_2.vtk'