text:
/*
 * Evaluation of the scheduling algorithms HET, RTA, RTA2 and RTA3.
 * v1.0 -- 12/01/2017 -- initial version.
 * v2.0 -- 16/05/2017 -- second version.
 */
#include <libxml/xmlreader.h>
#include <libxml/parser.h>
#include <libxml/xpath.h>
#include <math.h>
#include <unistd.h>
#include <stdlib.h>
#include <time.h>
#include <getopt.h>
#include <gsl/gsl_statistics.h>

/*
 * Ceil and floor operations without using the library math, when period
 * and wcet values are integers.
 */
#define U_CEIL( x, y )    ( ( x / y ) + ( x % y != 0 ) )
#define U_FLOOR( x, y )   ( x / y )

/*
 * Elementos en el archivo XML.
 */
#define ELEMENT             1                               // Tag end
#define END_ELEMENT         15                              // Tag start
#define SET_TAG             (const xmlChar*) "Set"          // <Set> tag -- set
#define S_TAG               (const xmlChar*) "S"            // <S> tag -- rts
#define I_TAG               (const xmlChar*) "i"            // <i> tag -- task
#define SET_SIZE_ATTR       (const xmlChar*) "size"         // "size" attribute in <Set> tag
#define SET_RTS_SIZE_ATTR   (const xmlChar*) "n"            // "n" attribute in <Set> tag
#define SET_UF_ATTR         (const xmlChar*) "u"            // "u" attribute in <Set> tag
#define RTS_ID_ATTR         (const xmlChar*) "count"        // "count" attribute in <S> tag
#define RTS_UF_ATTR         (const xmlChar*) "U"            // "U" attribute in <S> tag
#define ID_ATTR             (const xmlChar*) "nro"          // "nro" attribute in <i> tag
#define WCET_ATTR           (const xmlChar*) "C"            // "C" attribute in <i> tag
#define T_ATTR              (const xmlChar*) "T"            // "T" attribute in <i> tag
#define D_ATTR              (const xmlChar*) "D"            // "D" attribute in <i> tag

/*
 * Number of schedulability methods to test.
 */
#define NUM_SCHED_METHODS 4

/*
 * Name of the schedulability methods to evaluate.
 */
#define HET    "het"
#define RTA    "rta"
#define RTA2   "rta2"
#define RTA3   "rta3"

/*
 * Position of the method in the schedulable array.
 */
#define HET_ID    0
#define RTA_ID    1
#define RTA2_ID   2
#define RTA3_ID   3

/*
 * Return value for the schedulability methods.
 */
#define SCHED     1
#define NON_SCHED 0

/*
 * Global variables.
 */
int rts_founded = 0;        // Number of RTS in the XML file evaluated.
int verbose = 0;            // Print addtional info to stderr
FILE* out_file;             // Result file

// Tarea
struct task_t {
    int id;                         // task id
    int c;                          // wcet
    int t;                          // period
    int d;                          // deadline
    int tmc;                        // t - c
    int wcrt[NUM_SCHED_METHODS];    // wcrt
    int cc[NUM_SCHED_METHODS];      // cc
    int loops_w[NUM_SCHED_METHODS]; // number of while loops
    int loops_f[NUM_SCHED_METHODS]; // number of for loops
    int a_rta2;
    int b_rta2;
    int a_rta3;
    int b_rta3;
    int last_psi;           // het
    int last_workload;      // het
};

// rts
struct rts_t {
    int rts_id;
    int rts_uf;
    int rts_ntask;
    int *schedulable;
    struct task_t **tasks;
};

// set of rts
struct set_t {
    int set_size;
    int set_uf;
    int set_rts_ntask;
    struct rts_t **rts_list;
};

// prototipe for scheduling analysis methods
typedef int (*sched_test_method) (struct rts_t*);

// test method result
struct result_t {
    double *cc;
    double *loops;
    double cc_mean;
    double cc_std;
    double loops_mean;
    double loops_std;
};

struct method_t {
    char* method_name;
    int method_id;
    sched_test_method method;
    struct result_t *result;
};

/*
 * Prototipes
 */
int rta_wcrt(struct rts_t*);
int rta2_wcrt(struct rts_t*);
int rta3_wcrt(struct rts_t*);
int het_workload(int i, int b, int n, struct task_t**);
int het_wcrt(struct rts_t*);

int het_workload(int i, int b, int n, struct task_t **tasks)
{
    tasks[n]->loops_w[HET_ID] += 1;

    double tmp = (double) b / (double) tasks[i]->t;
    int f = (int) floor(tmp);
    int c = (int) ceil(tmp);

    tasks[n]->cc[HET_ID] += 2;

    int branch0 = b - f * (tasks[i]->t - tasks[i]->c);
    int branch1 = c * tasks[i]->c;

    if (i > 0) {
        int l_w = tasks[i - 1]->last_workload;
        int tmp = f * tasks[i]->t;
        if (tmp > tasks[i - 1]->last_psi) {
            l_w = het_workload(i - 1, tmp, n, tasks);
        }

        branch0 += l_w;
        branch1 += het_workload(i - 1, b, n, tasks);
    }

    tasks[i]->last_psi = b;

    if (branch0 <= branch1) {
        tasks[i]->last_workload = branch0;
    } else {
        tasks[i]->last_workload = branch1;
    }

    return tasks[i]->last_workload;
}

/*
 * HET
 * "Schedulability Analysis of Periodic Fixed Priority Systems"
 * http://ieeexplore.ieee.org/document/1336766/
 * --
 * See also:
 * "Efficient Exact Schedulability Tests for Fixed Priority Real-Time Systems"
 * http://ieeexplore.ieee.org/document/4487061/
 */
int het_wcrt(struct rts_t *rts)
{
    struct task_t **tasks = rts->tasks;

    int i;
    for (i = 1; i < rts->rts_ntask; i++) {
        tasks[i]->loops_f[HET_ID] += 1;

        int w = het_workload(i - 1, tasks[i]->d, i, tasks);

        if ((w + tasks[i]->c) > tasks[i]->d) {
            rts->schedulable[HET_ID] = NON_SCHED;
            return NON_SCHED;
        }

        tasks[i]->wcrt[HET_ID] = w + tasks[i]->c;
    }

    rts->schedulable[HET_ID] = SCHED;
    return SCHED;
}

/*
 * RTA
 * "Improved Response-Time Analysis Calculations"
 * http://doi.ieeecomputersociety.org/10.1109/REAL.1998.739773
 */
int rta_wcrt(struct rts_t *rts)
{
    struct task_t **tasks = rts->tasks;

    int w = 0;
    int tr = 0;
    int t = tasks[0]->c;
    tasks[0]->wcrt[RTA_ID] = tasks[0]->c;

    int i, j;
    for (i = 1; i < rts->rts_ntask; i++) {
        tr = t + tasks[i]->c;
        tasks[i]->loops_f[RTA_ID] += 1;

        do {
            tasks[i]->loops_w[RTA_ID] += 1;
            t = tr;
            w = tasks[i]->c;

            for (j = 0; j < i; j++) {
                tasks[i]->loops_f[RTA_ID] += 1;

                int c_j = tasks[j]->c;
                int t_j = tasks[j]->t;
                int a = (int) ceil( ((double) tr) / ((double) t_j) );
                tasks[i]->cc[RTA_ID] += 1;

                w = w + (a * c_j);

                if (w > tasks[i]->d) {
                    rts->schedulable[RTA_ID] = NON_SCHED;
                    return NON_SCHED;
                }
            }

            tr = w;

        } while (t != tr);

        tasks[i]->wcrt[RTA_ID] = t;
    }

    rts->schedulable[RTA_ID] = SCHED;
    return SCHED;
}

/*
 * RTA2
 * "Reduced computational cost in the calculation of worst case response time for real time systems"
 * http://sedici.unlp.edu.ar/handle/10915/9654
 */
int rta2_wcrt(struct rts_t *rts)
{
    struct task_t **tasks = rts->tasks;

    int tr = 0;
    int t = tasks[0]->c;
    tasks[0]->wcrt[RTA2_ID] = tasks[0]->c;

    int i, j;
    for (i = 1; i < rts->rts_ntask; i++) {
        tr = t + tasks[i]->c;
        tasks[i]->loops_f[RTA2_ID] += 1;

        do {
            tasks[i]->loops_w[RTA2_ID] += 1;
            t = tr;

            for (j = 0; j < i; j++) {
                tasks[i]->loops_f[RTA2_ID] += 1;

                int a = (int) ceil( ((double) tr) / ((double) tasks[j]->t) );
                tasks[i]->cc[RTA2_ID] += 1;
                a = a * tasks[j]->c;

                if (a > tasks[j]->a_rta2) {
                    tr = tr + a - tasks[j]->a_rta2;
                    tasks[j]->a_rta2 = a;

                    if (tr > tasks[i]->d) {
                        rts->schedulable[RTA2_ID] = NON_SCHED;
                        return NON_SCHED;
                    }
                }
            }
        } while (t != tr);

        tasks[i]->wcrt[RTA2_ID] = t;
    }

    rts->schedulable[RTA2_ID] = SCHED;
    return SCHED;
}

/*
 * RTA3
 * "Computational Cost Reduction for Real-Time Schedulability Tests Algorithms"
 * http://ieeexplore.ieee.org/document/7404899/
 */
int rta3_wcrt(struct rts_t *rts)
{
    struct task_t **tasks = rts->tasks;

    int tr = 0;
    int t = tasks[0]->c;
    tasks[0]->wcrt[RTA3_ID] = tasks[0]->c;

    int i, j;
    for (i = 1; i < rts->rts_ntask; i++) {
        tr = t + tasks[i]->c;
        tasks[i]->loops_f[RTA3_ID] += 1;

        do {
            tasks[i]->loops_w[RTA3_ID] += 1;
            t = tr;

            for (j = i - 1; j >= 0; j--) {
                tasks[i]->loops_f[RTA3_ID] += 1;

                if (tr > tasks[j]->b_rta3) {
                    int a_t = (int) ceil( ((double)tr) / ((double)tasks[j]->t) );
                    tasks[i]->cc[RTA3_ID] += 1;

                    int a = a_t * tasks[j]->c;
                    tr = tr + a - tasks[j]->a_rta3;

                    tasks[j]->a_rta3 = a;
                    tasks[j]->b_rta3 = a_t * tasks[j]->t;

                    // verifica vencimiento
                    if (tr > tasks[i]->d) {
                        rts->schedulable[RTA3_ID] = NON_SCHED;
                        return NON_SCHED;
                    }
                }
            }
        } while (t != tr);

        tasks[i]->wcrt[RTA3_ID] = t;
    }

    rts->schedulable[RTA3_ID] = SCHED;
    return SCHED;
}

void reset_rts(struct rts_t *rts)
{
    int i, j;
    for (i = 0; i < rts->rts_ntask; i++) {
        struct task_t *task = rts->tasks[i];

        task->a_rta2 = task->c;
        task->b_rta2 = task->t;
        task->a_rta3 = task->c;
        task->b_rta3 = task->t;
        task->last_psi = 0;
        task->last_workload = 0;

        for (j = 0; j < NUM_SCHED_METHODS; j++) {
            task->cc[j] = 0;
            task->loops_w[j] = 0;
            task->loops_f[j] = 0;
            task->wcrt[j] = 0;
        }
    }
}

/*
 * Parse the XML file. If a new RTS is found, it is evalutad with the methods in method array.
 */
void processXmlFile(xmlTextReaderPtr reader, struct set_t *rts_set, struct method_t *methods)
{
    xmlChar *name = xmlTextReaderLocalName(reader);

    // Tag <Set> -- initial tag
    if (xmlStrcasecmp(name, SET_TAG) == 0) {
        if (xmlTextReaderNodeType(reader) == ELEMENT) {
            xmlChar *c_set_rts_uf = xmlTextReaderGetAttribute(reader, SET_UF_ATTR);
            xmlChar *c_set_rts_size = xmlTextReaderGetAttribute(reader, SET_SIZE_ATTR);
            xmlChar *c_set_rts_ntask = xmlTextReaderGetAttribute(reader, SET_RTS_SIZE_ATTR);

            rts_set->set_uf = atoi((char*) c_set_rts_uf);
            rts_set->set_size = atoi((char*) c_set_rts_size);
            rts_set->set_rts_ntask = atoi((char*) c_set_rts_ntask);

            xmlFree(c_set_rts_uf);
            xmlFree(c_set_rts_size);
            xmlFree(c_set_rts_ntask);
        }
    }

    // Tag <S> -- RTS
    if (xmlStrcasecmp(name, S_TAG) == 0) {
        if (xmlTextReaderNodeType(reader) == ELEMENT) {
            xmlChar *c_rts_id = xmlTextReaderGetAttribute(reader, RTS_ID_ATTR);
            xmlChar *c_rts_uf = xmlTextReaderGetAttribute(reader, RTS_UF_ATTR);

            // reserve memory for the rts and add it to rts_list
            struct rts_t *new_rts = malloc(sizeof(struct rts_t));

            // reserve memory for the methods results
            new_rts->schedulable = malloc(sizeof(int) * NUM_SCHED_METHODS);
            // reserve memory for the rts tasks
            new_rts->tasks = malloc(sizeof(struct task_t*) * rts_set->set_rts_ntask);

            // complete data about this rts
            new_rts->rts_id = atoi((char*) c_rts_id);
            new_rts->rts_uf = atoi((char*) c_rts_uf);
            new_rts->rts_ntask = rts_set->set_rts_ntask;

            // add rts to set
            rts_set->rts_list[rts_founded] = new_rts;

            // free memory
            xmlFree(c_rts_id);
            xmlFree(c_rts_uf);
        }

        if (xmlTextReaderNodeType(reader) == END_ELEMENT) {
            struct rts_t *rts = rts_set->rts_list[rts_founded];
            reset_rts(rts);

            // evaluate the methods
            int i, j;
            for (i = 0; i < NUM_SCHED_METHODS; i++) {
                int method_id = methods[i].method_id;
                rts->schedulable[method_id] = (*methods[i].method)(rts);

                // store totals
                for (j = 0; j < rts->rts_ntask; j++) {
                    struct task_t *task = rts->tasks[j];
                    methods[i].result->cc[rts_founded] += task->cc[method_id];
                    methods[i].result->loops[rts_founded] += task->loops_w[method_id] + task->loops_f[method_id];
                }
            }

            rts_founded = rts_founded + 1;
        }
    }

    // Tag <i> -- a real-time task
    if (xmlStrcasecmp(name, I_TAG) == 0) {
        xmlChar *c_id = xmlTextReaderGetAttribute(reader, ID_ATTR);
        xmlChar *wcet = xmlTextReaderGetAttribute(reader, WCET_ATTR);
        xmlChar *t = xmlTextReaderGetAttribute(reader, T_ATTR);
        xmlChar *d = xmlTextReaderGetAttribute(reader, D_ATTR);

        int id = atoi((char*) c_id) - 1;

        struct rts_t *rts = rts_set->rts_list[rts_founded];

        // reserve memory for the task
        struct task_t *task = malloc(sizeof(struct task_t));

        // complete the basic task data
        task->id = id + 1;
        task->c = atoi((char*) wcet);
        task->t = atoi((char*) t);
        task->d = atoi((char*) d);
        task->tmc = task->t - task->c;

        // add the task to rts
        rts->tasks[id] = task;

        xmlFree(c_id);
        xmlFree(wcet);
        xmlFree(t);
        xmlFree(d);
    }

    xmlFree(name);
}

void testRtsInXml(char *file, struct set_t* rts_set, struct method_t *methods, int limit)
{
    // get read pointer
    xmlTextReaderPtr reader = xmlNewTextReaderFilename(file);
    if (reader == NULL) {
        fprintf(stderr, "Unable to open %s\n", file);
        exit(EXIT_FAILURE);
    }

    // parse xml file and evaluate schedulability methods
    int ret = xmlTextReaderRead(reader);
    while (ret == 1) {
        if (limit > 0 && rts_founded == limit) {
            break;
        }
        processXmlFile(reader, rts_set, methods);
        ret = xmlTextReaderRead(reader);
    }

    if (rts_set->set_size < limit) {
        fprintf(stderr, "Warning: %d str in file according to XML info, but %d to be tested.\n", rts_set->set_size, limit);
    }

    xmlFreeTextReader(reader);
}

/*
 * Print method results to out_file.
 */
void save_result(char* method, struct result_t* result, int use_csv, char *csv_sep)
{
    if (use_csv == 0) {
        fprintf(out_file, "%10s%15f%15f%15f%15f\n", method, result->cc_mean, result->cc_std,
                                                    result->loops_mean, result->loops_std);
    } else {
        fprintf(out_file, "%2$s%1$s%3$f%1$s%4$f%1$s%5$f%1$s%6$f\n", csv_sep, method, result->cc_mean, result->cc_std,
                                                                                result->loops_mean, result->loops_std);
    }
}

/*
 * Print help and usage information.
 */
void printUsage(char* progName, int exitCode)
{
    fprintf(stderr, "Usage: %s [options] file\n", progName);
    fprintf(stderr,
            "\t-v  --verbose\tDisplay additional information about the clock used.\n"
            "\t-h  --help\tDisplay this information.\n"
            "\t-l  --limit\tTest first n RTS in file.\n"
            "\t-c  --csv\tCSV output with specified line separator.\n");
    exit(exitCode);
}

int main(int argc, char **argv)
{
    int i, j, k;

    if (argc <= 1) {
        printUsage(argv[0], EXIT_FAILURE);
    }

    // options -- short format
    const char *shortOpts = "hvl:c:";
    // options -- long format
    const struct option longOpts[] = {
        {"help",    no_argument,        NULL, 'h'},
        {"verbose", no_argument,        NULL, 'v'},
        {"limit",   required_argument,  NULL, 'l'},
        {"csv",     required_argument,  NULL, 'c'},
        {0, 0, 0, 0}
    };

    verbose = 0;
    rts_founded = 0;
    int limit = 0;

    int use_csv = 0;
    char* csv_sep;

    int nextOption;

    do {
        nextOption = getopt_long(argc, argv, shortOpts, longOpts, NULL);
        switch (nextOption) {
            case 'h': // -h or --help
                printUsage(argv[0], EXIT_SUCCESS);
            case 'v': // -v or --verbose
                verbose = 1;
                break;
            case 'l': // -l or --limit
                limit = atoi(optarg);
                break;
            case 'c': // -c or --csv
                use_csv = 1;
                csv_sep = optarg;
                break;
            case '?': // invalid option
                printUsage(argv[0], EXIT_FAILURE);
            case -1: // no more options
                break;
            default:
                abort();
        }
    } while (nextOption != -1);

    // print info to stderr if requested
    if (verbose == 1) {
        fprintf(stderr, "Testing %s...\n", argv[optind]);
    }

    // stdout as default output file
    out_file = stdout;

    // reserve memory for the set
    struct set_t *rts_set = malloc(sizeof(struct set_t));
    rts_set->set_uf = 0;
    rts_set->set_size = 0;
    rts_set->set_rts_ntask = 0;
    // reserve memory to store the rts to test
    rts_set->rts_list = malloc(sizeof(struct rts_t*) * limit);

    // reserve memory for test results
    struct result_t* het_results = malloc(sizeof(struct result_t));
    struct result_t* rta_results = malloc(sizeof(struct result_t));
    struct result_t* rta2_results = malloc(sizeof(struct result_t));
    struct result_t* rta3_results = malloc(sizeof(struct result_t));

    // methods to test
    struct method_t methods[] = {[RTA_ID]  {RTA,  RTA_ID,  rta_wcrt,  rta_results},
                                 [RTA2_ID] {RTA2, RTA2_ID, rta2_wcrt, rta2_results},
                                 [RTA3_ID] {RTA3, RTA3_ID, rta3_wcrt, rta3_results},
                                 [HET_ID]  {HET,  HET_ID,  het_wcrt,  het_results}
                                 };

    for (i = 0; i < NUM_SCHED_METHODS; i++) {
        methods[i].result->cc = calloc(limit, sizeof(double));
        methods[i].result->loops = calloc(limit, sizeof(double));
    }

    if (verbose == 1) {
        fprintf(stderr, "Testing %d rts.\n", limit);
    }

    // read rts from xml file into rts_set
    char *filename = argv[optind];
    testRtsInXml(filename, rts_set, methods, limit);

    // compute means and stdev
    for (i = 0; i < NUM_SCHED_METHODS; i++) {
        methods[i].result->cc_mean = gsl_stats_mean(methods[i].result->cc, 1, rts_founded);
        methods[i].result->cc_std = gsl_stats_sd_m(methods[i].result->cc, 1, rts_founded, methods[i].result->cc_mean);
        methods[i].result->loops_mean = gsl_stats_mean(methods[i].result->loops, 1, rts_founded);
        methods[i].result->loops_std = gsl_stats_sd_m(methods[i].result->loops, 1, rts_founded, methods[i].result->loops_mean);
    }

    int rts_sched_cnt = 0;
    int rts_nonsched_cnt = 0;

    // verify that all methods give the same result
    for (i = 0; i < rts_founded; i++) {
        struct rts_t *rts = rts_set->rts_list[i];

        int sum = 0;
        for (j = 0; j < NUM_SCHED_METHODS; j++) {
            sum += rts->schedulable[j];
        }

        if (sum > 0 && sum < NUM_SCHED_METHODS) {
            fprintf(stderr, "Error! Method results are not the same. RTS %d\n", i);
            for (j = 0; j < NUM_SCHED_METHODS; j++) {
                fprintf(stderr, "%s: %d\n", methods[j].method_name, rts->schedulable[j]);
            }
            break;
        }

        if (rts->schedulable[RTA_ID] == SCHED) {
            rts_sched_cnt += 1;
        } else {
            rts_nonsched_cnt += 1;
        }
    }

    // verify that all wcrt are the same (only RTA methods)
    for (j = 0; j < rts_founded; j++) {
        struct rts_t *rts = rts_set->rts_list[j];
        for (k = 0; k < rts->rts_ntask; k++) {
            struct task_t *task = rts->tasks[k];
            int ref_wcrt = task->wcrt[RTA_ID];
            if (ref_wcrt != task->wcrt[RTA2_ID] || ref_wcrt != task->wcrt[RTA3_ID])
            {
                fprintf(stderr, "Error! WCRT are not the same. RTS %d, task %d\n", j, k);

                int i = 0;
                fprintf(stderr, "%13s%10s%10s%10s%10s\n", "RTA", "RTA2", "RTA3", "C_i", "D_i");
                for (i = 0; i < rts->rts_ntask; i++) {
                    fprintf(stderr, "%3d%10d%10d%10d%10d%10d\n", i, rts->tasks[i]->wcrt[RTA_ID], rts->tasks[i]->wcrt[RTA2_ID],
                                                                    rts->tasks[i]->wcrt[RTA3_ID], rts->tasks[i]->c, rts->tasks[i]->d );
                }

                return(EXIT_FAILURE);
            }
        }
    }

    // get timestamp for the test
    char test_date[50];
    time_t current_time = time(NULL);
    strftime(test_date, 50, "%R %d/%m/%Y", localtime(&current_time));

    // print header
    fprintf(out_file, "%s\n", filename);
    fprintf(out_file, "%s\n", test_date);
    fprintf(out_file, "Total: %d\n", rts_founded);
    fprintf(out_file, "Sched: %d\n", rts_sched_cnt);
    fprintf(out_file, "Non sched: %d\n", rts_nonsched_cnt);

    // print column names
    if (use_csv == 0) {
        fprintf(out_file, "%10s%15s%15s%15s%15s\n", "method", "cc_mean", "cc_mean_std", "loops_mean", "loops_mean_std");
    } else {
        fprintf(out_file, "method%1$scc_mean%1$scc_mean_std%1$sloops_mean%1$sloops_mean_std\n", csv_sep);
    }

    // print the results
    for (i = 0; i < NUM_SCHED_METHODS; i++) {
        save_result(methods[i].method_name, methods[i].result, use_csv, csv_sep);
    }

    return(EXIT_SUCCESS);
}

max_stars_repo_name: unpsjb-rtsg/wcrt

max_issues_repo_name: unpsjb-rtsg/wcrt

max_forks_repo_name: unpsjb-rtsg/wcrt


text:
#include <stdio.h>
#include <unistd.h>
#include <assert.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include "dlib.h"
#include "svec.h"
#include "rng.h"

const char *usage = "Usage: scode [OPTIONS] < file\n"
  "file should have columns of arbitrary tokens\n"
  "-r RESTART: number of restarts (default 1)\n"
  "-i NITER: number of iterations over data (default UINT32_MAX)\n"
  "-t THRESHOLD: quit if logL increase for iter <= this (default .001)\n"
  "-d NDIM: number of dimensions (default 25)\n"
  "-z Z: partition function approximation (default 0.166)\n"
  "-p PHI0: learning rate parameter (default 50.0)\n"
  "-u ETA0: learning rate parameter (default 0.2)\n"
  "-s SEED: random seed (default 0)\n"
  "-c calculate real Z (default false)\n"
  "-w The first line of the input is weights (default false)\n"
  "-v verbose messages (default false)\n";

//typedef uint32_t u32;
//typedef uint64_t u64;
u32 RESTART = 1;
u32 NITER = UINT32_MAX;
double THRESHOLD = 0.001;
u32 NDIM = 25;
double Z = 0.166;
double PHI0 = 50.0;
double ETA0 = 0.2;
unsigned long int SEED = 0;
bool CALCZ = false;
bool WEIGHT = false;
bool VERBOSE = false;

u32 NTOK = 0;
u64 NTUPLE = 0;

const gsl_rng_type *rng_T;
gsl_rng *rng_R = NULL;
darr_t data;
u64 **update_cnt;
double * weight = NULL;
double * uweight = NULL; /*Updated weights*/
u64 **cnt;
#define frq(i,j) ((double)cnt[i][j]*NTOK/len(data))
svec **vec;
svec **best_vec;
svec dummy_vec;
sym_t qmax;
sym_t NULLFEATID;
#define NULLFEATMARKER "/XX/"

int main(int argc, char **argv);
void init_rng();
void free_rng();
u64 init_data();
u32 init_weight();
void free_weight();
void randomize_vectors();
void copy_best_vec();
void free_data();
void update_tuple(sym_t *t);
double logL();
double calcZ();

#define vmsg(...) if(VERBOSE)msg(__VA_ARGS__)

int main(int argc, char **argv) {
  int opt;
  while((opt = getopt(argc, argv, "r:i:t:d:z:p:u:s:cwv")) != -1) {
    switch(opt) {
    case 'r': RESTART = atoi(optarg); break;
    case 'i': NITER = atoi(optarg); break;
    case 't': THRESHOLD = atof(optarg); break;
    case 'd': NDIM = atoi(optarg); break;
    case 'z': Z = atof(optarg); break;
    case 'p': PHI0 = atof(optarg); break;
    case 'u': ETA0 = atof(optarg); break;
    case 's': SEED = atoi(optarg); break;
    case 'c': CALCZ = true; break;
    case 'w': WEIGHT = true; break;
    case 'v': VERBOSE = true; break;
    default: die("%s",usage);
    }
  }

  vmsg("scode -r %u -i %u -t %g -d %u -z %g -p %g -u %g -s %lu %s%s%s",
       RESTART, NITER, THRESHOLD, NDIM, Z, PHI0, ETA0, SEED,
       (CALCZ ? "-c " : ""), (WEIGHT ? "-w " : ""), (VERBOSE ? "-v " : ""));

  init_rng();
  if (SEED) gsl_rng_set(rng_R, SEED);
  if (WEIGHT) NTOK = init_weight();
  NTUPLE = init_data();
  vmsg("Read %zu tuples %u uniq tokens", NTUPLE, qmax);

  double best_logL = 0;
  for (u32 start = 0; start < RESTART; start++) {
    randomize_vectors();
    double ll = logL();
    vmsg("Restart %u/%u logL0=%g best=%g", 1+start, RESTART, ll, best_logL);
    if (CALCZ) vmsg("Z=%g (approx %g)", calcZ(), Z);
    for (u32 iter = 0; iter < NITER; iter++) {
      for (u64 di = 0; di < NTUPLE; di++) {
	update_tuple(&val(data, di * NTOK, sym_t));
      }
      double ll0 = ll;
      ll = logL();
      vmsg("Iteration %u/%u logL=%g", 1+iter, NITER, ll);
      if (ll - ll0 <= THRESHOLD) break;
    }
    if (start == 0 || ll > best_logL) {
      vmsg("Updating best_vec with logL=%g", ll);
      best_logL = ll;
      copy_best_vec();
    }
    vmsg("Restart %u/%u logL1=%g best=%g", 1+start, RESTART, ll, best_logL);
    if (CALCZ) vmsg("Z=%g (approx %g)", calcZ(), Z);
  }
  for (u32 t = 0; t < NTOK; t++) {
    for (sym_t q = 1; q <= qmax; q++) {
      if (best_vec[t][q] == NULL) continue;
      printf("%u:%s\t%zu\t", t, sym2str(q), cnt[t][q]);
      svec_print(best_vec[t][q]);
      putchar('\n');
    }
  }
  fflush(stdout);
  free_data();
  free_rng();
  if (WEIGHT) free_weight();
  symtable_free();
  dfreeall();
  fprintf(stderr, "%f\n", best_logL);
  vmsg("bye");
}

double logL() {
  double l = 0;
  for (u64 i = 0; i < NTUPLE; i++) {
    sym_t *t = &val(data, i * NTOK, sym_t);
    sym_t x = t[0];
    sym_t y = t[1];
    float px = frq(0, x);
    float py = frq(1, y);
    svec vx = vec[0][x];
    svec vy = vec[1][y];
    float xy = svec_sqdist(vx, vy);
    l += log(px * py) - xy;
  }
  return (l / NTUPLE - log(Z));
}

double calcZ() {
  double z = 0;
  for (sym_t x = 1; x <= qmax; x++) {
    if (VERBOSE && (x % 1000 == 0)) fputc('.', stderr);
    if (cnt[0][x] == 0) continue;
    float px = frq(0, x);
    svec vx = vec[0][x];
    for (sym_t y = 1; y <= qmax; y++) {
      if (cnt[1][y] == 0) continue;
      float py = frq(1, y);
      svec vy = vec[1][y];
      float xy = svec_sqdist(vx, vy);
      z += px * py * exp(-xy);
    }
  }
  if (VERBOSE) fputc('\n', stderr);
  return z;
}

void update_tuple(sym_t *t) {
  /*weighted update*/
  static svec *u = NULL;
  static svec *v = NULL;
  static svec dx = NULL;
  if (u == NULL) u = _d_malloc(NTOK * sizeof(svec));
  if (v == NULL) v = _d_malloc(NTOK * sizeof(svec));
  if (dx == NULL) dx = svec_alloc(NDIM);
  for (u32 i = 0; i < NTOK; i++) u[i] = vec[i][t[i]];
  for (u32 i = 0; i < NTOK; i++) {
    /* Sampling values from the marginal distributions. */
    /* Can this be done once, or do we have to resample for every x? */
    if(i > 0 && t[i] == NULLFEATID) continue;
    for (u32 j = 0; j < NTOK; j++) {
      if (j==i) { v[j] = u[i]; continue;}
      u64 r = gsl_rng_get(rng_R);
      r = (r << 32) | gsl_rng_get(rng_R);
      r = r % NTUPLE;
      sym_t y = val(data, r * NTOK + j, sym_t);
      v[j] = vec[j][y];
      if(i > 0) break;
    }
    /* Compute the move for u[i] */
    svec_set_zero(dx);
    double ww;
    for (u32 j = 0; j < NTOK; j++) {
      if (j == i) continue;
      ww = weight == NULL ? 1 : (i > 0 ? weight[i] : weight[j]);
      double push = 0, pull = 0;
      if (v[j] == NULL) v[j] = dummy_vec;
      else push = exp(-svec_sqdist(u[i], v[j])) / Z;
      if(u[j] == NULL)  u[j] = dummy_vec;
      else pull = 1;
      if(push != 0 || pull != 0){
	for (u32 d = 0; d < NDIM; d++) {
	  float dxd = svec_get(dx, d);
	  float x = svec_get(u[i], d);
	  float y = svec_get(u[j], d);
	  float z = svec_get(v[j], d);
	  svec_set(dx, d, dxd + ww * ( pull * (y - x) + push * (x - z)));
	}
      }
      /*restore the vectors to original forms*/
      if(push == 0) v[j] = NULL;
      if(pull == 0) u[j] = NULL;
      if(i > 0) break;
      }
    /* Apply the move scaled by learning parameter */
    u64 cx = update_cnt[i][t[i]]++;
    float nx = ETA0 * (PHI0 / (PHI0 + cx));
    svec_scale(dx, nx);
    svec_add(u[i], dx);
    svec_normalize(u[i]);
  }
}

u32 init_weight(){
  u32 size = 100, i = 0;
  weight = _d_malloc(size * sizeof(double));
  forline (buf, NULL) {
    fortok (tok, buf) {
      weight[i] = atof(tok);
      assert(weight[i++] >= 0);
      if(i >= 100) {
	size *= 2;
	weight = _d_realloc(weight, size);
      }
    }
    assert(i > 0);
    break;
  }
  return i;
}

void free_weight() {
  if (weight != NULL) _d_free(weight);
}

u64 init_data() {
  qmax = 0;
  data = darr(0, sym_t);
  forline (buf, NULL) {
    u32 ntok = 0;
    fortok (tok, buf) {
      sym_t q = str2sym(tok, true);
      if (q > qmax) qmax = q;
      size_t lendata = len(data);
      val(data, lendata, sym_t) = q;
      if(strcmp(tok, NULLFEATMARKER) == 0) NULLFEATID = q;
      ntok++;
    }
    if(NTOK == 0) NTOK = ntok;
    assert(ntok == NTOK); //Each line has equal number of tokens
  }
  assert(NTOK > 0);
  update_cnt = _d_malloc(NTOK * sizeof(ptr_t));
  cnt = _d_malloc(NTOK * sizeof(ptr_t));
  vec = _d_malloc(NTOK * sizeof(ptr_t));
  best_vec = _d_malloc(NTOK * sizeof(ptr_t));
  dummy_vec = svec_alloc(NDIM);
  svec_zero(dummy_vec);
  uweight = _d_calloc(NTOK, sizeof(double));
  for (u32 i = 0; i < NTOK; i++) {
    update_cnt[i] = _d_calloc(qmax+1, sizeof(u64));
    cnt[i] = _d_calloc(qmax+1, sizeof(u64));
    vec[i] = _d_calloc(qmax+1, sizeof(svec));
    best_vec[i] = _d_calloc(qmax+1, sizeof(svec));
  }
  u64 N = len(data) / NTOK;
  for (u64 i = 0; i < N; i++) {
    sym_t *p = &val(data, i * NTOK, sym_t);
    for (u32 j = 0; j < NTOK; j++) {
      sym_t k = p[j];
      assert(k <= qmax);
      cnt[j][k]++;
      if(k == NULLFEATID){
	vec[j][k] = best_vec[j][k] = NULL;
      }
      else if (vec[j][k] == NULL) {
	vec[j][k] = svec_alloc(NDIM);
	best_vec[j][k] = svec_alloc(NDIM);
      }
    }
  }
  return N;
}

void free_data() {
  for (u32 i = 0; i < NTOK; i++) {
    for (sym_t j = 0; j <= qmax; j++) {
      if (vec[i][j] != NULL) {
	svec_free(vec[i][j]);
	svec_free(best_vec[i][j]);
      }
    }
    _d_free(best_vec[i]);
    _d_free(vec[i]);
    _d_free(cnt[i]);
    _d_free(update_cnt[i]);
  }
  _d_free(uweight);
  svec_free(dummy_vec);
  _d_free(best_vec);
  _d_free(vec);
  _d_free(cnt);
  _d_free(update_cnt);
  darr_free(data);
}

void randomize_vectors() {
  for (u32 j = 0; j < NTOK; j++) {
    for (sym_t q = 1; q <= qmax; q++) {
      if (vec[j][q] != NULL) {
	svec_randomize(vec[j][q]);
	update_cnt[j][q] = 0;
      }
    }
  }
}

void copy_best_vec() {
  for (u32 j = 0; j < NTOK; j++) {
    for (sym_t q = 1; q <= qmax; q++) {
      if (vec[j][q] != NULL) {
	svec_memcpy(best_vec[j][q], vec[j][q]);
      }
    }
  }
}

void init_rng() {
  gsl_rng_env_setup();
  rng_T = gsl_rng_mt19937;
  rng_R = gsl_rng_alloc(rng_T);
}

void free_rng() {
  gsl_rng_free(rng_R);
}

ext: c

max_stars_repo_name: ai-ku/scode

max_issues_repo_name: ai-ku/scode

max_forks_repo_name: ai-ku/scode


text:
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include "utils.h"

int main() {
    int success_flag = 1;
    char * data_path = "../Test_Data/linear.csv";

    int cols = column_count(data_path);
    int rows = row_count(data_path) - 1;
    double data[rows][cols];

    int success = read_data(data_path, rows, cols, data);
    if(success != 0){
        return 1;
    }

    double *x = (double *) malloc(sizeof(double)*rows);
    double *y = (double *) malloc(sizeof(double)*rows);

    fetch_column(x, data, 1, rows, cols);
    fetch_column(y, data, 2, rows, cols);

    double xTx[rows];

    /* optimization: reuse x as read from data, dont read for each vector used */
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, 1, 1, rows, 1, x, rows, x, rows, 0.0, xTx, 1 );

    double xTy[rows];

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, 1, 1, rows, 1, x, rows, y, rows, 0.0, xTy, 1 );

    printf("C: %.30f\n", xTx[0]);
    printf("C: %.30f\n", xTy[0]);

    free(x);
    free(y);
    return success_flag;
}