text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void compute_MCMC(MCMC_info *MCMC) {
    char          filename_output_dir[SID_MAX_FILENAME_LENGTH];
    char          filename_chain_dir[SID_MAX_FILENAME_LENGTH];
    char          filename_results_dir[SID_MAX_FILENAME_LENGTH];
    char          filename_plots_dir[SID_MAX_FILENAME_LENGTH];
    char          filename_run[SID_MAX_FILENAME_LENGTH];
    char          filename_chain[SID_MAX_FILENAME_LENGTH];
    char          filename_chain_config[SID_MAX_FILENAME_LENGTH];
    char          filename_chain_covariance[SID_MAX_FILENAME_LENGTH];
    char          filename_stats[SID_MAX_FILENAME_LENGTH];
    char          filename_coverage[SID_MAX_FILENAME_LENGTH];
    char          filename_histograms[SID_MAX_FILENAME_LENGTH];
    char          filename_results[SID_MAX_FILENAME_LENGTH];
    char          filename_stop[SID_MAX_FILENAME_LENGTH];
    char          column_txt[SID_MAX_FILENAME_LENGTH];
    char          problem_name_test[MCMC_NAME_SIZE];
    char          P_name_test[MCMC_NAME_SIZE];
    char          format_string[64];
    char          name_test[MCMC_NAME_SIZE];
    char          array_name_test[MCMC_NAME_SIZE];
    int           n_avg_test;
    int           flag_autocor_on_test;
    int           n_P_test;
    int           n_DS_test;
    int           n_arrays_test;
    int           n_M_test;
    double        P_init_test;
    double        M_target_test;
    double        dM_target_test;
    double        array_test;
    double        P_min_test;
    double        P_max_test;
    double        ln_Pr_min;
    double        ln_Pr_max;
    double        ln_Pr_avg;
    int           i_tune;
    int           i_P, j_P, k_P;
    int           i_M, j_M;
    int           ii_P, jj_P;
    int           i_avg;
    int           i_C;
    int           i_DS, n_DS;
    int           i_iteration, j_iteration;
    int           flag_continue;
    char          flag_success;
    int           flag_stop = GBP_FALSE;
    int           i_coverage;
    int           j_coverage;
    int           bin_x;
    int           bin_y;
    int           n_avg;
    int           n_avg_burn;
    int           n_avg_integrate;
    int           n_P, n_C;
    int           n_gals;
    double *      x_P;
    double **     x_M;
    double *      M_target;
    double *      dM_target;
    double *      P_min;
    double *      P_max;
    double *      P_new;
    double *      P_last;
    double *      P_init;
    double *      P_chain;
    double *      P_avg;
    double *      dP_avg;
    double **     M_best;
    double **     M_min;
    double **     M_max;
    double **     M_new;
    double **     M_lo_68;
    double **     M_hi_68;
    double **     M_lo_95;
    double **     M_hi_95;
    size_t ***    M_histogram;
    double **     M_last;
    double **     M_avg;
    double **     dM_avg;
    double *      slopes;
    double *      dP_sub;
    double *      ln_Pr_chain;
    double *      ln_Pr_proposals;
    double *      drift;
    double **     auto_cor;
    double **     P_chain_stats;
    double **     P_prop_stats;
    size_t **     P_histogram;
    size_t **     coverage_true;
    size_t **     coverage_false;
    size_t **     coverage_keep;
    double *      V_read;
    double        L_x, L_y, L_z;
    int           n_x, n_y, n_z;
    int           n_C_x;
    FILE *        fp_run;
    FILE *        fp_chain;
    FILE *        fp_chain_config;
    FILE *        fp_chain_covariance;
    FILE *        fp_stats;
    FILE *        fp_coverage;
    FILE *        fp_histograms;
    FILE *        fp_results;
    FILE *        fp_stop;
    int           seed = 182743;
    int           i_report;
    time_t        time_start, time_stop;
    double        time_diff;
    int           i_iteration_start;
    int           i_iteration_next_report;
    int           n_accepted;
    int           n_active;
    int           n_iterations;
    int           n_iterations_burn;
    int           n_iterations_integrate;
    int           n_iterations_phase;
    int           n_thin;
    int           n_thin_burn;
    int           n_thin_integrate;
    int           i_thin;
    int           n_chains_test;
    int           flag_autocor_on;
    int           n_coverage;
    int           coverage_size;
    int           i_phase;
    int           i_array;
    int *         n_M;
    int           flag_initialized;
    int           flag_no_map_write;
    int           flag_no_map_write_test;
    int           n_used;
    double        RN;
    double ***    P_arrays;
    int *         n_P_arrays;
    double        constant;
    size_t *      histogram_index;
    size_t *      coverage_keep_index;
    size_t        P_best_index;
    size_t        M_best_index;
    double        accumulator;
    size_t        P_lo_68_index;
    size_t        P_hi_68_index;
    size_t        P_lo_95_index;
    size_t        P_hi_95_index;
    size_t        M_lo_68_index;
    size_t        M_hi_68_index;
    size_t        M_lo_95_index;
    size_t        M_hi_95_index;
    double *      P_lo_68;
    double *      P_hi_68;
    double *      P_lo_95;
    double *      P_hi_95;
    double *      P_contour_68;
    double *      P_contour_95;
    double *      P_best;
    int           my_chain;
    int           flag_init;
    int           i_column;
    int           dummy_t;
    int           n_iterations_file_total;
    int           n_iterations_file_burn;
    int           flag_restart     = GBP_FALSE;
    int           flag_minimize_IO = GBP_FALSE;
    MCMC_DS_info *current_DS;
    MCMC_DS_info *next_DS;

    size_t i_iteration_buffer;
    size_t i_P_buffer;
    size_t i_M_buffer;

    SID_log("Performing MCMC...", SID_LOG_OPEN | SID_LOG_TIMER);

    SID_log("Initializing...", SID_LOG_OPEN);

    // Initialize the covariance matrix
    set_MCMC_covariance(MCMC, NULL); // First call with NULL must be after P_init is set!

    // Initialize arrays
    init_MCMC_arrays(MCMC);

    // Grab some stuff from the MCMC structure
    n_P                    = MCMC->n_P;
    n_DS                   = MCMC->n_DS;
    n_avg                  = MCMC->n_avg;
    n_thin                 = MCMC->n_thin;
    coverage_size          = MCMC->coverage_size;
    flag_autocor_on        = MCMC->flag_autocor_on;
    my_chain               = MCMC->my_chain;
    n_iterations           = MCMC->n_iterations;
    n_iterations_burn      = MCMC->n_iterations_burn;
    n_iterations_integrate = n_iterations - n_iterations_burn;
    flag_no_map_write      = SID_CHECK_BITFIELD_SWITCH(MCMC->mode, MCMC_MODE_NO_MAP_WRITE);
    flag_minimize_IO       = SID_CHECK_BITFIELD_SWITCH(MCMC->mode, MCMC_MODE_MINIMIZE_IO);

    SID_log("temperature            = %lf", SID_LOG_COMMENT, MCMC->temperature);
    SID_log("n_avg                  = %d", SID_LOG_COMMENT, n_avg);
    SID_log("n_thin                 = %d", SID_LOG_COMMENT, n_thin);
    SID_log("coverage_size          = %d", SID_LOG_COMMENT, coverage_size);
    SID_log("n_iterations           = %d", SID_LOG_COMMENT, n_iterations);
    SID_log("n_iterations_burn      = %d", SID_LOG_COMMENT, n_iterations_burn);
    SID_log("n_iterations_integrate = %d", SID_LOG_COMMENT, n_iterations_integrate);

    // Initialize random number generator
    SID_log("random seed            = %d", SID_LOG_COMMENT, MCMC->seed);

    if(MCMC->V != NULL)
        SID_log("Covariance matrix is initialized.", SID_LOG_COMMENT);
    else
        SID_log("Covariance matrix is NOT initialized.", SID_LOG_COMMENT);
    if(flag_autocor_on)
        SID_log("Auto-correlation  is on.", SID_LOG_COMMENT);
    else
        SID_log("Auto-correlation  is off.", SID_LOG_COMMENT);
    if(flag_minimize_IO)
        SID_log("Minimize I/O      is on.", SID_LOG_COMMENT);
    else
        SID_log("Minimize I/O      is off.", SID_LOG_COMMENT);

    // Initialize chain arrays
    n_M             = MCMC->n_M;
    M_new           = MCMC->M_new;
    M_last          = MCMC->M_last;
    P_new           = MCMC->P_new;
    P_last          = MCMC->P_last;
    P_init          = MCMC->P_init;
    P_chain         = MCMC->P_chain;
    P_min           = (double *)SID_malloc(sizeof(double) * n_P);
    P_max           = (double *)SID_malloc(sizeof(double) * n_P);
    P_avg           = (double *)SID_malloc(sizeof(double) * n_P);
    dP_avg          = (double *)SID_malloc(sizeof(double) * n_P);
    ln_Pr_chain     = (double *)SID_malloc(sizeof(double) * n_avg);
    ln_Pr_proposals = (double *)SID_malloc(sizeof(double) * n_avg);
    drift           = (double *)SID_malloc(sizeof(double) * n_P);
    slopes          = (double *)SID_malloc(sizeof(double) * n_P);
    dP_sub          = (double *)SID_malloc(sizeof(double) * n_P);
    P_chain_stats   = (double **)SID_malloc(sizeof(double *) * n_P);
    P_prop_stats    = (double **)SID_malloc(sizeof(double *) * n_P);
    if(flag_autocor_on)
        auto_cor = (double **)SID_malloc(sizeof(double *) * n_P);
    else
        auto_cor = NULL;
    for(i_P = 0; i_P < n_P; i_P++) {
        drift[i_P]         = 0.;
        slopes[i_P]        = 0.;
        dP_sub[i_P]        = 0.;
        P_chain_stats[i_P] = (double *)SID_malloc(sizeof(double) * n_avg);
        P_prop_stats[i_P]  = (double *)SID_malloc(sizeof(double) * n_avg);
        if(flag_autocor_on)
            auto_cor[i_P] = (double *)SID_malloc(sizeof(double) * n_avg);
    }

    // Set directories
    sprintf(filename_output_dir, "%s/", MCMC->filename_output_dir);
    sprintf(filename_chain_dir, "%s/chains/", filename_output_dir);
    sprintf(filename_results_dir, "%s/results/", filename_output_dir);
    sprintf(filename_plots_dir, "%s/plots/", filename_output_dir);

    // Set filenames
    sprintf(filename_run, "%s/run.dat", filename_output_dir);
    sprintf(filename_chain, "%s/chain_trace_%06d.dat", filename_chain_dir, my_chain);
    sprintf(filename_chain_config, "%s/chain_config_%06d.dat", filename_chain_dir, my_chain);
    sprintf(filename_chain_covariance, "%s/chain_covariance_%06d.dat", filename_chain_dir, my_chain);
    sprintf(filename_stats, "%s/chain_stats_%06d.dat", filename_chain_dir, my_chain);
    sprintf(filename_coverage, "%s/coverage.dat", filename_results_dir);
    sprintf(filename_histograms, "%s/histograms.dat", filename_results_dir);
    sprintf(filename_stop, "%s/stop", filename_output_dir);

    SID_log("Done.", SID_LOG_CLOSE);

    // Perform integration
    SID_log("Initializing integration...", SID_LOG_OPEN);
    MCMC->flag_integrate_on = GBP_TRUE;

    // Make sure the needed directories exist
    mkdir(filename_output_dir, 02755);
    mkdir(filename_chain_dir, 02755);
    mkdir(filename_results_dir, 02755);
    mkdir(filename_plots_dir, 02755);

    // Read/Write Header file
    if(SID.I_am_Master) {
        // If run.dat already exists, this is a restart.  Check that it is consistant with the current run.
        if((fp_run = fopen(filename_run, "rb")) != NULL) {
            flag_restart = GBP_TRUE;
            SID_log("Checking the consistancy of this run with the previous run...", SID_LOG_OPEN);
            SID_fread_verify(problem_name_test, sizeof(char), MCMC_NAME_SIZE, fp_run);
            if(strcmp(problem_name_test, MCMC->problem_name))
                SID_exit_error("Problem names are inconsistant (i.e. {%s}!={%s}).", SID_ERROR_LOGIC, MCMC->problem_name,
                               problem_name_test);
            SID_fread_verify(&n_chains_test, sizeof(int), 1, fp_run);
            MCMC->n_chains = GBP_MAX(MCMC->n_chains, n_chains_test);
            SID_fread_verify(&n_avg_test, sizeof(int), 1, fp_run);
            if(n_avg_test != n_avg)
                SID_exit_error("Integration averaging intervals are inconsistant (i.e. %d!=%d).", SID_ERROR_LOGIC,
                               n_avg, n_avg_test);
            SID_fread_verify(&flag_autocor_on_test, sizeof(int), 1, fp_run);
            if(flag_autocor_on_test != flag_autocor_on)
                SID_exit_error("Autocorrelation flags are inconsistant (i.e. %d!=%d).", SID_ERROR_LOGIC,
                               flag_autocor_on, flag_autocor_on_test);
            SID_fread_verify(&flag_no_map_write_test, sizeof(int), 1, fp_run);
            if(flag_no_map_write_test != flag_no_map_write)
                SID_exit_error("Map write flags are inconsistant (i.e. %d!=%d).", SID_ERROR_LOGIC, flag_no_map_write,
                               flag_no_map_write_test);
            SID_fread_verify(&n_P_test, sizeof(int), 1, fp_run);
            if(n_P_test != n_P)
                SID_exit_error("The number of paramaters is inconsistant (i.e. %d!=%d).", SID_ERROR_LOGIC, n_P,
                               n_P_test);
            MCMC->P_name_length = 0;
            for(i_P = 0; i_P < n_P; i_P++) {
                SID_fread_verify(P_name_test, sizeof(char), MCMC_NAME_SIZE, fp_run);
                if(strcmp(P_name_test, MCMC->P_names[i_P]))
                    SID_exit_error("Parameter #%d's names are inconsistant (i.e. {%s}!={%s}).", SID_ERROR_LOGIC, i_P,
                                   MCMC->P_names[i_P], P_name_test);
                SID_fread_verify(&P_init_test, sizeof(double), 1, fp_run);
                // if(P_init_test!=MCMC->P_init[i_P])
                //   SID_exit_error("Parameter #%d's initial values are inconsistant (i.e.
                //   %le!=%le).",SID_ERROR_LOGIC,i_P,MCMC->P_init[i_P],P_init_test);
                SID_fread_verify(&P_min_test, sizeof(double), 1, fp_run);
                if(P_min_test != MCMC->P_limit_min[i_P])
                    SID_exit_error("Parameter #%d's minimum values are inconsistant (i.e. %le!=%le).", SID_ERROR_LOGIC,
                                   i_P, MCMC->P_limit_min[i_P], P_min_test);
                SID_fread_verify(&P_max_test, sizeof(double), 1, fp_run);
                if(P_max_test != MCMC->P_limit_max[i_P])
                    SID_exit_error("Parameter #%d's maximum values are inconsistant (i.e. %le!=%le).", SID_ERROR_LOGIC,
                                   i_P, MCMC->P_limit_max[i_P], P_max_test);
                MCMC->P_name_length = GBP_MAX((size_t)(MCMC->P_name_length), strlen(MCMC->P_names[i_P]));
            }
            sprintf(MCMC->P_name_format, "%%-%ds", MCMC->P_name_length);
            SID_fread_verify(&n_arrays_test, sizeof(int), 1, fp_run);
            if(n_arrays_test != MCMC->n_arrays)
                SID_exit_error("Numbers of project arrays are inconsistant (i.e. %d!=%d).", SID_ERROR_LOGIC,
                               MCMC->n_arrays, n_arrays_test);
            for(i_array = 0; i_array < MCMC->n_arrays; i_array++) {
                SID_fread_verify(&array_name_test, sizeof(char), MCMC_NAME_SIZE, fp_run);
                if(strcmp(array_name_test, MCMC->array_name[i_array]))
                    SID_exit_error("Project array names are inconsisitant (i.e. {%s}!={%s}).", SID_ERROR_LOGIC,
                                   MCMC->array_name[i_array], array_name_test);
                for(i_P = 0; i_P < n_P; i_P++) {
                    SID_fread_verify(&array_test, sizeof(double), 1, fp_run);
                    if(array_test != MCMC->array[i_array][i_P])
                        SID_exit_error("Project array #%d element #%d is inconsistant (i.e. %le!=%le).",
                                       SID_ERROR_LOGIC, i_array, i_P, MCMC->array[i_array][i_P], array_test);
                }
            }
            SID_fread_verify(&n_DS_test, sizeof(int), 1, fp_run);
            if(n_DS_test != n_DS)
                SID_exit_error("The number of datasets is inconsistant (i.e. %d!=%d).", SID_ERROR_LOGIC, n_DS,
                               n_DS_test);
            current_DS = MCMC->DS;
            i_DS       = 0;
            while(current_DS != NULL) {
                next_DS = current_DS->next;
                SID_fread_verify(name_test, sizeof(char), MCMC_NAME_SIZE, fp_run);
                if(strcmp(name_test, current_DS->name))
                    SID_exit_error("Dataset #%d's names are inconsistant (i.e. {%s}!={%s}).", SID_ERROR_LOGIC, i_DS,
                                   current_DS->name, name_test);
                SID_fread_verify(&n_M_test, sizeof(int), 1, fp_run);
                if(n_M_test != n_M[i_DS])
                    SID_exit_error("The sizes of dataset #%d are inconsistant (i.e. %d!=%d).", SID_ERROR_LOGIC, i_DS,
                                   n_M_test, n_M[i_DS]);
                for(i_M = 0; i_M < current_DS->n_M; i_M++) {
                    SID_fread_verify(&M_target_test, sizeof(double), 1, fp_run);
                    if(M_target_test != current_DS->M_target[i_M])
                        SID_exit_error("Dataset #%d, element #%d is inconsistant (i.e. %le!=%le).", SID_ERROR_LOGIC,
                                       i_DS, i_M, current_DS->M_target[i_M], M_target_test);
                }
                for(i_M = 0; i_M < current_DS->n_M; i_M++) {
                    SID_fread_verify(&dM_target_test, sizeof(double), 1, fp_run);
                    if(dM_target_test != current_DS->dM_target[i_M])
                        SID_exit_error("Dataset #%d, uncertainty element #%d is inconsistant (i.e. %le!=%le).",
                                       SID_ERROR_LOGIC, i_DS, i_M, current_DS->M_target[i_M], M_target_test);
                }
                SID_fread_verify(&n_arrays_test, sizeof(int), 1, fp_run);
                if(n_arrays_test != current_DS->n_arrays)
                    SID_exit_error("The number of arrays in dataset #%d is inconsistant (i.e. %d!=%d).",
                                   SID_ERROR_LOGIC, i_DS, n_arrays_test, current_DS->n_arrays);
                for(i_array = 0; i_array < current_DS->n_arrays; i_array++) {
                    SID_fread_verify(array_name_test, sizeof(char), MCMC_NAME_SIZE, fp_run);
                    if(strcmp(array_name_test, current_DS->array_name[i_array]))
                        SID_exit_error("Array name #%d for dataset #%d is inconsisitant (i.e. {%s}!={%s}).",
                                       SID_ERROR_LOGIC, i_array, i_DS, current_DS->array_name[i_array], array_name_test);
                    for(i_M = 0; i_M < current_DS->n_M; i_M++) {
                        SID_fread_verify(&array_test, sizeof(double), 1, fp_run);
                        if(array_test != current_DS->array[i_array][i_M])
                            SID_exit_error("Array #%d, element #%d for dataset #%d is inconsistant (i.e. %le!=%le).",
                                           SID_ERROR_LOGIC, i_array, i_M, i_DS, current_DS->array[i_array][i_M],
                                           array_test);
                    }
                }
                current_DS = next_DS;
                i_DS++;
            }
            fclose(fp_run);
            SID_log("Done.", SID_LOG_CLOSE);
        }
        // ... else if run.dat does not exist, create it
        else {
            SID_log("Write header file...", SID_LOG_OPEN);
            if((fp_run = fopen(filename_run, "wb")) == NULL)
                SID_exit_error("Could not open file for writing {%s}.", SID_ERROR_IO_OPEN, filename_run);
            // Stuff relating to this MCMC project
            fwrite(MCMC->problem_name, sizeof(char), MCMC_NAME_SIZE, fp_run);
            fwrite(&(MCMC->n_chains), sizeof(int), 1, fp_run);
            fwrite(&n_avg, sizeof(int), 1, fp_run);
            fwrite(&flag_autocor_on, sizeof(int), 1, fp_run);
            fwrite(&flag_no_map_write, sizeof(int), 1, fp_run);
            fwrite(&n_P, sizeof(int), 1, fp_run);
            for(i_P = 0; i_P < n_P; i_P++) {
                fwrite(MCMC->P_names[i_P], sizeof(char), MCMC_NAME_SIZE, fp_run);
                fwrite(&(MCMC->P_init[i_P]), sizeof(double), 1, fp_run);
                fwrite(&(MCMC->P_limit_min[i_P]), sizeof(double), 1, fp_run);
                fwrite(&(MCMC->P_limit_max[i_P]), sizeof(double), 1, fp_run);
            }
            fwrite(&(MCMC->n_arrays), sizeof(int), 1, fp_run);
            for(i_array = 0; i_array < MCMC->n_arrays; i_array++) {
                fwrite(MCMC->array_name[i_array], sizeof(char), MCMC_NAME_SIZE, fp_run);
                fwrite(MCMC->array[i_array], sizeof(double), n_P, fp_run);
            }
            // Stuff relating to the constraining datasets
            fwrite(&n_DS, sizeof(int), 1, fp_run);
            current_DS = MCMC->DS;
            while(current_DS != NULL) {
                next_DS = current_DS->next;
                fwrite(current_DS->name, sizeof(char), MCMC_NAME_SIZE, fp_run);
                fwrite(&(current_DS->n_M), sizeof(int), 1, fp_run);
                fwrite(current_DS->M_target, sizeof(double), current_DS->n_M, fp_run);
                fwrite(current_DS->dM_target, sizeof(double), current_DS->n_M, fp_run);
                fwrite(&(current_DS->n_arrays), sizeof(int), 1, fp_run);
                for(i_array = 0; i_array < current_DS->n_arrays; i_array++) {
                    fwrite(current_DS->array_name[i_array], sizeof(char), MCMC_NAME_SIZE, fp_run);
                    fwrite(current_DS->array[i_array], sizeof(double), current_DS->n_M, fp_run);
                }
                current_DS = next_DS;
            }
            fclose(fp_run);
            SID_log("Done.", SID_LOG_CLOSE);
        }
    }
    SID_Bcast(&flag_restart, 1, SID_INT, SID_MASTER_RANK, MCMC->comm);
    SID_Bcast(&(MCMC->n_chains), 1, SID_INT, SID_MASTER_RANK, MCMC->comm);
    SID_Bcast(&(MCMC->P_name_length), 1, SID_INT, SID_MASTER_RANK, MCMC->comm);
    sprintf(MCMC->P_name_format, "%%-%ds", MCMC->P_name_length);

    // If this is NOT a restart, start from scratch ...
    n_iterations_file_total = 0;
    n_iterations_file_burn  = 0;
    if(!flag_restart) {
        SID_log("Initializing state...", SID_LOG_OPEN);

        // Report the starting conditions
        if(my_chain == SID.My_rank) {
            SID_log("Initial parameters & priors:", SID_LOG_ALLRANKS | SID_LOG_OPEN);
            sprintf(format_string, "%s = %%13.6le (%%13.6le ->%%13.6le)", MCMC->P_name_format);
            for(i_P = 0; i_P < n_P; i_P++)
                SID_log(format_string,
                        SID_LOG_ALLRANKS | SID_LOG_COMMENT,
                        MCMC->P_names[i_P],
                        P_last[i_P],
                        MCMC->P_limit_min[i_P],
                        MCMC->P_limit_max[i_P]);
            SID_log("", SID_LOG_ALLRANKS | SID_LOG_NOPRINT | SID_LOG_CLOSE);
        }

        // Perform autotuning (if requested)
        if(SID_CHECK_BITFIELD_SWITCH(MCMC->mode, MCMC_MODE_AUTOTUNE))
            autotune_MCMC(MCMC);
        MCMC->flag_init_chain = GBP_TRUE;

        // Set the initial state
        SID_log("Writing chain config file...", SID_LOG_OPEN);
        if((fp_chain = fopen(filename_chain, "wb")) == NULL)
            SID_exit_error("Could not open file for writing {%s}.", SID_ERROR_IO_OPEN, filename_chain);
        if((fp_stats = fopen(filename_stats, "wb")) == NULL)
            SID_exit_error("Could not open file for writing {%s}.", SID_ERROR_IO_OPEN, filename_stats);
        if((fp_chain_config = fopen(filename_chain_config, "wb")) == NULL)
            SID_exit_error("Could not open file for writing {%s}.", SID_ERROR_IO_OPEN, filename_chain_config);
        fwrite(&n_iterations_file_total, sizeof(int), 1, fp_chain_config);
        fwrite(&n_iterations_file_burn, sizeof(int), 1, fp_chain_config);
        fwrite(&(MCMC->temperature), sizeof(double), 1, fp_chain_config);
        fwrite(MCMC->V, sizeof(double), n_P * n_P, fp_chain_config);
        fclose(fp_chain_config);
        SID_log("Done.", SID_LOG_CLOSE);

        SID_log("Done.", SID_LOG_CLOSE);
    }
    // ... else read the state we left-off from ...
    else if(my_chain == SID.My_rank) {
        SID_log("Loading previous state...", SID_LOG_OPEN);
        MCMC->flag_init_chain = GBP_FALSE;
        // ... fetch the number of intervals that have already been computed ...
        n_iterations_file_total = 0;
        if((fp_chain_config = fopen(filename_chain_config, "rb")) != NULL) {
            SID_log("Reading the existant number of iterations...", SID_LOG_OPEN);
            V_read = (double *)SID_malloc(sizeof(double) * n_P * n_P);
            SID_fread_verify(&n_iterations_file_total, sizeof(int), 1, fp_chain_config);
            SID_fread_verify(&n_iterations_file_burn, sizeof(int), 1, fp_chain_config);
            SID_fread_verify(&(MCMC->temperature), sizeof(double), 1, fp_chain_config);
            SID_fread_verify(V_read, sizeof(double), n_P * n_P, fp_chain_config);
            set_MCMC_covariance(MCMC, V_read);
            SID_free(SID_FARG V_read);
            SID_log("# burn  iterations = %d (%d requested)", SID_LOG_COMMENT, n_iterations_file_burn, n_iterations_burn);
            SID_log("# total iterations = %d (%d requested)", SID_LOG_COMMENT, n_iterations_file_total, n_iterations);
            SID_log("Temperature        = %le", SID_LOG_COMMENT, MCMC->temperature);
            fclose(fp_chain_config);
            SID_log("Done.", SID_LOG_CLOSE);
        }

        // ... scan to the end of the chain and read the last-used parameter set
        if(n_iterations_file_total > 0) {
            SID_log("Reading the last-used parameter set...", SID_LOG_OPEN);
            fp_chain = fopen(filename_chain, "rb");
            for(i_iteration = 0; i_iteration < n_iterations_file_total; i_iteration++) {
                for(i_avg = 0; i_avg < n_avg; i_avg++) {
                    SID_fread_verify(&flag_success, sizeof(char), 1, fp_chain);
                    if((i_iteration + i_avg) != 0) {
                        MCMC->ln_likelihood_last = MCMC->ln_likelihood_new;
                        memcpy(P_last, P_new, n_P * sizeof(double));
                        if(!flag_no_map_write) {
                            for(i_DS = 0; i_DS < n_DS; i_DS++)
                                memcpy(M_last[i_DS], M_new[i_DS], n_M[i_DS] * sizeof(double));
                        }
                    }
                    SID_fread_verify(&(MCMC->ln_likelihood_new), sizeof(double), 1, fp_chain);
                    SID_fread_verify(P_new, sizeof(double), n_P, fp_chain);
                    if(!flag_no_map_write) {
                        for(i_DS = 0; i_DS < n_DS; i_DS++)
                            SID_fread_verify(M_new[i_DS], sizeof(double), n_M[i_DS], fp_chain);
                    }
                    if((i_iteration + i_avg) == 0) {
                        MCMC->ln_likelihood_chain = MCMC->ln_likelihood_new;
                        MCMC->ln_likelihood_last  = MCMC->ln_likelihood_new;
                    }
                    if(flag_success) {
                        MCMC->ln_likelihood_chain = MCMC->ln_likelihood_new;
                        memcpy(P_chain, P_new, n_P * sizeof(double));
                    }
                }
            }
            fclose(fp_chain);
            fp_chain = fopen(filename_chain, "ab");
            SID_log("Done.", SID_LOG_CLOSE);

            // ... scan to the end of the chain stats and read the last-generated statistics (particularly the drift)
            SID_log("Reading the last-generated statistics...", SID_LOG_OPEN);
            fp_stats = fopen(filename_stats, "rb");
            for(i_iteration = 0; i_iteration < n_iterations_file_total; i_iteration++) {
                SID_fread_verify(P_min, sizeof(double), n_P, fp_stats); // Read proposal stats
                SID_fread_verify(P_avg, sizeof(double), n_P, fp_stats);
                SID_fread_verify(P_max, sizeof(double), n_P, fp_stats);
                SID_fread_verify(dP_avg, sizeof(double), n_P, fp_stats);
                SID_fread_verify(dP_sub, sizeof(double), n_P, fp_stats);
                SID_fread_verify(&ln_Pr_min, sizeof(double), 1, fp_stats);
                SID_fread_verify(&ln_Pr_avg, sizeof(double), 1, fp_stats);
                SID_fread_verify(&ln_Pr_max, sizeof(double), 1, fp_stats);
                if(flag_autocor_on)
                    SID_fread_verify(auto_cor, sizeof(double), n_avg - 1, fp_stats);
                SID_fread_verify(P_min, sizeof(double), n_P, fp_stats); // Read chain stats
                SID_fread_verify(P_avg, sizeof(double), n_P, fp_stats);
                SID_fread_verify(P_max, sizeof(double), n_P, fp_stats);
                SID_fread_verify(dP_avg, sizeof(double), n_P, fp_stats);
                SID_fread_verify(dP_sub, sizeof(double), n_P, fp_stats);
                SID_fread_verify(&ln_Pr_min, sizeof(double), 1, fp_stats);
                SID_fread_verify(&ln_Pr_avg, sizeof(double), 1, fp_stats);
                SID_fread_verify(&ln_Pr_max, sizeof(double), 1, fp_stats);
                if(flag_autocor_on)
                    SID_fread_verify(auto_cor, sizeof(double), n_avg - 1, fp_stats);
                SID_fread_verify(slopes, sizeof(double), n_P, fp_stats);
                SID_fread_verify(drift, sizeof(double), n_P, fp_stats);
            }
            fclose(fp_stats);
            fp_stats = fopen(filename_stats, "ab");
            SID_log("Done.", SID_LOG_CLOSE);
        } else {
            memcpy(P_last, P_init, n_P * sizeof(double));
            MCMC->flag_init_chain = GBP_TRUE;
            fp_chain              = fopen(filename_chain, "wb");
            fp_stats              = fopen(filename_stats, "wb");
        }

        // Report the starting conditions
        if(my_chain == SID.My_rank) {
            SID_log("Resume from parameters:", SID_LOG_ALLRANKS | SID_LOG_OPEN);
            sprintf(format_string, "%s = %%13.6le", MCMC->P_name_format);
            for(i_P = 0; i_P < n_P; i_P++)
                SID_log(format_string, SID_LOG_ALLRANKS | SID_LOG_COMMENT, MCMC->P_names[i_P], P_last[i_P]);
            SID_log("", SID_LOG_ALLRANKS | SID_LOG_NOPRINT | SID_LOG_CLOSE);
        }

        SID_log("Done.", SID_LOG_CLOSE);

    } // End of restart stuff
    SID_Bcast(&n_iterations_file_total, 1, SID_INT, my_chain, MCMC->comm);
    SID_Bcast(&n_iterations_file_burn, 1, SID_INT, my_chain, MCMC->comm);
    SID_Bcast(&n_iterations, 1, SID_INT, my_chain, MCMC->comm);
    SID_Bcast(&(MCMC->flag_init_chain), 1, SID_INT, my_chain, MCMC->comm);
    SID_Bcast(&(MCMC->ln_likelihood_chain), 1, SID_DOUBLE, my_chain, MCMC->comm);
    SID_Bcast(&(MCMC->ln_likelihood_new), 1, SID_DOUBLE, my_chain, MCMC->comm);

    // Remove the existant iterations from the totals we need to perform still
    if(n_iterations_file_total < n_iterations_burn) {
        i_phase     = 0;
        i_iteration = n_iterations_file_total;
    } else if(n_iterations_file_total < n_iterations) {
        i_phase     = 1;
        i_iteration = n_iterations_file_total - n_iterations_burn;
    } else {
        i_phase     = 2; // ... in other words, we're done; skip the integration.
        i_iteration = n_iterations;
    }
    i_iteration_start = i_iteration;

    // This is the end of the initialization
    SID_log("Done.", SID_LOG_CLOSE);

    // Create the chain in 2 phases: a burn-in phase and an integration phase
    SID_log("Performing integration...", SID_LOG_OPEN | SID_LOG_TIMER);
    i_iteration_buffer = 0;
    i_P_buffer         = 0;
    i_M_buffer         = 0;
    for(; i_phase < 2; i_phase++) {
        switch(i_phase) {
            case 0:
                SID_log("Performing burn-in phase...", SID_LOG_OPEN | SID_LOG_TIMER);
                n_iterations_phase = n_iterations_burn;
                break;
            case 1:
                SID_log("Performing integration phase...", SID_LOG_OPEN | SID_LOG_TIMER);
                n_iterations_phase = n_iterations - n_iterations_burn;
                break;
        }
        // Initialize progress reporting
        int n_report = 10; // This is the number of progress reports that will be creates. eg. 10 -> every 10%
        i_report     = 0;
        if(n_iterations_phase > 20)
            i_iteration_next_report = i_iteration_start + (n_iterations_phase - i_iteration_start) / n_report;
        else
            i_iteration_next_report = 20;

        // Loop until this phase (burn or integrate) is done
        flag_continue = GBP_TRUE;
        while(flag_continue) {
            // Process one averaging interval at a time

            // Start timer
            time(&time_start);

            for(i_avg = 0, i_thin = 1; i_avg < n_avg; i_thin++) {
                // Generate new proposal and determine it's chi^2
                flag_success = generate_MCMC_chain(MCMC);

                // If this rank is processing chain information ...
                if(my_chain == SID.My_rank) {
                    // ... and this proposal is kept by the thining then ...
                    if(i_thin == n_thin) {
                        // ... reformat parameter arrays for statistics generation
                        for(i_P = 0; i_P < n_P; i_P++) {
                            P_chain_stats[i_P][i_avg % n_avg] = P_chain[i_P];
                            P_prop_stats[i_P][i_avg % n_avg]  = P_new[i_P];
                        }
                        ln_Pr_chain[i_avg % n_avg]     = MCMC->ln_Pr_chain;
                        ln_Pr_proposals[i_avg % n_avg] = MCMC->ln_Pr_new;

                        // ... and write to the chain file
                        switch(flag_minimize_IO) {
                            case GBP_TRUE:
                                MCMC->flag_success_buffer[i_iteration_buffer]      = flag_success;
                                MCMC->ln_likelihood_new_buffer[i_iteration_buffer] = MCMC->ln_likelihood_new;
                                i_iteration_buffer++;
                                memcpy(&(MCMC->P_new_buffer[i_P_buffer]), P_new, (size_t)n_P * sizeof(double));
                                i_P_buffer += n_P;
                                if(!flag_no_map_write) {
                                    for(i_DS = 0; i_DS < n_DS; i_DS++) {
                                        memcpy(&(MCMC->M_new_buffer[i_M_buffer]), M_new[i_DS], (size_t)n_M[i_DS] * sizeof(double));
                                        i_M_buffer += n_M[i_DS];
                                    }
                                }
                                break;
                            default:
                                fwrite(&flag_success, sizeof(char), 1, fp_chain);
                                fwrite(&(MCMC->ln_likelihood_new), sizeof(double), 1, fp_chain);
                                fwrite(P_new, sizeof(double), n_P, fp_chain);
                                if(!flag_no_map_write) {
                                    for(i_DS = 0; i_DS < n_DS; i_DS++)
                                        fwrite(M_new[i_DS], sizeof(double), n_M[i_DS], fp_chain);
                                }
                                break;
                        }
                    }
                }

                // Change counters at the end of a thining interval
                if(i_thin == n_thin) {
                    i_thin = 0;
                    i_avg++;
                }
            } // End of averaging interval loop
            if(i_phase == 0)
                n_iterations_file_burn++;
            n_iterations_file_total++;
            i_iteration++;

            // Stop timer
            time(&time_stop);
            time_diff = difftime(time_stop, time_start);

            // Generate statistics for the averaging interval we just completed
            if(my_chain == SID.My_rank) {
                // Write the statistics to the chain stats file
                compute_MCMC_chain_stats(P_prop_stats,
                                         n_P,
                                         n_avg,
                                         P_min,
                                         P_avg,
                                         P_max,
                                         dP_avg,
                                         auto_cor,
                                         slopes,
                                         dP_sub,
                                         ln_Pr_proposals,
                                         &ln_Pr_min,
                                         &ln_Pr_avg,
                                         &ln_Pr_max);
                fwrite(P_min, sizeof(double), n_P, fp_stats);
                fwrite(P_avg, sizeof(double), n_P, fp_stats);
                fwrite(P_max, sizeof(double), n_P, fp_stats);
                fwrite(dP_avg, sizeof(double), n_P, fp_stats);
                fwrite(dP_sub, sizeof(double), n_P, fp_stats);
                fwrite(&ln_Pr_min, sizeof(double), 1, fp_stats);
                fwrite(&ln_Pr_avg, sizeof(double), 1, fp_stats);
                fwrite(&ln_Pr_max, sizeof(double), 1, fp_stats);
                if(flag_autocor_on)
                    fwrite(auto_cor, sizeof(double), n_avg - 1, fp_stats);
                compute_MCMC_chain_stats(P_chain_stats,
                                         n_P,
                                         n_avg,
                                         P_min,
                                         P_avg,
                                         P_max,
                                         dP_avg,
                                         auto_cor,
                                         slopes,
                                         dP_sub,
                                         ln_Pr_chain,
                                         &ln_Pr_min,
                                         &ln_Pr_avg,
                                         &ln_Pr_max);
                fwrite(P_min, sizeof(double), n_P, fp_stats);
                fwrite(P_avg, sizeof(double), n_P, fp_stats);
                fwrite(P_max, sizeof(double), n_P, fp_stats);
                fwrite(dP_avg, sizeof(double), n_P, fp_stats);
                fwrite(dP_sub, sizeof(double), n_P, fp_stats);
                fwrite(&ln_Pr_min, sizeof(double), 1, fp_stats);
                fwrite(&ln_Pr_avg, sizeof(double), 1, fp_stats);
                fwrite(&ln_Pr_max, sizeof(double), 1, fp_stats);
                if(flag_autocor_on)
                    fwrite(auto_cor, sizeof(double), n_avg - 1, fp_stats);
                for(i_P = 0; i_P < n_P; i_P++)
                    drift[i_P] += slopes[i_P];
                fwrite(slopes, sizeof(double), n_P, fp_stats);
                fwrite(drift, sizeof(double), n_P, fp_stats);
            }

            // Report progress
            if(i_iteration == i_iteration_next_report) {
                i_report++;
                SID_log("%3d%% complete.", SID_LOG_COMMENT | SID_LOG_TIMER, (100 / n_report) * (i_report));
                time_diff /= (double)(n_avg);
                if(time_diff > 0.1)
                    SID_log("\tMean time for single model call: %.1f", SID_LOG_COMMENT, time_diff);
                i_iteration_next_report =
                    GBP_MIN(n_iterations_phase, i_iteration_start + (n_iterations_phase - i_iteration_start) * (i_report + 1) / n_report);
            }

            // Check to see if a stop has been called to the run ...
            if(SID.I_am_Master) {
                if((fp_stop = fopen(filename_stop, "r")) != NULL) {
                    fclose(fp_stop);
                    // remove(filename_stop);
                    flag_stop = GBP_TRUE;
                }
            }
            SID_Bcast(&flag_stop, 1, SID_INT, SID_MASTER_RANK, MCMC->comm);

            // ... if so, stop all ranks and cancel the subsequent analysis stage
            if(flag_stop) {
                flag_continue          = GBP_FALSE;
                MCMC->flag_analysis_on = GBP_FALSE;
                i_phase                = 3;
                SID_log("*** Stop file found.  Terminating run. ***", SID_LOG_COMMENT);
            }

            // Check to see if this phase's iterations are complete
            if(i_iteration >= n_iterations_phase)
                flag_continue = GBP_FALSE;

        } // while flag_continue=GBP_TRUE

        // Report progress
        SID_log("Completed iterations in this phase: %d (%d requested)", SID_LOG_COMMENT, i_iteration, n_iterations_phase);

        // Report success rate
        SID_log("Proposal success: %5.3f%% (%lld of %lld)",
                SID_LOG_COMMENT,
                1e2 * ((float)(MCMC->n_success) / (float)(MCMC->n_propositions)),
                MCMC->n_success,
                MCMC->n_propositions);

        i_iteration       = 0;
        i_iteration_start = i_iteration;
        SID_log("Done.", SID_LOG_CLOSE);
    }
    if(my_chain == SID.My_rank) {
        fp_chain_config = fopen(filename_chain_config, "wb");
        fwrite(&n_iterations_file_total, sizeof(int), 1, fp_chain_config);
        fwrite(&n_iterations_file_burn, sizeof(int), 1, fp_chain_config);
        fwrite(&(MCMC->temperature), sizeof(double), 1, fp_chain_config);
        fwrite(MCMC->V, sizeof(double), n_P * n_P, fp_chain_config);
        fclose(fp_chain_config);
        fclose(fp_chain);
        fclose(fp_stats);
    }
    SID_log("Done.", SID_LOG_CLOSE);

    // Perform analysis on chain(s)
    if(MCMC->flag_analysis_on)
        analyze_MCMC(MCMC);

    // Clean-up
    SID_log("Cleaning-up...", SID_LOG_OPEN);
    SID_free(SID_FARG P_min);
    SID_free(SID_FARG P_max);
    SID_free(SID_FARG P_avg);
    SID_free(SID_FARG dP_avg);
    SID_free(SID_FARG ln_Pr_chain);
    SID_free(SID_FARG ln_Pr_proposals);
    SID_free(SID_FARG drift);
    SID_free(SID_FARG slopes);
    SID_free(SID_FARG dP_sub);
    for(i_P = 0; i_P < n_P; i_P++) {
        SID_free(SID_FARG P_chain_stats[i_P]);
        SID_free(SID_FARG P_prop_stats[i_P]);
        if(flag_autocor_on)
            SID_free(SID_FARG auto_cor[i_P]);
    }
    SID_free(SID_FARG P_chain_stats);
    SID_free(SID_FARG P_prop_stats);
    if(flag_autocor_on)
        SID_free(SID_FARG auto_cor);
    SID_log("Done.", SID_LOG_CLOSE);

    SID_log("Done.", SID_LOG_CLOSE);
}

ext: c

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/*
main.h: contains all forward declarations and necessary includes
*/

#ifndef MAIN_H
#define MAIN_H

#define DEBUG false

#include <cstdlib>
#include <iostream>
#include <string>
#include <cmath>
#include <list>
#include <ctime>
#include <unistd.h>
#include <gsl/gsl_rng.h>

using namespace std;

extern gsl_rng *r; // global random number generator:

class Agent; // a thinking human agent
class Link; // a link between two agents

bool validLink(Agent* agent1, Agent* agent2, list<Link> &relation);
// returns false if the Agents are the same or already linked

void removeLink(Agent* agent1, Agent* agent2, list<Link> &relation);
// searches for and removes the link between two agents

void printXML(int nr_of_agents, list<Link> &relation);
// prints XML output

void printDOT(list<Link> &relation);
// prints GraphViz DOT output

double connectivity();
// TODO: returns some statistic indicating the fat-tailness of the connectivity distribution

double assortativity(list<Link> &relation);
// returns the assortativity coefficient of the network

double clustering(int nr_of_agents, Agent population[], list<Link> &relation);
// returns the global clustering coefficient

bool exists(Agent* target, list<Agent*> &agents);
// returns true if the target is in the list of agents

double avgpath(int nr_of_agents, Agent population[]);
// returns the average path length

#endif
// MAIN_H

ext: h


text:
/*
 * author: Achim Gaedke
 * created: January 2002
 * file: pygsl/src/multiminmodule.c
 * $Id: multiminmodule.c,v 1.5 2004/03/08 08:42:16 schnizer Exp $
 */
#include <setjmp.h>
#include <pygsl/block_helpers.h>
#include <pygsl/error_helpers.h>
#include <pygsl/function_helpers.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_multimin.h>
#include "multiminmodule_doc.h"
/*
 * I have two different types to minimize and I am too lazy to implement the same
 * twice. Therefore I use unions and add a flag to the minimizer for the type
 */

union pygsl_multimin_minimizer_type{
     const gsl_multimin_fminimizer_type   *f;
     const gsl_multimin_fdfminimizer_type *fdf;
};

union pygsl_multimin_minimizer{
     gsl_multimin_fminimizer   *f;
     gsl_multimin_fdfminimizer *fdf;
};

union pygsl_multimin_func{
     gsl_multimin_function     *f;
     gsl_multimin_function_fdf *fdf;
};

struct pygsl_solver_types{
     const char * f_minimizer;
     const char * fdf_minimizer;
};
static struct pygsl_solver_types my_solvers = {"F-Minimizer", "Fdf-Minimizer"};

/*
 *  type:
 *	 == 0 -> f
 *	 != 0 -> fdf
 */
typedef struct {
     PyObject_HEAD
     PyObject* py_f;
     PyObject* py_df;
     PyObject* py_fdf;
     PyObject* trailing_params;
     union pygsl_multimin_func      func;
     union pygsl_multimin_minimizer min;
     size_t n;
     const char *mytype;
     int isset; /* Used as a flag if the jmp_buf is set */
     jmp_buf buffer;
} PyGSL_multimin;

#define PyGSL_multimin_check(op)  ((op)->ob_type == &PyGSL_multimin_pytype)
#define PyGSL_multimin_isf(op)    ((op)->mytype  == my_solvers.f_minimizer)

static void
PyGSL_multimin_dealloc(PyGSL_multimin* self);
static PyObject*
PyGSL_multimin_getattr(PyGSL_multimin * obj, char *name);

PyTypeObject PyGSL_multimin_pytype = {
  PyObject_HEAD_INIT(NULL)	 /* fix up the type slot in initcrng */
  0,				 /* ob_size */
  "PyGSL_multimin",		 /* tp_name */
  sizeof(PyGSL_multimin),	 /* tp_basicsize */
  0,				 /* tp_itemsize */

  /* standard methods */
  (destructor)  PyGSL_multimin_dealloc, /* tp_dealloc  ref-count==0  */
  (printfunc)   0,		 /* tp_print    "print x"     */
  (getattrfunc) PyGSL_multimin_getattr,/* tp_getattr  "x.attr"      */
  (setattrfunc) 0,		 /* tp_setattr  "x.attr=v"    */
  (cmpfunc)     0,		   /* tp_compare  "x > y"       */
  (reprfunc)    0,                 /* tp_repr     `x`, print x  */

  /* type categories */
  0,				/* tp_as_number   +,-,*,/,%,&,>>,pow...*/
  0,				/* tp_as_sequence +,[i],[i:j],len, ...*/
  0,				/* tp_as_mapping  [key], len, ...*/

  /* more methods */
  (hashfunc)     0,		/* tp_hash    "dict[x]" */
  (ternaryfunc)  0,      /* tp_call    "x()"     */
  (reprfunc)     0,             /* tp_str     "str(x)"  */
  (getattrofunc) 0,		/* tp_getattro */
  (setattrofunc) 0,		/* tp_setattro */
  0,				/* tp_as_buffer */
  0L,				/* tp_flags */
  (char *) PyGSL_multimin_type_doc		/* tp_doc */
};

/* Reference to this module */
PyObject *module = NULL;
static const char filename[] = __FILE__;

/* The Callbacks */
double
PyGSL_multimin_function_f(const gsl_vector* x, void* params)
{
     double result;
     int flag;
     int i;

     FUNC_MESS_BEGIN();
     PyGSL_multimin *min_o;
     min_o = (PyGSL_multimin *) params;
     assert(PyGSL_multimin_check(min_o));
     for(i = 0; i<x->size; i++){
	  DEBUG_MESS(2, "Got a x[%d] of %f", i, gsl_vector_get(x, i));
     }
     flag = PyGSL_function_wrap_On_O(x, min_o->py_f, min_o->trailing_params, &result,
				     NULL, x->size, __FUNCTION__);
     if (flag!= GSL_SUCCESS){
	  result = gsl_nan();
	  if(min_o->isset == 1) longjmp(min_o->buffer,flag);
     }
     DEBUG_MESS(2, "Got a result of %f", result);
     FUNC_MESS_END();
     return result;
}

void
PyGSL_multimin_function_df(const gsl_vector* x, void* params, gsl_vector *df)
{
     int flag, i;
     PyGSL_multimin *min_o;
     min_o = (PyGSL_multimin *) params;

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(min_o));
     for(i = 0; i<x->size; i++){
	  DEBUG_MESS(2, "Got a x[%d] of %f", i, gsl_vector_get(x, i));
     }
     flag = PyGSL_function_wrap_Op_On(x, df, min_o->py_df, min_o->trailing_params,
				      x->size, x->size, __FUNCTION__);
     for(i = 0; i<df->size; i++){
	  DEBUG_MESS(2, "Got df x[%d] of %f", i, gsl_vector_get(df, i));
     }
     if(flag!=GSL_SUCCESS){
	  if(min_o->isset == 1) longjmp(min_o->buffer,flag);
     }
     FUNC_MESS_END();
}
void
PyGSL_multimin_function_fdf(const gsl_vector* x, void* params, double *f, gsl_vector *df)
{
     int flag, i;
     PyGSL_multimin *min_o;

     FUNC_MESS_BEGIN();
     min_o = (PyGSL_multimin *) params;
     assert(PyGSL_multimin_check(min_o));
     for(i = 0; i<x->size; i++){
	  DEBUG_MESS(2, "Got a x[%d] of %f", i, gsl_vector_get(x, i));
     }
     flag = PyGSL_function_wrap_On_O(x, min_o->py_fdf, min_o->trailing_params, f,
				     df, x->size, __FUNCTION__);
     DEBUG_MESS(2, "Got a result of %f", *f);
     for(i = 0; i<df->size; i++){
	  DEBUG_MESS(2, "Got df x[%d] of %f", i, gsl_vector_get(df, i));
     }
     if (flag!= GSL_SUCCESS){
	  *f = gsl_nan();
	  if(min_o->isset == 1) longjmp(min_o->buffer,flag);
     }
     FUNC_MESS_END();
     return;
}

static PyObject*
PyGSL_multimin_set_f(PyGSL_multimin *self, PyObject *args, PyObject *kw)
{

     PyObject* func = NULL, * params = NULL, * x = NULL, * steps = NULL;
     PyArrayObject * xa = NULL, * stepsa = NULL;
     int n=0, flag=GSL_EFAILED;
     int stride_recalc;
     gsl_vector_view gsl_x;
     gsl_vector_view gsl_steps;
     static const char *kwlist[] = {"f", "x0", "args", "steps", NULL};

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if (self->min.f == NULL) {
	  gsl_error("Got a NULL Pointer of min.f",  filename, __LINE__ - 3, GSL_EFAULT);
	  return NULL;
     }

     assert(args);
     /* arguments PyFunction, Parameters, start Vector, step Vector */
     if (0==PyArg_ParseTupleAndKeywords(args,kw,"OOOO", (char **) kwlist, &func,&x,&params,&steps))
	  return NULL;

     if(!PyCallable_Check(func)){
	  gsl_error("First argument must be callable",  filename, __LINE__ - 3, GSL_EBADFUNC);
	  return NULL;
     }
     n=self->n;
     xa  = PyGSL_PyArray_PREPARE_gsl_vector_view(x, PyArray_DOUBLE, 0, n, 4, NULL);
     if (xa == NULL){
	  PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__ - 1);
	  goto fail;
     }
     if(PyGSL_STRIDE_RECALC(xa->strides[0],sizeof(double), &stride_recalc) != GSL_SUCCESS)
	  goto fail;

     gsl_x = gsl_vector_view_array_with_stride((double *)(xa->data), stride_recalc, xa->dimensions[0]);
     stepsa =  PyGSL_PyArray_PREPARE_gsl_vector_view(steps, PyArray_DOUBLE, 0, n, 5, NULL);
     if (stepsa == NULL){
	  PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__ - 1);
	  goto fail;
     }
     if(PyGSL_STRIDE_RECALC(stepsa->strides[0],sizeof(double), &stride_recalc) != GSL_SUCCESS)
	  goto fail;

     gsl_steps = gsl_vector_view_array_with_stride((double *)(stepsa->data), stride_recalc, stepsa->dimensions[0]);

     if (self->func.f != NULL) {
	  /* free the previous function and params */
	  Py_XDECREF(self->trailing_params);
	  Py_XDECREF(self->py_f);
     } else {
	  /* allocate function space */
	  self->func.f=calloc(1, sizeof(gsl_multimin_function));
	  if (self->func.f==NULL) {
	       gsl_error("Could not allocate the object for the minimizer function",
			 filename, __LINE__ - 3, GSL_ENOMEM);
	       goto fail;
	  }
     }
     /* add new function  and parameters */
     self->trailing_params=params;
     Py_INCREF(params);
     self->py_f=func;
     Py_INCREF(func);

     /* initialize the function struct */
     self->func.f->n=n;
     self->func.f->f=PyGSL_multimin_function_f;
     self->func.f->params=(void*)self;

     if((flag = setjmp(self->buffer)) == 0){
	  self->isset = 1;
	  flag = gsl_multimin_fminimizer_set(self->min.f,self->func.f, &gsl_x.vector, &gsl_steps.vector);
	  if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS){
	       goto fail;
	  }
     } else {
	  goto fail;
     }
     Py_DECREF(xa);
     Py_DECREF(stepsa);

     Py_INCREF(Py_None);
     self->isset = 0;
     FUNC_MESS_END();
     return Py_None;

 fail:
     FUNC_MESS("Fail");
     PyGSL_ERROR_FLAG(flag);
     self->isset = 0;
     Py_XDECREF(xa);
     Py_XDECREF(stepsa);
     return NULL;

}
static PyObject*
PyGSL_multimin_set_fdf(PyGSL_multimin *self, PyObject *args, PyObject *kw)
{

     PyObject * f = NULL, * df = NULL, * fdf = NULL, * params = NULL,
	      * x = NULL;
     PyArrayObject * xa = NULL;
     int n=0, flag=GSL_EFAILED;
     int stride_recalc=-1;
     double step=0.01, tol=1e-4;
     gsl_vector_view gsl_x;
     static const char *kwlist[] = {"f", "df", "fdf",  "x0", "args", "step", "tol", NULL};

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if (self->min.fdf == NULL) {
	  gsl_error("Got a NULL Pointer of min.fdf", filename, __LINE__ - 3, GSL_EFAULT);
	       return NULL;
	  }

     /* arguments PyFunction, Parameters, start Vector, step Vector */
     if (0==PyArg_ParseTupleAndKeywords(args, kw, "OOOO|Odd", (char **)kwlist,
					&f, &df, &fdf, &x, &params, &step, &tol))
	  return NULL;

     n=self->n;
     xa  = PyGSL_PyArray_PREPARE_gsl_vector_view(x, PyArray_DOUBLE, 0, n, 4, NULL);
     if (xa == NULL){
	  PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__ - 1);
	  goto fail;
     }
     if(params == NULL){
          /* Reference counter increased later, when the parameters are set */
	  params = Py_None;
     }
     if(PyGSL_STRIDE_RECALC(xa->strides[0],sizeof(double), &stride_recalc) != GSL_SUCCESS)
	  goto fail;
     gsl_x = gsl_vector_view_array_with_stride((double *)(xa->data), stride_recalc, xa->dimensions[0]);

     if (self->func.fdf != NULL) {
	  /* free the previous function and params */
	  Py_XDECREF(self->trailing_params);
	  Py_XDECREF(self->py_f);
	  Py_XDECREF(self->py_df);
	  Py_XDECREF(self->py_fdf);
     } else {
	  /* allocate function space */
	  self->func.fdf=malloc(sizeof(gsl_multimin_function_fdf));
	  if (self->func.fdf==NULL) {
	       gsl_error("Could not allocate the object for the minimizer function",
			 filename, __LINE__ - 3, GSL_ENOMEM);
	       goto fail;
	  }
     }
     /* add new function  and parameters */
     self->trailing_params=params;     Py_INCREF(params);
     self->py_f=f;                     Py_INCREF(f);
     self->py_df=df;                   Py_INCREF(df);
     self->py_fdf=fdf;                 Py_INCREF(fdf);

     /* initialize the function struct */
     self->func.fdf->n=n;
     self->func.fdf->f  =PyGSL_multimin_function_f;
     self->func.fdf->df =PyGSL_multimin_function_df;
     self->func.fdf->fdf=PyGSL_multimin_function_fdf;
     self->func.fdf->params=(void*)self;

     if((flag = setjmp(self->buffer)) == 0){
	  self->isset = 1;
	  flag = gsl_multimin_fdfminimizer_set(self->min.fdf,self->func.fdf, &gsl_x.vector, step, tol);
	  if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS){
	       goto fail;
	  }
     }else{
	  goto fail;
     }
     self->isset = 0;
     Py_DECREF(xa);

     Py_INCREF(Py_None);
     FUNC_MESS_END();
     return Py_None;

 fail:
     PyGSL_ERROR_FLAG(flag);
     self->isset = 0;
     Py_XDECREF(xa);
     return NULL;

}
/*
static PyObject*
PyGSL_multimin_set(PyGSL_multimin *self, PyObject *args)
{

     if(PyGSL_multimin_isf(self)){
	  return PyGSL_multimin_set_f(self, args);
     }else {
	  return PyGSL_multimin_set_fdf(self, args);
     }
}
*/
static PyObject*
PyGSL_multimin_iterate(PyGSL_multimin *self, PyObject *args)
{
     int result, flag;

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if ( self->min.f ==NULL || self->func.f ==NULL) {
	  gsl_error("Got a NULL Pointer of min.f", filename, __LINE__ - 3, GSL_EFAULT);
	  return NULL;
     }

     if((flag = setjmp(self->buffer)) == 0){
	  self->isset = 1;
	  if(PyGSL_multimin_isf(self)){
	       result = gsl_multimin_fminimizer_iterate(self->min.f);
	  } else {
	       result = gsl_multimin_fdfminimizer_iterate(self->min.fdf);
	  }
     } else {
	  PyGSL_ERROR_FLAG(flag);
	  self->isset = 0;
	  return NULL;
     }
     self->isset = 0;
     FUNC_MESS_END();
     return PyGSL_error_flag_to_pyint(result);
}

static PyObject*
PyGSL_multimin_x(PyGSL_multimin *self, PyObject *args)
{
     gsl_vector* result;

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if (self->min.f==NULL || self->func.f==NULL) {
	  gsl_error("Got a NULL Pointer for min.f", filename, __LINE__ - 3, GSL_EFAULT);
	  return NULL;
     }
     if(PyGSL_multimin_isf(self)){
	  result=gsl_multimin_fminimizer_x(self->min.f);
     } else {
	  result=gsl_multimin_fdfminimizer_x(self->min.fdf);
     }
     if (result==NULL) {
	  gsl_error("How could that happen?",  filename, __LINE__ - 3, GSL_ESANITY);
	  return NULL;
     }
     FUNC_MESS_END();
     return (PyObject *) PyGSL_copy_gslvector_to_pyarray(result);
}

static PyObject*
PyGSL_multimin_minimum(PyGSL_multimin *self, PyObject *args)
{
     double min;

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if (self->min.f == NULL || self->func.f ==NULL) {
	  gsl_error("Got a NULL Pointer of min.f", filename, __LINE__ - 3, GSL_EFAULT);
	  return NULL;
     }
     if(PyGSL_multimin_isf(self)){
	  min = gsl_multimin_fminimizer_minimum(self->min.f);
     } else {
	  min = gsl_multimin_fdfminimizer_minimum(self->min.fdf);
     }
     FUNC_MESS_END();
     return PyFloat_FromDouble(min);
}

static PyObject*
PyGSL_multimin_name(PyGSL_multimin *self, PyObject *args)
{
     const char * name;
     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if(PyGSL_multimin_isf(self)){
	  name = gsl_multimin_fminimizer_name(self->min.f);
     } else {
	  name = gsl_multimin_fdfminimizer_name(self->min.fdf);
     }
     FUNC_MESS_END();
     return PyString_FromString(name);
}

static PyObject*
PyGSL_multimin_size(PyGSL_multimin *self, PyObject *args)
{
     double size;
     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if ( self->min.f ==NULL || self->func.f ==NULL) {
	  PyErr_SetString(PyExc_RuntimeError,"no function specified!");
	  return NULL;
     }
     if(PyGSL_multimin_isf(self)){
	  size = gsl_multimin_fminimizer_size(self->min.f);
     } else {
	  gsl_error("Can not calculate size for a FDF", filename, __LINE__, GSL_ESANITY);
	  return NULL;
     }
     FUNC_MESS_END();
     return PyFloat_FromDouble(size);
}

static PyObject*
PyGSL_multimin_test_size_method(PyGSL_multimin *self, PyObject *args)
{
     int flag=GSL_EFAILED;
     double epsabs;

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if ( self->min.f ==NULL || self->func.f ==NULL) {
	  PyErr_SetString(PyExc_RuntimeError,"no function specified!");
	  return NULL;
     }

     if(PyGSL_multimin_isf(self)){
	  if (0==PyArg_ParseTuple(args,"d", &epsabs))
	       return NULL;
	  flag = gsl_multimin_test_size(gsl_multimin_fminimizer_size(self->min.f), epsabs);
     } else {
	  gsl_error("Can not calculate size for a FDF", filename, __LINE__, GSL_ESANITY);
	  return NULL;
     }

     FUNC_MESS_END();
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);
}

static PyObject*
PyGSL_multimin_restart(PyGSL_multimin *self, PyObject *args)
{
     int flag;

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if(PyGSL_multimin_isf(self)){
	  gsl_error("Can not restart for a F type solver", filename, __LINE__, GSL_ESANITY);
	  return NULL;
     }
     flag = gsl_multimin_fdfminimizer_restart(self->min.fdf);

     if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS){
	  return NULL;
     }
     FUNC_MESS_END();
     Py_INCREF(Py_None);
     return Py_None;
}
static PyObject*
PyGSL_multimin_vec_fdf(PyGSL_multimin *self, PyObject *args,
		       gsl_vector *(*func)(gsl_multimin_fdfminimizer * s))
{
     gsl_vector *result;
     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if(PyGSL_multimin_isf(self)){
	  gsl_error("Can not retrieve this information for a F type solver!", filename, __LINE__, GSL_ESANITY);
	  return NULL;
     } else {
	  result=func(self->min.fdf);
     }
     FUNC_MESS_END();
     return (PyObject *) PyGSL_copy_gslvector_to_pyarray(result);
}

static PyObject*
PyGSL_multimin_istype(PyGSL_multimin *self, PyObject *args)
{
     static const char *f = "F-Minimizer",  *fdf = "Fdf-Minimizer";
     const char *p;

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if(PyGSL_multimin_isf(self)){
	  p = f;
     } else {
	  p = fdf;
     }
     FUNC_MESS_END();
     return PyString_FromString(p);
}

static PyObject*
PyGSL_multimin_dx(PyGSL_multimin *self, PyObject *args)
{
     return PyGSL_multimin_vec_fdf(self, args, gsl_multimin_fdfminimizer_dx);
}

static PyObject*
PyGSL_multimin_gradient(PyGSL_multimin *self, PyObject *args)
{
     return PyGSL_multimin_vec_fdf(self, args, gsl_multimin_fdfminimizer_gradient);
}

static PyObject*
PyGSL_multimin_test_gradient_method(PyGSL_multimin * self, PyObject *args)
{

     double epsabs;
     int flag;
     FUNC_MESS_BEGIN();

     assert(PyGSL_multimin_check(self));
     if (0==PyArg_ParseTuple(args,"d", &epsabs))
	  return NULL;

     if(PyGSL_multimin_isf(self)){
	  gsl_error("Can not retrieve this information for a F type solver!", filename, __LINE__, GSL_ESANITY);
	  return NULL;
     }
     flag = gsl_multimin_test_gradient(gsl_multimin_fdfminimizer_gradient(self->min.fdf), epsabs);
     FUNC_MESS_END();
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);

}

#define PyGSL_MULTIMIN_COMMON_METHODS  \
  {"iterate",     (PyCFunction)PyGSL_multimin_iterate,METH_NOARGS,(char *)multimin_iterate_doc,}, \
  {"x",           (PyCFunction)PyGSL_multimin_x,      METH_NOARGS,(char *)multimin_x_doc,      }, \
  {"minimum",     (PyCFunction)PyGSL_multimin_minimum,METH_NOARGS,(char *)multimin_minimum_doc,}, \
  {"name",        (PyCFunction)PyGSL_multimin_name,   METH_NOARGS,(char *)multimin_name_doc,   }, \
  {"type",        (PyCFunction)PyGSL_multimin_istype, METH_NOARGS,(char *)multimin_istype_doc, },

static PyMethodDef PyGSL_multimin_fmethods[] = {
     PyGSL_MULTIMIN_COMMON_METHODS
     {"set",      (PyCFunction)PyGSL_multimin_set_f,           METH_VARARGS|METH_KEYWORDS, (char *)multimin_set_f_doc    },
     {"size",     (PyCFunction)PyGSL_multimin_size,            METH_NOARGS,  (char *)multimin_size_doc     },
     {"test_size",(PyCFunction)PyGSL_multimin_test_size_method,METH_VARARGS, (char *)multimin_test_size_doc},
     {NULL, NULL, 0, NULL}           /* sentinel */
};

static PyMethodDef PyGSL_multimin_fdfmethods[] = {
     PyGSL_MULTIMIN_COMMON_METHODS
     {"set",      (PyCFunction)PyGSL_multimin_set_fdf,                 METH_VARARGS|METH_KEYWORDS, (char *)multimin_set_fdf_doc      },
     {"restart",  (PyCFunction)PyGSL_multimin_restart,                 METH_NOARGS,  (char *)multimin_restart_doc      },
     {"dx",       (PyCFunction)PyGSL_multimin_dx,                      METH_NOARGS,  (char *)multimin_dx_doc           },
     {"gradient", (PyCFunction)PyGSL_multimin_gradient,                METH_NOARGS,  (char *)multimin_gradient_doc     },
     {"test_gradient",(PyCFunction)PyGSL_multimin_test_gradient_method,METH_VARARGS, (char *)multimin_test_gradient_doc},
     {NULL, NULL, 0, NULL}           /* sentinel */
};

static PyObject*
PyGSL_multimin_init(PyObject *self, PyObject *args,
		    union pygsl_multimin_minimizer_type type, int t)
{

     PyGSL_multimin *min_o=NULL;
     size_t n;
     /* static const char functionname [] = __FUNCTION__; */

     FUNC_MESS_BEGIN();
     min_o =  (PyGSL_multimin *) PyObject_NEW(PyGSL_multimin, &PyGSL_multimin_pytype);
     if(min_o == NULL){
	  return NULL;
     }

     if (0==PyArg_ParseTuple(args,"l", &n))
	  return NULL;

     if (n<=0) {
	  PyErr_SetString(PyExc_RuntimeError, "dimension must be >0");
	  return NULL;
     }
     min_o->n=n;
     min_o->min.f=NULL;
     min_o->min.fdf=NULL;
     min_o->func.f=NULL;
     min_o->func.fdf=NULL;
     min_o->py_f=NULL;
     min_o->py_df=NULL;
     min_o->py_fdf=NULL;
     min_o->trailing_params=NULL;
     min_o->mytype=NULL;

     if(t == 0){
	  min_o->min.f = gsl_multimin_fminimizer_alloc(type.f,n);
	  if (min_o->min.f == NULL) {
		gsl_error("Could not allocate the object for the minimizer",
			  filename, __LINE__ - 3, GSL_ENOMEM);
	       goto fail;
	  }
	  min_o->mytype = my_solvers.f_minimizer;
	  assert(PyGSL_multimin_isf(min_o));
     } else {
	  min_o->min.fdf = gsl_multimin_fdfminimizer_alloc(type.fdf,n);
	  if (min_o->min.fdf == NULL) {
		gsl_error("Could not allocate the object for the fdfminimizer",
			  filename, __LINE__ - 3, GSL_ENOMEM);
	       goto fail;
	  }
	  min_o->mytype = my_solvers.fdf_minimizer;
	  assert(!PyGSL_multimin_isf(min_o));
     }
     FUNC_MESS_END();
     return (PyObject *) min_o;
 fail:
     Py_XDECREF(min_o);
     return NULL;
}

#define AMINIMIZER(name)                                                      \
static PyObject* PyGSL_multimin_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     union pygsl_multimin_minimizer_type type;                                \
     FUNC_MESS_BEGIN();                                                       \
     type.f  = gsl_multimin_fminimizer_ ## name;                              \
     tmp = PyGSL_multimin_init(self, args, type, 0);                          \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, (char *) filename, __FUNCTION__, __LINE__);      \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}
#define AMINIMIZER_FDF(name)                                                  \
static PyObject* PyGSL_multimin_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     union pygsl_multimin_minimizer_type type;                                \
     FUNC_MESS_BEGIN();                                                       \
     type.fdf  = gsl_multimin_fdfminimizer_ ## name;                          \
     tmp = PyGSL_multimin_init(self, args,  type, 1);                         \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, (char *) filename, __FUNCTION__, __LINE__);      \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}

AMINIMIZER(nmsimplex)
AMINIMIZER_FDF(steepest_descent)
AMINIMIZER_FDF(vector_bfgs)
AMINIMIZER_FDF(conjugate_pr)
AMINIMIZER_FDF(conjugate_fr)

static void
PyGSL_multimin_dealloc(PyGSL_multimin *self)
{
     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if(PyGSL_multimin_isf(self)){
	  if (self->min.f  != NULL)   gsl_multimin_fminimizer_free(self->min.f);
	  if (self->func.f != NULL)  free(self->func.f);
     }else{
	  if (self->min.fdf  != NULL)   gsl_multimin_fdfminimizer_free(self->min.fdf);
	  if (self->func.fdf != NULL)  free(self->func.fdf);
     }

     Py_XDECREF(self->trailing_params);
     Py_XDECREF(self->py_f);
     Py_XDECREF(self->py_df);
     Py_XDECREF(self->py_fdf);
     PyMem_Free(self);
     FUNC_MESS_END();
}

static PyObject*
PyGSL_multimin_getattr(PyGSL_multimin *self, char *name)
{

     PyObject *tmp = NULL;

     FUNC_MESS_BEGIN();
     assert(PyGSL_multimin_check(self));
     if(PyGSL_multimin_isf(self)){
	  tmp = Py_FindMethod(PyGSL_multimin_fmethods, (PyObject *) self, name);
     } else {
	  tmp = Py_FindMethod(PyGSL_multimin_fdfmethods, (PyObject *) self, name);
     }
     FUNC_MESS_END();
     return tmp;
}

static PyObject*
PyGSL_multimin_test_size(PyObject * self, PyObject *args)
{
     double size, epsabs;
     int flag = GSL_EFAILED;
     FUNC_MESS_BEGIN();
     if (0==PyArg_ParseTuple(args,"dd", &size, &epsabs))
	  return NULL;
     flag = gsl_multimin_test_size(size, epsabs);
     FUNC_MESS_END();
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);

}

static PyObject*
PyGSL_multimin_test_gradient(PyObject * self, PyObject *args)
{
     PyObject *g=NULL;
     PyArrayObject *ga=NULL;
     gsl_vector_view gradient;

     double epsabs;
     int flag = GSL_EFAILED, stride_recalc=-1;

     FUNC_MESS_BEGIN();
     if (0==PyArg_ParseTuple(args,"Od", &g, &epsabs))
	  return NULL;

     ga  = PyGSL_PyArray_PREPARE_gsl_vector_view(g, PyArray_DOUBLE, 0, -1, 1, NULL);
     if (ga == NULL){
	  PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__ - 1);
	  return NULL;
     }
     if((PyGSL_STRIDE_RECALC(ga->strides[0],sizeof(double), &stride_recalc)) != GSL_SUCCESS){
	  Py_XDECREF(ga);
	  return NULL;
     }
     gradient = gsl_vector_view_array_with_stride((double *)(ga->data), stride_recalc, ga->dimensions[0]);

     flag = gsl_multimin_test_gradient(&gradient.vector, epsabs);
     FUNC_MESS_END();
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);

}

static PyMethodDef multiminMethods[] = {
     {"nmsimplex",        PyGSL_multimin_init_nmsimplex,        METH_VARARGS, (char *)nmsimplex_doc       },
     {"steepest_descent", PyGSL_multimin_init_steepest_descent, METH_VARARGS, (char *)steepest_descent_doc},
     {"vector_bfgs",      PyGSL_multimin_init_vector_bfgs,      METH_VARARGS, (char *)vector_bfgs_doc     },
     {"conjugate_pr",     PyGSL_multimin_init_conjugate_pr,     METH_VARARGS, (char *)conjugate_pr_doc    },
     {"conjugate_fr",     PyGSL_multimin_init_conjugate_fr,     METH_VARARGS, (char *)conjugate_fr_doc    },
     {"test_size",        PyGSL_multimin_test_size,             METH_VARARGS, (char *)test_size_doc       },
     {"test_gradient",    PyGSL_multimin_test_gradient,         METH_VARARGS, (char *)test_gradient_doc   },
     {NULL, NULL, 0, NULL}        /* Sentinel */
};

void
initmultimin(void)
{
  PyObject* m, *dict, *item;
  m=Py_InitModule("multimin", multiminMethods);
  import_array();
  init_pygsl();
  /* init multimin type */
  PyGSL_multimin_pytype.ob_type  = &PyType_Type;

  module = m;

  Py_INCREF((PyObject*)&PyGSL_multimin_pytype);

  dict = PyModule_GetDict(m);
  if(!dict)
       goto fail;

  if (!(item = PyString_FromString((char*)PyGSL_multimin_module_doc))){
       PyErr_SetString(PyExc_ImportError,
		       "I could not generate module doc string!");
       goto fail;
  }
  if (PyDict_SetItemString(dict, "__doc__", item) != 0){
       PyErr_SetString(PyExc_ImportError,
		       "I could not init doc string!");
       goto fail;
  }

 fail:
  return;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
//
// Created by L. Jonathan Feldstein
//

#ifndef CIMPLE_CIMPLE_CONTROLLER_H
#define CIMPLE_CIMPLE_CONTROLLER_H

#include <stddef.h>
#include <math.h>
#include "cimple_system.h"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include "cimple_polytope_library.h"
#include <pthread.h>
#include "cimple_mpc_computation.h"

/**
 * @brief Action to get plant from current abstract state to target_abs_state.
 *
 * @param target target region the plant is supposed to reach
 * @param now current state of the plant
 * @param d_dyn discrete abstraction of the system
 * @param s_dyn system dynamics including auxiliary matrices
 * @param f_cost cost function to be minimized on the path
 */
void ACT(int target,
         current_state * now,
         discrete_dynamics * d_dyn,
         system_dynamics * s_dyn,
         cost_function * f_cost,
         double sec);
/**
 * @brief Apply the calculated control to the current state using system dynamics
 * @param x current state at time [0]
 * @param u matrix with next N inputs calculated by the MPC controller
 * @param A system dynamics
 * @param B input dynamics
 */
void apply_control(gsl_vector *x,
                   gsl_vector *u,
                   gsl_matrix *A,
                   gsl_matrix *B,
                   gsl_matrix *E,
                   gsl_vector *w,
                   size_t current_time);

/**
 * Fill a vector with gaussian distributed noise
 * @param w
 */
void simulate_disturbance(gsl_vector *w,
                          double mu,
                          double sigma);

void * main_computation(void *arg);

#endif //CIMPLE_CIMPLE_CONTROLLER_H

max_stars_repo_name: shaesaert/TuLiPXML

max_issues_repo_name: shaesaert/TuLiPXML

max_forks_repo_name: shaesaert/TuLiPXML


text:
#pragma once

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

gsl_rng *RandomNumberGenerator; /* Random number generator. */
unsigned RandomSeed; /* Random seed. */

unsigned GetRandomSeed()
{
  FILE *random_file;
  unsigned seed;
  int error;

  seed = 0;
  random_file = fopen("/dev/urandom", "rb");
  if (!random_file)
    return 0;
  do
  {
    error = fread(&seed, 1, sizeof(unsigned int), random_file);
    if (error != sizeof(unsigned int))
      return 0;
  }
  while (!seed);
  fclose(random_file);
  return seed;
}

/** InitialiseRandom.

    - Initialises the random number generator. **/

int InitialiseRandom(const gsl_rng_type *random_generator)
{
  RandomSeed = GetRandomSeed();
  if (!RandomSeed)
    return EXIT_FAILURE;
  RandomNumberGenerator = gsl_rng_alloc(random_generator);
  if (!RandomNumberGenerator)
    return EXIT_FAILURE;
  gsl_rng_set(RandomNumberGenerator, RandomSeed);
  return EXIT_SUCCESS;
}

/** FinaliseRandom.

    - Finalises the random number generator. **/

int FinaliseRandom()
{
  /* Release random number generator memory. */

  if (RandomNumberGenerator)
  {
    gsl_rng_free(RandomNumberGenerator);
    return EXIT_SUCCESS; // EXIT_SUCCESS = 0
  }
  return EXIT_FAILURE; // EXIT_FAILURE = 1
}

/////////////


text:
// -*- C++ -*-
//
// michael a.g. aïvázis
// orthologue
// (c) 1998-2021 all rights reserved
//

// get the headers
#include <cblas.h>

// smallest possible driver
int main() {
    // allocate the scala
    double alpha = 1;
    // allocate a vector
    double x[] = {0.0, 0.0, 0.0};
    double y[] = {0.0, 0.0, 0.0};

    // do it
    cblas_daxpy(3, alpha, x, 1, y, 1);

    // all done
    return 0;
}

// end of file

max_stars_repo_name: gmgunter/pyre

max_issues_repo_name: gmgunter/pyre

max_forks_repo_name: gmgunter/pyre


text:
#pragma once

#include <gsl.h>

#include <boost/optional.hpp>

#include <mitkDataStorage.h>
#include <mitkNodePredicateBase.h>

#include <mitkOperationActor.h>

#include "uk_ac_kcl_HierarchyManager_Export.h"

#define DECLARE_HIERARCHY_MANAGER_NODE_TYPE(nodeTypeName) static const crimson::HierarchyManager::NodeType& nodeTypeName();

#define DEFINE_HIERARCHY_MANAGER_NODE_TYPE(containerClass, nodeTypeName)                                                       \
    const crimson::HierarchyManager::NodeType& containerClass::nodeTypeName()                                                  \
    {                                                                                \
        static const crimson::HierarchyManager::NodeType type = crimson::HierarchyManager::getInstance()->createNodeType();    \
        return type;                                                                                \
    }

namespace crimson
{

/*! \brief   A singleton class for managing the hierarchical relations of CRIMSON data objects.
 *
 * The hierarchical relations define the dependencies between various data objects.
 * This mostly affects how the nodes are presented to the user in the MITK's Data Manager view.
 * For example, the vessel paths belonging to the same vessel tree must be grouped under the node
 * representing this vessel tree. Thus, deleting the vessel tree node will also delete the vessel
 * paths belonging to it.
 *
 * Hierarchical grouping also facilitates the discovery of th nodes using the mitk::DataStorage.
 * For example, to find all the contours defining a vessel, it is sufficient to iterate the
 * child nodes of the vessel path and select the nodes containing contours.
 *
 * HierarchyManager allows adding new node types and relationships between them. The most
 * convenient way to do so is to use the DECLARE_HIERARCHY_MANAGER_NODE_TYPE and
 * DEFINE_HIERARCHY_MANAGER_NODE_TYPE macro pair. This will correctly register the node types in
 * the hierarchy mananger.
 *
 * Example usage:
 *
 *     // MyNodeTypes.h
 *     #pragma once
 *     #include <HierarchyManager.h>
 *
 *     struct MyNodeTypes {
 *         DECLARE_HIERARCHY_MANAGER_NODE_TYPE(MyParentNode)
 *         DECLARE_HIERARCHY_MANAGER_NODE_TYPE(MyChildNode)
 *     };
 *
 *     // MyNodeTypes.cpp
 *     #include "MyNodeTypes.h"
 *
 *     DEFINE_HIERARCHY_MANAGER_NODE_TYPE(MyNodeTypes, MyParentNode)
 *     DEFINE_HIERARCHY_MANAGER_NODE_TYPE(MyNodeTypes, MyChildNode)
 *
 *     // UsingHierarchyManager.cpp
 *     void initialize()
 *     {
 *          auto hm = crimson::HierarchyManager::getInstance();
 *          hm->addNodeType(MyNodeTypes::MyParentNode(), ParentNodePredicate::New());
 *          hm->addNodeType(MyNodeTypes::MyChildNode(),  ChildNodePredicate::New());
 *          hm->addRelation(MyNodeTypes::MyParentNode(), MyNodeTypes::MyChildNode(), crimson::HierarchyManager::rtOneToMany);
 *     }
 *
 *     void addNode(mitk::DataNode* parent, mitkDataNode* child)
 *     {
 *          auto hm = crimson::HierarchyManager::getInstance();
 *          hm->addNodeToHierarchy(parent, MyNodeTypes::MyParentNode(), child, MyNodeTypes::MyChildNode());
 *          // ...
 *     }
 */

class HIERARCHYMANAGER_EXPORT HierarchyManager : public itk::Object, public mitk::OperationActor
{
public:
    /*! \name Singleton interface */
    ///@{
    static bool init();
    static HierarchyManager* getInstance();
    static void term();
    ///@}

    using NodeType = int;

    /*!
     * \brief   Creates a unique id for a new node type.
     */
    int createNodeType() { return _lastAssignedNodeType++; }

    /*! \brief   Values that represent relation types between parent and child nodes. */
    enum RelationType {
        rtOneToOne = 0, ///< Only one child of a particular type is allowed for the parent
        rtOneToMany, ///< Multiple children of a particular type are allowed for the parent
        rtUnknown   ///< No information about the relations between types
    };

    /*! \brief   Values that represent node type flags. */
    enum NodeTypeFlags {
        ntfNone = 0x00, ///< No flags
        ntfRecursiveDeletion = 0x01,    ///< Deleting a parent node also deletes all its children
        ntfUndoableDeletion = 0x02, ///< If the node is deleted, it should be put onto undo stack to allow deletion undo
        ntfPickable = 0x04  ///< The node should be allowed to be picked with mouse clicks in 3D rendering window
    };

    ///@{
    /*!
     * \brief   Adds a new node type.
     *
     * \param   type        The new node type identifier.
     * \param   predicate   The predicate defining the node.
     * \param   flags       The flags to be applied for nodes of this type.
     */
    bool addNodeType(const NodeType& type, const mitk::NodePredicateBase::Pointer& predicate, int flags = ntfNone);

    /*!
     * \brief   Adds a relation between the parent and child node types.
     *
     * \param   parentType  Type of the parent.
     * \param   childType   Type of the child.
     * \param   relation    The type of relation.
     */
    bool addRelation(const NodeType& parentType, const NodeType& childType, const RelationType& relation);

    /*!
     * \brief   Gets a relation between the parent and child node types.
     *
     * \param   parentNodeType  Type of the parent node.
     * \param   childNodeType   Type of the child node.
     */
    RelationType getRelation(const NodeType& parentNodeType, const NodeType& childNodeType) const;
    ///@}

    /*!
     * \brief   Gets data storage that the HierarchyManager operates upon.
     */
    mitk::DataStorage::Pointer getDataStorage() const;

    /*!
     * \brief   Set the flag whether to start a new undo group when adding nodes to hierarchy.
     *
     *  By default, addNodeToHierarchy() will start a new undo group meaning that undoing the
     *  addition will remove nodes one by one. However, when multiple nodes are to be added as a
     *  single unit, set this flag to false (and reset it to true after all nodes are added).
     */
    void setStartNewUndoGroup(bool start);

    /*!
     * \brief   Determine if a child node of a particular type can be added to the parent node of a different type.
     *
     * \param   parentNode          The desired parent node.
     * \param   parentNodeType      Type of the parent node.
     * \param   childNodeType       Type of the child node.
     * \param   allowReplacement    true if replacing a node is allowed (useful for rtOneToOne relationship).
     */
    bool canAddNode(mitk::DataNode* parentNode, const NodeType& parentNodeType, const NodeType& childNodeType,
                    bool allowReplacement) const;

    /*!
     * \brief   Adds a node to hierarchy.
     *
     * \param   parentNode      The parent node.
     * \param   parentNodeType  Type of the parent node.
     * \param   newNode         The new node.
     * \param   newNodeType     Type of the new node.
     *
     * \return  false if the relationship between the node types is unknown.
     */
    bool addNodeToHierarchy(const mitk::DataNode* parentNode, const NodeType& parentNodeType, const mitk::DataNode* newNode,
                            const NodeType& newNodeType);

    /*!
     * \brief   Searches for the first node type that the node satisfies.
     */
    boost::optional<NodeType> findFittingNodeType(mitk::DataNode* node) const;

    /*!
     * \brief   Searches for the potential parents of the node with a particular nodeType.
     */
    std::unordered_set<mitk::DataNode*> findPotentialParents(mitk::DataNode* node, NodeType nodeType) const;

    /*!
     * \brief   Change the parent of the node.
     */
    bool reparentNode(const mitk::DataNode::Pointer& node, const mitk::DataNode::Pointer& newParent);

    /*!
     * \brief   Gets a predicate associated with the node type.
     */
    mitk::NodePredicateBase::Pointer getPredicate(const NodeType& type) const;

    ///@{
    /*!
     * \brief   Gets an ancestor of the node of a particular type.
     *
     * \param   node        The child node.
     * \param   parentType  Type of the parent.
     * \param   directOnly  false to search parents of parents.
     */
    mitk::DataNode::Pointer getAncestor(const mitk::DataNode* node, const NodeType& parentType, bool directOnly = false) const;

    /*!
     * \brief   Gets the first descendant node of a particular type.
     *
     * \param   node            The parent node.
     * \param   descendantType  Type of the descendant.
     * \param   directOnly      false to search children of children.
     */
    mitk::DataNode::Pointer getFirstDescendant(const mitk::DataNode* node, const NodeType& descendantType,
                                               bool directOnly = false) const;

    /*!
     * \brief   Gets all the descendant nodes of a particular type.
     *
     * \param   node            The parent node.
     * \param   descendantsType Type of the descendants.
     * \param   directOnly      false to search children of children.
     */
    mitk::DataStorage::SetOfObjects::ConstPointer getDescendants(const mitk::DataNode* node, const NodeType& descendantsType,
                                                                 bool directOnly = false) const;

    ///@}

    static const char* nodeUIDPropertyName; ///< The name of the DataNode property where the nodeUID is stored

    /*!
     * \brief   Searches for the node with a particular UID.
     */
    mitk::DataNode::Pointer findNodeByUID(gsl::cstring_span<> nodeUID) const;

    /*!
     * \brief   Enables/disables the undo for adding and deleting nodes.
     */
    void enableUndo(bool enable);

    /*!
     * \brief   Overriden from mitk::Operation for undo/redo support.
     */
    void ExecuteOperation(mitk::Operation* operation) override;

private:
    HierarchyManager();
    ~HierarchyManager();

    HierarchyManager(const HierarchyManager&) = delete;
    HierarchyManager& operator=(const HierarchyManager&) = delete;

    void _connectDataStorageEvents();
    void _disconnectDataStorageEvents();

    void nodeAdded(const mitk::DataNode*);
    void nodeRemoved(const mitk::DataNode*);

private:
    static HierarchyManager* _instance;

    struct HierarchyManagerImpl;

    std::unique_ptr<HierarchyManagerImpl> _impl;

    int _lastAssignedNodeType = 0;

    bool _testFlags(const mitk::DataNode* node, NodeTypeFlags flags) const;
    std::vector<std::pair<std::vector<const mitk::DataNode*>, mitk::DataNode::ConstPointer>>
    _getUndoableRemoveNodes(const mitk::DataNode* node) const;
};

} // namespace crimson

max_stars_repo_name: carthurs/CRIMSONGUI

max_issues_repo_name: carthurs/CRIMSONGUI

max_forks_repo_name: carthurs/CRIMSONGUI


text:
/*
 * Copyright (c) 2015, Aleksas Mazeliauskas and Derek Teaney
 * All rights reserved.
 *
 * rnavier is distributed under MIT license;
 * see the LICENSE file that should be present in the root
 * of the source distribution, or alternately available at:
 * https://github.com/rnavier/rnavier/
 */
#ifndef RNAVIER_TBjFluctuation_h
#define RNAVIER_TBjFluctuation_h

#include <gsl/gsl_odeiv.h>
#include <memory>
#include "THModel3D.h"
#include "TIC3D.h"

class THydro3DBj ;
class TRNavier3DBj ;

class TBjFluctuation{

   private:

      THModel3D *fModel ;

      double fTau ;
      double fH ;
      double fKx;
      double fKy;
      double fKeta;

      static constexpr double fEpsAbs = 1.e-10;
      static constexpr double fEpsRel = 1.e-4  ;

      static const size_t fSize= 9 ;
      double fY[fSize] ;

      const gsl_odeiv_step_type *fType ;
      gsl_odeiv_step      *fStep ;
      gsl_odeiv_control   *fControl ;
      gsl_odeiv_evolve    *fEvolve ;
      gsl_odeiv_system    fSystem  ;

      void setY(const double &e0,
              const std::complex<double> &de,
              const std::complex<double> &gx,
              const std::complex<double> &gy,
              const std::complex<double> &gz, double y[]) ;

      void getY(const double y[], double &e0,
              std::complex<double> &de,
              std::complex<double> &gx,
              std::complex<double> &gy,
              std::complex<double> &gz) ;

      int step(double &t , const double &tf, double &h, double y[])  {
          return gsl_odeiv_evolve_apply(fEvolve, fControl, fStep,
                  &fSystem,  &t, tf, &h, y) ;
      }

      void eos_params(const double &e0,
             double &p0,
             double &cs0,
             double &s0,
             double &T0,
             double &eta0) ;

   public:

      TBjFluctuation(THModel3D *hmodel) ;
      ~TBjFluctuation() ;

      void setic(const double &t0,
              const double &kx, const double &ky, const double &keta,
              const double &e0,
              const std::complex<double> &de,
              const std::complex<double> &gx,
              const std::complex<double> &gy,
              const std::complex<double> &gz) ;

      void get_hydro_ic(const TICPoint3D &pnt, double &e, double &n, double *uI, double *piIJ_ns);

      int run(const double &t1, const bool &verbose = false) ;

      void get_timescales_cartesian(double &tauR, double &tau_sound, double &tau_damp) ;

      friend int tbjfluctuation_func (double t, const double y[], double f[],
         void *params) ;
} ;

int tbjfluctuation_func (double t, const double y[], double f[],
      void *params) ;

//////////////////////////////////////////////////////////////////////////////

class TICBjFluctuation : public TIC3D
{
   private:

      THModel3D *fModel ;

      double fTau0 ;
      double fLx;
      double fLy ;
      double fLeta;

      double fInitialS0;
      double fInitialE0;
      double fNKx;
      double fNKy;
      double fNKeta;
      double fDeOverE0;
      double fPhaseE;

      std::unique_ptr<TBjFluctuation> fFluctuation;

   public:

      TICBjFluctuation(TRNavier3DBj *rn, TGrid3D *gr) ;
      virtual void nextEvent(const double &bgen) override ;
      virtual void IC(const TICPoint3D &pnt, const TNDArray1D &auxi) override ;
      using TIC3D::fillGrid ;

      void evolve_to_time(const double &t)
      {
          fFluctuation->run(t) ;
      }
      void get_timescales_cartesian(double &tauR, double &tau_sound, double &tau_damp)
      {
          fFluctuation->get_timescales_cartesian(tauR, tau_sound, tau_damp) ;
      }

};

std::unique_ptr<TIC3D> make_ic_icbjfluctuation(THydro3DBj *hy, const std::string &icname) ;

#endif //RNAVIER_TBjFluctuation_h

max_stars_repo_name: rnavier/rnavier

max_issues_repo_name: rnavier/rnavier

max_forks_repo_name: rnavier/rnavier


text:
/**
 *
 * @file qwrapper_claswp.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:58 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_claswp(Quark *quark, Quark_Task_Flags *task_flags,
                       int n, PLASMA_Complex32_t *A, int lda,
                       int i1,  int i2, const int *ipiv, int inc)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_claswp_quark, task_flags,
        sizeof(int),                      &n,    VALUE,
        sizeof(PLASMA_Complex32_t)*lda*n,  A,        INOUT | LOCALITY,
        sizeof(int),                      &lda,  VALUE,
        sizeof(int),                      &i1,   VALUE,
        sizeof(int),                      &i2,   VALUE,
        sizeof(int)*n,                     ipiv,     INPUT,
        sizeof(int),                      &inc,  VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_claswp_quark = PCORE_claswp_quark
#define CORE_claswp_quark PCORE_claswp_quark
#endif
void CORE_claswp_quark(Quark *quark)
{
    int n, lda, i1, i2, inc;
    int *ipiv;
    PLASMA_Complex32_t *A;

    quark_unpack_args_7(quark, n, A, lda, i1, i2, ipiv, inc);
    LAPACKE_claswp_work(LAPACK_COL_MAJOR, n, A, lda, i1, i2, ipiv, inc );
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_claswp_f2(Quark *quark, Quark_Task_Flags *task_flags,
                          int n, PLASMA_Complex32_t *A, int lda,
                          int i1,  int i2, const int *ipiv, int inc,
                          PLASMA_Complex32_t *fake1, int szefake1, int flag1,
                          PLASMA_Complex32_t *fake2, int szefake2, int flag2)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_claswp_f2_quark, task_flags,
        sizeof(int),                        &n,     VALUE,
        sizeof(PLASMA_Complex32_t)*lda*n,    A,         INOUT | LOCALITY,
        sizeof(int),                        &lda,   VALUE,
        sizeof(int),                        &i1,    VALUE,
        sizeof(int),                        &i2,    VALUE,
        sizeof(int)*n,                       ipiv,      INPUT,
        sizeof(int),                        &inc,   VALUE,
        sizeof(PLASMA_Complex32_t)*szefake1, fake1,     flag1,
        sizeof(PLASMA_Complex32_t)*szefake2, fake2,     flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_claswp_f2_quark = PCORE_claswp_f2_quark
#define CORE_claswp_f2_quark PCORE_claswp_f2_quark
#endif
void CORE_claswp_f2_quark(Quark* quark)
{
    int n, lda, i1, i2, inc;
    int *ipiv;
    PLASMA_Complex32_t *A;
    void *fake1, *fake2;

    quark_unpack_args_9(quark, n, A, lda, i1, i2, ipiv, inc, fake1, fake2);
    LAPACKE_claswp_work(LAPACK_COL_MAJOR, n, A, lda, i1, i2, ipiv, inc );
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_claswp_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, PLASMA_Complex32_t *Aij,
                              int i1,  int i2, const int *ipiv, int inc, PLASMA_Complex32_t *fakepanel)
{
    DAG_CORE_LASWP;
    if (fakepanel == Aij) {
        QUARK_Insert_Task(
            quark, CORE_claswp_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(PLASMA_Complex32_t)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(PLASMA_Complex32_t)*1,      fakepanel,     SCRATCH,
            0);
    } else {
        QUARK_Insert_Task(
            quark, CORE_claswp_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(PLASMA_Complex32_t)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(PLASMA_Complex32_t)*1,      fakepanel,     INOUT,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_claswp_ontile_quark = PCORE_claswp_ontile_quark
#define CORE_claswp_ontile_quark PCORE_claswp_ontile_quark
#endif
void CORE_claswp_ontile_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    PLASMA_Complex32_t *A, *fake;
    PLASMA_desc descA;

    quark_unpack_args_7(quark, descA, A, i1, i2, ipiv, inc, fake);
    CORE_claswp_ontile(descA, i1, i2, ipiv, inc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_claswp_ontile_f2(Quark *quark, Quark_Task_Flags *task_flags,
                                 PLASMA_desc descA, PLASMA_Complex32_t *Aij,
                                 int i1,  int i2, const int *ipiv, int inc,
                                 PLASMA_Complex32_t *fake1, int szefake1, int flag1,
                                 PLASMA_Complex32_t *fake2, int szefake2, int flag2)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_claswp_ontile_f2_quark, task_flags,
        sizeof(PLASMA_desc),                &descA, VALUE,
        sizeof(PLASMA_Complex32_t)*1,        Aij,       INOUT | LOCALITY,
        sizeof(int),                        &i1,    VALUE,
        sizeof(int),                        &i2,    VALUE,
        sizeof(int)*(i2-i1+1)*abs(inc),      ipiv,      INPUT,
        sizeof(int),                        &inc,   VALUE,
        sizeof(PLASMA_Complex32_t)*szefake1, fake1, flag1,
        sizeof(PLASMA_Complex32_t)*szefake2, fake2, flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_claswp_ontile_f2_quark = PCORE_claswp_ontile_f2_quark
#define CORE_claswp_ontile_f2_quark PCORE_claswp_ontile_f2_quark
#endif
void CORE_claswp_ontile_f2_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    PLASMA_Complex32_t *A;
    PLASMA_desc descA;
    void *fake1, *fake2;

    quark_unpack_args_8(quark, descA, A, i1, i2, ipiv, inc, fake1, fake2);
    CORE_claswp_ontile(descA, i1, i2, ipiv, inc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_cswptr_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, PLASMA_Complex32_t *Aij,
                              int i1,  int i2, const int *ipiv, int inc,
                              const PLASMA_Complex32_t *Akk, int ldak)
{
    DAG_CORE_TRSM;
    QUARK_Insert_Task(
        quark, CORE_cswptr_ontile_quark, task_flags,
        sizeof(PLASMA_desc),              &descA, VALUE,
        sizeof(PLASMA_Complex32_t)*1,      Aij,       INOUT | LOCALITY,
        sizeof(int),                      &i1,    VALUE,
        sizeof(int),                      &i2,    VALUE,
        sizeof(int)*(i2-i1+1)*abs(inc),    ipiv,      INPUT,
        sizeof(int),                      &inc,   VALUE,
        sizeof(PLASMA_Complex32_t)*ldak,   Akk,       INPUT,
        sizeof(int),                      &ldak,  VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cswptr_ontile_quark = PCORE_cswptr_ontile_quark
#define CORE_cswptr_ontile_quark PCORE_cswptr_ontile_quark
#endif
void CORE_cswptr_ontile_quark(Quark *quark)
{
    int i1, i2, inc, ldak;
    int *ipiv;
    PLASMA_Complex32_t *A, *Akk;
    PLASMA_desc descA;

    quark_unpack_args_8(quark, descA, A, i1, i2, ipiv, inc, Akk, ldak);
    CORE_cswptr_ontile(descA, i1, i2, ipiv, inc, Akk, ldak);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_claswpc_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, PLASMA_Complex32_t *Aij,
                              int i1,  int i2, const int *ipiv, int inc, PLASMA_Complex32_t *fakepanel)
{
    DAG_CORE_LASWP;
    if (fakepanel == Aij) {
        QUARK_Insert_Task(
            quark, CORE_claswpc_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(PLASMA_Complex32_t)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(PLASMA_Complex32_t)*1,      fakepanel,     SCRATCH,
            0);
    } else {
        QUARK_Insert_Task(
            quark, CORE_claswpc_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(PLASMA_Complex32_t)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(PLASMA_Complex32_t)*1,      fakepanel,     INOUT,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_claswpc_ontile_quark = PCORE_claswpc_ontile_quark
#define CORE_claswpc_ontile_quark PCORE_claswpc_ontile_quark
#endif
void CORE_claswpc_ontile_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    PLASMA_Complex32_t *A, *fake;
    PLASMA_desc descA;

    quark_unpack_args_7(quark, descA, A, i1, i2, ipiv, inc, fake);
    CORE_claswpc_ontile(descA, i1, i2, ipiv, inc);
}