text:
#include "np.h"
#include "main.h"
#include "argv.h"
#include "memory.h"
#include "utf8.h"

#include "module_categorical.h"
#include "module_lde.h"

#include <stdlib.h>
#include <string.h>

#include <math.h>
#include <float.h>

#include <gsl/gsl_sf.h>
#include "gslsupp.h"

#include "sort.h"

#ifndef TEST_DEACTIVATE

#   include "test.h"
#   include "test_ll.h"
#   include "test_np.h"
#   include "test_sort.h"
#   include "test_utf8.h"

static bool test_main(struct log *log)
{
    const struct test_group group_arr[] = {
        {
            NULL,
            sizeof(struct test_ll_a),
            CLII((test_generator_callback[]) {
                test_ll_generator_a,
            }),
            CLII((test_callback[]) {
                test_ll_a_1,
                test_ll_a_2,
                test_ll_a_3
            })
        },
        {
            NULL,
            sizeof(struct test_ll_a),
            CLII((test_generator_callback[]) {
                test_ll_generator_b,
            }),
            CLII((test_callback[]) {
                test_ll_b,
            })
        },
        {
            test_np_disposer_a,
            sizeof(struct test_np_a),
            CLII((test_generator_callback[]) {
                test_np_generator_a,
            }),
            CLII((test_callback[]) {
                test_np_a,
            })
        },
        {
            test_sort_disposer_a,
            sizeof(struct test_sort_a),
            CLII((test_generator_callback[]) {
                test_sort_generator_a_1,
                test_sort_generator_a_2,
                test_sort_generator_a_3
            }),
            CLII((test_callback[]) {
                test_sort_a,
            })
        },
        {
            test_sort_disposer_b,
            sizeof(struct test_sort_b),
            CLII((test_generator_callback[]) {
                test_sort_generator_b_1
            }),
            CLII((test_callback[]) {
                test_sort_b_1,
                test_sort_b_2
            })
        },
        {
            test_sort_disposer_c,
            sizeof(struct test_sort_c),
            CLII((test_generator_callback[]) {
                test_sort_generator_c_1
            }),
            CLII((test_callback[]) {
                test_sort_c_1,
                test_sort_c_2
            })
        },
        {
            NULL,
            sizeof(struct test_utf8),
            CLII((test_generator_callback[]) {
                test_utf8_generator,
            }),
            CLII((test_callback[]) {
                test_utf8_len,
                test_utf8_encode,
                test_utf8_decode,
                test_utf16_encode,
                test_utf16_decode
            })
        }
    };
    log_message_generic(log, CODE_METRIC, MESSAGE_NOTE, "Test mode triggered!\n");
    return test(group_arr, countof(group_arr), log);
}

#else

static bool test_main(struct log *log)
{
    log_message_generic(log, CODE_METRIC, MESSAGE_NOTE, "Test mode deactivated!\n");
    return 1;
}

#endif

// Main routine arguments management
struct main_args main_args_default()
{
    struct main_args res = { .thread_cnt = get_processor_count() };
    uint8_bit_set(res.bits, MAIN_ARGS_BIT_POS_THREAD_CNT);
    return res;
}

struct main_args main_args_override(struct main_args args_hi, struct main_args args_lo)
{
    struct main_args res = { .log_path = args_hi.log_path ? args_hi.log_path : args_lo.log_path };
    if (res.log_path && !strcmp(res.log_path, "stderr")) res.log_path = NULL;
    memcpy(res.bits, args_hi.bits, UINT8_CNT(MAIN_ARGS_BIT_CNT));
    if (uint8_bit_test(args_hi.bits, MAIN_ARGS_BIT_POS_THREAD_CNT)) res.thread_cnt = args_hi.thread_cnt;
    else
    {
        res.thread_cnt = args_lo.thread_cnt;
        uint8_bit_set(res.bits, MAIN_ARGS_BIT_POS_THREAD_CNT);
    }
    return res;
}

// Main routine (see below for actual entry point)
static int Main(int argc, char **argv)
{
    /*const char *strings[] = {
        __FUNCTION__,
        "NOTE (%s): No input data specified.\n",
        "INFO (%s): Help mode triggered.\n",
        "WARNING (%s): Unable to compile XML document \"%s\"!\n",
        "WARNING (%s): Unable to execute one or more branches of the XML document \"%s\"!\n",
        "Overall program execution"
    };

    enum
    {
        STR_FN = 0,
        STR_FR_ND,
        STR_FR_IH,
        STR_FR_WC,
        STR_FR_WE,
        STR_M_EXE
    };
    */

    // All names should be sorted according to 'strncmp'!!!
    struct argv_par_sch argv_par_sch =
    {
        CLII((struct tag[]) { { STRI("help"), 0 }, { STRI("log"), 1 }, { STRI("test"), 2 }, { STRI("threads"), 3 }}),
        CLII((struct tag[]) { { STRI("C"), 4 }, { STRI("L"), 5 }, { STRI("T"), 2 }, { STRI("h"), 0 }, { STRI("l"), 1 }, { STRI("t"), 3 } }),
        CLII((struct par_sch[])
        {
            { 0, &(struct handler_context) { offsetof(struct main_args, bits), MAIN_ARGS_BIT_POS_HELP }, empty_handler, 1 },
            { offsetof(struct main_args, log_path), NULL, p_str_handler, 0 },
            { 0, &(struct handler_context) { offsetof(struct main_args, bits), MAIN_ARGS_BIT_POS_TEST }, empty_handler, 1 },
            { offsetof(struct main_args, thread_cnt), &(struct handler_context) { offsetof(struct main_args, bits), MAIN_ARGS_BIT_POS_THREAD_CNT }, size_handler, 0 },
            { 0, &(struct handler_context) { offsetof(struct main_args, bits), MAIN_ARGS_BIT_POS_CAT }, empty_handler, 1 },
            { 0, &(struct handler_context) { offsetof(struct main_args, bits), MAIN_ARGS_BIT_POS_LDE }, empty_handler, 1 },
        })
    };

    /*

    // All names on every sub-level should be sorted according to 'strncmp'!!!
    struct { xmlNode *node; size_t cnt; } dscsch = CLII((xmlNode[])
    {
        {
            .name = STRI("Density"),
            .sz = sizeof(densityContext),
            .prologue = (prologueCallback) densityPrologue,
            .epilogue = (epilogueCallback) densityEpilogue,
            .dispose = (disposeCallback) densityContextDispose,
            CLII((struct att[])
            {
                { STRI("pos"), 0, &(handlerContext) { offsetof(densityContext, bits), DENSITYCONTEXT_BIT_POS_POS }, (readHandlerCallback) boolHandler },
                { STRI("radius"), offsetof(densityContext, radius), &(handlerContext) { offsetof(densityContext, bits) - offsetof(densityContext, radius), DENSITYCONTEXT_BIT_POS_RADIUS }, (readHandlerCallback) uint32Handler },
                { STRI("type"), offsetof(densityContext, type), NULL, (readHandlerCallback) densityTypeHandler }
            }),
            CLII((xmlNode[])
            {
                {
                    .name = STRI("Fold"),
                    .sz = sizeof(densityFoldContext),
                    .prologue = (prologueCallback) densityFoldPrologue,
                    .epilogue = (epilogueCallback) densityFoldEpilogue,
                    .dispose = (disposeCallback) densityFoldContextDispose,
                    CLII((struct att[])
                    {
                        { STRI("group"), offsetof(densityFoldContext, group), NULL, (readHandlerCallback) strHandler },
                        { STRI("type"), offsetof(densityFoldContext, type), NULL, (readHandlerCallback) densityFoldHandler }
                    }),
                    CLII((xmlNode[])
                    {
                        {
                            .name = STRI("Report"),
                            .sz = sizeof(dfReportContext),
                            .prologue = (prologueCallback) dfReportPrologue,
                            .epilogue = (epilogueCallback) dfReportEpilogue,
                            .dispose = (disposeCallback) dfReportContextDispose,
                            CLII((struct att[])
                            {
                                { STRI("header"), 0, &(handlerContext) { offsetof(dfReportContext, bits), DFREPORTCONTEXT_BIT_POS_HEADER }, (readHandlerCallback) boolHandler },
                                { STRI("limit"), offsetof(dfReportContext, limit), &(handlerContext) { offsetof(dfReportContext, bits) - offsetof(dfReportContext, limit), DFREPORTCONTEXT_BIT_POS_LIMIT }, (readHandlerCallback) uint32Handler },
                                { STRI("path"), offsetof(dfReportContext, path), NULL, (readHandlerCallback) strHandler },
                                { STRI("semicolon"), 0, &(handlerContext) { offsetof(dfReportContext, bits), DFREPORTCONTEXT_BIT_POS_SEMICOLON }, (readHandlerCallback) boolHandler },
                                { STRI("threshold"), offsetof(dfReportContext, threshold), &(handlerContext) { offsetof(dfReportContext, bits) - offsetof(dfReportContext, threshold), DFREPORTCONTEXT_BIT_POS_THRESHOLD }, (readHandlerCallback) float64Handler },
                                { STRI("type"), offsetof(dfReportContext, type), NULL, (readHandlerCallback) dfReportHandler }
                            })
                        }
                    })
                },
                {
                    .name = STRI("Regions"),
                    .sz = sizeof(regionsContext),
                    .prologue = (prologueCallback) regionsPrologue,
                    .epilogue = (epilogueCallback) regionsEpilogue,
                    .dispose = (disposeCallback) regionsContextDispose,
                    CLII((struct att[])
                    {
                        { STRI("decay"), offsetof(regionsContext, decay), &(handlerContext) { offsetof(regionsContext, bits) - offsetof(regionsContext, decay), REGIONSCONTEXT_BIT_POS_DECAY }, (readHandlerCallback) uint32Handler },
                        { STRI("depth"), offsetof(regionsContext, depth), &(handlerContext) { offsetof(regionsContext, bits) - offsetof(regionsContext, depth), REGIONSCONTEXT_BIT_POS_DEPTH }, (readHandlerCallback) uint32Handler },
                        { STRI("length"), offsetof(regionsContext, length), &(handlerContext) { offsetof(regionsContext, bits) - offsetof(regionsContext, length), REGIONSCONTEXT_BIT_POS_LENGTH }, (readHandlerCallback) uint32Handler },
                        { STRI("threshold"), offsetof(regionsContext, threshold), &(handlerContext) { offsetof(regionsContext, bits) - offsetof(regionsContext, threshold), REGIONSCONTEXT_BIT_POS_THRESHOLD }, (readHandlerCallback) float64Handler },
                        { STRI("tolerance"), offsetof(regionsContext, tolerance), &(handlerContext) { offsetof(regionsContext, bits) - offsetof(regionsContext, tolerance), REGIONSCONTEXT_BIT_POS_TOLERANCE }, (readHandlerCallback) float64Handler },
                        { STRI("slope"), offsetof(regionsContext, slope), &(handlerContext) { offsetof(regionsContext, bits) - offsetof(regionsContext, slope), REGIONSCONTEXT_BIT_POS_SLOPE }, (readHandlerCallback) float64Handler }
                    })
                },
                {
                    .name = STRI("Report"),
                    .sz = sizeof(dReportContext),
                    .prologue = (prologueCallback) dReportPrologue,
                    .epilogue = (epilogueCallback) dReportEpilogue,
                    .dispose = (disposeCallback) dReportContextDispose,
                    CLII((struct att[])
                    {
                        { STRI("header"), 0, &(handlerContext) { offsetof(dReportContext, bits), DREPORTCONTEXT_BIT_POS_HEADER }, (readHandlerCallback) boolHandler },
                        { STRI("limit"), offsetof(dReportContext, limit), &(handlerContext) { offsetof(dReportContext, bits) - offsetof(dReportContext, limit), DREPORTCONTEXT_BIT_POS_LIMIT }, (readHandlerCallback) uint32Handler },
                        { STRI("path"), offsetof(dReportContext, path), NULL, (readHandlerCallback) strHandler },
                        { STRI("semicolon"), 0, &(handlerContext) { offsetof(dReportContext, bits), DREPORTCONTEXT_BIT_POS_SEMICOLON }, (readHandlerCallback) boolHandler },
                        { STRI("threshold"), offsetof(dReportContext, threshold), &(handlerContext) { offsetof(dReportContext, bits) - offsetof(dReportContext, threshold), DREPORTCONTEXT_BIT_POS_THRESHOLD }, (readHandlerCallback) float64Handler }
                    })
                }
            })
        },
        {
            .name = STRI("MySQL.Dispatch"),
            .sz = sizeof(mysqlDispatchContext),
            .prologue = (prologueCallback) mysqlDispatchPrologue,
            .epilogue = (epilogueCallback) mysqlDispatchEpilogue,
            .dispose = (disposeCallback) mysqlDispatchContextDispose,
            CLII((struct att[])
            {
                { STRI("host"), offsetof(mysqlDispatchContext, host), NULL, (readHandlerCallback) strHandler },
                { STRI("login"), offsetof(mysqlDispatchContext, login), NULL, (readHandlerCallback) strHandler },
                { STRI("password"), offsetof(mysqlDispatchContext, password), NULL, (readHandlerCallback) strHandler },
                { STRI("port"), offsetof(mysqlDispatchContext, port), NULL, (readHandlerCallback) uint32Handler },
                { STRI("schema"), offsetof(mysqlDispatchContext, schema), NULL, (readHandlerCallback) strHandler },
                { STRI("temp.prefix"), offsetof(mysqlDispatchContext, temppr), NULL, (readHandlerCallback) strHandler },
            })
        }
    });

    // All names on every sub-level should be sorted according to 'strncmp'!!!
    xmlNode xmlsch =
    {
        .name = STRI("RegionsMT"),
        .sz = sizeof(frameworkContext),
        .prologue = (prologueCallback) frameworkPrologue,
        .epilogue = (epilogueCallback) frameworkEpilogue,
        .dispose = (disposeCallback) frameworkContextDispose,
        CLII((struct att[])
        {
            { STRI("log"), offsetof(frameworkContext, logFile), NULL, (readHandlerCallback) strHandler },
            { STRI("threads"), offsetof(frameworkContext, threadCount), &(handlerContext) { offsetof(frameworkContext, bits) - offsetof(frameworkContext, threadCount), FRAMEWORKCONTEXT_BIT_POS_THREADCOUNT }, (readHandlerCallback) uint32Handler }
        }),
        CLII((xmlNode[])
        {
            {
                .name = STRI("Data.Load"),
                .sz = sizeof(loadDataContext),
                .prologue = (prologueCallback) loadDataPrologue,
                .epilogue = (epilogueCallback) loadDataEpilogue,
                .dispose = (disposeCallback) loadDataContextDispose,
                CLII((struct att[])
                {
                    { STRI("header"), 0, &(handlerContext) { offsetof(loadDataContext, bits), LOADDATACONTEXT_BIT_POS_HEADER }, (readHandlerCallback) boolHandler },
                    { STRI("logarithm"), 0, &(handlerContext) { offsetof(loadDataContext, bits), LOADDATACONTEXT_BIT_POS_LOGARITHM }, (readHandlerCallback) boolHandler },
                    { STRI("no.ranks"), 0, &(handlerContext) { offsetof(loadDataContext, bits), LOADDATACONTEXT_BIT_POS_NORANKS }, (readHandlerCallback) boolHandler },
                    { STRI("path.chr"), offsetof(loadDataContext, pathchr), NULL, (readHandlerCallback) strHandler },
                    { STRI("path.row"), offsetof(loadDataContext, pathrow), NULL, (readHandlerCallback) strHandler },
                    { STRI("path.test"), offsetof(loadDataContext, pathtest), NULL, (readHandlerCallback) strHandler },
                    { STRI("path.val"), offsetof(loadDataContext, pathval), NULL, (readHandlerCallback) strHandler },
                }),
                { dscsch.node, dscsch.cnt }
            },
            {
                .name = STRI("Genotypes"),
                .sz = sizeof(genotypesContext),
                .prologue = (prologueCallback) genotypesPrologue,
                .epilogue = (epilogueCallback) genotypesEpilogue,
                .dispose = (disposeCallback) genotypesContextDispose,
                CLII((struct att[])
                {
                    { STRI("option0"), 0, &(handlerContext) { offsetof(genotypesContext, bits), GENOTYPESCONTEXT_BIT_POS_OPTION0 }, (readHandlerCallback) boolHandler },
                    { STRI("option1"), 0, &(handlerContext) { offsetof(genotypesContext, bits), GENOTYPESCONTEXT_BIT_POS_OPTION1 }, (readHandlerCallback) boolHandler },
                    { STRI("option2"), 0, &(handlerContext) { offsetof(genotypesContext, bits), GENOTYPESCONTEXT_BIT_POS_OPTION2 }, (readHandlerCallback) boolHandler },
                    { STRI("option3"), 0, &(handlerContext) { offsetof(genotypesContext, bits), GENOTYPESCONTEXT_BIT_POS_OPTION3 }, (readHandlerCallback) boolHandler },
                    { STRI("path.bed"), offsetof(genotypesContext, path_bed), NULL, (readHandlerCallback) strHandler },
                    { STRI("path.bim"), offsetof(genotypesContext, path_bim), NULL, (readHandlerCallback) strHandler },
                    { STRI("path.fam"), offsetof(genotypesContext, path_fam), NULL, (readHandlerCallback) strHandler }
                }),
                { dscsch.node, dscsch.cnt }
            },
            {
                .name = STRI("MySQL.Fetch"),
                .sz = sizeof(mysqlFetchContext),
                .prologue = (prologueCallback) mysqlFetchPrologue,
                .epilogue = (epilogueCallback) mysqlFetchEpilogue,
                .dispose = (disposeCallback) mysqlFetchContextDispose,
                CLII((struct att[])
                {
                    { STRI("host"), offsetof(mysqlFetchContext, host), NULL, (readHandlerCallback) strHandler },
                    { STRI("login"), offsetof(mysqlFetchContext, login), NULL, (readHandlerCallback) strHandler },
                    { STRI("no.ranks"), 0, &(handlerContext) { offsetof(mysqlFetchContext, bits), MYSQLFETCHCONTEXT_BIT_POS_NORANKS }, (readHandlerCallback) boolHandler },
                    { STRI("password"), offsetof(mysqlFetchContext, password), NULL, (readHandlerCallback) strHandler },
                    { STRI("port"), offsetof(mysqlFetchContext, port), NULL, (readHandlerCallback) uint32Handler },
                    { STRI("schema"), offsetof(mysqlFetchContext, schema), NULL, (readHandlerCallback) strHandler }
                }),
                { dscsch.node, dscsch.cnt }
            }
        })
    };

    */

    //FILE *f = stderr; // fopen("f:\\test\\h", "w");
    //fputs("\xef\xbb\xbf", f);
    //for (size_t i = 0; i < (size_t) argc; i++) fprintf(stderr, "%s\n", argv[i]);
    //fclose(f);

    struct style style = {
        .ttl = { INIT_ENV_COL(FG_GREEN), INIT_ENV_COL(FG_RED), INIT_ENV_COL(FG_YELLOW), INIT_ENV_COL(FG_MAGENTA), INIT_ENV_COL(FG_CYAN) },
        .ts = INIT_ENV_COL(FG_BR_BLACK),
        .src = INIT_ENV_COL(FG_BR_BLACK),
        .str = INIT_ENV_COL(FG_BR_MAGENTA),
        .num = INIT_ENV_COL(FG_BR_CYAN),
        .time = INIT_ENV_COL(FG_BR_YELLOW)
    };

    struct log log;
    if (log_init(&log, NULL, 1, 0, style, NULL))
    {
        size_t pos_cnt;
        char **pos_arr;
        struct main_args main_args = { 0 };
        if (argv_parse(argv_par_selector, &argv_par_sch, &main_args, argv, argc, &pos_arr, &pos_cnt, &log))
        {
            main_args = main_args_override(main_args, main_args_default());
            if (uint8_bit_test(main_args.bits, MAIN_ARGS_BIT_POS_HELP))
            {
                // Help mode
                log_message_generic(&log, CODE_METRIC, MESSAGE_INFO, "Help mode triggered!\n");
            }
            else if (uint8_bit_test(main_args.bits, MAIN_ARGS_BIT_POS_TEST))
            {
                test_main(&log);
            }
            else if (uint8_bit_test(main_args.bits, MAIN_ARGS_BIT_POS_CAT))
            {
                if (pos_cnt >= 6)
                {
                    size_t rpl = (size_t) strtoull(pos_arr[4], NULL, 10);
                    uint64_t seed = (uint64_t) strtoull(pos_arr[5], NULL, 10);
                    categorical_run(pos_arr[0], pos_arr[1], pos_arr[2], pos_arr[3], rpl, seed, &log);
                }
            }
            else if (uint8_bit_test(main_args.bits, MAIN_ARGS_BIT_POS_LDE))
            {
                if (pos_cnt >= 2) lde_run(pos_arr[0], pos_arr[1], &log);
            }
            else
            {
                if (!pos_cnt) log_message_generic(&log, CODE_METRIC, MESSAGE_NOTE, "No input data specified.\n");
                else
                {

                }
            }
            free(pos_arr);
        }
        log_close(&log);
    }

    /*
    if (bitTest(args.bits, CMDARGS_BIT_POS_HELP))
    {
        logMsg(&loginfo, strings[STR_FR_IH], strings[STR_FN]);

        // Help message should be here
    }
    else
    {
        if (!inputcnt) logMsg(&loginfo, strings[STR_FR_ND], strings[STR_FN]);

        for (size_t i = 0; i < inputcnt; i++) // Processing input files
        {
            programObject *obj = programObjectFromXML(&xmlsch, input[i], &loginfo);

            if (obj)
            {
                if (!programObjectExecute(obj, &args.in)) logMsg(&loginfo, strings[STR_FR_WE], strings[STR_FN], input[i]);
                programObjectDispose(obj);
            }
            else logMsg(&loginfo, strings[STR_FR_WC], strings[STR_FN], input[i]);
        }
    }
    */

    return EXIT_SUCCESS;
}

#if _WIN32
#   include <windows.h>

int wmain(int argc, wchar_t **wargv)
{
    // Memory leaks will be reported at the program exit
    _CrtSetDbgFlag(_CRTDBG_ALLOC_MEM_DF | _CRTDBG_LEAK_CHECK_DF);
    _CrtSetReportMode(_CRT_WARN, _CRTDBG_MODE_FILE);
    _CrtSetReportFile(_CRT_WARN, _CRTDBG_FILE_STDERR);
    _CrtSetReportMode(_CRT_ERROR, _CRTDBG_MODE_FILE);
    _CrtSetReportFile(_CRT_ERROR, _CRTDBG_FILE_STDERR);
    _CrtSetReportMode(_CRT_ASSERT, _CRTDBG_MODE_FILE);
    _CrtSetReportFile(_CRT_ASSERT, _CRTDBG_FILE_STDERR);

    // Making console output UTF-8 friendly
    if (!SetConsoleOutputCP(CP_UTF8)) return EXIT_FAILURE;

    int main_res = EXIT_FAILURE;

    // Translating UTF-16 command-line parameters to UTF-8
    char **argv;
    if (array_init(&argv, NULL, argc, sizeof(*argv), 0, ARRAY_STRICT))
    {
        size_t base_cnt = 0, i;
        for (i = 0; i < (size_t) argc; i++) // Determining total length and performing error checking
        {
            wchar_t *word = wargv[i];
            uint32_t val;
            uint8_t context = 0;
            for (; *word; word++)
            {
                if (!utf16_decode((uint16_t) *word, &val, NULL, NULL, &context)) break;
                if (!context) base_cnt += utf8_len(val);
            }
            if (*word) break;
            else base_cnt++;
        }
        if (i == (size_t) argc)
        {
            char *base;
            if (array_init(&base, NULL, base_cnt, sizeof(*base), 0, ARRAY_STRICT))
            {
                char *byte = base;
                for (i = 0; i < (size_t) argc; i++) // Performing translation
                {
                    argv[i] = byte;
                    wchar_t *word = wargv[i];
                    uint32_t val;
                    uint8_t context = 0;
                    for (; *word; word++)
                    {
                        if (!utf16_decode((uint16_t) *word, &val, NULL, NULL, &context)) break;
                        if (!context)
                        {
                            uint8_t len;
                            utf8_encode(val, (uint8_t *) byte, &len);
                            byte += len;
                        }
                    }
                    if (*word) break;
                    else *(byte++) = '\0';
                }
                if (i == (size_t) argc) main_res = Main(argc, argv);
                free(base);
            }
        }
        free(argv);
    }
    return main_res;
}

#else

int main(int argc, char **argv)
{
    return Main(argc, argv);
}

#endif


text:
/* vector/gsl_vector_complex_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_VECTOR_COMPLEX_DOUBLE_H__
#define __GSL_VECTOR_COMPLEX_DOUBLE_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_block_complex_double.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  double *data;
  gsl_block_complex *block;
  int owner;
} gsl_vector_complex;

typedef struct
{
  gsl_vector_complex vector;
} _gsl_vector_complex_view;

typedef _gsl_vector_complex_view gsl_vector_complex_view;

typedef struct
{
  gsl_vector_complex vector;
} _gsl_vector_complex_const_view;

typedef const _gsl_vector_complex_const_view gsl_vector_complex_const_view;

/* Allocation */

gsl_vector_complex *gsl_vector_complex_alloc (const size_t n);
gsl_vector_complex *gsl_vector_complex_calloc (const size_t n);

gsl_vector_complex *
gsl_vector_complex_alloc_from_block (gsl_block_complex * b,
                                           const size_t offset,
                                           const size_t n,
                                           const size_t stride);

gsl_vector_complex *
gsl_vector_complex_alloc_from_vector (gsl_vector_complex * v,
                                             const size_t offset,
                                             const size_t n,
                                             const size_t stride);

void gsl_vector_complex_free (gsl_vector_complex * v);

/* Views */

_gsl_vector_complex_view
gsl_vector_complex_view_array (double *base,
                                     size_t n);

_gsl_vector_complex_view
gsl_vector_complex_view_array_with_stride (double *base,
                                                 size_t stride,
                                                 size_t n);

_gsl_vector_complex_const_view
gsl_vector_complex_const_view_array (const double *base,
                                           size_t n);

_gsl_vector_complex_const_view
gsl_vector_complex_const_view_array_with_stride (const double *base,
                                                       size_t stride,
                                                       size_t n);

_gsl_vector_complex_view
gsl_vector_complex_subvector (gsl_vector_complex *base,
                                         size_t i,
                                         size_t n);

_gsl_vector_complex_view
gsl_vector_complex_subvector_with_stride (gsl_vector_complex *v,
                                                size_t i,
                                                size_t stride,
                                                size_t n);

_gsl_vector_complex_const_view
gsl_vector_complex_const_subvector (const gsl_vector_complex *base,
                                               size_t i,
                                               size_t n);

_gsl_vector_complex_const_view
gsl_vector_complex_const_subvector_with_stride (const gsl_vector_complex *v,
                                                      size_t i,
                                                      size_t stride,
                                                      size_t n);

_gsl_vector_view
gsl_vector_complex_real (gsl_vector_complex *v);

_gsl_vector_view
gsl_vector_complex_imag (gsl_vector_complex *v);

_gsl_vector_const_view
gsl_vector_complex_const_real (const gsl_vector_complex *v);

_gsl_vector_const_view
gsl_vector_complex_const_imag (const gsl_vector_complex *v);

/* Operations */

gsl_complex
gsl_vector_complex_get (const gsl_vector_complex * v, const size_t i);

void gsl_vector_complex_set (gsl_vector_complex * v, const size_t i,
                                   gsl_complex z);

gsl_complex
*gsl_vector_complex_ptr (gsl_vector_complex * v, const size_t i);

const gsl_complex
*gsl_vector_complex_const_ptr (const gsl_vector_complex * v, const size_t i);

void gsl_vector_complex_set_zero (gsl_vector_complex * v);
void gsl_vector_complex_set_all (gsl_vector_complex * v,
                                       gsl_complex z);
int gsl_vector_complex_set_basis (gsl_vector_complex * v, size_t i);

int gsl_vector_complex_fread (FILE * stream,
                                    gsl_vector_complex * v);
int gsl_vector_complex_fwrite (FILE * stream,
                                     const gsl_vector_complex * v);
int gsl_vector_complex_fscanf (FILE * stream,
                                     gsl_vector_complex * v);
int gsl_vector_complex_fprintf (FILE * stream,
                                      const gsl_vector_complex * v,
                                      const char *format);

int gsl_vector_complex_memcpy (gsl_vector_complex * dest, const gsl_vector_complex * src);

int gsl_vector_complex_reverse (gsl_vector_complex * v);

int gsl_vector_complex_swap (gsl_vector_complex * v, gsl_vector_complex * w);
int gsl_vector_complex_swap_elements (gsl_vector_complex * v, const size_t i, const size_t j);

int gsl_vector_complex_isnull (const gsl_vector_complex * v);
int gsl_vector_complex_ispos (const gsl_vector_complex * v);
int gsl_vector_complex_isneg (const gsl_vector_complex * v);

#ifdef HAVE_INLINE

extern inline
gsl_complex
gsl_vector_complex_get (const gsl_vector_complex * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      gsl_complex zero = {{0, 0}};
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, zero);
    }
#endif
  return *GSL_COMPLEX_AT (v, i);
}

extern inline
void
gsl_vector_complex_set (gsl_vector_complex * v,
                              const size_t i, gsl_complex z)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  *GSL_COMPLEX_AT (v, i) = z;
}

extern inline
gsl_complex *
gsl_vector_complex_ptr (gsl_vector_complex * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_AT (v, i);
}

extern inline
const gsl_complex *
gsl_vector_complex_const_ptr (const gsl_vector_complex * v,
                                    const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_AT (v, i);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_COMPLEX_DOUBLE_H__ */

ext: h

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef __BLAS_H__
#define __BLAS_H__

#define sgemm sgemm_
#define dgemm dgemm_
#define ssteqr ssteqr_
#define dsteqr dsteqr_
#define dgemv dgemv_
#define ddot ddot_
#define daxpy daxpy_
#define dscal dscal_
#define dasum dasum_

#ifdef __cplusplus
extern "C"{
#endif
    //#include <cblas.h>
    //#include <clapack.h>

    void sgemm_(const char* TRANSA, const char* TRANSB,
        const int* M, const int* N, const int* K,
        const float* ALPHA, const float* A, const int* LDA,
        const float* B, const int* LDB,
        const float* BETA, float* C, const int* LDC);

    void dgemm_(const char* TRANSA, const char* TRANSB,
        const int* M, const int* N, const int* K,
        const double* ALPHA, const double* A, const int* LDA,
        const double* B, const int* LDB,
        const double* BETA, double* C, const int* LDC);
/*
    void ssteqr_(char *compz, const int *n, float *d, float *e,
        float *z, const int *ldz, float *work, const int *info);

    void dsteqr_(char *compz, const int *n, double *d, double *e,
        double *z, const int *ldz, double *work, const int *info);
 */

     void dgemv_(const char *trans, const int *m, const int *n,
          const double *alpha, const double *a, const int *lda,
          const double *x, const int *incx, const double *beta,
          double *y, const int *incy);

     double ddot_(const int* n, const double* x, const int* incx,
          const double* y, const int* incy);

     void daxpy_(const int* n, const double *alpha, const double* x,
	const int* incx, double* y, const int* incy);

	 void dscal_(const int*n, const double *alpha, const double* x, const int* incx);

	 double dasum_(const int*n, const double *x, const int* incx);

#ifdef __cplusplus
}
#endif

#endif

max_stars_repo_name: maumueller/rehashing

max_issues_repo_name: maumueller/rehashing

max_forks_repo_name: maumueller/rehashing


text:
/****************************************************************************/
/*                                                                          */
/* File:      typedefs.h                                                    */
/*                                                                          */
/* Purpose:   header file for all necessary inculdes and typedefs for files */
/*            associated with the NeuroMorph or CellCount projects          */
/*                                                                          */
/* Author:    Marcel Oberlaender                                            */
/*            Max-Planck-Institute for Neurobiologie                        */
/*            Am Kolpferspitz 18                                            */
/*            D-82152 Martinsried (Munich)                                  */
/*                                                                          */
/* Co-Author: Robert Egger                                                  */
/*            Max-Planck-Institute for Medical Research                     */
/*            Jahnstrasse 19                                                */
/*            D-69120 Heidelberg                                            */
/*                                                                          */
/* EMail:     regger@mpimf-heidelberg.mpg.de                                */
/*                                                                          */
/* History:   17.01.2008                                                    */
/*                                                                          */
/* Remarks:   All rights are reserved by the Max-Planck-Society             */
/*                                                                          */
/* Revised by: 	Daniel Udvary												*/
/* 				In Silico Brain Sciences									*/
/*				Max Planck Institute for Neurobiology of Behavior – caesar	*/
/*				Ludwig-Erhard-Allee 2										*/
/*				53175 Bonn, Germany											*/
/*                                                                          */
/* e-mail: 		daniel.udvary@mpinb.mpg.de									*/
/*                                                                          */
/* Date: 		21.11.2019                                                  */
/*                                                                          */
/****************************************************************************/

#ifndef TYPEDEF
#define TYPEDEF

//#define DEBUG

#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#ifdef _OPENMP
#include <omp.h>
#endif

#define PI 3.1415926

#define _CRT_SECURE_NO_DEPRECATE
#define _SCL_SECURE_NO_DEPRECATE
#define _CRT_SECURE_NO_WARNINGS

//STL includes
#include "string"
#include <iostream>
#include <cstdio>
#include <fstream>
#include <string>
#include <list>
#include <vector>
#include <map>
#include <cmath>
#include <algorithm>
#include <complex>
#include <utility>
#include <ctime>

//GSL includes
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

//ITK includes
#include "itkImage.h"
#include "itkImageFileWriter.h"
#include "itkImageFileReader.h"
#include "itkVTKImageImport.h"
#include "itkVTKImageExport.h"
#include "itkNeighborhoodIterator.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkShapedNeighborhoodIterator.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSignedDanielssonDistanceMapImageFilter.h"
#include "itkDanielssonDistanceMapImageFilter.h"
#include "itkShrinkImageFilter.h"
// // for Voronoi diagram creation & use
// #include "itkPointSet.h"
// #include "itkVoronoiDiagram2D.h"
// #include "itkVoronoiDiagram2DGenerator.h"
// #include "itkMeshSpatialObject.h"
// #include "itkSpatialObjectToImageFilter.h"

//VTK includes
//basics
#include "vtkSmartPointer.h"
#include "vtkMath.h"
#include "vtkPoints.h"
#include "vtkImageData.h"
#include "vtkPolyData.h"
#include "vtkUnstructuredGrid.h"
#include "vtkPointData.h"
#include "vtkDataArray.h"
#include "vtkFloatArray.h"
#include "vtkDoubleArray.h"
#include "vtkIdTypeArray.h"
#include "vtkTIFFWriter.h"
#include "vtkPolyDataWriter.h"
#include "vtkGenericCell.h"
#include "vtkPolygon.h"
#include "vtkTriangle.h"
#include "vtkTetra.h"
#include "vtkPlane.h"
#include "vtkLine.h"
#include "vtkCylinder.h"
#include "vtkCutter.h"
#include "vtkClipPolyData.h"
#include "vtkTransform.h"
#include "vtkParametricSpline.h"
#include "vtkKochanekSpline.h"
#include "vtkBoundingBox.h"
#include "vtkImageImport.h"
#include "vtkImageExport.h"
#include "vtkBox.h"
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>

//algorithms
#include "vtkContourFilter.h"
#include "vtkMarchingCubes.h"
#include "vtkSmoothPolyDataFilter.h"
#include "vtkWindowedSincPolyDataFilter.h"
#include "vtkSurfaceReconstructionFilter.h"
#include "vtkCellLocator.h"
#include "vtkSelectPolyData.h"
#include "vtkExtractPolyDataGeometry.h"
#include "vtkTransformPolyDataFilter.h"
#include "vtkAppendFilter.h"
#include "vtkAppendPolyData.h"
#include "vtkPolyDataNormals.h"
#include "vtkHull.h"
#include "vtkPointsProjectedHull.h"
#include "vtkDelaunay2D.h"
#include "vtkDelaunay3D.h"
#include "vtkDataSetSurfaceFilter.h"
#include "vtkButterflySubdivisionFilter.h"
#include "vtkLoopSubdivisionFilter.h"
#include "vtkSelectEnclosedPoints.h"

  typedef float 									CalcPixelType;
  typedef unsigned char									PixelType;
  typedef itk::Image< PixelType, 3 >							ImageType;
  typedef itk::Image< PixelType, 2 >							Image2DType;
  typedef itk::Image< CalcPixelType, 3 >						CalcImageType;
  typedef itk::Image< CalcPixelType, 2 >						CalcImage2DType;
  typedef itk::NeighborhoodIterator< ImageType >					SegNeighborhoodIteratorType;
  typedef itk::NeighborhoodIterator< CalcImageType >					CalcNeighborhoodIteratorType;
  typedef itk::ShapedNeighborhoodIterator< ImageType >					ShapedNeighborhoodIteratorType;
  typedef itk::ShapedNeighborhoodIterator< CalcImageType >				ShapedCalcNeighborhoodIteratorType;
  typedef itk::ImageRegionIterator< ImageType >						IteratorType2;
  typedef itk::ImageRegionConstIterator< ImageType >					ConstIteratorType;
  typedef itk::ImageRegionIterator< Image2DType >					Iterator2DType;
  typedef itk::ImageRegionIterator< CalcImageType >					CalcIteratorType;
  typedef itk::ImageRegionIterator< CalcImage2DType >					Calc2DIteratorType;
  typedef itk::ImageRegionConstIterator< CalcImageType >				ConstCalcIteratorType;
  typedef itk::ImageRegionIteratorWithIndex< ImageType >				IndexIteratorType;
  typedef itk::RescaleIntensityImageFilter< CalcImage2DType, Image2DType >	CalcImage2DToImage2DFilterType;
  typedef itk::SignedDanielssonDistanceMapImageFilter< ImageType, CalcImageType >	DistanceMapImageFilterType;
  typedef itk::DanielssonDistanceMapImageFilter< ImageType, CalcImageType >		DistanceMapImageFilterType2;
  typedef itk::ImageFileWriter< Image2DType >					Image2DWriterType;
  typedef itk::ImageFileReader< Image2DType > 					Image2DReaderType;
  typedef itk::VTKImageImport< ImageType >						VTK2ITKImageImportType;
  typedef itk::VTKImageExport< CalcImageType >						ITK2VTKCalcImageExportType;
  typedef std::vector< SegNeighborhoodIteratorType::OffsetType >			NeighborhoodOffsetVectorType;
  typedef std::vector< ShapedNeighborhoodIteratorType::OffsetType >			ShapedNeighborhoodOffsetVectorType;
  typedef itk::ShrinkImageFilter<Image2DType, Image2DType> 	ShrinkImageFilterType;

//   typedef itk::PointSet< PixelType, 2 >							itkPointSetType;
//   typedef itkPointSetType::PointType							itkPointType;
//   typedef itk::VoronoiDiagram2D< float >						VoronoiDiagramType;
//   typedef itk::VoronoiDiagram2DGenerator< float >					VoronoiDiagramGeneratorType;
//   typedef itk::MeshSpatialObject< VoronoiDiagramType >					VoronoiDiagramSpatialObjectType;
//   typedef itk::SpatialObjectToImageFilter< VoronoiDiagramSpatialObjectType, Image2DType >	VoronoiDiagramToImageFilterType;

  typedef vtkSmartPointer< vtkPoints >		PointsPointerType;
  typedef vtkSmartPointer< vtkDataArray >	DataArrayPointerType;
  typedef vtkSmartPointer< vtkFloatArray >	FloatArrayPointerType;
  typedef vtkSmartPointer< vtkDoubleArray >	DoubleArrayPointerType;
  typedef vtkSmartPointer< vtkIdTypeArray >	IdTypeArrayPointerType;
  typedef vtkSmartPointer< vtkPolygon >		PolygonPointerType;
  typedef vtkSmartPointer< vtkTriangle >	TrianglePointerType;
  typedef vtkSmartPointer< vtkTetra >		TetraPointerType;
  typedef vtkSmartPointer< vtkIdList >		IdListPointerType;
  typedef vtkSmartPointer< vtkPlane >		PlanePointerType;
  typedef vtkSmartPointer< vtkLine >		LinePointerType;
  typedef vtkSmartPointer< vtkCylinder >	CylinderPointerType;
  typedef vtkSmartPointer< vtkCutter >		CutterPointerType;
  typedef vtkSmartPointer< vtkClipPolyData >	ClipPolyDataPointerType;
  typedef vtkSmartPointer< vtkTransform >	TransformPointerType;
  typedef vtkSmartPointer< vtkMatrix4x4 >	HomogeneousMatrixPointerType;
  typedef vtkSmartPointer< vtkParametricSpline >	ParametricSplinePointerType;
  typedef vtkSmartPointer< vtkKochanekSpline >	KochanekSplinePointerType;
  typedef vtkSmartPointer< vtkBoundingBox >	BoundingBoxPointerType;
  typedef vtkSmartPointer< vtkImageData >	ImageDataPointerType;
  typedef vtkSmartPointer< vtkPolyData >	PolyDataPointerType;
  typedef vtkSmartPointer< vtkUnstructuredGrid >	UnstructuredGridPointerType;
  typedef vtkSmartPointer< vtkTIFFWriter >	TiffWriterPointerType;
  typedef vtkSmartPointer< vtkPolyDataWriter >	PolyDataWriterPointerType;
  typedef vtkSmartPointer< vtkGenericCell >	GenericCellPointerType;
  typedef vtkSmartPointer< vtkCell >		CellPointerType;
  typedef vtkSmartPointer< vtkImageImport >	ITK2VTKImageImportPointerType;
  typedef vtkSmartPointer< vtkImageExport >	VTK2ITKImageExportPointerType;

  typedef vtkSmartPointer< vtkMarchingCubes >			MarchingCubesPointerType;
  typedef vtkSmartPointer< vtkContourFilter >			ContourFilterPointerType;
  typedef vtkSmartPointer< vtkSmoothPolyDataFilter >		AveragePolyDataFilterType;
  typedef vtkSmartPointer< vtkWindowedSincPolyDataFilter >	LowpassPolyDataFilterType;
  typedef vtkSmartPointer< vtkSurfaceReconstructionFilter >	SurfaceReconstructionFilterType;
  typedef vtkSmartPointer< vtkCellLocator >			CellLocatorPointerType;
  typedef vtkSmartPointer< vtkSelectPolyData >			SelectPolyDataPointerType;
  typedef vtkSmartPointer< vtkExtractPolyDataGeometry >		ExtractPolyDataGeometryPointerType;
  typedef vtkSmartPointer< vtkTransformPolyDataFilter >		TransformFilterType;
  typedef vtkSmartPointer< vtkAppendFilter >			AppendFilterPointerType;
  typedef vtkSmartPointer< vtkAppendPolyData >			AppendPolyDataPointerType;
  typedef vtkSmartPointer< vtkPolyDataNormals >			PolyDataNormalsPointerType;
  typedef vtkSmartPointer< vtkHull >				ConvexHullFilterPointerType;
  typedef vtkSmartPointer< vtkPointsProjectedHull >		ConvexHull2DFilterPointerType;
  typedef vtkSmartPointer< vtkDelaunay2D >			Delaunay2DFilterPointerType;
  typedef vtkSmartPointer< vtkDelaunay3D >			Delaunay3DFilterPointerType;
  typedef vtkSmartPointer< vtkDataSetSurfaceFilter >		DataSetSurfaceFilterPointerType;
  typedef vtkSmartPointer< vtkButterflySubdivisionFilter >	MeshRefinementFilterPointerType;
  typedef vtkSmartPointer< vtkLoopSubdivisionFilter >		MeshRefinementFilter2PointerType;
  typedef vtkSmartPointer< vtkSelectEnclosedPoints >		SelectEnclosedPointsFilterType;

  typedef std::pair< unsigned int, unsigned int > ColumnCellTypePair;
  typedef std::pair< unsigned int, unsigned int > MatrixIndexType;
  typedef std::vector< unsigned int > SelectionType;

  extern float XYSAMPLING;
  extern float ZSAMPLING;
  extern float averageSomaRadius;
  extern float zScale;
  extern unsigned long BINSIZE;
  extern unsigned long BOXSIZE;

  //VoxelFeatures of the Measurement Vector Type
  #define X_COORD		0
  #define Y_COORD		1
  #define Z_COORD		2

  //celltypes for automatic dendrite detection
  #define SUPRA			1
  #define GRAN			2
  #define INFRA			3

  //Label IDs for Amira Spatial Graphs
  #define Neuron		2
  #define Dendrite		3
  #define ApicalDendrite	4
  #define BasalDendrite		5
  #define Axon			6
  #define Soma			7
  #define Landmark		8
  #define Pia			9
  #define WhiteMatter		48
  #define Vessel		10
  #define Barrel		11
  #define ZAxis			50
  #define aRow			12
  #define A1			13
  #define A2			14
  #define A3			15
  #define A4			16
  #define bRow			17
  #define B1			18
  #define B2			19
  #define B3			20
  #define B4			21
  #define cRow			22
  #define C1			23
  #define C2			24
  #define C3			25
  #define C4			26
  #define C5			27
  #define C6			28
  #define dRow			29
  #define D1			30
  #define D2			31
  #define D3			32
  #define D4			33
  #define D5			34
  #define D6			35
  #define eRow			36
  #define E1			37
  #define E2			38
  #define E3			39
  #define E4			40
  #define E5			41
  #define E6			42
  #define greekRow		43
  #define Alpha			44
  #define Beta			45
  #define Gamma			46
  #define Delta			47
  #define Septum		0

  // cell types for NeuroNet Connectome analysis
  #define L2			1
  #define L34			2
  #define L4py			3
  #define L4sp			4
  #define L4ss			5
  #define L5st			6
  #define L5tt			7
  #define L6cc			8
  #define L6ccinv		9
  #define L6ct			10
  #define VPM			11
  #define L2axon		12
  #define L34axon		13
  #define L4pyaxon		14
  #define L4spaxon		15
  #define L4ssaxon		16
  #define L5staxon		17
  #define L5ttaxon		18
  #define L6ccaxon		19
  #define L6ccinvaxon	20
  #define L6ctaxon		21
  #define SymLocal		22
  #define SymLocal1		23
  #define SymLocal2		24
  #define SymLocal3		25
  #define SymLocal4		26
  #define SymLocal5		27
  #define SymLocal6		28
  #define L1			29
  #define L23Trans		30
  #define L45Sym		31
  #define L45Peak		32
  #define L56Trans		33
  #define SymLocalaxon	34
  #define SymLocal1axon	35
  #define SymLocal2axon	36
  #define SymLocal3axon	37
  #define SymLocal4axon	38
  #define SymLocal5axon	39
  #define SymLocal6axon	40
  #define L1axon		41
  #define L23Transaxon	42
  #define L45Symaxon	43
  #define L45Peakaxon	44
  #define L56Transaxon	45
//168fines
#endif


text:
#include <stdlib.h>
#include <gsl/gsl_vector.h>
#include "../include/type.h"
#include "../include/global.h"
#include "../include/cthreadpool.h"

void *pool_thread(void *arg) {
    threadpool_t *pool = (threadpool_t*)arg;
    task_t task;
    int head;

    while(1) {
        pthread_mutex_lock(&(pool->pool_lock));

        while(pool->task_count == 0 && pool->flag != POOL_RECLAIM) {
            pthread_cond_wait(&(pool->notify), &(pool->pool_lock));
        }

        if (pool->task_count == 0 && pool->flag == POOL_RECLAIM) {
            break;
        }

        head = pool->task_head;
        task.func = pool->tasks[head].func;
        task.arg = pool->tasks[head].arg;
        pool->task_head++;
        pool->task_count--;

        pthread_mutex_unlock(&(pool->pool_lock));

        *(pool->tasks[head].ret) = (task.func)(task.arg);
    }

    pthread_mutex_unlock(&(pool->pool_lock));

    return NULL;
}

threadpool_t* threadpool_create(const int nthread, const int queue_size) {
    int i = 0;
    threadpool_t *pool = NULL;

    if ((pool = (threadpool_t *)malloc(sizeof(threadpool_t))) == NULL) {
        return NULL;
    }

    pool->nthread = nthread;
    pool->queue_size = queue_size;
    pool->task_head = 0;
    pool->task_tail = -1;
    pool->task_count = 0;

    pthread_mutex_init(&pool->pool_lock, NULL);
    pthread_cond_init(&pool->notify, NULL);

    pool->working = (int*)malloc(sizeof(int)*nthread);
    pool->threads = (pthread_t*)malloc(sizeof(pthread_t)*nthread);
    pool->tasks = (task_t*)malloc(sizeof(task_t)*queue_size);
    for (i = 0; i < nthread; ++i) {
        pool->working[i] = IDLE;
    }

    for(i = 0; i < nthread; i++) {
        if(pthread_create(&(pool->threads[i]), NULL, pool_thread, (void*)pool) != 0) {
            threadpool_reclaim(pool);
            return NULL;
        }
        pool->working[i] = INCOMPLETED;
    }

    return pool;
}

int threadpool_add(threadpool_t *pool, void* (*func)(void *), void *arg, void **ret) {
    int tail;

    if (pool == NULL || func == NULL) {
        return POOL_ADD_ERR;
    }

    if(pthread_mutex_lock(&(pool->pool_lock)) != 0) {
        return POOL_LOCK_FAILURE;
    }

    tail = pool->task_tail + 1;
    tail = tail == pool->queue_size? 0:tail;
    if (pool->task_count == pool->queue_size) {
        // TODO erroneous.
        if ((realloc(pool->tasks, 2 * pool->queue_size)) != 0) {
            return POOL_ADD_ERR;
        }
        tail = pool->queue_size;
        pool->queue_size *= 2;
    }

    pool->tasks[tail].func = func;
    pool->tasks[tail].arg = arg;
    pool->tasks[tail].ret = ret;
    pool->task_tail = tail;
    pool->task_count++;

    if(pthread_cond_signal(&(pool->notify)) != 0 ||
       pthread_mutex_unlock(&pool->pool_lock) != 0) {
        return POOL_LOCK_FAILURE;
    }

    return 0;
}

int threadpool_reclaim(threadpool_t *pool) {
    int i = 0;

    if(pthread_mutex_lock(&(pool->pool_lock)) != 0) {
        return POOL_LOCK_FAILURE;
    }

    pool->flag = POOL_RECLAIM;
    pthread_mutex_unlock(&(pool->pool_lock));

    for (i = 0; i < pool->nthread; ++i) {
        if (pool->working[i] != IDLE) {
            pthread_join(pool->threads[i], NULL);
        }
    }

    free(pool->threads);
    free(pool->tasks);
    free(pool);

    return 0;
}

ext: c

max_stars_repo_name: chunjie-sam-liu/rmats-turbo

max_issues_repo_name: chunjie-sam-liu/rmats-turbo

max_forks_repo_name: chunjie-sam-liu/rmats-turbo


text:
/***************************************************************************
Copyright (c) 2014, The OpenBLAS Project
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the
distribution.
3. Neither the name of the OpenBLAS project nor the names of
its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*****************************************************************************/

// #include "bench.h"

#include <cblas.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>

// mlir-clang gemm.c -I .. -I /home/lchelini/scratch/polygeist/llvm-project/llvm/../clang/lib/Headers -o gemm --emit-llvm

int main(int argc, char *argv[]) {

  float *a, *b;
  float *c;
  float alpha = 1.0;
  float beta = 1.0;
  char transa = 'N';
  char transb = 'N';
  int m, n, k, i, j, lda, ldb, ldc;

  m = 500;
  k = 500;
  n = 500;

  a = (float *)malloc(sizeof(float) * m * k);
  b = (float *)malloc(sizeof(float) * k * n);
  c = (float *)malloc(sizeof(float) * m * n);

  for (i = 0; i < m * k; i++) {
    a[i] = ((float)rand() / (float)RAND_MAX) - 0.5;
  }
  for (i = 0; i < k * n; i++) {
    b[i] = ((float)rand() / (float)RAND_MAX) - 0.5;
  }
  for (i = 0; i < m * n; i++) {
    c[i] = ((float)rand() / (float)RAND_MAX) - 0.5;
  }

  if (transa == 'N') {
    lda = m;
  } else {
    lda = k;
  }
  if (transb == 'N') {
    ldb = k;
  } else {
    ldb = n;
  }
  ldc = m;

  sgemm_(&transa, &transb, &m, &n, &k, &alpha, &(a[5]), &lda, b, &ldb, &beta, c, &ldc);
  fprintf(stderr, "%.6f\n", c[100]);

  return 0;
}

max_stars_repo_name: chelini/openBlas

max_issues_repo_name: chelini/openBlas

max_forks_repo_name: chelini/openBlas


text:
//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2008 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/

#ifndef __Pulsar_Interpolation_h
#define __Pulsar_Interpolation_h

#include "Reference.h"

#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>

#include <vector>

//! Interface to GSL interpolation routines
class Interpolation : public Reference::Able
{

public:

  //! Default constructor
  Interpolation ();

  //! Destructor
  virtual ~Interpolation ();

  //! Initialize interpolation object
  void init (const std::vector<double>& x, const std::vector<double>& y);

  //! Evaluate at the given abscissa
  double eval (double x);

protected:

  const double* xa;
  const double* ya;
  size_t size;

  gsl_interp* interp;
  gsl_interp_accel* acc;

  void destroy ();
};

#endif


text:
static char help[] =
"ODE system solver example using TS, but with Jacobian.  Sets TS type to\n"
"implicit Crank-Nicolson.  Compare ode.c.\n\n";

#include <petsc.h>

extern PetscErrorCode ExactSolution(PetscReal, Vec);
extern PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void*);
extern PetscErrorCode FormRHSJacobian(TS, PetscReal, Vec, Mat, Mat, void*);

int main(int argc,char **argv) {
  PetscErrorCode ierr;
  PetscInt   steps;
  PetscReal  t0 = 0.0, tf = 20.0, dt = 0.1, err;
  Vec        y, yexact;
  Mat        J;
  TS         ts;

  ierr = PetscInitialize(&argc,&argv,NULL,help); if (ierr) return ierr;

  ierr = VecCreate(PETSC_COMM_WORLD,&y); CHKERRQ(ierr);
  ierr = VecSetSizes(y,PETSC_DECIDE,2); CHKERRQ(ierr);
  ierr = VecSetFromOptions(y); CHKERRQ(ierr);
  ierr = VecDuplicate(y,&yexact); CHKERRQ(ierr);

  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
  ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,NULL); CHKERRQ(ierr);

//STARTMATJ
  ierr = MatCreate(PETSC_COMM_WORLD,&J); CHKERRQ(ierr);
  ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2); CHKERRQ(ierr);
  ierr = MatSetFromOptions(J); CHKERRQ(ierr);
  ierr = MatSetUp(J); CHKERRQ(ierr);
  ierr = TSSetRHSJacobian(ts,J,J,FormRHSJacobian,NULL); CHKERRQ(ierr);
  ierr = TSSetType(ts,TSCN); CHKERRQ(ierr);
//ENDMATJ

  // set time axis
  ierr = TSSetTime(ts,t0); CHKERRQ(ierr);
  ierr = TSSetMaxTime(ts,tf); CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,dt); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts); CHKERRQ(ierr);

  // set initial values and solve
  ierr = TSGetTime(ts,&t0); CHKERRQ(ierr);
  ierr = ExactSolution(t0,y); CHKERRQ(ierr);
  ierr = TSSolve(ts,y); CHKERRQ(ierr);

  // compute error and report
  ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr);
  ierr = TSGetTime(ts,&tf); CHKERRQ(ierr);
  ierr = ExactSolution(tf,yexact); CHKERRQ(ierr);
  ierr = VecAXPY(y,-1.0,yexact); CHKERRQ(ierr);    // y <- y - yexact
  ierr = VecNorm(y,NORM_INFINITY,&err); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
              "error at tf = %.3f with %d steps:  |y-y_exact|_inf = %g\n",
              tf,steps,err); CHKERRQ(ierr);

  MatDestroy(&J);
  VecDestroy(&y);  VecDestroy(&yexact);  TSDestroy(&ts);
  return PetscFinalize();
}

PetscErrorCode ExactSolution(PetscReal t, Vec y) {
    PetscReal *ay;
    VecGetArray(y,&ay);
    ay[0] = t - PetscSinReal(t);
    ay[1] = 1.0 - PetscCosReal(t);
    VecRestoreArray(y,&ay);
    return 0;
}

PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec y, Vec g,
                               void *ptr) {
    const PetscReal *ay;
    PetscReal       *ag;
    VecGetArrayRead(y,&ay);
    VecGetArray(g,&ag);
    ag[0] = ay[1];            // = g_1(t,y)
    ag[1] = - ay[0] + t;      // = g_2(t,y)
    VecRestoreArrayRead(y,&ay);
    VecRestoreArray(g,&ag);
    return 0;
}

//STARTJACOBIAN
PetscErrorCode FormRHSJacobian(TS ts, PetscReal t, Vec y, Mat J, Mat P,
                               void *ptr) {
    PetscErrorCode ierr;
    PetscInt   row[2] = {0, 1},  col[2] = {0, 1};
    PetscReal  v[4] = { 0.0, 1.0,
                       -1.0, 0.0};
    ierr = MatSetValues(P,2,row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}
//ENDJACOBIAN


text:
#ifndef INCLUDED_TELEPATH_MATRIX_H
#define INCLUDED_TELEPATH_MATRIX_H

#include <cblas.h>
#include "telepath/blas/traits.h"

namespace telepath{

    enum ORDER {
        ROW_MAJOR = CblasRowMajor,
        COL_MAJOR = CblasColMajor
    };
    enum TRANSPOSE {
        NONE = CblasNoTrans,
        TRANS = CblasTrans,
        CONJ_TRANS = CblasConjTrans
    };
    enum UPLO{
        UPPER = CblasUpper,
        LOWER = CblasLower,
    };

    /*
    Struct for storing BLAS matrix data
    */
    template< typename T >
    struct matrix{
        ORDER layout;
        TRANSPOSE trans;
        T* array;
        std::size_t rows;
        std::size_t cols;
        std::size_t ld;
    };

    template< typename T >
    struct sym_matrix{
        ORDER layout;
        TRANSPOSE trans;
        UPLO uplo;
        T* array;
        std::size_t rows;
        std::size_t cols;
        std::size_t ld;
    };

} //namespace telepath

namespace blas{

    template< typename T >
    struct Traits< const telepath::matrix<T> >{
        using mat_t = typename telepath::matrix<T>;
        const mat_t& mat;
        constexpr Traits( const mat_t& m ) : mat( m ) { }
        constexpr telepath::ORDER layout() const{ return mat.layout; }
        constexpr telepath::TRANSPOSE trans() const{ return mat.trans; }
        constexpr const T* array() const{ return mat.array; }
        constexpr auto nrows() const{ return mat.rows; }
        constexpr auto ncols() const{ return mat.cols; }
        constexpr auto ld() const{ return mat.ld; }
    };

    template< typename T >
    struct Traits< telepath::matrix<T> >{
        using mat_t = typename telepath::matrix<T>;
        mat_t& mat;
        constexpr Traits( mat_t& m ) : mat( m ) { }
        constexpr telepath::ORDER layout() const{ return mat.layout; }
        constexpr telepath::TRANSPOSE trans() const{ return mat.trans; }
        inline T* array() { return mat.array; }
        constexpr const T* array() const{ return mat.array; }
        constexpr auto nrows() const{ return mat.rows; }
        constexpr auto ncols() const{ return mat.cols; }
        constexpr auto ld() const{ return mat.ld; }
    };

    template< typename T >
    struct Traits< const telepath::sym_matrix<T> >{
        using mat_t = telepath::sym_matrix<T>;
        const mat_t& mat;
        constexpr Traits( const mat_t& m ) : mat( m ) { }
        constexpr telepath::ORDER layout() const{ return mat.layout; }
        constexpr telepath::TRANSPOSE trans() const{ return mat.trans; }
        constexpr telepath::UPLO uplo() const{ return mat.uplo; }
        constexpr const T* array() const{ return mat.array; }
        constexpr auto nrows() const{ return mat.rows; }
        constexpr auto ncols() const{ return mat.cols; }
        constexpr auto ld() const{ return mat.ld; }
    };

    template< typename T >
    struct Traits< telepath::sym_matrix<T> >{
        using mat_t = typename telepath::sym_matrix<T>;
        mat_t& mat;
        constexpr Traits( mat_t& m ) : mat( m ) { }
        constexpr telepath::ORDER layout() const{ return mat.layout; }
        constexpr telepath::TRANSPOSE trans() const{ return mat.trans; }
        constexpr telepath::UPLO uplo() const{ return mat.uplo; }
        constexpr const T* array() const{ return mat.array; }
        inline T* array() { return mat.array; }
        constexpr auto nrows() const{ return mat.rows; }
        constexpr auto ncols() const{ return mat.cols; }
        constexpr auto ld() const{ return mat.ld; }
    };

    template< typename T >
    struct MatrixTraits< telepath::matrix<T> >{
        using mat = typename telepath::matrix<T>;
        using scalar_t = T;
        constexpr static CBLAS_ORDER layout( const mat& m ){
            return (CBLAS_ORDER) m.layout;
        }
        constexpr static CBLAS_TRANSPOSE trans( const mat& m ){
            return (CBLAS_TRANSPOSE) m.trans;
        }
        constexpr static auto nrows( const mat& m ){ return m.rows; }
        constexpr static auto ncols( const mat& m ){ return m.cols; }
        constexpr static T* array( mat& m ){ return m.array; }
        constexpr static const T* array( const mat& m ){ return m.array; }
        constexpr static auto ld( const mat& m ){ return m.ld; }
    };

    template< typename T >
    struct MatrixTraits< telepath::sym_matrix<T> >{
        using mat = typename telepath::sym_matrix<T>;
        using scalar_t = T;
        constexpr static CBLAS_ORDER layout( const mat& m ){
            return (CBLAS_ORDER) m.layout;
        }
        constexpr static CBLAS_TRANSPOSE trans( const mat& m ){
            return (CBLAS_TRANSPOSE) m.trans;
        }
        constexpr static CBLAS_UPLO uplo( const mat& m ){
            return (CBLAS_UPLO) m.uplo;
        }
        constexpr static auto nrows( const mat& m ){ return m.rows; }
        constexpr static auto ncols( const mat& m ){ return m.cols; }
        constexpr static T* array( mat& m ){ return m.array; }
        constexpr static const T* array( const mat& m ){ return m.array; }
        constexpr static auto ld( const mat& m ){ return m.ld; }
    };

} //namespace blas

#endif

max_stars_repo_name: tbepler/telepath

max_issues_repo_name: tbepler/telepath

max_forks_repo_name: tbepler/telepath


text:
#include <stdio.h>
#include <gsl/gsl_fit.h>

int main (void)
{
  int n = 11;
  double x[11] = {10.0, 8.0, 13.0, 9.0, 11.0,
                  14.0 ,6.0, 4.0 , 12.0,7.0,5.0};
  double y[11] = {8.04, 6.95,7.68, 8.81, 8.33,
                  9.96, 7.24,4.26,10.84, 4.82,5.68 };

  double c0, c1, cov00, cov01, cov11, sumsq;

  gsl_fit_linear(x, 1, y, 1, n,
                 &c0, &c1, &cov00, &cov01, &cov11,
                 &sumsq);

  printf ("best fit: Y = %g + %g X\n", c0, c1);
  printf ("covariance matrix:\n");
  printf ("[ %g, %g\n  %g, %g]\n", cov00, cov01, cov01, cov11);
  printf ("sumsq = %g\n", sumsq);
  printf ("\n");

  // plot
  FILE *pipe = popen("gnuplot -persist", "w"); // Open a pipe to gnuplot
  if (pipe) // If gnuplot is found
  {
     fprintf(pipe, "set term wx\n");         // set the terminal
     fprintf(pipe, "set xlabel 'X'\n");
     fprintf(pipe, "set ylabel 'Y'\n");
     fprintf(pipe, "set xrange [0:20]\n");
     fprintf(pipe, "set yrange [2:14]\n");
     fprintf(pipe, "set title '<X,Y> and Linear fit:y=%.4f*x+%.4f'\n",c1,c0);

     /* In this case, the datafile is written directly to the gnuplot pipe with no need for a temporary file.
           The special filename '-' specifies that the data are inline; i.e., they follow the command.
        1 sending gnuplot the plot '-' command
        2 followed by data points
        3 followed by the letter "e"
     */

     // 1 sending gnuplot the plot '-' command
     fprintf(pipe, "plot '-' title '<x,y>' with points  pt 7 lc rgb 'blue',\
                         '-' title 'Line' with  linespoints  pt  6 lc rgb 'red'\n");

     // 2 followed by data points: <x,y>
     for (int i = 0; i < n; i++)
     {
        fprintf(pipe, "%lf %lf\n", x[i], y[i]);
     }
     // 3 followed by the letter "e"
     fprintf(pipe, "e");

     // linear fit
     fprintf(pipe,"\n"); // start a new draw item
     fprintf(pipe, "%lf %lf\n", 0.0, c0+c1*0,0);
     for (int i = 0; i < n; i++)
     {
        fprintf(pipe, "%lf %lf\n", x[i], c0+c1*x[i]);
     }
     fprintf(pipe, "%lf %lf\n", 20.0,c0+c1*20,0);
     fprintf(pipe, "e");

     fflush(pipe);
     fprintf(pipe, "exit \n"); // exit gnuplot
     pclose(pipe);             //close pipe
  }

  return 0;
}


text:
#include <stdlib.h>
#include <stdio.h>

//#define VERY_VERBOSE

#define MSG(x) { printf(x); fflush(stdout); }
#define MSGF(x) { printf x ; fflush(stdout); }

#include "buffer.h"

static void test_buffer()
{
}

#include "array.h"

static void test_array()
{
	int k;
	const int * p;
	const int * q;
	OLIO_ARRAY_STACK(st, int, 64);
	olio_array * hp;

	OLIO_ARRAY_ALLOC(hp, int, 0);

	printf("[testing array...]\n"); fflush(stdout);

	MSGF(("st %u\n", st->base.allocated));
	MSGF(("hp %u\n", hp->base.allocated));

#define ARRAY_COUNT 1000

	MSG("appending to arrays...\n");
	for (k = 0; k < ARRAY_COUNT; k++) {
		olio_array_append(st, &k, 1);
		olio_array_append(hp, &k, 1);
	}

	MSG("prepending to arrays...\n");
	for (k = 0; k < ARRAY_COUNT; k++) {
		olio_array_prepend(st, &k, 1);
		olio_array_prepend(hp, &k, 1);
	}

	MSG("inserting to arrays...\n");
	for (k = 0; k <= ARRAY_COUNT * 2; k++) {
		olio_array_insert(st, k * 2, &k, 1);
		olio_array_insert(hp, k * 2, &k, 1);
	}

	MSG("checking array contents...\n");
	p = olio_array_contents(st);
	q = olio_array_contents(hp);
	for (k = 0; k < ARRAY_COUNT; k++) {
		if (p[k * 2 + 1] != ARRAY_COUNT - 1 - k ||
			p[k * 2 + 1 + ARRAY_COUNT * 2] != k)
			MSG("ERROR: odd numbered stack item not as expected\n");

		if (p[k * 2] != k ||
			p[k * 2 + ARRAY_COUNT * 2] != k + ARRAY_COUNT)
			MSG("ERROR: even numbered stack item not as expected\n");

		if (q[k * 2 + 1] != ARRAY_COUNT - 1 - k ||
			q[k * 2 + 1 + ARRAY_COUNT * 2] != k)
			MSG("ERROR: odd numbered heap item not as expected\n");

		if (q[k * 2] != k ||
			q[k * 2 + ARRAY_COUNT * 2] != k + ARRAY_COUNT)
			MSG("ERROR: even numbered heap item not as expected\n");
	}
	if (p[ARRAY_COUNT * 4] != ARRAY_COUNT * 2)
		MSG("ERROR: last item in stack array not as expected\n");
	if (q[ARRAY_COUNT * 4] != ARRAY_COUNT * 2)
		MSG("ERROR: last item in heap array not as expected\n");

#ifdef VERY_VERBOSE
	for (k = 0; k < olio_array_length(st) ||
		k < olio_array_length(hp); k++) {
		MSGF(("%i: st % 4i hp % 4i\n", k,
			(k < olio_array_length(st)) ? p[k] : -1,
			(k < olio_array_length(hp)) ? q[k] : -1));
	}
#endif

	MSGF(("st %p = %p\n", st->base.data, st->stack));
	MSGF(("hp %p = %p\n", hp->base.data, hp->stack));

	MSG("freeing arrays...\n");
	olio_array_free(st);
	olio_array_free(hp);
}

#include "string.h"

static void test_string()
{
}

#include "graph.h"

static void test_graph()
{
	int rc;
	olio_graph g;

	printf("[testing graph...]\n"); fflush(stdout);

	printf("graph initialization...\n"); fflush(stdout);
	olio_graph_init(&g, 10);

	printf("graph setting edges...\n"); fflush(stdout);
	olio_graph_set_edge(&g, 0, 1, 10);
	olio_graph_set_edge(&g, 1, 2, 10);

	printf("graph testing cyclic...\n"); fflush(stdout);
	rc = olio_graph_acyclic(&g);
	if (rc == 0) {
		printf("ERROR: graph NOT acyclic\n"); fflush(stdout);
	}

	printf("graph setting edges...\n"); fflush(stdout);
	olio_graph_set_edge(&g, 2, 0, 10);

	MSG("graph testing cyclic...\n");
	rc = olio_graph_acyclic(&g);
	if (rc != 0) MSG("ERROR: graph NOT cyclic\n");

	MSG("freeing graph...\n");
	olio_graph_free(&g);
}

#include "skiplist.h"

static void test_skiplist()
{
	uint32_t test_data[] = {0xf, 0x436, 0x53547, 0x65, 0x54, 0x5654, 0x8756, 0x8767};

	olio_skiplist sk;
	int i;

	printf("sizeof(entry) = %li\n", sizeof(olio_skiplist_entry));
	printf("sizeof(block) = %li\n", sizeof(olio_skiplist_block));
	printf("sizeof(skiplist) = %li\n", sizeof(olio_skiplist));

	olio_skiplist_init(&sk);
	printf("adding...\n");
	for (i = 0; i < sizeof(test_data) / sizeof(uint32_t); i++) {
		int rv = olio_skiplist_add(&sk, test_data[i], (void*) test_data[i]);
		printf("k: 0x%x, rv: %i\n", test_data[i], rv);
	}

	printf("searching...\n");
	for (i = 0; i < sizeof(test_data) / sizeof(uint32_t); i++) {
		void * d = olio_skiplist_find(&sk, test_data[i]);
		printf("k: 0x%x, data: %p\n", test_data[i], d);
	}

	olio_skiplist_display(&sk);

	printf("freeing...\n");
	olio_skiplist_free(&sk);

	olio_skiplist_display(&sk);

	olio_skiplist_init(&sk);
	printf("adding...\n");
	for (i = 0; i < sizeof(test_data) / sizeof(uint32_t); i++) {
		int rv = olio_skiplist_add(&sk, test_data[i], (void*) test_data[i]);
		printf("k: 0x%x, rv: %i\n", test_data[i], rv);
	}

	printf("deleting...\n");
	for (i = 0; i < sizeof(test_data) / sizeof(uint32_t); i++) {
		void * d;
		int rv = olio_skiplist_remove(&sk, test_data[i], &d);
		printf("k: 0x%x, rv: %i, data: %p\n", test_data[i], rv, d);
		olio_skiplist_display(&sk);
	}

	printf("freeing...\n");
	olio_skiplist_free(&sk);

	olio_skiplist_display(&sk);
}

#include "error.h"

static void test_error()
{
}

#ifdef HAVE_GSL
#include <gsl/gsl_rng.h>
#include "random.h"

gsl_rng * gsl;
olio_random olio;

static void test_random()
{
    uint32_t seed;
    uint32_t a, b;
    long i;

    printf("[testing random...]\n"); fflush(stdout);

    gsl = gsl_rng_alloc(gsl_rng_mt19937);

    //for (seed = 0; seed < 1000; seed++)
    {
	seed = 0;
	gsl_rng_set(gsl, 0);
	olio_random_set_seed(&olio, 4357);
	for (i = 0; i < 10000; i++) {
	    a = gsl_rng_uniform_int(gsl, 0xffffffff);
	    b = olio_random_integer(&olio);
	    if (a != b) {
		printf("mismatch\n");
	    }
	}
    }
    gsl_rng_free(gsl);
}
#else
static void test_random() {}
#endif

#include "phash.h"
#include "hash.h"

static void test_hashes()
{
    unsigned long i;
    const char * s[] = {
	"alpha", "beta", "gamma", "delta", "epsilon",
	"lamda", "mu", "nu", "omicron", "pi", "phi", "psi",
	"tau", "theta", "zeta", NULL
    };

    olio_phash_build m;
    int rc;

    printf("[testing hashes...]\n"); fflush(stdout);

    rc = olio_phash_init(&m, NULL);
    if (rc == -1) {
	printf("error, aborting\n");
	exit(1);
    }

    for (i = 0; s[i] != NULL; i++) {
	rc = olio_phash_add_entry(&m, s[i], strlen(s[i]), i * 1000);
	if (rc == -1) {
	    printf("error, aborting\n");
	    exit(1);
	}
    }

    rc = olio_phash_generate(&m, 1000, 1);
    if (rc < 0) {
	printf("error, aborting\n");
	exit(1);
    }
    if (rc > 0) {
	printf("failed to generate with attempt/extra limits\n");
	exit(1);
    }

    printf("phash generated:\n");
    printf("A hash: %08x\nB hash: %08x\nG size: %u\n",
	   m.phash->a_hash_seed, m.phash->b_hash_seed,
	   m.phash->g_table_size);

    printf("{ ");
    for (i = 0; i < m.phash->g_table_size; i++) {
	int16_t * gv = (int16_t *) m.phash->data;
	printf("%i, ", gv[i]);
    }
    printf("}\n");

    for (i = 0; s[i] != NULL; i++) {
	int32_t v = olio_phash_value(m.phash,
				     s[i], strlen(s[i]));
	printf(" %s -> %li\n",
	       s[i], v);
    }

    olio_phash_free(&m);
}

#include "varint.h"

static void test_varint_sub(uint32_t x, uint8_t bytes_expected)
{
    uint32_t a;
    uint8_t storage[5];
    uint8_t bytes_set, bytes_get;
    uint8_t i;

    bytes_set = olio_varint_set(x, storage);
    a = olio_varint_get(storage, &bytes_get);

    if (a != x || bytes_set != bytes_get || bytes_set != bytes_expected)
    {
	printf("x = %lu\n(%u) 0x", x, bytes_set);
	for (i = 0; i < bytes_set; i++)
	    printf("%02x", storage[i]);
	printf("\na = %lu\n(%u)\n\n", a, bytes_get);
    }
}

static void test_varint()
{
    printf("[testing varint...]\n"); fflush(stdout);

    test_varint_sub(0, 1);
    test_varint_sub(1,1);
    test_varint_sub(127, 1);
    test_varint_sub(200, 2);
    test_varint_sub(32000, 3);
    test_varint_sub(32000000, 3);
    test_varint_sub(100000000, 4);
    test_varint_sub(400000000, 5);
    test_varint_sub(0, 1);
    test_varint_sub(1, 1);
}

int main(int argv, char ** argc)
{
	test_buffer();
	test_array();
	test_string();
	test_graph();
	test_skiplist();
	test_error();
	test_random();
	test_hashes();
	test_varint();
}

max_stars_repo_name: jabr/olio

max_issues_repo_name: jabr/olio

max_forks_repo_name: jabr/olio


text:
#pragma once

#include "AsyncSystem.h"
#include "IAssetRequest.h"
#include "Library.h"

#include <gsl/span>

#include <cstddef>
#include <memory>
#include <string>
#include <vector>

namespace CesiumAsync {
class AsyncSystem;

/**
 * @brief Provides asynchronous access to assets, usually files downloaded via
 * HTTP.
 */
class CESIUMASYNC_API IAssetAccessor {
public:
  /**
   * @brief An HTTP header represented as a key/value pair.
   */
  typedef std::pair<std::string, std::string> THeader;

  virtual ~IAssetAccessor() = default;

  /**
   * @brief Starts a new request for the asset with the given URL.
   * The request proceeds asynchronously without blocking the calling thread.
   *
   * @param asyncSystem The async system used to do work in threads.
   * @param url The URL of the asset.
   * @param headers The headers to include in the request.
   * @return The in-progress asset request.
   */
  virtual CesiumAsync::Future<std::shared_ptr<IAssetRequest>> requestAsset(
      const AsyncSystem& asyncSystem,
      const std::string& url,
      const std::vector<THeader>& headers = {}) = 0;

  /**
   * @brief Starts a new POST request to the given URL.
   *
   * The request proceeds asynchronously without blocking the calling thread.
   *
   * @param asyncSystem The async system used to do work in threads.
   * @param url The URL of the asset.
   * @param headers The headers to include in the request.
   * @param contentPayload The payload data of the POST.
   * @return The in-progress asset request.
   */
  virtual CesiumAsync::Future<std::shared_ptr<IAssetRequest>> post(
      const AsyncSystem& asyncSystem,
      const std::string& url,
      const std::vector<THeader>& headers = std::vector<THeader>(),
      const gsl::span<const std::byte>& contentPayload = {}) = 0;

  /**
   * @brief Ticks the asset accessor system while the main thread is blocked.
   *
   * If the asset accessor is not dependent on the main thread to
   * dispatch requests, this method does not need to do anything.
   */
  virtual void tick() noexcept = 0;
};

} // namespace CesiumAsync

max_stars_repo_name: JiangMuWen/cesium-native

max_issues_repo_name: JiangMuWen/cesium-native

max_forks_repo_name: JiangMuWen/cesium-native


text:
/* -*- Mode:C++; c-file-style:"gnu"; indent-tabs-mode:nil; -*- */
/*
 * Copyright (c) 2010 The Boeing Company
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License version 2 as
 * published by the Free Software Foundation;
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * Author: Gary Pei <guangyu.pei@boeing.com>
 */

/*
 * This program is used to generate plots found in the paper
 * G. Pei and Tom Henderson, "Validation of ns-3 802.11b PHY model",
 * available online at http://www.nsnam.org/~pei/80211b.pdf
 *
 * It can be compiled as a C program and relies on a library installation of
 * the GNU Scientific Library (gsl).  To compile:
 *   gcc 80211b.c -o 80211b -lm -lgsl -lgslcblas
 *
 * The executable output should be redirected into a text file 80211b.txt
 *   ./80211b > 80211b.txt
 *
 * Then gnuplot can load the associated plot file which references 80211b.txt:
 *   gnuplot 80211b.plt
 */

#include "math.h"
#include "stdlib.h"
#include "stdio.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_sf_bessel.h>
#define min(a,b) ((a) < (b) ? (a) : (b))
#define max(a,b) ((a) > (b) ? (a) : (b))
#define WLAN_SIR_perfect 10.0 // if SIR > 10dB, perfect reception
#define WLAN_SIR_impossible 0.1 // if SIR < -10dB, impossible to receive

/**
 * \ingroup wifi
 * \defgroup wifi-test wifi module tests
 */

/**
 * \ingroup wifi-test
 * \ingroup tests
 *
 * \brief fn_parameter_t structure
 */
typedef struct fn_parameter_t
{
  double beta; ///< beta
  double n; ///< n
} fn_parameters;

double QFunction (double x)
{
  return 0.5 * erfc (x / sqrt (2.0));
}

double f (double x, void * params)
{
  double beta = ((fn_parameters *) params)->beta;
  double n = ((fn_parameters *) params)->n;
  double f = pow ( 2 * gsl_cdf_ugaussian_P (x + beta) - 1, n - 1) * exp (-x * x / 2.0) / sqrt (2.0 * M_PI);
  return f;
}

double p_e2 (double e2)
{
  double sep;
  double error;
  fn_parameters params;
  params.beta = sqrt (2.0 * e2);
  params.n = 8.0;
  gsl_integration_workspace* w = gsl_integration_workspace_alloc (1000);
  gsl_function F;
  F.function = &f;
  F.params = &params;
  gsl_integration_qagiu (&F, -params.beta, 0, 1e-7, 1000, w, &sep, &error);
  gsl_integration_workspace_free (w);
  if (error == 0.0)
    {
      sep = 1.0;
    }
  return 1.0 - sep;
}

double p_e1 (double e1)
{
  return 1.0 - pow ( 1.0 - p_e2 (e1 / 2.0), 2.0);
}

double DbToNoneDb (double x)
{
  return pow (10.0, x / 10.0);
}

double NoneDbToDb (double x)
{
  return 10.0 * log10 (x);
}

double DQPSKFunction (double x)
{
  double pi = acos (-1.0);
  return ( (sqrt (2.0) + 1.0) / sqrt (8.0 * pi * sqrt (2.0))) * (1.0 / sqrt (x)) * exp ( -(2.0 - sqrt (2.0)) * x);
}

double Get80211bDsssDbpskBerIeee (double EcNc)
{
  double ber;
  if (EcNc > WLAN_SIR_perfect)
    {
      ber = 0;
    }
  else if (EcNc < WLAN_SIR_impossible)
    {
      ber =  0.5;
    }
  else
    {
      ber = min (QFunction (sqrt (11.0 * EcNc)),0.5);
    }
  return ber;
}

double Get80211bDsssDbpskBer (double sinr)
{
  double EbN0 = sinr * 22000000.0 / 1000000.0;
  double ber = 0.5 * exp (-EbN0);
  return ber;
}

double Get80211bDsssDqpskBerIeee (double EcNc)
{
  double ber;
  if (EcNc > WLAN_SIR_perfect)
    {
      ber = 0;
    }
  else if (EcNc < WLAN_SIR_impossible)
    {
      ber = 0.5;
    }
  else
    {
      ber = min (QFunction (sqrt (5.5 * EcNc)),0.5);
    }
  return ber;
}

double Get80211bDsssDqpskBer (double sinr)
{
// 2 bits per symbol, 1 MSPS
  double EbN0 = sinr * 22000000.0 / 1000000.0 / 2.0;
  double ber = DQPSKFunction (EbN0);
  return ber;
}

double Get80211bDsssDqpskCCK5_5BerIeee (double EcNc)
{
  double ber;
  if (EcNc > WLAN_SIR_perfect)
    {
      ber = 0.0;
    }
  else if (EcNc < WLAN_SIR_impossible)
    {
      ber = 0.5;
    }
  else
    {
      double pew = 14.0 * QFunction (sqrt (EcNc * 8.0)) + QFunction (sqrt (EcNc * 16.0));
      pew = min (pew, 0.99999);
      ber = 8.0 / 15.0 * pew;
    }
  return ber;
}

double Get80211bDsssDqpskCCK11BerIeee (double EcNc)
{
  double ber;
  if (EcNc > WLAN_SIR_perfect)
    {
      ber = 0.0;
    }
  else if (EcNc < WLAN_SIR_impossible)
    {
      ber = 0.5;
    }
  else
    {
      double pew = 24.0 * QFunction (sqrt (EcNc * 4.0))  + 16.0 * QFunction (sqrt (EcNc * 6.0))  + 174.0 * QFunction (sqrt (EcNc * 8.0)) + 16.0 * QFunction (sqrt (EcNc * 10.0)) + 24.0 * QFunction (sqrt (EcNc * 12.0)) + QFunction (sqrt (EcNc * 16.0));
      pew = min (pew, 0.99999);
      ber = 128.0 / 255.0 * pew;
    }
  return ber;
}

int main (int argc, char * argv[])
{
  double rss, sinr;
  double totalPkt = 200.0;
//double noise = 1.552058;  // (dB) this noise figure value corresponds to
  // -99 dBm noise floor reported in CMU paper
  double noise = 7;  // (dB) this noise figure value corresponds to the
                     // default in YansWifiPhy, and matches CMU testbed results
  double EcNc, EbN01, EbN02, EbN05, EbN011;
  double ieee1,ieee2,ieee5,ieee11;
  double numBits = (1024. + 40. + 14.) * 8.;
  double dbpsk,dqpsk,cck16,cck256,sepcck16,sepcck256;
  noise = DbToNoneDb (noise) * 1.3803e-23 * 290.0 * 22000000;
  for (rss = -102.0; rss <= -80.0; rss += 0.1)
    {
      sinr   = DbToNoneDb (rss) / 1000.0 / noise;
      EcNc   = sinr * 22000000.0 / 11000000.0; // IEEE sir
      EbN01  = sinr * 22000000.0 / 1000000.0;
      // 2 bits per symbol, 1 MSPS
      EbN02  = sinr * 22000000.0 / 1000000.0 / 2.0;
      EbN05  = sinr * 22000000.0 / 1375000.0 / 4.0;
      EbN011 = sinr * 22000000.0 / 1375000.0 / 8.0;
      // 1=rss, 2=EcNc, 3=EbN01, 4=EbN02, 5=EBN05, 6=EbN011
      printf ("%g %g %g %g %g %g ", rss, NoneDbToDb (EcNc),
              NoneDbToDb (EbN01),NoneDbToDb (EbN02),
              NoneDbToDb (EbN05),NoneDbToDb (EbN011));
      ieee1  = Get80211bDsssDbpskBerIeee (EcNc);
      ieee2  = Get80211bDsssDqpskBerIeee (EcNc);
      ieee5  = Get80211bDsssDqpskCCK5_5BerIeee (EcNc);
      ieee11 = Get80211bDsssDqpskCCK11BerIeee (EcNc);
      // 7=ber_ieee1, 8=ber_ieee2, 9=ber_ieee5, 10=ber_ieee11
      printf (" %g %g %g %g ", ieee1, ieee2,ieee5,ieee11);
      ieee1  = totalPkt * pow (1 - ieee1,  numBits);
      ieee2  = totalPkt * pow (1 - ieee2,  numBits);
      ieee5  = totalPkt * pow (1 - ieee5,  numBits);
      ieee11 = totalPkt * pow (1 - ieee11, numBits);
      // 11=pkt_ieee1, 12=pkt_ieee2, 13=pkt_ieee5, 14=pkt_ieee11
      printf (" %g %g %g %g ", ieee1, ieee2,ieee5,ieee11);
      dbpsk  = Get80211bDsssDbpskBer (sinr);
      dqpsk  = Get80211bDsssDqpskBer (sinr);
      cck16  = max (0, 8.0 / 15.0 * p_e2 (4.0 * EbN05 / 2.0));
      cck256 = max (0, 128.0 / 255.0 * p_e1 (8.0 * EbN011 / 2.0));
      // 15=ber_dbpsk, 16=ber_dqpsk, 17=ber_cck16, 18=ber_cck256
      printf (" %g %g %g %g ", dbpsk, dqpsk,cck16,cck256);
      dbpsk     = totalPkt * pow (1 - dbpsk,numBits);
      dqpsk     = totalPkt * pow (1 - dqpsk,numBits);
      sepcck16  = p_e2 (4.0 * EbN05 / 2.0);
      sepcck256 = p_e1 (8.0 * EbN011 / 2.0);
      cck16     = totalPkt * pow (1.0 - sepcck16,numBits / 4.0);
      cck256    = totalPkt * pow (1.0 - sepcck256,numBits / 8.0);
      // 19=pkt_dbpsk, 20=pkt_dqpsk, 21=pkt_cck16, 22=pkt_cck256
      printf (" %g %g %g %g ", dbpsk, dqpsk,cck16,cck256);
      // 23=sinr
      printf (" %g \n",NoneDbToDb (sinr));
    }
  return 0;
}


text:
/*
 * This file is part of the Visual Computing Library (VCL) release under the
 * license.
 *
 * Copyright (c) 2017 Basil Fierz
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */
#pragma once

// VCL configuration
#include <vcl/config/global.h>

// C++ Standard library
#include <memory>
#include <string>
#include <tuple>
#include <vector>

// GSL
#include <gsl/gsl>

// Windows API
#define WIN32_LEAN_AND_MEAN
#include <Windows.h>
extern "C" {
#    include <hidsdi.h>
}

// VCL
#include <vcl/hid/device.h>

namespace Vcl { namespace HID { namespace Windows
{
	//! \note Structure definition taken from
	//! https://zfx.info/viewtopic.php?f=11&t=2977
	//! https://www.codeproject.com/Articles/297312/Minimal-Key-Logger-using-RAWINPUT
	//! https://www.codeproject.com/Articles/185522/Using-the-Raw-Input-API-to-Process-Joystick-Input
	struct Axis
	{
		//! Usage page as defined in the standard (e.g. "generic (0001)")
		USAGE usagePage;

		//! Usage of the axis as defined in the standard (e.g. "slider (0036)")
		USAGE usage;

		//! Index as defined through Hidp_GetData()
		USHORT index;

		//! Minimum value defined by the HID device
		int32_t                logicalMinimum;

		//! Maximum value defined by the HID device
		int32_t                logicalMaximum;

		//! Indicate wheter DirectInput calibration data was applied
		bool                   isCalibrated;

		//! Minimum value after calibration
		int32_t                logicalCalibratedMinimum;

		//! Maximum value after calibration
		int32_t                logicalCalibratedMaximum;

		//! Through calibration defined center value of the axis
		int32_t                logicalCalibratedCenter;

		//! Physical minimum value
		int32_t physicalMinimum;

		//! Physical maxiumum value
		int32_t physicalMaximum;

		//! Name as given by the driver
		std::wstring name;
	};

	struct Button
	{
		//! Usage page as defined in the standard (e.g. "buttons (0009)")
		USAGE usagePage;

		//! Usage of the axis as defined in the standard (e.g. "secondary (0002)")
		USAGE usage;

		//! Index as defined through Hidp_GetData()
		USHORT index;

		//! Name as given by the driver
		std::wstring name;
	};

	//! Normalize the axix value using the device configuration data
	inline float normalizeAxis(ULONG value, const Axis& axis)
	{
		if (static_cast<LONG>(value) < axis.logicalCalibratedCenter)
		{
			float range = static_cast<float>(axis.logicalCalibratedCenter - axis.logicalCalibratedMinimum);
			return (static_cast<LONG>(value) - axis.logicalCalibratedCenter) / range;
		}
		else
		{
			float range = static_cast<float>(axis.logicalCalibratedMaximum - axis.logicalCalibratedCenter);
			return (static_cast<LONG>(value) - axis.logicalCalibratedCenter) / range;
		}
	}

	class GenericHID
	{
	public:
		GenericHID(HANDLE raw_handle);

		//! Access the raw-input API handle
		//! \returns The input device handle
		HANDLE rawHandle() const { return _rawInputHandle; }

		//! Access the windows internal handle
		//! \returns The file handle to the device
		HANDLE fileHandle() const { return _fileHandle; }

		//! Access the vendor ID
		//! \returns The vendor ID
		DWORD vendorId() const { return _vendorId; }

		//! Access the product ID
		//! \returns The product ID
		DWORD productId() const { return _productId; }

		//! Read the device name from the hardware
		//! \returns The vendor defined names (vendor, product)
		auto readDeviceName() const -> std::pair<std::wstring, std::wstring>;

		const std::vector<Axis>& axes() const { return _axes; }
		const std::vector<Button>& buttons() const { return _buttons; }

		const std::vector<HIDP_VALUE_CAPS>& axisCaps() const { return _axesCaps; }
		const std::vector<HIDP_BUTTON_CAPS>& buttonCaps() const { return _buttonCaps; }

	private:
		//! Read the device capabilities
		auto readDeviceCaps() const -> std::tuple<std::vector<HIDP_BUTTON_CAPS>, std::vector<HIDP_VALUE_CAPS>>;

		//! Map buttons and fetch calibration data
		//! \param mapping Direct Input related remapping of buttons
		void readButtonCalibration(gsl::span<struct DirectInputButtonMapping> mapping) const;

		//! Map axes and fetch calibration data
		//! \param axes_caps
		void readAxisCalibration(const std::vector<HIDP_VALUE_CAPS>& axes_caps, gsl::span<struct DirectInputAxisMapping> mapping) const;

		//! Convert and store the button caps
		//! \param button_caps
		//! \param mapping Direct Input related remapping of buttons
		void storeButtons(std::vector<HIDP_BUTTON_CAPS>&& button_caps, gsl::span<struct DirectInputButtonMapping> mapping);

		//! Convert and store the axes caps
		//! \param axes_caps
		void storeAxes(std::vector<HIDP_VALUE_CAPS>&& axes_caps, gsl::span<struct DirectInputAxisMapping> mapping);

	private:
		//! Handle provided by the raw input API
		HANDLE _rawInputHandle{ nullptr };

		//! Handle from the file API
		HANDLE _fileHandle{ nullptr };

		//! Vendor ID
		DWORD _vendorId;

		//! Product ID
		DWORD _productId;

		//! Buttons associated with the device
		std::vector<Button> _buttons;

		//! Axes associated with the device
		std::vector<Axis> _axes;

		//! HID button representation
		std::vector<HIDP_BUTTON_CAPS> _buttonCaps;

		//! HID axis representation
		std::vector<HIDP_VALUE_CAPS> _axesCaps;
	};

	class AbstractHID
	{
	public:
		AbstractHID(std::unique_ptr<GenericHID> device) : _device{ std::move(device) } {}

		const GenericHID* device() const { return _device.get(); }

		virtual bool processInput(HWND window_handle, UINT input_code, PRAWINPUT raw_input) = 0;

	private:
		//! Actual hardware device implementation
		std::unique_ptr<GenericHID> _device;
	};

	template<typename JoystickType>
	class JoystickHID : public AbstractHID, public JoystickType
	{
	public:
		JoystickHID(std::unique_ptr<GenericHID> device);

		bool processInput(HWND window_handle, UINT input_code, PRAWINPUT raw_input) override;
	};

	template<typename GamepadType>
	class GamepadHID : public AbstractHID, public GamepadType
	{
	public:
		GamepadHID(std::unique_ptr<GenericHID> device);

		bool processInput(HWND window_handle, UINT input_code, PRAWINPUT raw_input) override;
	};

	template<typename ControllerType>
	class MultiAxisControllerHID : public AbstractHID, public ControllerType
	{
	public:
		MultiAxisControllerHID(std::unique_ptr<GenericHID> device);

		bool processInput(HWND window_handle, UINT input_code, PRAWINPUT raw_input) override;
	};

	class DeviceManager
	{
	public:
		DeviceManager();

		gsl::span<Device const* const> devices() const;

		//! Register devices with a specific window
		//! \param device_types Types of devices for which input should
		//!                     be processed.
		//! \param window_handle Handle to the window to which devices should be
		//!                      registered.
		void registerDevices(Flags<DeviceType> device_types, HWND window_handle);

		//! Polls all the devices instead of processing only a single device
		//! \param window_handle Handle of the window calling this method
		//! \returns True, if any device was successfully polled.
		//! \note This method is implemented according to the documenation
		//!       on MSDN. However, it seems not possible to call it.
		bool poll(HWND window_handle, UINT input_code);

		//! Process the input of a specific device
		bool processInput(HWND window_handle, UINT message, WPARAM wide_param, LPARAM low_param);

	private:
		//! List of Windows HID
		std::vector<std::unique_ptr<AbstractHID>> _devices;

		//! List of device pointers (links to `_devices`)
		std::vector<Device*> _deviceLinks;
	};
}}}


text:
/*
# This file is part of the Astrometry.net suite.
# Licensed under a 3-clause BSD style license - see LICENSE
*/

#ifndef AN_GSL_UTILS_H
#define AN_GSL_UTILS_H

#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_vector_double.h>

void gslutils_use_error_system();

/**
 Solves a least-squares matrix equation
 A X_i = B_i

 For NB pairs of X_i, B_i.

 NOTE: THIS DESTROYS A!

 A: MxN matrix
 B: array of NB x length-M vectors
 X: must be an array big enough to hold NB vectors.
   (they will be length-N).
 resids: if non-NULL, must be an array big enough to hold NB vectors
   (they will be length-M).

 The result vectors are freshly allocated and should be freed with gsl_vector_free().
 */
int gslutils_solve_leastsquares(gsl_matrix* A, gsl_vector** B,
                                gsl_vector** X, gsl_vector** resids,
                                int NB);

/**
 Same as above, but using varargs.  There must be exactly 3 * NB additional
 arguments, in the order:

 B0, &X0, &resid0, B1, &X1, &resid1, ...

 ie, the types must be repeating triples of:

 gsl_vector* b, gsl_vector** x, gsl_vector** resid

 */
int gslutils_solve_leastsquares_v(gsl_matrix* A, int NB, ...);

// C = A B
void gslutils_matrix_multiply(gsl_matrix* C, const gsl_matrix* A, const gsl_matrix* B);

int gslutils_invert_3x3(const double* A, double* B);

#endif