/*!
   \file vector/neta/flow.c

   \brief Network Analysis library - flow in graph

   Computes the length of the shortest path between all pairs of nodes
   in the network.

   (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team

   This program is free software under the GNU General Public License
   (>=v2). Read the file COPYING that comes with GRASS for details.

   \author Daniel Bundala (Google Summer of Code 2009)
 */

#include <stdio.h>
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/vector.h>
#include <grass/glocale.h>
#include <grass/dgl/graph.h>
#include <grass/neta.h>

dglInt32_t sign(dglInt32_t x)
{
    if (x >= 0)
        return 1;
    return -1;
}

/*!
   \brief Get max flow from source to sink.

   Array flow stores flow for each edge. Negative flow corresponds to a
   flow in opposite direction. The function assumes that the edge costs
   correspond to edge capacities.

   \param graph input graph
   \param source_list list of sources
   \param sink_list list of sinks
   \param[out] flow max flows

   \return number of flows
   \return -1 on failure
 */
int NetA_flow(dglGraph_s *graph, struct ilist *source_list,
              struct ilist *sink_list, int *flow)
{
    int nnodes, nlines, i;
    dglEdgesetTraverser_s et;
    dglInt32_t *queue;
    dglInt32_t **prev;
    char *is_source, *is_sink;
    int begin, end, total_flow;
    int have_node_costs;
    dglInt32_t ncost;

    nnodes = dglGet_NodeCount(graph);
    nlines = dglGet_EdgeCount(graph) /
             2; /*each line corresponds to two edges. One in each direction */
    queue = (dglInt32_t *)G_calloc(nnodes + 3, sizeof(dglInt32_t));
    prev = (dglInt32_t **)G_calloc(nnodes + 3, sizeof(dglInt32_t *));
    is_source = (char *)G_calloc(nnodes + 3, sizeof(char));
    is_sink = (char *)G_calloc(nnodes + 3, sizeof(char));
    if (!queue || !prev || !is_source || !is_sink) {
        G_fatal_error(_("Out of memory"));
        return -1;
    }

    for (i = 0; i < source_list->n_values; i++)
        is_source[source_list->value[i]] = 1;
    for (i = 0; i < sink_list->n_values; i++)
        is_sink[sink_list->value[i]] = 1;

    for (i = 0; i <= nlines; i++)
        flow[i] = 0;

    ncost = 0;
    have_node_costs = dglGet_NodeAttrSize(graph);

    total_flow = 0;
    while (1) {
        dglInt32_t node, edge_id, min_residue;
        int found = -1;

        begin = end = 0;
        for (i = 0; i < source_list->n_values; i++)
            queue[end++] = source_list->value[i];

        for (i = 1; i <= nnodes; i++) {
            prev[i] = NULL;
        }
        while (begin != end && found == -1) {
            dglInt32_t vertex = queue[begin++];
            dglInt32_t *edge, *node = dglGetNode(graph, vertex);

            dglEdgeset_T_Initialize(&et, graph,
                                    dglNodeGet_OutEdgeset(graph, node));
            for (edge = dglEdgeset_T_First(&et); edge;
                 edge = dglEdgeset_T_Next(&et)) {
                dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
                dglInt32_t id = dglEdgeGet_Id(graph, edge);
                dglInt32_t to =
                    dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
                if (!is_source[to] && prev[to] == NULL &&
                    cap > sign(id) * flow[labs(id)]) {
                    prev[to] = edge;
                    if (is_sink[to]) {
                        found = to;
                        break;
                    }
                    /* do not go through closed nodes */
                    if (have_node_costs) {
                        memcpy(&ncost,
                               dglNodeGet_Attr(graph,
                                               dglEdgeGet_Tail(graph, edge)),
                               sizeof(ncost));
                    }
                    if (ncost >= 0)
                        queue[end++] = to;
                }
            }
            dglEdgeset_T_Release(&et);
        }
        if (found == -1)
            break; /*no augmenting path */
        /*find minimum residual capacity along the augmenting path */
        node = found;
        edge_id = dglEdgeGet_Id(graph, prev[node]);
        min_residue = dglEdgeGet_Cost(graph, prev[node]) -
                      sign(edge_id) * flow[labs(edge_id)];
        while (!is_source[node]) {
            dglInt32_t residue;

            edge_id = dglEdgeGet_Id(graph, prev[node]);
            residue = dglEdgeGet_Cost(graph, prev[node]) -
                      sign(edge_id) * flow[labs(edge_id)];
            if (residue < min_residue)
                min_residue = residue;
            node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
        }
        total_flow += min_residue;
        /*update flow along the augmenting path */
        node = found;
        while (!is_source[node]) {
            edge_id = dglEdgeGet_Id(graph, prev[node]);
            flow[labs(edge_id)] += sign(edge_id) * min_residue;
            node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
        }
    }

    G_free(queue);
    G_free(prev);
    G_free(is_source);
    G_free(is_sink);

    return total_flow;
}

/*!
   \brief Calculates minimum cut between source(s) and sink(s).

   Flow is the array produced by NetA_flow() method when called with
   source_list and sink_list as the input. The output of this and
   NetA_flow() method should be the same.

   \param graph input graph
   \param source_list list of sources
   \param sink_list list of sinks (unused)
   \param flow
   \param[out] cut list of edges (cut)

   \return number of edges
   \return -1 on failure
 */
int NetA_min_cut(dglGraph_s *graph, struct ilist *source_list,
                 struct ilist *sink_list UNUSED, int *flow, struct ilist *cut)
{
    int nnodes, i;
    dglEdgesetTraverser_s et;
    dglInt32_t *queue;
    char *visited;
    int begin, end, total_flow;

    nnodes = dglGet_NodeCount(graph);
    queue = (dglInt32_t *)G_calloc(nnodes + 3, sizeof(dglInt32_t));
    visited = (char *)G_calloc(nnodes + 3, sizeof(char));
    if (!queue || !visited) {
        G_fatal_error(_("Out of memory"));
        return -1;
    }

    total_flow = begin = end = 0;

    for (i = 1; i <= nnodes; i++)
        visited[i] = 0;

    for (i = 0; i < source_list->n_values; i++) {
        queue[end++] = source_list->value[i];
        visited[source_list->value[i]] = 1;
    }

    /* find vertices reachable from source(s) using only non-saturated edges */
    while (begin != end) {
        dglInt32_t vertex = queue[begin++];
        dglInt32_t *edge, *node = dglGetNode(graph, vertex);

        dglEdgeset_T_Initialize(&et, graph, dglNodeGet_OutEdgeset(graph, node));
        for (edge = dglEdgeset_T_First(&et); edge;
             edge = dglEdgeset_T_Next(&et)) {
            dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
            dglInt32_t id = dglEdgeGet_Id(graph, edge);
            dglInt32_t to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
            if (!visited[to] && cap > sign(id) * flow[labs(id)]) {
                visited[to] = 1;
                queue[end++] = to;
            }
        }
        dglEdgeset_T_Release(&et);
    }
    /*saturated edges from reachable vertices to non-reachable ones form a
     * minimum cost */
    Vect_reset_list(cut);
    for (i = 1; i <= nnodes; i++) {
        if (!visited[i])
            continue;
        dglInt32_t *node, *edgeset, *edge;

        node = dglGetNode(graph, i);
        edgeset = dglNodeGet_OutEdgeset(graph, node);
        dglEdgeset_T_Initialize(&et, graph, edgeset);
        for (edge = dglEdgeset_T_First(&et); edge;
             edge = dglEdgeset_T_Next(&et)) {
            dglInt32_t to, edge_id;

            to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
            edge_id = labs(dglEdgeGet_Id(graph, edge));
            if (!visited[to] && flow[edge_id] != 0) {
                Vect_list_append(cut, edge_id);
                total_flow += abs(flow[labs(edge_id)]);
            }
        }
        dglEdgeset_T_Release(&et);
    }

    G_free(visited);
    G_free(queue);
    return total_flow;
}

/*!
   \brief Splits each vertex of in graph into two vertices

   The method splits each vertex of in graph into two vertices: in
   vertex and out vertex. Also, it adds an edge from an in vertex to
   the corresponding out vertex (capacity=2) and it adds an edge from
   out vertex to in vertex for each edge present in the in graph
   (forward capacity=1, backward capacity=0). If the id of a vertex is
   v then id of in vertex is 2*v-1 and of out vertex 2*v.

   \param in from graph
   \param out to graph
   \param node_costs list of node costs

   \return number of undirected edges in the graph
   \return -1 on failure
 */
int NetA_split_vertices(dglGraph_s *in, dglGraph_s *out, int *node_costs)
{
    dglInt32_t opaqueset[16] = {360000, 0, 0, 0, 0, 0, 0, 0,
                                0,      0, 0, 0, 0, 0, 0, 0};
    dglNodeTraverser_s nt;
    dglInt32_t edge_cnt;
    dglInt32_t *cur_node;

    dglInitialize(out, (dglByte_t)1, (dglInt32_t)0, (dglInt32_t)0, opaqueset);
    dglNode_T_Initialize(&nt, in);
    edge_cnt = 0;
    dglInt32_t max_node_cost = 0;

    for (cur_node = dglNode_T_First(&nt); cur_node;
         cur_node = dglNode_T_Next(&nt)) {
        dglInt32_t v = dglNodeGet_Id(in, cur_node);

        edge_cnt++;
        dglInt32_t cost = 1;

        if (node_costs)
            cost = node_costs[v];
        /* skip closed nodes */
        if (cost < 0)
            continue;
        if (cost > max_node_cost)
            max_node_cost = cost;
        dglAddEdge(out, 2 * v - 1, 2 * v, cost, edge_cnt);
        dglAddEdge(out, 2 * v, 2 * v - 1, (dglInt32_t)0, -edge_cnt);
    }
    dglNode_T_Release(&nt);
    dglNode_T_Initialize(&nt, in);
    for (cur_node = dglNode_T_First(&nt); cur_node;
         cur_node = dglNode_T_Next(&nt)) {
        dglEdgesetTraverser_s et;
        dglInt32_t *edge;
        dglInt32_t v = dglNodeGet_Id(in, cur_node);
        dglInt32_t cost = 1;

        if (node_costs)
            cost = node_costs[v];
        /* skip closed nodes */
        if (cost < 0)
            continue;

        dglEdgeset_T_Initialize(&et, in, dglNodeGet_OutEdgeset(in, cur_node));
        for (edge = dglEdgeset_T_First(&et); edge;
             edge = dglEdgeset_T_Next(&et)) {
            dglInt32_t to;

            to = dglNodeGet_Id(in, dglEdgeGet_Tail(in, edge));
            edge_cnt++;
            dglAddEdge(out, 2 * v, 2 * to - 1, max_node_cost + 1, edge_cnt);
            dglAddEdge(out, 2 * to - 1, 2 * v, (dglInt32_t)0, -edge_cnt);
        }
        dglEdgeset_T_Release(&et);
    }
    dglNode_T_Release(&nt);
    if (dglFlatten(out) < 0)
        G_fatal_error(_("GngFlatten error"));
    return edge_cnt;
}