polygen.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <config.h>

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <math.h>

#include <jwsc/base/error.h>
#include <jwsc/base/option.h>
#include <jwsc/vector/vector.h>
#include <jwsc/vector/vector_io.h>
#include <jwsc/matrix/matrix.h>
#include <jwsc/matrix/matrix_io.h>
#include <jwsc/matrix/matrix_math.h>
#include <jwsc/prob/pdf.h>
#include <jwsc/prob/mv_pdf.h>
#include <jwsc/stat/regress.h>


#define  NUM_OPTS_NO_ARG    2
#define  NUM_OPTS_WITH_ARG  9

#define  DEFAULT_SEED             289375
#define  DEFAULT_SIGMA            1.0f
#define  DEFAULT_DEGREE           1
#define  DEFAULT_NUM_PTS          100
#define  DEFAULT_W_MU             0.0
#define  DEFAULT_W_SIGMA          0.1
#define  DEFAULT_W_MIN           -2
#define  DEFAULT_W_MAX            2
#define  DEFAULT_X_MIN           -10
#define  DEFAULT_X_MAX            10
#define  DEFAULT_OUTLIER          0
#define  DEFAULT_OUTLIER_SIGMA    6.0f


static Option_no_arg opts_no_arg[NUM_OPTS_NO_ARG];

static Option_with_arg opts_with_arg[NUM_OPTS_WITH_ARG];

static uint32_t seed = DEFAULT_SEED;

static float sigma = DEFAULT_SIGMA;

static uint32_t degree = DEFAULT_DEGREE;

static uint32_t num_pts = DEFAULT_NUM_PTS;

static float w_mu = DEFAULT_W_MU;

static float w_sigma = DEFAULT_W_SIGMA;

static float w_min = DEFAULT_W_MIN;

static float w_max = DEFAULT_W_MAX;

static float x_min = DEFAULT_X_MIN;

static float x_max = DEFAULT_X_MAX;

static float outlier = DEFAULT_OUTLIER;

static float outlier_sigma = DEFAULT_OUTLIER_SIGMA;

static enum { GAUSS=0, UNIFORM } w_sample_t;


static void print_usage(void)
{
    fprintf(stderr, "usage: polygen OPTIONS\n");
    fprintf(stderr, "where OPTIONS is one or more of the following:\n");
    print_options(stderr, 25, NUM_OPTS_NO_ARG, opts_no_arg, NUM_OPTS_WITH_ARG, 
            opts_with_arg);
}

static Error* process_help_opt(void)
{
    print_usage();
    exit(EXIT_SUCCESS);
}

static Error* process_version_opt(void)
{
    fprintf(stderr, "%s\n", POLYFIT_PACKAGE_STRING);
    exit(EXIT_SUCCESS);
}

static Error* process_seed_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'seed' requires an argument");
    }
    if (sscanf(arg, "%u", &seed) != 1)
    {
        return JWSC_EARG("Option 'seed' requires an argument");
    }
    return NULL;
}

static Error* process_sigma_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'sigma' requires an argument");
    }
    if (sscanf(arg, "%f", &sigma) != 1 || sigma <= 1.0e-6)
    {
        return JWSC_EARG("Option 'sigma' must be > 0");
    }
    return NULL;
}

static Error* process_degree_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'degree' requires an argument");
    }
    if (sscanf(arg, "%u", &degree) != 1)
    {
        return JWSC_EARG("Option 'degree' requires an argument");
    }
    return NULL;
}

static Error* process_num_pts_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'num-pts' requires an argument");
    }
    if (sscanf(arg, "%u", &num_pts) != 1)
    {
        return JWSC_EARG("Option 'num-pts' requires an argument");
    }
    return NULL;
}

static Error* process_w_gauss_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'w-gauss' requires an argument");
    }
    if (sscanf(arg, "%f,%f", &w_mu, &w_sigma) != 2)
    {
        return JWSC_EARG("Option 'w-gauss' has format mu,sigma");
    }
    if (w_sigma < 0.01)
    {
        return JWSC_EARG("Option 'w-uni': min > max");
    }
    w_sample_t = GAUSS;
    return NULL;
}

static Error* process_w_uni_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'w-uni' requires an argument");
    }
    if (sscanf(arg, "%f,%f", &w_min, &w_max) != 2)
    {
        return JWSC_EARG("Option 'w-uni' has format min,max");
    }
    if (w_min > w_max)
    {
        return JWSC_EARG("Option 'w-uni': min > max");
    }
    w_sample_t = UNIFORM;
    return NULL;
}

static Error* process_x_range_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'x-range' requires an argument");
    }
    if (sscanf(arg, "%f,%f", &x_min, &x_max) != 2)
    {
        return JWSC_EARG("Option 'x-range' has format min,max");
    }
    if (x_min > x_max)
    {
        return JWSC_EARG("Option 'x-range': min > max");
    }
    return NULL;
}

static Error* process_outlier_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'outlier' requires an argument");
    }
    if (sscanf(arg, "%f", &outlier) != 1 || outlier < 0 || outlier > 1)
    {
        return JWSC_EARG("Option 'outlier' must be in [0,1]");
    }
    return NULL;
}

static Error* process_outlier_sigma_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'outlier-sigma' requires an argument");
    }
    if (sscanf(arg, "%f", &outlier_sigma) != 1 || outlier_sigma <= 1.0e-6)
    {
        return JWSC_EARG("Option 'outlier-sigma' must be > 0");
    }
    return NULL;
}


static void sample_w(Vector_d** w_out)
{
    uint32_t d;
    Vector_d* w_mu_v = 0;
    Matrix_d* w_sigma_m = 0;

    if (w_sample_t == GAUSS)
    {
        create_init_vector_d(&w_mu_v, degree+1, w_mu);
        create_identity_matrix_d(&w_sigma_m, degree+1);
        multiply_matrix_by_scalar_d(&w_sigma_m, w_sigma_m, w_sigma);
        assert(sample_mv_gaussian_pdf_d(w_out, w_mu_v, w_sigma_m) == 0);
    }
    else if (w_sample_t == UNIFORM)
    {
        create_vector_d(w_out, degree+1);
        for (d = 0; d <= degree; d++)
        {
            (*w_out)->elts[ d ] = sample_uniform_pdf_d(w_min, w_max);
        }
    }
    else
    {
        abort();
    }

    free_vector_d(w_mu_v);
    free_matrix_d(w_sigma_m);
}


int main(int argc, const char** argv)
{
    int      argi;
    uint32_t i, d;

    Option_string        l_name;
    Option_flag          s_name;
    Option_desc          desc;
    Option_func_no_arg   fnoarg;
    Option_func_with_arg farg;

    Vector_d* w = NULL;
    Vector_d* x = NULL;
    Vector_d* y = NULL;
    Error*    err;

    i      = 0;
    l_name = "help";
    s_name = 'h';
    desc   = "Prints program usage.";
    fnoarg = process_help_opt;
    init_option_no_arg(&(opts_no_arg[i++]), l_name, s_name, desc, fnoarg);

    l_name = "version";
    s_name = 'v';
    desc   = "Prints program version.";
    fnoarg = process_version_opt;
    init_option_no_arg(&(opts_no_arg[i++]), l_name, s_name, desc, fnoarg);

    i      = 0;
    l_name = "seed";
    s_name = 0;
    desc   = "Random seed.";
    farg   = process_seed_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "w-gauss";
    s_name = 0;
    desc   = "Gaussian over w to generate points from. Has format mu,sigma.";
    farg   = process_w_gauss_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "w-uni";
    s_name = 0;
    desc   = "Uniform range over w to generate points from. Has format min,max.";
    farg   = process_w_uni_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "x-range";
    s_name = 0;
    desc   = "Range over x to generate points from.";
    farg   = process_x_range_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "sigma";
    s_name = 0;
    desc   = "Gaussian sigma for polynomial point generation.";
    farg   = process_sigma_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "outlier";
    s_name = 0;
    desc   = "Percentage of points to generate as outliers.";
    farg   = process_outlier_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "outlier-sigma";
    s_name = 0;
    desc   = "Gaussian sigma for polynomial outlier point generation.";
    farg   = process_outlier_sigma_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "degree";
    s_name = 'd';
    desc   = "Polynomial degree.";
    farg   = process_degree_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    l_name = "num-pts";
    s_name = 'n';
    desc   = "Number of points to generate.";
    farg   = process_num_pts_opt;
    init_option_with_arg(&(opts_with_arg[i++]), l_name, s_name, desc, farg);

    if ((err = process_options(argc, argv, &argi, NUM_OPTS_NO_ARG, opts_no_arg, 
                    NUM_OPTS_WITH_ARG, opts_with_arg)) != NULL)
    {
        print_error_msg_exit("polygen", err->msg);
    }

    srand(seed); rand(); rand();

    sample_w(&w);
    create_vector_d(&x, num_pts);
    create_vector_d(&y, num_pts);

    for (i = 0; i < num_pts; i++)
    {
        x->elts[ i ] = sample_uniform_pdf_d(x_min, x_max);

        if (outlier > sample_uniform_pdf_d(0, 1))
        {
            y->elts[ i ] = sample_gaussian_pdf_d(0, outlier_sigma, 
                    -6*outlier_sigma, 6*outlier_sigma);
            fprintf(stderr, "outlier[ %2d ]: %f, %f\n", i, x->elts[ i ], y->elts[ i ]);
        }
        else
        {
            y->elts[ i ] = sample_gaussian_pdf_d(0, sigma, -6*sigma, 6*sigma);
        }


        for (d = 0; d <= degree; d++)
        {
            y->elts[ i ] += w->elts[ d ] * pow(x->elts[ i ], d);
        }

        fprintf(stdout, "%13.6e  %13.6e\n", x->elts[ i ], y->elts[ i ]);
    }

    fprintf(stderr, "y(x) = %.6e", w->elts[ 0 ]);
    for (d = 1; d <= degree; d++)
    {
        fprintf(stderr, " + %.6e*x**%d", w->elts[ d ], d);
    }
    fprintf(stderr, "\n");


    free_vector_d(w);
    free_vector_d(x);
    free_vector_d(y);

    return EXIT_SUCCESS;
}