/****************************************************************************
 *
 * MODULE:       v.perturb
 * AUTHOR(S):    James Darrell McCauley darrell@mccauley-usa.com
 *                  http://mccauley-usa.com/
 * PURPOSE:      Random location perturbations of vector points
 *
 * COPYRIGHT:    (C) 1994-2009 by James Darrell McCauley and the GRASS
 *               Development Team
 *
 * Modification History:
 * 3/2006              added min and seed MN/SM ITC-irst
 * 2005                updated to GRASS 6 RB ITC-irst
 * 0.1B <18 Feb 1994>  first version (jdm)
 * 0.2B <02 Jan 1995>  clean man page, added html page (jdm)
 * 0.3B <25 Feb 1995>  cleaned up 'gcc -Wall' warnings (jdm)
 * <13 Sept 2000>      released under GPL
 *
 *               This program is free software under the GNU General
 *               Public License (>=v2). Read the file COPYING that
 *               comes with GRASS for details.
 *
 *****************************************************************************/

#include <stdlib.h>
#include <math.h>

#include <grass/gis.h>
#include <grass/dbmi.h>
#include <grass/vector.h>
#include <grass/glocale.h>

#include "perturb.h"
#include "zufall.h"

int main(int argc, char **argv)
{
    double p1, p2, numbers[1000], numbers2[1000];
    int (*rng)(int, double *);
    int i;
    int line, nlines, ttype, n, ret, seed, field;
    struct field_info *Fi, *Fin;
    int out_is_3d;
    double min = 0.;
    int debuglevel = 3;

    struct Map_info In, Out;

    struct line_pnts *Points;
    struct line_cats *Cats;

    struct Cell_head window;
    struct GModule *module;
    struct {
        struct Option *in, *out, *dist, *pars, *min, *seed, *field;
        struct Flag *gen_seed, *no_topo;
    } parm;

    G_gisinit(argv[0]);

    module = G_define_module();
    G_add_keyword(_("vector"));
    G_add_keyword(_("geometry"));
    G_add_keyword(_("statistics"));
    G_add_keyword(_("random"));
    G_add_keyword(_("point pattern"));
    G_add_keyword(_("level1"));

    module->description = _("Random location perturbations of vector points.");

    parm.in = G_define_standard_option(G_OPT_V_INPUT);

    parm.field = G_define_standard_option(G_OPT_V_FIELD_ALL);

    parm.out = G_define_standard_option(G_OPT_V_OUTPUT);

    parm.dist = G_define_option();
    parm.dist->key = "distribution";
    parm.dist->type = TYPE_STRING;
    parm.dist->required = NO;
    parm.dist->options = "uniform,normal";
    parm.dist->answer = "uniform";
    parm.dist->description = _("Distribution of perturbation");

    parm.pars = G_define_option();
    parm.pars->key = "parameters";
    parm.pars->type = TYPE_DOUBLE;
    parm.pars->required = YES;
    parm.pars->multiple = YES;
    parm.pars->label = _("Parameter(s) of distribution");
    parm.pars->description =
        _("If the distribution "
          "is uniform, only one parameter, the maximum, is needed. "
          "For a normal distribution, two parameters, the mean and "
          "standard deviation, are required.");

    parm.min = G_define_option();
    parm.min->key = "minimum";
    parm.min->type = TYPE_DOUBLE;
    parm.min->required = NO;
    parm.min->answer = "0.0";
    parm.min->description = _("Minimum deviation in map units");

    parm.seed = G_define_option();
    parm.seed->key = "seed";
    parm.seed->type = TYPE_INTEGER;
    parm.seed->required = NO;
    parm.seed->description = _("Seed for random number generation");

    parm.gen_seed = G_define_flag();
    parm.gen_seed->key = 's';
    parm.gen_seed->description =
        _("Generate random seed (result is non-deterministic)");

    parm.no_topo = G_define_standard_flag(G_FLG_V_TOPO);

    /* Either explicit seed or explicitly generate seed, but not both. */
    G_option_exclusive(parm.seed, parm.gen_seed, NULL);
    G_option_required(parm.seed, parm.gen_seed, NULL);

    if (G_parser(argc, argv))
        exit(EXIT_FAILURE);

    min = atof(parm.min->answer);
    if (parm.seed->answer) {
        seed = atol(parm.seed->answer);
        G_debug(3, "Read random seed from seed=: %d", seed);
    }
    else {
        seed = G_srand48_auto();
        G_debug(3, "Generated random seed (-s): %d", seed);
    }

    switch (parm.dist->answer[0]) {
    case 'u':
        rng = zufall;
        break;
    case 'n':
        rng = normalen;
        break;
    }
    if (rng == zufall) {
        i = sscanf(parm.pars->answer, "%lf", &p1);
        if (i != 1) {
            G_fatal_error(_("Error scanning arguments"));
        }
        else if (p1 <= 0)
            G_fatal_error(_("Maximum of uniform distribution must be >= zero"));
    }
    else {
        if ((i = sscanf(parm.pars->answer, "%lf,%lf", &p1, &p2)) != 2) {
            G_fatal_error(_("Error scanning arguments"));
        }
        if (p2 <= 0)
            G_fatal_error(
                _("Standard deviation of normal distribution must be >= zero"));
    }

    G_get_window(&window);

    /* Open input */
    Vect_set_open_level(2);
    if (Vect_open_old_head2(&In, parm.in->answer, "", parm.field->answer) < 0)
        G_fatal_error(_("Unable to open vector map <%s>"), parm.in->answer);

    field = Vect_get_field_number(&In, parm.field->answer);

    out_is_3d = WITHOUT_Z;
    if (Vect_is_3d(&In))
        out_is_3d = WITH_Z;

    /* Open output */
    if (Vect_open_new(&Out, parm.out->answer, out_is_3d) < 0)
        G_fatal_error(_("Unable to create vector map <%s>"), parm.out->answer);
    /* TODO: Add also optional z perturb support */

    Vect_hist_copy(&In, &Out);
    Vect_hist_command(&Out);

    /* Generate a bunch of random numbers */
    zufalli(&seed);
    myrng(numbers, 1000, rng, p1 - min, p2);
    myrng(numbers2, 1000, rng, p1, p2);

    Points = Vect_new_line_struct();
    Cats = Vect_new_cats_struct();

    nlines = Vect_get_num_lines(&In);

    /* Close input, re-open on level 1 */
    Vect_close(&In);
    Vect_set_open_level(1);
    if (Vect_open_old2(&In, parm.in->answer, "", parm.field->answer) < 0)
        G_fatal_error(_("Unable to open vector map <%s>"), parm.in->answer);

    i = 0;
    line = 0;
    while (1) {
        int type = Vect_read_next_line(&In, Points, Cats);

        /* Note: check for dead lines is not needed, because they are skipped by
         * V1_read_next_line_nat() */
        if (type == -1) {
            G_fatal_error(_("Unable to read vector map"));
        }
        else if (type == -2) {
            break;
        }
        G_percent(line++, nlines, 4);

        if (type & GV_POINT) {
            if (field != -1 && !Vect_cat_get(Cats, field, NULL))
                continue;

            if (i >= 800) {
                /* Generate some more random numbers */
                myrng(numbers, 1000, rng, p1 - min, p2);
                myrng(numbers2, 1000, rng, p1, p2);
                i = 0;
            }

            G_debug(debuglevel, "x:      %f y:      %f", Points->x[0],
                    Points->y[0]);

            /* perturb */
            /* TODO: tends to concentrate in box corners when min is used */
            if (numbers2[i] >= 0) {
                if (numbers[i] >= 0) {
                    G_debug(debuglevel, "deltax: %f", numbers[i] + min);
                    Points->x[0] += numbers[i++] + min;
                }
                else {
                    G_debug(debuglevel, "deltax: %f", numbers[i] - min);
                    Points->x[0] += numbers[i++] - min;
                }
                Points->y[0] += numbers2[i++];
            }
            else {
                if (numbers[i] >= 0) {
                    G_debug(debuglevel, "deltay: %f", numbers[i] + min);
                    Points->y[0] += numbers[i++] + min;
                }
                else {
                    G_debug(debuglevel, "deltay: %f", numbers[i] - min);
                    Points->y[0] += numbers[i++] - min;
                }
                Points->x[0] += numbers2[i++];
            }

            G_debug(debuglevel, "x_pert: %f y_pert: %f", Points->x[0],
                    Points->y[0]);
        }

        Vect_write_line(&Out, type, Points, Cats);
    }
    G_percent(1, 1, 1);

    /* Copy tables */
    n = Vect_get_num_dblinks(&In);
    ttype = GV_1TABLE;
    if (n > 1)
        ttype = GV_MTABLE;

    for (i = 0; i < n; i++) {
        Fi = Vect_get_dblink(&In, i);
        if (Fi == NULL) {
            G_fatal_error(_("Cannot get db link info"));
        }
        Fin = Vect_default_field_info(&Out, Fi->number, Fi->name, ttype);
        Vect_map_add_dblink(&Out, Fi->number, Fi->name, Fin->table, Fi->key,
                            Fin->database, Fin->driver);

        ret = db_copy_table(Fi->driver, Fi->database, Fi->table, Fin->driver,
                            Vect_subst_var(Fin->database, &Out), Fin->table);
        if (ret == DB_FAILED) {
            G_warning("Cannot copy table");
        }
    }

    Vect_close(&In);

    if (!parm.no_topo->answer)
        Vect_build(&Out);
    Vect_close(&Out);

    return (EXIT_SUCCESS);
}