#!/usr/bin/env python
# *-* coding: utf-8 *-*
"""Generate a decouplings.dat for CRTomo

From a given extra_lines.dat file (see triangular grid generation in manual),
generate a decouplings.dat for the decoupling of cell interfaces from the
regularisation.

As only the start and end points are given, the algorithm must find all nodes
generated by GMSH along those lines. This is done by comparing the orthogonal
distance of each node to the line. All nodes within a certain threshold
distance are assumed to lie on the line. The threshold is computed by taking
the minimal node distance of the grid (pdist function) and then multiplying
with 0.75. This probably is a little bit overkill and we could also use some
measure of machine precision.

For convenience this threshold can be overwritten using the --eps option.

Some of the for-loops can probably be optimized or replaced by some numpy
broadcasting rules.

Examples
--------

    # enter the output directory of cr_trig_create.py
    grid_extralines_gen_decouplings.py
    # this create the file decouplings.dat

    grid_extralines_gen_decouplings.py --plot_node_nrs --plot_elm_nrs

"""
import os
from optparse import OptionParser
import numpy as np
from crtomo.mpl_setup import *
import crtomo.grid as CRGrid


def handle_cmd_options():
    parser = OptionParser()
    parser.add_option(
        '-e',
        "--elem",
        dest="elem_file",
        type="string",
        help="elem.dat file (default: elem.dat)",
        default="elem.dat",
    )
    parser.add_option(
        "--eps",
        dest="eps",
        type="float",
        help="User override for distance eps",
        default=None,
    )
    parser.add_option(
        "--eta",
        action="store",
        dest="eta",
        type="float",
        help="User override for coupling coefficient",
        default=0.001,
    )
    parser.add_option(
        '-l',
        "--linefile",
        dest="linefile",
        help="Line file (default: extra_lines.dat)",
        default='extra_lines.dat',
    )
    parser.add_option(
        "-o",
        "--output",
        dest="output",
        help="Output file (default: decouplings.dat)",
        metavar="FILE", default="decouplings.dat",
    )

    parser.add_option(
        "--ln",
        action="store",
        dest="line_nr",
        type="int",
        help="Process only one line with index N (zero indexed)",
        default=None,
    )

    parser.add_option(
        '--debug_plot',
        action='store_true',
        dest='debug_plot',
        default=True,
        help='plot a debug plot',
    )

    parser.add_option(
        '--plot_node_nrs',
        action='store_true',
        dest='plot_node_nrs',
        default=False,
        help='plot node numbers in the debug plot. default: False',
    )

    parser.add_option(
        '--plot_elm_nrs',
        action='store_true',
        dest='plot_elm_nrs',
        default=False,
        help='plot elements numbers in the debug plot. default: False',
    )

    (options, args) = parser.parse_args()
    return options


def line_line_intersect(x, y):
    """Compute the intersection point of two lines

    Parameters
    ----------
    x = x4 array: x1, x2, x3, x4
    y = x4 array: y1, y2, y3, y4
    line 1 is defined by p1,p2
    line 2 is defined by p3,p4

    Returns
    -------
    Ix: x-coordinate of intersection
    Iy: y-coordinate of intersection
    """
    A = x[0] * y[1] - y[0] * x[1]
    B = x[2] * y[3] - y[2] * x[4]
    C = (x[0] - x[1]) * (y[2] - y[3]) - (y[0] - y[1]) * (x[2] - x[3])

    Ix = (A * (x[2] - x[3]) - (x[0] - x[1]) * B) / C
    Iy = (A * (y[2] - y[3]) - (y[0] - y[1]) * B) / C
    return Ix, Iy


def distances(x, y, px, py, plot=False):
    """
    Compute shortest distance of points (px, py) to the line described by the
    end points in (x, y).

    Parameters
    ----------
    x,y: x and y coordinates of line
    px: x coordinates of data points
    py: y coordiante of data points

    Returns
    -------
    distances: list of distances
    """

    length_line = np.sqrt((x[1] - x[0]) ** 2 + (y[1] - y[0]) ** 2)

    horizontal = (np.diff(y) == 0)
    vertical = (np.diff(x) == 0)

    dist = []
    dpl = []
    counter = 0
    for xp, yp in zip(px, py):
        # compute distances to the end points
        dp = np.sqrt((xp - x) ** 2 + (yp - y) ** 2)

        dpl.append(dp)
        if np.any(dp > length_line):
            point_dist = np.min(dp)
        elif horizontal:
            point_dist = np.abs(yp - y[0])
        elif vertical:
            point_dist = np.abs(xp - x[0])
        else:
            nominator = np.abs(
                (y[1] - y[0]) * xp -
                (x[1] - x[0]) * yp +
                x[1] * y[0] -
                y[1] * x[0])
            denominator = np.sqrt((y[1] - y[0]) ** 2 + (x[1] - x[0]) ** 2)
            point_dist = nominator / denominator

        dist.append(point_dist)
        counter += 1

    # debug plot
    if plot:
        fig, ax = plt.subplots(1, 1, figsize=(10, 8))
        ax.scatter(x, y)
        ax.scatter(px, py, color='r')

        index = 0
        for xp, yp in zip(px, py):
            ax.annotate('{0:.5}, {1:.5}, {2:.5}'.format(
                dist[index], dpl[index][0], dpl[index][1]),
                xy=(xp, yp))
            index += 1
        ax.plot(x, y)
        ax.set_aspect('equal')
    return dist


def get_decouplings_for_line(grid, line, settings, fids=None):
    """Compute the element pairs for regulazisation decoupling, given a grid
    object and a line, denoted by start and end coordinates, in the grid.

    Parameters
    ----------
    grid: a crtomo.grid.crt_grid object
    line: np.array/list with 4 entries: x1, y1, x2, x2, denoting start and end
          point
    fids: None, or: list with file ids for
        debug_nodes



    Returns
    -------
    neighbors: Mx2 array with the numbers of adjoining elements along the given
               line
    """
    key = 'presort'
    # x/z coordinates of 'macro'-line
    x = [line[0], line[2]]
    y = [line[1], line[3]]

    # x/z coordinates of all nodes in the grid
    nx = grid.nodes[key][:, 1]
    ny = grid.nodes[key][:, 2]
    nxy = np.vstack((nx, ny)).T

    # shortest distance of all nodes to this line
    dist = np.array(distances(x, y, nx, ny))

    # set the epsilon environment
    if settings.get('eps', None) is not None:
        eps = settings['eps']
        print('Taking user supplied eps value: {0}'.format(eps))
    else:
        import scipy.spatial.distance as spdist
        pair_distances = spdist.pdist(nxy)
        print('automatic eps determination')
        print('minimal distance: {0}'.format(pair_distances.min()))
        eps = np.sqrt(
            pair_distances.min() ** 2 - (pair_distances.min() ** 2) / 4
        )
        print('estimated smallest distance to next node: {0}'.format(eps))
        eps *= 0.5
        print('final eps (0.5 * estimate): {0}'.format(eps))
        if eps == 0:
            raise Exception('eps == 0')

    # only consider nodes that lie in the eps environment
    indices = np.where(dist < eps)[0]

    # extract these nodes
    nodes = grid.nodes[key][indices, 0]
    nodes = indices

    nodes_full = np.hstack((
        grid.nodes[key][indices, :].squeeze(),
        np.atleast_2d(dist[indices]).T,
    ))
    if fids is not None:
        np.savetxt(
            fids[0],
            nodes_full.squeeze()
        )
        np.savetxt(
            fids[1],
            np.atleast_2d(dist).T,
        )

    elm_indices = []
    elm_nodes = []
    # find elements with two nodes in it
    for elmnr, element in enumerate(grid.elements):
        if len(np.intersect1d(element, nodes)) == 2:
            elm_indices.append(elmnr)
            elm_nodes.append(element)

    elms = np.array(elm_nodes)

    # find neighboring elements
    neighbors = []
    eta = settings['eta']
    # for each element on the line
    for index, elm in enumerate(elms):
        found_it = False
        for index1, elm1 in enumerate(elms):
            ints = np.intersect1d(elm, elm1)
            # only two nodes can be on a line
            if len(ints) == 2:
                # this check ensures that we do not identify the
                # boundary between two adjacent cells on the line
                # as a decoupling line.
                if np.all(dist[ints] < eps):
                    found_it = True
                    break
                else:
                    pass
                    # print(
                    #     'distances of neighbor not correct:', dist[ints], eps
                    # )
            # probably the element itself
            elif len(ints) > 2:
                # print('found more than two common nodes!')
                # print('element nodes:', elm)
                # print('element1 nodes:', elm1)
                pass
        if found_it:
            nb = (elm_indices[index] + 1, elm_indices[index1] + 1, eta)
            # reversed neighbor
            # nb_rev = (nb[1], nb[0], nb[2])
            # if nb_rev not in neighbors:
            if True:
                neighbors.append(
                    nb
                )
        else:
            pass
            print('No neighbors found. Strange...')

    return np.array(neighbors)


def check_options(options):
    for filename in (options.elem_file, options.linefile):
        if not os.path.isfile(filename):
            raise IOError('File not found: {0}'.format(filename))


def debug_plot(grid, extra_lines, decoupling_elements, options,):
    print('Creating debug plot...')
    fig, ax = plt.subplots(1, 1, figsize=(10, 10))

    all_xz = []
    for x, z in zip(grid.grid['x'], grid.grid['z']):
        tmp = np.vstack((x, z)).T
        all_xz.append(tmp)

    collection = mpl.collections.PolyCollection(
        all_xz,
        edgecolor='k',
        facecolor='none',
        linewidth=0.4,
    )
    ax.add_collection(collection)
    if grid.electrodes is not None:
        ax.scatter(
            grid.electrodes[:, 1],
            grid.electrodes[:, 2],
            color='k',
            clip_on=False,
            label='electrodes',
        )

    # plot extra lines
    for line in np.atleast_2d(extra_lines):
        ax.plot(
            [line[0], line[2]],
            [line[1], line[3]],
            '.-',
            color='c',
            label='extra line',
        )
        # only line 0
        break

    # plot decouplings
    label = 'decoupling line'
    for (el1, el2, coef) in decoupling_elements:
        n1 = grid.elements[int(el1) - 1]
        n2 = grid.elements[int(el2) - 1]

        ints = np.intersect1d(n1, n2)

        x = grid.nodes['presort'][ints, 1]
        z = grid.nodes['presort'][ints, 2]

        ax.plot(
            x,
            z,
            '.-',
            color='r',
            linestyle='dashed',
            label=label,
        )
        label = ''

    # plot search nodes
    nodes = np.loadtxt('debug_search_nodes.dat')
    ax.scatter(
        nodes[:, 1],
        nodes[:, 2],
        color='g',
        s=40,
        clip_on=False,
        label='search nodes for decoupling',
    )
    for node in nodes:
        circle = plt.Circle(
            (node[1], node[2]),
            radius=node[3],
            fc='none',
        )
        ax.add_patch(circle)

    # plot grid numbers
    if options.plot_elm_nrs:
        for nr, quads in enumerate(grid.elements):
            # nr = grid.nodes['presort'][quads, 0],
            x = grid.nodes['presort'][quads, 1],
            z = grid.nodes['presort'][quads, 2],

            # print('NXZ', nr, x, z)
            ax.annotate(
                '{0}'.format(nr + 1),
                xy=(
                    np.mean(x),
                    np.mean(z),
                ),
                xycoords='data',
            )

    # # plot node numbers
    if options.plot_node_nrs:
        for nr, node in enumerate(grid.nodes['sorted']):
            ax.annotate(
                '{0}({1})'.format(nr, node[0]),
                xy=(
                    node[1],
                    node[2],
                ),
                xycoords='data',
                color='y',
            )

    ax.set_xlim(grid.grid['x'].min(), grid.grid['x'].max())
    ax.set_ylim(grid.grid['z'].min(), grid.grid['z'].max())
    # ax.autoscale_view()
    ax.set_aspect('equal')
    ax.set_xlabel('x')
    ax.set_ylabel('z')
    ax.set_title('Debug plot decoupling lines')
    fig.tight_layout()
    ax.legend(
        loc="lower center",
        ncol=4,
        bbox_to_anchor=(0, 0, 1, 1),
        bbox_transform=fig.transFigure,
    )
    fig.subplots_adjust(bottom=0.1)

    fig.savefig('debug_plot_decoupling.png', dpi=300, bbox_inches='tight')


def main():
    options = handle_cmd_options()
    check_options(options)
    grid = CRGrid.crt_grid()
    grid.load_elem_file(options.elem_file)

    extra_lines = np.atleast_2d(np.loadtxt(options.linefile))

    fids = []
    fids.append(
        open('debug_search_nodes.dat', 'wb')
    )
    fids.append(
        open('debug_dist.dat', 'wb')
    )
    neighbors = None
    if options.line_nr is None:
        lines = extra_lines
    else:
        print('processing only line: {0}'.format(options.line_nr))
        lines = (extra_lines[options.line_nr], )
    for line in lines:
        data = get_decouplings_for_line(
            grid,
            line,
            {
                'eps': options.eps,
                'eta': options.eta,
            },
            fids
        )
        if neighbors is None:
            neighbors = data
        else:
            try:
                neighbors = np.vstack(
                    (neighbors, data)
                )
            except:
                import IPython
                IPython.embed()

    for fid in fids:
        fid.close()

    with open(options.output, 'w') as fid:
        fid.write('{0}\n'.format(neighbors.shape[0]))
    with open(options.output, 'ab') as fid:
        np.savetxt(fid, np.array(neighbors), fmt='%i %i %.2f')

    if options.debug_plot:
        debug_plot(
            grid=grid,
            extra_lines=extra_lines,
            decoupling_elements=np.array(neighbors),
            options=options,
        )


if __name__ == '__main__':
    main()
