polyfit.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/matrix/matrix.h>
#include <jwsc/matrix/matrix_io.h>
#include <jwsc/matrix/matrix_math.h>
#include <jwsc/stat/regress.h>


#define  NUM_OPTS_NO_ARG    2
#define  NUM_OPTS_WITH_ARG  4

#define  DEFAULT_SEED    289375
#define  DEFAULT_SIGMA   1.0f
#define  DEFAULT_DEGREE  1

#define  DEFAULT_RANSAC_ENABLED  0
#define  DEFAULT_RANSAC_N        0
#define  DEFAULT_RANSAC_K        0
#define  DEFAULT_RANSAC_T        0
#define  DEFAULT_RANSAC_D        0


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 struct
{
    uint8_t enabled;
    uint32_t n;
    uint32_t k;
    float    t;
    uint32_t d;
}
ransac = { DEFAULT_RANSAC_ENABLED, DEFAULT_RANSAC_N, DEFAULT_RANSAC_K, DEFAULT_RANSAC_T, DEFAULT_RANSAC_D };


static void print_usage(void)
{
    fprintf(stderr, "usage: polyfit 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 <= 0)
    {
        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_ransac_opt(Option_arg arg)
{
    if (arg == NULL)
    {
        return JWSC_EARG("Option 'ransac' requires an argument");
    }
    if (sscanf(arg, "%u,%u,%f,%u", &(ransac.n), &(ransac.k), &(ransac.t), 
                &(ransac.d)) != 4)
    {
        return JWSC_EARG("Option 'ransac' has format n,k,t,d");
    }
    ransac.enabled = 1;
    return NULL;
}




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

    float sse;

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

    Matrix_f* M = NULL;
    Matrix_f* x = NULL;
    Vector_f* y = NULL;
    Vector_f* w = NULL;
    Vector_f* r = 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 = "sigma";
    s_name = 's';
    desc   = "Gaussian sigma.";
    farg   = process_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 = "ransac";
    s_name = 0;
    desc   = "RANSAC robust regression options. Has format n,k,t,d to represent: (n) Number of points to start with, (k) Iterations, (t) Goodness-of-fit threshold for adding a point, and (d) Min number of points needed for a good fit.";
    farg   = process_ransac_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("polyfit", err->msg);
    }

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


    if ((err = read_matrix_fp_f(&M, stdin)) != NULL)
    {
        print_error_msg_exit("polyfit", err->msg);
    }

    assert(M->num_cols == 2);

    create_vector_f(&y, M->num_rows);
    for (row = 0; row < M->num_rows; row++)
    {
        y->elts[ row ] = M->elts[ row ][ 1 ];
    }

    create_matrix_f(&x, M->num_rows, degree+1);
    for (row = 0; row < M->num_rows; row++)
    {
        for (i = 0; i <= degree; i++)
        {
            x->elts[ row ][ i ] = pow(M->elts[ row ][ 0 ], i);
        }
    }

    if (ransac.enabled)
    {
        if ((err = linear_least_squares_ransac_regress_f(&w, x, y, 
                        ransac.n, ransac.k, ransac.t, ransac.d)))
        {
            print_error_msg_exit("polyfit", err->msg);
        }
    }
    else
    {
        if ((err = linear_least_squares_regress_f(&w, x, y)))
        {
            print_error_msg_exit("polyfit", err->msg);
        }
    }

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

    // Compute the SSE of the residuals
    multiply_matrix_and_vector_f(&r, x, w);
    sse = 0;
    for (i = 0; i < y->num_elts; i++)
    {
        sse += (y->elts[ i ] - r->elts[ i ])*(y->elts[ i ] - r->elts[ i ]);
    }

    fprintf(stdout, "SSE: %f\n", sse);


    free_matrix_f(M);
    free_matrix_f(x);
    free_vector_f(y);
    free_vector_f(w);
    free_vector_f(r);

    return EXIT_SUCCESS;
}