text:
/**
 *
 * @file qwrapper_clantr.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:57 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

void
CORE_clantr_quark(Quark *quark);
void
CORE_clantr_f1_quark(Quark *quark);

/***************************************************************************//**
 *
 **/
void QUARK_CORE_clantr(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum norm, PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
                       const PLASMA_Complex32_t *A, int LDA, int szeA,
                       int szeW, float *result)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;
    QUARK_Insert_Task(quark, CORE_clantr_quark, task_flags,
        sizeof(PLASMA_enum),                &norm,  VALUE,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(PLASMA_enum),                &diag,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(PLASMA_Complex32_t)*szeA,     A,             INPUT,
        sizeof(int),                        &LDA,   VALUE,
        sizeof(float)*szeW,                 NULL,          SCRATCH,
        sizeof(float),                      result,        OUTPUT,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_clantr_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum norm, PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
                          const PLASMA_Complex32_t *A, int LDA, int szeA,
                          int szeW, float *result,
                          float *fake, int szeF)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;

    if ( result == fake ) {
        QUARK_Insert_Task(quark, CORE_clantr_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(PLASMA_enum),                &diag,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(PLASMA_Complex32_t)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(float)*szeW,                 NULL,          SCRATCH,
            sizeof(float),                      result,        OUTPUT | GATHERV,
            0);
    } else {
        QUARK_Insert_Task(quark, CORE_clantr_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(PLASMA_enum),                &diag,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(PLASMA_Complex32_t)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(float)*szeW,                 NULL,          SCRATCH,
            sizeof(float),                      result,        OUTPUT,
            sizeof(float)*szeF,                 fake,          OUTPUT | GATHERV,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_clantr_quark = PCORE_clantr_quark
#define CORE_clantr_quark PCORE_clantr_quark
#endif
void CORE_clantr_quark(Quark *quark)
{
    float *normA;
    PLASMA_enum norm, uplo, diag;
    int M;
    int N;
    PLASMA_Complex32_t *A;
    int LDA;
    float *work;

    quark_unpack_args_9(quark, norm, uplo, diag, M, N, A, LDA, work, normA);
    *normA = LAPACKE_clantr_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        lapack_const(diag),
        M, N, A, LDA, work);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_clantr_f1_quark = PCORE_clantr_f1_quark
#define CORE_clantr_f1_quark PCORE_clantr_f1_quark
#endif
void CORE_clantr_f1_quark(Quark *quark)
{
    float *normA;
    PLASMA_enum norm, uplo, diag;
    int M;
    int N;
    PLASMA_Complex32_t *A;
    int LDA;
    float *work;
    float *fake;

    quark_unpack_args_10(quark, norm, uplo, diag, M, N, A, LDA, work, normA, fake);
    *normA = LAPACKE_clantr_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        lapack_const(diag),
        M, N, A, LDA, work);
}


text:
/**
 *
 * Copyright (c) 2017-2020  King Abdullah University of Science and Technology
 * All rights reserved.
 *
 * ExaGeoStat is a software package provided by KAUST
 **/
/**
 *
 * @file MLE_misc.h
 *
 * Header file of  auxiliary functions that are needed by ExaGeoStat.
 *
 * @version 1.1.0
 *
 * @author Sameh Abdulah
 * @date 2020-06-06
 *
 **/
#ifndef _MLE_MISC_H_
#define _MLE_MISC_H_
#include <stdbool.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include <sys/time.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>
#include "common.h"
#include "descriptor.h"
#include "chameleon/morse_runtime.h"
#include <nlopt.h>
#include <math.h>
#include <morse.h>
#include <lapacke.h>
#include <starpu.h>
#include <cblas.h>
#include "../../include/flops.h"
#include <starpu_profiling.h>
#if defined(CHAMELEON_USE_MPI)
#include <starpu_mpi.h>
#else
#endif
#if defined(CHAMELEON_USE_CUDA) && !defined(CHAMELEON_SIMULATION)
#include <starpu_scheduler.h>
#include <starpu_cuda.h>
#endif
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include "../../include/morse_starpu.h"
/** ****************************************************************************
 * PI value
 **/
#define PI (3.141592653589793)
/** ****************************************************************************
 * The radius of the  eartch (used by Great Circle Distance (GCD)
 **/
#define earthRadiusKm 6371.0
/** ****************************************************************************
 * Start timing macro
 **/
#define START_TIMING(_t)     _t =- cWtime();
/** ****************************************************************************
 * Stop timing macro
 **/
#define STOP_TIMING(_t)      _t += cWtime();
/*******************************************************************************/
/**
 *  Internal function to return address of block (m,n) with m,n = block indices
 */

inline static void *chameleon_getaddr_null(const MORSE_desc_t *A, int m, int n)
{
    (void)A; (void)m; (void)n;
    return NULL;
}

/**
 *  Internal function to return the leading dimension of element A(m,*) with m,n = block indices
 */
inline static int chameleon_getblkldd_ccrb(const MORSE_desc_t *A, int m)
{
    int mm = m + A->i / A->mb;
    return ( ((mm+1) == A->lmt) && ((A->lm % A->mb) != 0)) ? A->lm % A->mb : A->mb;
}

/**
 *  Internal function to return MPI rank of element A(m,n) with m,n = block indices
 */
inline static int chameleon_getrankof_2d(const MORSE_desc_t *desc, int m, int n)
{

    if(m>=n)
        return (m % desc->p) * desc->q + (n % desc->q);
    else
        return (n % desc->p) * desc->q + (m % desc->q);
}

/** ****************************************************************************
 * Allocate matrix in different modes
 **/
#define EXAGEOSTAT_ALLOCATE_MATRIX_TILE(_desc_, _memspace_, _type2_, _mb_, _nb_, _mbXnb_ , _lda_, _n_, _smb_, _snb_, _m_, _n2_, _p_, _q_) \
    if (data->ooc && _memspace_ == NULL && _mb_ != 1  && _nb_ !=1)                                                         \
MORSE_Desc_Create_OOC(_desc_, _type2_, _mb_, _nb_, _mbXnb_, _lda_, _n_, _smb_, _snb_, _m_, _n2_, \
        _p_, _q_);             \
else                                                            \
MORSE_Desc_Create(_desc_, _memspace_, _type2_, _mb_, _nb_, _mbXnb_, _lda_, _n_,_smb_, _snb_ , _m_, _n2_, \
        _p_, _q_); 	\

#define EXAGEOSTAT_ALLOCATE_FULL_MATRIX_TILE(_desc_, _memspace_, _type2_, _mb_, _nb_, _mbXnb_ , _lda_, _n_, _smb_, _snb_, _m_, _n2_, _p_, _q_) \
    if (data->ooc && _memspace_ == NULL && _mb_ != 1  && _nb_ !=1)                                                         \
MORSE_Desc_Create_OOC(_desc_, _type2_, _mb_, _nb_, _mbXnb_, _lda_, _n_, _smb_, _snb_, _m_, _n2_, \
        _p_, _q_);             \
else                                                            \
MORSE_Desc_Create_User(_desc_, _memspace_, _type2_, _mb_, _nb_, _mbXnb_, _lda_, _n_,_smb_, _snb_ , _m_, _n2_, \
        _p_, _q_, morse_getaddr_null, morse_getblkldd_ccrb, chameleon_getrankof_2d );	\

#define EXAGEOSTAT_ALLOCATE_DIAG_MATRIX_TILE(_desc_, _memspace_, _type2_, _mb_, _nb_, _mbXnb_ , _lda_, _n_, _smb_, _snb_, _m_, _n2_, _p_, _q_) \
    if (data->ooc && _memspace_ == NULL && _mb_ != 1  && _nb_ !=1)                                                         \
MORSE_Desc_Create_OOC(_desc_, _type2_, _mb_, _nb_, _mbXnb_, _lda_, _n_, _smb_, _snb_, _m_, _n2_, \
        _p_, _q_);             \
else                                                            \
MORSE_Desc_Create_User(_desc_, _memspace_, _type2_, _mb_, _nb_, _mbXnb_, _lda_, _n_,_smb_, _snb_ , _m_, _n2_, \
        _p_, _q_, morse_getaddr_null, morse_getblkldd_ccrb, morse_getrankof_2d_diag );   \

/** ****************************************************************************
 *  Structure for  identifies two dimensions struct with two vectors (x and y).
 **/
typedef struct {
    double *x;				///< Values in X dimension.
    double *y;				///< Values in Y dimension.
    double *z;                              ///< Values in Z dimension.
} location;

/** ****************************************************************************
 *  Structure for reading real datasets through STARS-H
 **/
typedef struct {
    double xy;                              ///< Locations (x, y)
    double z;                               ///< Measurements.
} sdata;

/** ****************************************************************************
 *  Structure for output results
 **/
typedef struct {
    int problem_size;
    char* computation;
    char* kernel;
    char* ds_type;
    char* precision;
    int z_sample;
    int dense_ts;
    int lr_ts;
    int lr_acc;
    int lr_maxrank;
    int ncores;
    int ngpus;
    int p;
    int q;
    int num_params;
    double *initial_theta;
    double *starting_theta;
    double *estimated_theta;
    double final_loglik;
    double time_per_iteration;
    double flops_per_iteration;
    double total_mle_time;
    double mse_pred1;
    double mse_pred2;
    double mse_pred;
    double total_pred_time;
    double total_pred_flops;
    double mloe;
    double mmom;
    char* mloe_exec;
    double total_mloe_mmom_time;
    double matrix_gen_mloe_mmom_time;
    double cho_fact_mloe_mmom_time;
    double loop_mloe_mmom_time;
    double total_mloe_mmom_flops;
} output;

/** ****************************************************************************
 *  Structure for  uniquely identifies three different parameters for accuracy
 *  measurments accuracyDenseAppDiff, normDenseAppDiff, and normA.
 **/
/*typedef struct{
  double accuracyDenseAppDiff; 		///< Accuracy different between dense format and approx format.
  double normDenseAppDiff;		///< Norm difference between dense format and approx format.
  double normA;   			///< normA.
  } acc_struct;
  */

/** ****************************************************************************
 *  Structure for  uniquely identifies different variables that are needed
 *  by  EXAGEOSTAT to  ease arguments pass.
 **/
typedef struct
{
    double variance;            ///< Variance parameter.
    double variance1;           ///< Variance1 parameter.
    double variance2;           ///< Variance2 parameter.
    char *computation;		    ///< Exact or approx computation.
    char *c_fun;	            ///< Matern or pow-exp kernels.
    int test;			        ///< Synthetic or real dataset execution.
    int async;			        ///< Running mode: synchronous or asynchronous.
    int iter_count;			    ///< Number of iterations to converge.
    location  l1;			    ///< 2D locations for the first dataset.
    location lmiss;			    ///< 2D locations for the missing data (prediction stage).
    location lobs;			    ///< 2D locations for the observed data (prediction stage).
    location lm;  		        ///< 2D locations for the median data point.
    void *descC;			    ///< Covariance matrix C descriptor.
    void *descsubC11;           ///< Covariance sub matrix C11 descriptor.
    void *descsubC12;           ///< Covariance sub matrix C12 descriptor.
    void *descsubC21;           ///< Covariance sub matrix C21 descriptor.
    void *descsubC22;           ///< Covariance sub matrix C22 descriptor.
    void *descZ;			    ///< Measurements Z descriptor.
    void *descZ1;               ///< Measurements Z1 submatrix descriptor.
    void *descZ2;               ///< Measurements Z2 submatrix descriptor.
    double *Adense;			    ///< Dense matrix descriptor in the case of approximation mode - accuracy check.
    double *Adense2;		    ///< Dense matrix descriptor2 in the case of approximation mode - accuracy check.
    void *descZcpy;			    ///< A copy of Measurements Z descriptor.
    void *descdet;			    ///< Determinant descriptor.
    void *descproduct;		    ///< Dot product descriptor.
    void *descproduct1;         ///< Dot product descriptor.
    void *descproduct2;         ///< Dot product descriptor.
    void *descZmiss;		    ///< Missing measurements descriptor.
    void *descC12;			    ///< Covariance Matrix C12 descriptor.
    void *descC22;			    ///< Covariance Matrix C22 descriptor.
    void *descZactual;		    ///< Actual Measurements Z descriptor.
    void *descZobs;			    ///< observed Measurements Z descriptor.
    void *descmse1;        		///< Mean Square Error (MSE) descriptor.
    void *descmse2;             ///< Mean Square Error (MSE) descriptor.
    void *descmse;              ///< Mean Square Error (MSE) descriptor.
    void *sequence;			    ///< MORSE sequence.
    void *request;			    ///< MORSE request.
    int verbose;			    ///< Verbose indicator.
    int check;			        ///< Check indicator -- approximation mode.
    int log;			        ///< Log files generation indicator, 0--> no, 1--> yes.
    double avg_exec_time_per_iter;	///< Avergae execution time per iteration (only used in verbose mode).
    double total_exec_time;		///< Total execution time (only used in verbose mode).
    double avg_flops_per_iter;	///< Avergae flops per iteration (only used in verbose mode).
    double final_loglik;		///< Final log likelihood value.
    char *locsFPath;		    ///< Locations file path -- in the case of real dataset (real mode).
    char *obsFPath;			    ///< Observations file path --  in the case of real dataset (real mode).
    char *obsFPath2;            ///< Observations file path2 (bivariate case) --  in the case of real dataset (real mode).
    char *actualZFPath;		    ///< Actual observations file path -- in the case of prediction.
    char *actualZFPath2;        ///< Actual observations file path -- in the case of prediction.
    char *actualZLocFPath;		///< Actial locations file path -- in the case of prediction.
    double det;			        ///< determinant value.
    double  dotp;			    ///< double dot product value.
    double  dotp1;              ///< double dot2 product value.
    double  dotp2;              ///< double dot3 product value.
    float sdotp;			    ///< single dot product value.
    double mserror;			    ///< Mean Square Error (MSE) value.
    double mserror1;            ///< Mean Square Error (MSE) value, variable 1 in case of bivariate.
    double mserror2;            ///< Mean Square Error (MSE) value, variable 2 in case of bivariate.
    char *dm;			        ///< Distance metric to be used ed->Euclidian Distance -- gcd->Great Circle Distance.
    int diag_thick;		        ///< The thick of used diagonal in the case of diagonal approximation approach.
    char *nFileLog;			    ///< log file name (only used if log -->1).
    FILE *pFileLog;			    ///< log file path (only used if log -->1).
    int hicma_maxrank;		    ///< Max Rank in the case of LR-HiCMA approx
    int hicma_data_type;        ///< To define the type of the problem to HiCMA (HICMA_STARSH_PROB_GEOSTAT (Synthetic) or HICMA_STARSH_PROB_GEOSTAT_POINT (real))
    void *hicma_descC;		    ///< HiCMA descC descriptor (for accuracy check).
    void *hicma_descZ;		    ///< HiCMA descZ descriptor.
    void *hicma_descCD;		    ///< HiCMA descCD descriptor.
    void *hicma_descCUV;		///< HiCMA descCUV descriptor.
    void *hicma_descCrk;		///< HiCMA descCrk descriptor.
    void *hicma_descZcpy;       ///< A copy of Measurements Z descriptor.
    void *hicma_descdet;        ///< Determinant descriptor.
    void *hicma_descproduct;    ///< Dot product descriptor.
    void *hicma_descC12D;		///< HiCMA descCD descriptor.
    void *hicma_descC12UV;		///< HiCMA descCUV descriptor.
    void *hicma_descC12rk;		///< HiCMA descCrk descriptor.
    void *hicma_descC22D;		///< HiCMA descCD descriptor.
    void *hicma_descC22UV;		///< HiCMA descCUV descriptor.
    void *hicma_descC22rk;		///< HiCMA descCrk descriptor.
    double hicma_acc;		    ///< Accuracy in the case of LR-HiCMA approx.
    void *hsequence;            ///< HiCMA sequence.
    void *hrequest;             ///< HiCMA request.
    int opt_tol; 	            ///< The parameter tol is a tolerance that is used for the purpose of stopping criteria only.
    int opt_max_iters;	        ///< Maximum number of mle iterations.
    int ooc;                    ///< Support Out-Of-Core execution, 0-->no, 1-->yes.
    char* kernel_fun;           ///< stationary_matern, or non_stationary_matern.
    int precision;              ///< (0)Double, (1)Single, and (2)Mixed.
    //Mixed Precision
    void *desctemp;             ///< Temporary descriptor for mixed precision Cholesky factorization.
    void *desctemp22;           ///< Temporary descriptor for mixed precision Cholesky factorization.
    //MLOE and MMOM
    void *desck_t;
    void *desck_a;
    void *desck_ttmp;
    void *desck_atmp;
    void *descK_ttmp;
    void *descK_t;
    void *descK_a;
    void *descexpr1;
    void *descexpr2;
    void *descexpr3;
    void *descexpr4;
    void *descestimatedalpha;
    void *desctruthalpha;
    void *desc_mloe_mmom;
    double expr1;
    double expr2;
    double expr3;
    double expr4;
    double mloe;
    double mmom;
    int mloe_mmom;
    int mloe_mmom_async;
    int mspe;
    char* recovery_file;
    char* checkpoint_file;
    int time_slots;
    int idw;
} MLE_data;

/** ****************************************************************************
 * Verbose Macro.
 **/
#define VERBOSE(str)	\
    if (data->verbose == 1 && MORSE_My_Mpi_Rank() == 0){	\
        fprintf(stdout, "%s", str);	\
        fflush(stdout);\
    }

/** ****************************************************************************
 *  Success Macro.
 **/
#define SUCCESS(success, str) \
    if (success != MORSE_SUCCESS){ \
        fprintf(stdout, "%s", str);\
        fflush(stdout);\
        exit(EXIT_FAILURE);\
    }

output results;

void pick_random_points(MLE_data *data, double *Zobs, double *Zactual,
        int nZmiss, int nZobs, int N);

void generate_interior_points(MLE_data *data, double *Zobs, double *Zactual,
        int nZmiss, int nZobs, int N);

void split_data(MLE_data * data, location *locations, double *Z,
        double *Zactual,
        int *N, int nZmiss);

void init_data_values(MLE_data *data);

int locations_obs_zsort_inplace(int n, location *locations, double *z);

void zsort_locations_obs(int n, location *locations, double *z);

double uniform_distribution(double rangeLow, double rangeHigh);

void print_dmatrix(char* desc, int m, int n,
        double* a, int lda);

void print_diagonal(char* desc, int m, double* a,
        int lda);

void print_smatrix(char* desc, int m, int n,
        float* a, int lda);

location* GenerateXYLoc(int n, int seed);

int countlines(char *filename);

void checkpointing(char *path, int iter_count, double* theta,
        double loglik, int num_params);

bool recover(char *path, int iter_count, double* theta,
        double* loglik, int num_params);

void write_to_file(char * path, int matrix_size,int ncores,
        int tile_size, int test, char *computation,
        int async, char *obsFPath,double total_exec_time,
        double avg_exec_time_per_iter, double avg_flops_per_iter,
        int p_grid, int q_grid, double final_loglik, int n);

void theta_parser2(double *theta_vec, char * kern, int num_params);

void write_vectors(double * zvec, MLE_data * data, int n);

void write_to_thetafile(char * path, double *theta, int num_params,
        int n, double time_per_iter, int total_no_iters,
        double prediction_error, double mloe,  double mmom);
//void readObsFile(char *obsfile, int n, double * streamdata);

void shuffle(double *array, location* locations, size_t n);

void theta_parser(double *initial_theta, double *target_theta,
        double *starting_theta, char *ikernel, char *kernel,
        double *lb, double *up, int test, int num_params);

void init_optimizer(nlopt_opt* opt, double *lb, double *up,
        double tol);

void print_summary(int test, int N, int ncores, int gpus,
        int ts, int lts, char *computation,
        int zvecs, int p_grid, int q_grid, int precision);

int print_result(MLE_data *data, double *starting_theta, int N,
        int zvecs, int ncores, int ts, int test,
        double *initial_theta, char *computation,
        int p_grid, int q_grid, double final_loglik,
        double prediction_error);

double cWtime(void);

void readlocfile(char* loc_file, int n,  location* l1);

void write_prediction_result(char *path, int matrix_size, int no_missing,
        double MSE1, double MSE2, double MSE,
        double solve_time, double flops);

int doesFileExist(const char *filename);

void init_log (MLE_data * data);

void finalize_log (MLE_data * data);
//acc_struct check_acc(MLE_data * HICMA_data, int n, int ts);
double core_matern_vector (double x0, double y0, double x1,
        double y1, double *localtheta, int distance_metric);

void pick_random_points2(MLE_data *data, double *Zobs, double *Zactual,
        int nZmiss, int nZobs, int N);

location* GenerateXYLoc_ST(int n, int t_slots, int seed);
void pick_random_points_noshuffle(MLE_data *data, double *Zobs, double *Zactual,
        int nZmiss, int nZobs, int N);

double* pred_idw(MLE_data *data, double *z_miss, double *z_actual,
        double*z_obs, int nZmiss, int nZobs);

void write_to_estimatedtheta(char * path, double *theta, int num_params,
        int n, double prediction_time,
        double mloe_mmom_time,
        double prediction_error1,
        double prediction_error2,
        double prediction_error,
        double mloe,
        double mmom,
        int zvecs);
#endif

ext: h


text:
/*	$Id$ */
/*
 * Copyright (c) 2014, 2015 Kristaps Dzonsons <kristaps@kcons.eu>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <assert.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#ifdef MAC_INTEGRATION
#include <gtkosxapplication.h>
#endif
#include <gtk/gtk.h>
#include <gdk/gdk.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_histogram.h>
#include <kplot.h>

#include "extern.h"

static void
swin_init(struct swin *c, enum view view, GtkBuilder *b)
{

	c->window = win_init_window(b, "window1");
	c->menu = win_init_menubar(b, "menubar1");
	c->draw = win_init_draw(b, "drawingarea1");
	c->boxconfig = win_init_box(b, "box3");
	c->menufile = win_init_menuitem(b, "menuitem1");
	c->menuview = win_init_menuitem(b, "menuitem2");
	c->menutools = win_init_menuitem(b, "menuitem3");
	c->viewclone = win_init_menuitem(b, "menuitem15");
	c->viewpause = win_init_menuitem(b, "menuitem20");
	c->viewunpause = win_init_menuitem(b, "menuitem21");
	c->views[VIEW_ISLANDMEAN] = win_init_menucheck(b, "menuitem45");
	c->views[VIEW_ISLANDERMEAN] = win_init_menucheck(b, "menuitem46");
	c->views[VIEW_MEAN] = win_init_menucheck(b, "menuitem8");
	c->views[VIEW_SMEAN] = win_init_menucheck(b, "menuitem37");
	c->views[VIEW_SEXTM] = win_init_menucheck(b, "menuitem43");
	c->views[VIEW_SEXTI] = win_init_menucheck(b, "menuitem44");
	c->views[VIEW_EXTM] = win_init_menucheck(b, "menuitem25");
	c->views[VIEW_TIMESPDF] = win_init_menucheck(b, "menuitem38");
	c->views[VIEW_TIMESCDF] = win_init_menucheck(b, "menuitem39");
	c->views[VIEW_EXTMMAXPDF] = win_init_menucheck(b, "menuitem28");
	c->views[VIEW_EXTMMAXCDF] = win_init_menucheck(b, "menuitem29");
	c->views[VIEW_EXTI] = win_init_menucheck(b, "menuitem26");
	c->views[VIEW_EXTIMINPDF] = win_init_menucheck(b, "menuitem27");
	c->views[VIEW_EXTIMINCDF] = win_init_menucheck(b, "menuitem30");
	c->views[VIEW_EXTIMINS] = win_init_menucheck(b, "menuitem35");
	c->views[VIEW_DEV] = win_init_menucheck(b, "menuitem6");
	c->views[VIEW_POLY] = win_init_menucheck(b, "menuitem7");
	c->views[VIEW_POLYMINPDF] = win_init_menucheck(b, "menuitem9");
	c->views[VIEW_POLYMINCDF] = win_init_menucheck(b, "menuitem11");
	c->views[VIEW_MEANMINPDF] = win_init_menucheck(b, "menuitem10");
	c->views[VIEW_MEANMINCDF] = win_init_menucheck(b, "menuitem12");
	c->views[VIEW_MEANMINQ] = win_init_menucheck(b, "menuitem13");
	c->views[VIEW_MEANMINS] = win_init_menucheck(b, "menuitem22");
	c->views[VIEW_POLYMINS] = win_init_menucheck(b, "menuitem31");
	c->views[VIEW_EXTMMAXS] = win_init_menucheck(b, "menuitem33");
	c->views[VIEW_POLYMINQ] = win_init_menucheck(b, "menuitem14");
	c->menuquit = win_init_menuitem(b, "menuitem5");
	c->menuautoexport = win_init_menuitem(b, "menuitem49");
	c->menuunautoexport = win_init_menuitem(b, "menuitem50");
	c->menuclose = win_init_menuitem(b, "menuitem24");
	c->menusave = win_init_menuitem(b, "menuitem34");
	c->menusaveall = win_init_menuitem(b, "menuitem47");

	gtk_widget_show_all(GTK_WIDGET(c->window));

	gtk_window_set_title(GTK_WINDOW(c->window),
		gtk_menu_item_get_label
		(GTK_MENU_ITEM(c->views[view])));

	gtk_widget_hide(GTK_WIDGET(c->menuunautoexport));

#ifdef MAC_INTEGRATION
	gtk_widget_hide(GTK_WIDGET(c->menu));
	gtk_widget_hide(GTK_WIDGET(c->menuquit));
	gtkosx_application_set_menu_bar
		(gtkosx_application_get(),
		 GTK_MENU_SHELL(c->menu));
	gtkosx_application_sync_menubar
		(gtkosx_application_get());
#endif
}

/*
 * Dereference a simulation.
 * This means that we've closed an output window that's looking at a
 * particular simulation.
 * When the simulation has no more references (i.e., no more windows are
 * painting that simulation), then the simulation destroys itself.
 */
static void
on_sim_deref(gpointer dat)
{
	struct sim	*sim = dat;

	g_debug("%p: Simulation deref (now %zu)", sim, sim->refs - 1);
	if (0 != --sim->refs)
		return;
	g_debug("%p: Simulation terminating", sim);
	sim_stop(sim, NULL);
}

/*
 * Dereference all simulations owned by a given window.
 * This happens when a window is closing.
 */
static void
on_sims_deref(gpointer dat)
{

	g_list_free_full(dat, on_sim_deref);
}

static void
curwin_free(gpointer dat)
{
	struct curwin	*cur = dat;
	size_t		 i;

	g_debug("%p: Simwin freeing", cur);
	kdata_destroy(cur->winmean);
	kdata_destroy(cur->winstddev);
	kdata_destroy(cur->winfitmean);
	kdata_destroy(cur->winfitstddev);
	kdata_destroy(cur->winmextinctmean);
	kdata_destroy(cur->winmextinctstddev);
	kdata_destroy(cur->winiextinctmean);
	kdata_destroy(cur->winiextinctstddev);
	for (i = 0; i < VIEW__MAX; i++)
		kplot_free(cur->views[i]);
	cur->b->windows = g_list_remove(cur->b->windows, cur);
	on_sims_deref(cur->sims);
	g_free(cur->autosave);
	g_free(cur);
}

static void
window_add_configmarkup(GtkWidget *box, const gchar *fmt, ...)
{
	gchar		 buf[1024];
	GtkWidget	*w;
	va_list		 ap;

	va_start(ap, fmt);
	g_vsnprintf(buf, sizeof(buf), fmt, ap);
	va_end(ap);

	w = gtk_label_new(NULL);
	gtk_label_set_markup(GTK_LABEL(w), buf);
	gtk_misc_set_alignment(GTK_MISC(w), 0.0, 0.5);
	gtk_container_add(GTK_CONTAINER(box), w);
}

static void
window_add_config(GtkWidget *box, const gchar *fmt, ...)
{
	gchar		 buf[1024];
	GtkWidget	*w;
	va_list		 ap;

	va_start(ap, fmt);
	g_vsnprintf(buf, sizeof(buf), fmt, ap);
	va_end(ap);

	w = gtk_label_new(buf);
	gtk_misc_set_alignment(GTK_MISC(w), 0.0, 0.5);
	gtk_container_add(GTK_CONTAINER(box), w);
}

static void
window_add_sim(struct curwin *cur, struct sim *sim)
{
	GtkWidget	*box, *outbox, *leftbox;
	struct kdata	*stats[2];
	enum kplottype	 ts[2];
	struct kdatacfg	 solidcfg, transcfg;
	struct kdatacfg	*cfgs[2];
	double		 solid[4], trans[4];
	gchar		 label[64];
	struct ksmthcfg	 smth;

	/* Get our colours. */

	memcpy(solid, cur->b->clrs[sim->colour %
		cur->b->clrsz].rgba, sizeof(solid));
	memcpy(trans, solid, sizeof(trans));
	trans[3] = 0.4;

	/* Append to the per-window simulation views. */
	kdata_vector_append(cur->winmean,
		g_list_length(cur->sims), 0.0);
	kdata_vector_append(cur->winstddev,
		g_list_length(cur->sims), 0.0);
	kdata_vector_append(cur->winfitmean,
		g_list_length(cur->sims), 0.0);
	kdata_vector_append(cur->winfitstddev,
		g_list_length(cur->sims), 0.0);
	kdata_vector_append(cur->winmextinctmean,
		g_list_length(cur->sims), 0.0);
	kdata_vector_append(cur->winmextinctstddev,
		g_list_length(cur->sims), 0.0);
	kdata_vector_append(cur->winiextinctmean,
		g_list_length(cur->sims), 0.0);
	kdata_vector_append(cur->winiextinctstddev,
		g_list_length(cur->sims), 0.0);

	/* Append our configuration. */
	outbox = gtk_box_new(GTK_ORIENTATION_HORIZONTAL, 2);

	g_snprintf(label, sizeof(label),
		""
		"    "
		"",
		(unsigned int)(solid[0] * 255),
		(unsigned int)(solid[1] * 255),
		(unsigned int)(solid[2] * 255));
	leftbox = gtk_label_new(NULL);
	gtk_misc_set_alignment(GTK_MISC(leftbox), 0.0, 0.0);
	gtk_label_set_markup(GTK_LABEL(leftbox), label);
	gtk_container_add(GTK_CONTAINER(outbox), leftbox);
	box = gtk_box_new(GTK_ORIENTATION_VERTICAL, 2);
	gtk_container_add(GTK_CONTAINER(outbox), box);
	window_add_config(box, "Name: %s", sim->name);
	window_add_configmarkup(box, "Payoffs: π = %s; "
		"T = %zu", sim->func, sim->stop);
	window_add_configmarkup(box, "Poisson offspring: "
		"λ = %g(1 + %g * π)",
		sim->alpha, sim->delta);
	window_add_config(box, "Incumbents: "
		"x = [%g, %g), %zu slices",
		sim->xmin, sim->xmax, sim->dims);
	if (MUTANTS_DISCRETE == sim->mutants)
		window_add_config(box, "Mutants: y = [%g, %g), "
			"%zu slices", sim->ymin, sim->ymax, sim->dims);
	else
		window_add_configmarkup(box, "Mutants: "
			"Gaussian σ = %g, Y = [%g, %g)",
			sim->mutantsigma, sim->ymin, sim->ymax);

	switch (sim->input) {
	case (INPUT_UNIFORM):
		window_add_config(box, "Population: uniform %zu "
			"islands, %zu islanders (%zu total), m = %g",
			sim->pop, sim->islands, sim->totalpop, sim->m);
		break;
	case (INPUT_VARIABLE):
		window_add_config(box, "Population: variable %zu "
			"islands (%zu total, %suniform), m = %g",
			sim->islands, sim->totalpop,
			NULL != sim->ms ? "non-" : "", sim->m);
		if (sim->ideathmean)
			window_add_config(box,
				"Island death Poisson mean: %zu, "
				"coefficient %g",
				sim->ideathmean, sim->ideathcoef);
		break;
	case (INPUT_MAPPED):
		window_add_config(box, "Population: mapped %zu "
			"islands (%zu total, %suniform), m = %g",
			sim->islands, sim->totalpop,
			NULL != sim->ms ? "non-" : "", sim->m);
		switch (sim->maptop) {
		case (MAPTOP_RECORD):
			window_add_config(box, "Map topology: "
				"KML records");
			break;
		case (MAPTOP_RAND):
			window_add_config(box, "Map topology: "
				"random records");
			break;
		case (MAPTOP_TORUS):
			window_add_config(box, "Map topology: "
				"toroidal records");
			break;
		default:
			abort();
		}
		switch (sim->migrant) {
		case (MAPMIGRANT_UNIFORM):
			window_add_config(box, "Map migration: "
				"uniform random");
			break;
		case (MAPMIGRANT_DISTANCE):
			window_add_config(box, "Map migration: "
				"inverse square distance");
			break;
		case (MAPMIGRANT_NEAREST):
			window_add_config(box, "Map migration: "
				"first nearest neighbour");
			break;
		case (MAPMIGRANT_TWONEAREST):
			window_add_config(box, "Map migration: "
				"two-nearest neighbours");
			break;
		default:
			abort();
		}
		break;
	default:
		abort();
	}

	if (0 == sim->fitpoly)
		window_add_config(box, "Polynomial fitting: disabled");
	else
		window_add_config(box, "Polynomial fitting: order "
			"%zu (%sweighted)", sim->fitpoly,
			sim->weighted ?  "" : "un");

	gtk_container_add(GTK_CONTAINER(cur->wins.boxconfig), outbox);
	gtk_widget_show_all(GTK_WIDGET(cur->wins.boxconfig));

	/* Configure our line and point style. */
	ksmthcfg_defaults(&smth);
	smth.movsamples = sim->smoothing;
	kdatacfg_defaults(&solidcfg);
	solidcfg.point.clr.type = KPLOTCTYPE_RGBA;
	memcpy(solidcfg.point.clr.rgba, solid, sizeof(solid));
	solidcfg.line.clr.type = KPLOTCTYPE_RGBA;
	memcpy(solidcfg.line.clr.rgba, solid, sizeof(solid));

	kdatacfg_defaults(&transcfg);
	transcfg.point.clr.type = KPLOTCTYPE_RGBA;
	memcpy(transcfg.point.clr.rgba, trans, sizeof(trans));
	transcfg.line.clr.type = KPLOTCTYPE_RGBA;
	memcpy(transcfg.line.clr.rgba, trans, sizeof(trans));

	/* Mean view. */
	kplot_attach_data(cur->views[VIEW_MEAN],
		sim->bufs.means->cold, KPLOT_LINES, &solidcfg);

	/* Mutant mean view. */
	kplot_attach_data(cur->views[VIEW_EXTM],
		sim->bufs.mextinct->cold, KPLOT_LINES, &solidcfg);

	/* Incumbent mean view. */
	kplot_attach_data(cur->views[VIEW_EXTI],
		sim->bufs.iextinct->cold, KPLOT_LINES, &solidcfg);

	/* Mean and poly-fitted line. */
	kplot_attach_data(cur->views[VIEW_POLY],
		sim->bufs.fitpolybuf, KPLOT_LINES, &solidcfg);
	kplot_attach_data(cur->views[VIEW_POLY],
		sim->bufs.means->cold, KPLOT_LINES, &transcfg);

	/* Mutant mean and smoothed line. */
	kplot_attach_smooth(cur->views[VIEW_SEXTM],
		sim->bufs.mextinct->cold, KPLOT_LINES, &solidcfg,
		KSMOOTH_MOVAVG, &smth);
	kplot_attach_data(cur->views[VIEW_SEXTM],
		sim->bufs.mextinct->cold, KPLOT_LINES, &transcfg);

	/* Mutant mean and smoothed line. */
	kplot_attach_smooth(cur->views[VIEW_SEXTI],
		sim->bufs.iextinct->cold, KPLOT_LINES, &solidcfg,
		KSMOOTH_MOVAVG, &smth);
	kplot_attach_data(cur->views[VIEW_SEXTI],
		sim->bufs.iextinct->cold, KPLOT_LINES, &transcfg);

	/* Mean and smoothed lines. */
	kplot_attach_smooth(cur->views[VIEW_SMEAN],
		sim->bufs.means->cold, KPLOT_LINES, &solidcfg,
		KSMOOTH_MOVAVG, &smth);
	kplot_attach_data(cur->views[VIEW_SMEAN],
		sim->bufs.means->cold, KPLOT_LINES, &transcfg);

	/* Mean and stddev.  */
	ts[0] = ts[1] = KPLOT_LINES;
	stats[0] = sim->bufs.means->cold;
	stats[1] = sim->bufs.stddevs->cold;
	cfgs[0] = &solidcfg;
	cfgs[1] = &transcfg;
	kplot_attach_datas(cur->views[VIEW_DEV], 2, stats, ts,
		(const struct kdatacfg *const *)cfgs,
		KPLOTS_YERRORLINE);

	/* Island mean and stddev. */
	ts[0] = ts[1] = KPLOT_POINTS;
	stats[0] = sim->bufs.islandmeans->cold;
	stats[1] = sim->bufs.islandstddevs->cold;
	cfgs[0] = &solidcfg;
	cfgs[1] = &transcfg;
	kplot_attach_datas(cur->views[VIEW_ISLANDERMEAN], 2, stats, ts,
		(const struct kdatacfg *const *)cfgs,
		KPLOTS_YERRORBAR);

	/* Island mean and stddev. */
	ts[0] = ts[1] = KPLOT_POINTS;
	stats[0] = sim->bufs.imeans->cold;
	stats[1] = sim->bufs.istddevs->cold;
	cfgs[0] = &solidcfg;
	cfgs[1] = &transcfg;
	kplot_attach_datas(cur->views[VIEW_ISLANDMEAN], 2, stats, ts,
		(const struct kdatacfg *const *)cfgs,
		KPLOTS_YERRORBAR);

	/* Time PMF views. */
	kplot_attach_smooth(cur->views[VIEW_TIMESPDF],
		sim->bufs.times->cold, KPLOT_LINES, &solidcfg,
		KSMOOTH_PMF, NULL);
	kplot_attach_smooth(cur->views[VIEW_TIMESCDF],
		sim->bufs.times->cold, KPLOT_LINES, &solidcfg,
		KSMOOTH_CDF, NULL);

	/* Mean PMF views. */
	kplot_attach_smooth(cur->views[VIEW_MEANMINPDF],
		sim->bufs.meanmins, KPLOT_LINES, &solidcfg,
		KSMOOTH_PMF, NULL);
	kplot_attach_smooth(cur->views[VIEW_MEANMINCDF],
		sim->bufs.meanmins, KPLOT_LINES, &solidcfg,
		KSMOOTH_CDF, NULL);

	/* Mutant PMF views. */
	kplot_attach_smooth(cur->views[VIEW_EXTMMAXPDF],
		sim->bufs.mextinctmaxs, KPLOT_LINES, &solidcfg,
		KSMOOTH_PMF, NULL);
	kplot_attach_smooth(cur->views[VIEW_EXTMMAXCDF],
		sim->bufs.mextinctmaxs, KPLOT_LINES, &solidcfg,
		KSMOOTH_CDF, NULL);

	/* Incumbent PMF views. */
	kplot_attach_smooth(cur->views[VIEW_EXTIMINPDF],
		sim->bufs.iextinctmins, KPLOT_LINES, &solidcfg,
		KSMOOTH_PMF, NULL);
	kplot_attach_smooth(cur->views[VIEW_EXTIMINCDF],
		sim->bufs.iextinctmins, KPLOT_LINES, &solidcfg,
		KSMOOTH_CDF, NULL);

	/* Fit poly PMF views. */
	kplot_attach_smooth(cur->views[VIEW_POLYMINPDF],
		sim->bufs.fitpolymins, KPLOT_LINES, &solidcfg,
		KSMOOTH_PMF, NULL);
	kplot_attach_smooth(cur->views[VIEW_POLYMINCDF],
		sim->bufs.fitpolymins, KPLOT_LINES, &solidcfg,
		KSMOOTH_CDF, NULL);

	kplot_attach_data(cur->views[VIEW_MEANMINQ],
		sim->bufs.meanminqbuf, KPLOT_LINES, &solidcfg);

	kplot_attach_data(cur->views[VIEW_POLYMINQ],
		sim->bufs.fitminqbuf, KPLOT_LINES, &solidcfg);
}

/*
 * Indicate that a given simulation is now being referenced by a new
 * window.
 */
static void
sim_ref(gpointer dat, gpointer unused)
{
	struct sim	*sim = dat;

	++sim->refs;
	g_debug("%p: Simulation ref (now %zu)", sim, sim->refs);
}

/*
 * Transfer the other widget's data into our own.
 */
void
on_drag_recv(GtkWidget *widget, GdkDragContext *ctx,
	gint x, gint y, GtkSelectionData *sel,
	guint target, guint time, gpointer dat)
{
	GObject		*srcptr, *dstptr;
	struct curwin	*cur;
	GList		*srcsims, *l, *ll;

	/* Get pointers to our and the other's window. */
	assert(NULL != sel);
	dstptr = G_OBJECT(gtk_widget_get_toplevel(widget));
	srcptr = G_OBJECT(*(void **)
		gtk_selection_data_get_data(sel));
	assert(NULL != srcptr);
	gtk_drag_finish(ctx, TRUE, FALSE, time);

	/* Don't copy into ourselves. */
	if (dstptr == srcptr)
		return;

	/* Get the simulation lists. */
	cur = g_object_get_data(srcptr, "cfg");
	srcsims = cur->sims;
	assert(NULL != srcsims);
	cur = g_object_get_data(dstptr, "cfg");

	/* Concatenate the simulation lists. */
	/* XXX: use g_list_concat? */
	for (l = srcsims; NULL != l; l = l->next) {
		g_debug("%p: Copying simulation", l->data);
		for (ll = cur->sims; NULL != ll; ll = ll->next)
			if (ll->data == l->data)
				break;

		if (NULL != ll) {
			g_debug("Simulation %p duplicate", l->data);
			continue;
		}

		sim_ref(l->data, NULL);
		cur->sims = g_list_append(cur->sims, l->data);
		window_add_sim(cur, l->data);
	}
}

/*
 * Signifity sight-unseen that we can transfer data.
 * We're only using DnD for one particular thing, so there's no need for
 * elaborate security measures.
 */
gboolean
on_drag_drop(GtkWidget *widget, GdkDragContext *ctx,
	gint x, gint y, guint time, gpointer dat)
{

	gtk_drag_get_data(widget, ctx, 0, time);
	return(TRUE);
}

/*
 * Send our identifier to the destination of the DnD.
 * They'll query our data separately.
 */
void
on_drag_get(GtkWidget *widget, GdkDragContext *ctx,
	GtkSelectionData *sel, guint targ, guint time, gpointer dat)
{
	void	*ptr;

	ptr = gtk_widget_get_toplevel(widget);
	gtk_selection_data_set(sel,
		gtk_selection_data_get_target(sel),
		sizeof(intptr_t) * 8,
		(const guchar *)&ptr, sizeof(intptr_t));
}

/*
 * Initialise a simulation (or simulations) window.
 * This is either called when we've just made a new simulation
 */
static void
window_init(struct bmigrate *b, struct curwin *cur, GList *sims)
{
	GtkBuilder	*builder;
	GtkTargetEntry   target;
	GList		*l;
	struct kdata	*stats[2];
	enum kplottype	 ts[2];
	size_t		 i;
	struct kplotcfg	 cfg;

	builder = builder_get("simwin.glade");
	g_assert(NULL != builder);
	swin_init(&cur->wins, cur->view, builder);
	gtk_builder_connect_signals(builder, cur);
	g_object_unref(G_OBJECT(builder));

	cur->winmean = kdata_vector_alloc(1);
	cur->winstddev = kdata_vector_alloc(1);
	cur->winfitmean = kdata_vector_alloc(1);
	cur->winfitstddev = kdata_vector_alloc(1);
	cur->winmextinctmean = kdata_vector_alloc(1);
	cur->winmextinctstddev = kdata_vector_alloc(1);
	cur->winiextinctmean = kdata_vector_alloc(1);
	cur->winiextinctstddev = kdata_vector_alloc(1);

	for (i = 0; i < VIEW__MAX; i++) {
		kplotcfg_defaults(&cfg);
		switch (i) {
		case (VIEW_DEV):
		case (VIEW_POLY):
		case (VIEW_EXTIMINS):
		case (VIEW_EXTMMAXS):
		case (VIEW_MEANMINS):
		case (VIEW_POLYMINS):
		case (VIEW_ISLANDMEAN):
		case (VIEW_ISLANDERMEAN):
			cfg.extrema_ymin = 0.0;
			cfg.extrema = EXTREMA_YMIN;
			break;
		default:
			break;
		}

		switch (i) {
		case (VIEW_DEV):
		case (VIEW_EXTI):
		case (VIEW_EXTIMINCDF):
		case (VIEW_EXTIMINPDF):
		case (VIEW_EXTM):
		case (VIEW_EXTMMAXCDF):
		case (VIEW_EXTMMAXPDF):
		case (VIEW_MEAN):
		case (VIEW_MEANMINCDF):
		case (VIEW_MEANMINPDF):
		case (VIEW_POLY):
		case (VIEW_POLYMINCDF):
		case (VIEW_POLYMINPDF):
		case (VIEW_SEXTM):
		case (VIEW_SEXTI):
		case (VIEW_SMEAN):
			cfg.xaxislabel = "incumbent";
			break;
		case (VIEW_ISLANDMEAN):
		case (VIEW_ISLANDERMEAN):
			cfg.xaxislabel = "island";
			break;
		case (VIEW_TIMESCDF):
		case (VIEW_TIMESPDF):
			cfg.xaxislabel = "exit time";
			break;
		case (VIEW_EXTIMINS):
		case (VIEW_EXTMMAXS):
		case (VIEW_MEANMINS):
		case (VIEW_POLYMINS):
			cfg.xaxislabel = "simulation";
			break;
		case (VIEW_MEANMINQ):
		case (VIEW_POLYMINQ):
			cfg.xaxislabel = "relative time";
			break;
		default:
			abort();
		}

		cfg.yaxislabelrot = M_PI / 2;

		switch (i) {
		case (VIEW_POLY):
		case (VIEW_DEV):
		case (VIEW_MEAN):
		case (VIEW_SMEAN):
			cfg.y2axislabel = "mutant fraction";
			break;
		case (VIEW_SEXTI):
		case (VIEW_EXTI):
			cfg.y2axislabel = "incumbent extinction";
			break;
		case (VIEW_SEXTM):
		case (VIEW_EXTM):
			cfg.y2axislabel = "mutant extinction";
			break;
		case (VIEW_EXTIMINCDF):
		case (VIEW_EXTIMINPDF):
		case (VIEW_EXTMMAXCDF):
		case (VIEW_EXTMMAXPDF):
		case (VIEW_MEANMINCDF):
		case (VIEW_MEANMINPDF):
		case (VIEW_POLYMINCDF):
		case (VIEW_POLYMINPDF):
		case (VIEW_TIMESCDF):
		case (VIEW_TIMESPDF):
			break;
		case (VIEW_ISLANDMEAN):
		case (VIEW_ISLANDERMEAN):
			cfg.y2axislabel = "mutant fraction mean";
			break;
		case (VIEW_EXTMMAXS):
		case (VIEW_EXTIMINS):
		case (VIEW_MEANMINS):
		case (VIEW_POLYMINS):
		case (VIEW_MEANMINQ):
		case (VIEW_POLYMINQ):
			cfg.y2axislabel = "incumbent";
			break;
		default:
			abort();
		}
		cur->views[i] = kplot_alloc(&cfg);
		g_assert(NULL != cur->views[i]);
	}

	ts[0] = ts[1] = KPLOT_POINTS;
	stats[0] = cur->winmean;
	stats[1] = cur->winstddev;
	kplot_attach_datas(cur->views[VIEW_MEANMINS], 2,
		stats, ts, NULL, KPLOTS_YERRORBAR);

	ts[0] = ts[1] = KPLOT_POINTS;
	stats[0] = cur->winfitmean;
	stats[1] = cur->winfitstddev;
	kplot_attach_datas(cur->views[VIEW_POLYMINS], 2,
		stats, ts, NULL, KPLOTS_YERRORBAR);

	ts[0] = ts[1] = KPLOT_POINTS;
	stats[0] = cur->winmextinctmean;
	stats[1] = cur->winmextinctstddev;
	kplot_attach_datas(cur->views[VIEW_EXTMMAXS], 2,
		stats, ts, NULL, KPLOTS_YERRORBAR);

	ts[0] = ts[1] = KPLOT_POINTS;
	stats[0] = cur->winiextinctmean;
	stats[1] = cur->winiextinctstddev;
	kplot_attach_datas(cur->views[VIEW_EXTIMINS], 2,
		stats, ts, NULL, KPLOTS_YERRORBAR);

	cur->redraw = 1;
	cur->sims = sims;
	cur->b = b;

	g_object_set_data_full(G_OBJECT
		(cur->wins.window), "cfg", cur, curwin_free);
	b->windows = g_list_append(b->windows, cur);

	/*
	 * Coordinate drag-and-drop.
	 * We only allow drag-and-drop between us and other windows of
	 * this same simulation.
	 */
	target.target = g_strdup("integer");
	target.flags = GTK_TARGET_SAME_APP|GTK_TARGET_OTHER_WIDGET;
	target.info = 0;

	gtk_drag_dest_set(GTK_WIDGET(cur->wins.draw),
		GTK_DEST_DEFAULT_ALL, &target, 1, GDK_ACTION_COPY);
	gtk_drag_source_set(GTK_WIDGET(cur->wins.draw),
		GDK_BUTTON1_MASK, &target, 1, GDK_ACTION_COPY);
	g_free(target.target);

	for (l = cur->sims; NULL != l; l = g_list_next(l))
		window_add_sim(cur, l->data);
}

/*
 * Clone the current window.
 * We create a new window, intialised in the same way as for a new
 * single simulation, but give it an existing list of simulations.
 */
void
onclone(GtkMenuItem *menuitem, gpointer dat)
{
	struct curwin	*oldcur = dat, *newcur;
	GList		*oldsims, *newsims;

	oldsims = oldcur->sims;
	g_list_foreach(oldsims, sim_ref, NULL);
	newsims = g_list_copy(oldsims);
	newcur = g_malloc0(sizeof(struct curwin));
	newcur->view = oldcur->view;
	window_init(oldcur->b, newcur, newsims);
}

/*
 * Brute-force scan all possible pi values (and Poisson means) by
 * scanning through the strategy space.
 */
static gboolean
on_rangefind_idle(gpointer dat)
{

	return(rangefind(dat));
}

static int
entry2func(GtkEntry *entry, struct hnode ***exp, GtkLabel *error)
{
	const gchar	*txt;
	GdkRGBA	 	 bad = { 1.0, 0.0, 0.0, 0.5 };

	txt = gtk_entry_get_text(entry);
	*exp = hnode_parse((const gchar **)&txt);
	if (NULL != *exp) {
		gtk_widget_override_background_color
			(GTK_WIDGET(entry),
			 GTK_STATE_FLAG_NORMAL, NULL);
		return(1);
	}

	gtk_label_set_text(error, "Error: not a function.");
	gtk_widget_show_all(GTK_WIDGET(error));
	gtk_widget_override_background_color
		(GTK_WIDGET(entry), GTK_STATE_FLAG_NORMAL, &bad);
	return(0);
}

static int
entryworder(GtkEntry *mine, GtkEntry *maxe,
	double min, double max, GtkLabel *error)
{
	GdkRGBA	 bad = { 1.0, 0.0, 0.0, 0.5 };

	if (min < max) {
		gtk_widget_override_background_color
			(GTK_WIDGET(mine),
			 GTK_STATE_FLAG_NORMAL, NULL);
		gtk_widget_override_background_color
			(GTK_WIDGET(maxe),
			 GTK_STATE_FLAG_NORMAL, NULL);
		return(1);
	}

	gtk_label_set_text(error, "Error: bad weak ordering.");
	gtk_widget_show_all(GTK_WIDGET(error));
	gtk_widget_override_background_color
		(GTK_WIDGET(mine), GTK_STATE_FLAG_NORMAL, &bad);
	gtk_widget_override_background_color
		(GTK_WIDGET(maxe), GTK_STATE_FLAG_NORMAL, &bad);
	return(0);
}

static int
entryorder(GtkEntry *mine, GtkEntry *maxe,
	double min, double max, GtkLabel *error)
{
	GdkRGBA	 bad = { 1.0, 0.0, 0.0, 0.5 };

	if (min <= max) {
		gtk_widget_override_background_color
			(GTK_WIDGET(mine),
			 GTK_STATE_FLAG_NORMAL, NULL);
		gtk_widget_override_background_color
			(GTK_WIDGET(maxe),
			 GTK_STATE_FLAG_NORMAL, NULL);
		return(1);
	}

	gtk_label_set_text(error, "Error: bad ordering.");
	gtk_widget_show_all(GTK_WIDGET(error));
	gtk_widget_override_background_color
		(GTK_WIDGET(mine), GTK_STATE_FLAG_NORMAL, &bad);
	gtk_widget_override_background_color
		(GTK_WIDGET(maxe), GTK_STATE_FLAG_NORMAL, &bad);
	return(0);
}

static int
entry2double(GtkEntry *entry, gdouble *sz, GtkLabel *error)
{
	gchar	*ep;
	GdkRGBA	 bad = { 1.0, 0.0, 0.0, 0.5 };

	*sz = g_ascii_strtod(gtk_entry_get_text(entry), &ep);

	if (ERANGE != errno && '\0' == *ep) {
		gtk_widget_override_background_color
			(GTK_WIDGET(entry),
			 GTK_STATE_FLAG_NORMAL, NULL);
		return(1);
	}

	gtk_label_set_text(error, "Error: not a decimal number.");
	gtk_widget_show_all(GTK_WIDGET(error));
	gtk_widget_override_background_color
		(GTK_WIDGET(entry), GTK_STATE_FLAG_NORMAL, &bad);
	return(0);
}

static int
entry2size(GtkEntry *entry, size_t *sz, GtkLabel *error, size_t min)
{
	guint64	 v;
	gchar	*ep;
	GdkRGBA	 bad = { 1.0, 0.0, 0.0, 0.5 };

	v = g_ascii_strtoull(gtk_entry_get_text(entry), &ep, 10);
	if (ERANGE == errno || '\0' != *ep || v >= SIZE_MAX) {
		gtk_label_set_text
			(error, "Error: not a natural number.");
		gtk_widget_show_all(GTK_WIDGET(error));
		gtk_widget_override_background_color
			(GTK_WIDGET(entry),
			 GTK_STATE_FLAG_NORMAL, &bad);
		return(0);
	} else if (v < min) {
		gtk_label_set_text
			(error, "Error: number too small.");
		gtk_widget_show_all(GTK_WIDGET(error));
		gtk_widget_override_background_color
			(GTK_WIDGET(entry),
			 GTK_STATE_FLAG_NORMAL, &bad);
		return(0);
	}

	*sz = (size_t)v;
	gtk_widget_override_background_color
		(GTK_WIDGET(entry),
		 GTK_STATE_FLAG_NORMAL, NULL);
	return(1);
}

/*
 * We have various ways of auto-setting the name of the simulation.
 * By default, it's set to the current date-time.
 */
static void
donamefill(struct hwin *c)
{
	enum namefill	 v;
	GTimeVal	 gt;
	enum input	 input;
	gchar		 buf[1024];
	gchar		*bufp;
	enum mutants	 mutants;

	input = gtk_notebook_get_current_page(c->inputs);
	g_assert(input < INPUT__MAX);

	for (mutants = 0; mutants < MUTANTS__MAX; mutants++)
		if (gtk_toggle_button_get_active
			(GTK_TOGGLE_BUTTON(c->mutants[mutants])))
			break;

	for (v = 0; v < NAMEFILL__MAX; v++)
		if (gtk_toggle_button_get_active(c->namefill[v]))
			break;

	bufp = buf;
	switch (v) {
	case (NAMEFILL_M):
		g_snprintf(buf, sizeof(buf), "m=%s",
			gtk_entry_get_text(c->migrate[input]));
		break;
	case (NAMEFILL_T):
		g_snprintf(buf, sizeof(buf), "T=%s",
			gtk_entry_get_text(c->stop));
		break;
	case (NAMEFILL_MUTANTS):
		if (MUTANTS_DISCRETE == mutants) {
			g_snprintf(buf, sizeof(buf),
				"discrete [%s,%s)",
				gtk_entry_get_text(c->xmin),
				gtk_entry_get_text(c->xmax));
			break;
		}
		g_snprintf(buf, sizeof(buf),
			"Gaussian s=%s, [%s,%s)",
			gtk_entry_get_text(c->mutantsigma),
			gtk_entry_get_text(c->ymin),
			gtk_entry_get_text(c->ymax));
		break;
	case (NAMEFILL_DATE):
		g_get_current_time(&gt);
		bufp = g_time_val_to_iso8601(&gt);
		break;
	case (NAMEFILL_NONE):
		return;
	default:
		abort();
	}

	gtk_entry_set_text(c->name, bufp);
	if (bufp != buf)
		g_free(bufp);
}

void
onnametoggle(GtkToggleButton *editable, gpointer dat)
{
	struct bmigrate	*b = dat;

	donamefill(&b->wins);
}

void
onnameupdate(GtkEditable *editable, gpointer dat)
{
	struct bmigrate	*b = dat;

	donamefill(&b->wins);
}

/*
 * We want to run the given simulation.
 * First, verify all data; then, start the simulation; lastly, open a
 * window assigned specifically to that simulation.
 */
void
onactivate(GtkButton *button, gpointer dat)
{
	gint	 	  input;
	struct bmigrate	 *b = dat;
	struct hnode	**exp;
	enum mapmigrant	  migrants;
	GtkWidget	 *w;
	struct kml	 *kml;
	GList		 *l, *cl;
	GtkLabel	 *err = b->wins.error;
	const gchar	 *name, *func;
	gchar	  	 *file;
	gdouble		**ms;
	gdouble		  xmin, xmax, delta, alpha, m, sigma,
			  ymin, ymax, idcoef, strat;
	enum mutants	  mutants;
	size_t		  i, totalpop, islands, stop, ideathmean,
			  slices, islandpop, mapindexfix;
	size_t		 *islandpops;
	struct sim	 *sim;
	struct curwin	 *cur;
	struct kmlplace	 *kmlp;
	GError		 *er;
	enum maptop	  maptop;
	enum mapindex	  mapindex;

	sigma = idcoef = 0.0;
	islandpops = NULL;
	islandpop = 0;
	ms = NULL;
	kml = NULL;
	exp = NULL;

	if ( ! entry2size(b->wins.stop, &stop, err, 1))
		goto cleanup;

	for (migrants = 0; migrants < MAPMIGRANT__MAX; migrants++)
		if (gtk_toggle_button_get_active
			 (b->wins.mapmigrants[migrants]))
			break;
	g_assert(migrants < MAPMIGRANT__MAX);

	for (mapindex = 0; mapindex < MAPINDEX__MAX; mapindex++)
		if (gtk_toggle_button_get_active
			 (b->wins.mapindices[mapindex]))
			break;
	g_assert(mapindex < MAPINDEX__MAX);
	mapindexfix = MAPINDEX_FIXED != mapindex ? 0 :
		gtk_adjustment_get_value(b->wins.mapindexfix);

	for (maptop = 0; maptop < MAPTOP__MAX; maptop++)
		if (gtk_toggle_button_get_active
			 (b->wins.maptop[maptop]))
			break;
	g_assert(maptop < MAPTOP__MAX);

	input = gtk_notebook_get_current_page(b->wins.inputs);
	switch (input) {
	case (INPUT_UNIFORM):
		/*
		 * In the simplest possible case, we use uniform values
		 * for both migration (no inter-island migration) and
		 * populations.
		 */
		islands = gtk_adjustment_get_value(b->wins.islands);
		islandpop = gtk_adjustment_get_value(b->wins.pop);
		break;
	case (INPUT_VARIABLE):
		/*
		 * Variable number of islands.
		 * We also add a check to see if we're really running
		 * with different island sizes or not.
		 */
		if ( ! entry2double(b->wins.ideathcoef, &idcoef, err))
			goto cleanup;
		if (idcoef < 0.0 || idcoef > 1.0) {
			gtk_label_set_text(err, "Death coefficient "
				"must be in the unit interval.");
			gtk_widget_show_all(GTK_WIDGET(err));
			goto cleanup;
		}
		l = gtk_container_get_children
			(GTK_CONTAINER(b->wins.mapbox));
		islands = g_list_length(l);
		islandpops = g_malloc0_n(islands, sizeof(size_t));
		for (i = 0; i < islands; i++) {
			w = GTK_WIDGET(g_list_nth_data(l, i));
			cl = gtk_container_get_children(GTK_CONTAINER(w));
			islandpops[i] = gtk_spin_button_get_value_as_int
				(GTK_SPIN_BUTTON(g_list_next(cl)->data));
			g_list_free(cl);
		}
		g_list_free(l);
		break;
	case (INPUT_MAPPED):
		/*
		 * Possibly variable island sizes, possibly variable
		 * inter-island migration.
		 */
		switch (maptop) {
		case (MAPTOP_RECORD):
			/* Try to read KML from file. */
			file = gtk_file_chooser_get_filename
				(b->wins.mapfile);
			if (NULL == file) {
				gtk_label_set_text(err, "Error: "
					"map file not specified.");
				gtk_widget_show_all(GTK_WIDGET(err));
				goto cleanup;
			}
			er = NULL;
			kml = kml_parse(file, &er);
			g_free(file);
			if (NULL == kml) {
				file = g_strdup_printf("Error: "
					"bad map file: %s",
					NULL != er ? er->message :
					"cannot load file");
				gtk_label_set_text(err, file);
				gtk_widget_show_all(GTK_WIDGET(err));
				g_free(file);
				if (NULL != er)
					g_error_free(er);
				goto cleanup;
			}
			break;
		case (MAPTOP_RAND):
			islands = gtk_adjustment_get_value
				(b->wins.maprandislands);
			islandpop = gtk_adjustment_get_value
				(b->wins.maprandislanders);
			kml = kml_rand(islands, islandpop);
			islandpop = islands = 0;
			break;
		case (MAPTOP_TORUS):
			islands = gtk_adjustment_get_value
				(b->wins.maptorusislands);
			islandpop = gtk_adjustment_get_value
				(b->wins.maptorusislanders);
			kml = kml_torus(islands, islandpop);
			islandpop = islands = 0;
			break;
		default:
			abort();
		}

		/*
		 * Grok our island populations from the input file.
		 * This will have a reasonable default, but make sure
		 * anyway with some assertions.
		 */
		islands = (size_t)g_list_length(kml->kmls);
		islandpops = g_malloc0_n(islands, sizeof(size_t));
		for (i = 0; i < islands; i++) {
			kmlp = g_list_nth_data(kml->kmls, i);
			islandpops[i] = kmlp->pop;
		}

		/*
		 * If we're uniformly migrating, then stop processing
		 * right now.
		 * If we're distance-migrating, then have the kml file
		 * set the inter-island migration probabilities.
		 */
		switch (migrants) {
		case (MAPMIGRANT_UNIFORM):
			ms = kml_migration_distance(kml->kmls, maptop);
			break;
		case (MAPMIGRANT_NEAREST):
			ms = kml_migration_nearest(kml->kmls, maptop);
			break;
		case (MAPMIGRANT_TWONEAREST):
			ms = kml_migration_twonearest(kml->kmls, maptop);
			break;
		default:
			break;
		}
		break;
	default:
		abort();
	}

	/*
	 * Base check: we need at least two islands.
	 */
	if (islands < 2) {
		gtk_label_set_text(err,
			"Error: need at least two islands.");
		gtk_widget_show_all(GTK_WIDGET(err));
		goto cleanup;
	}

	ideathmean = gtk_adjustment_get_value(b->wins.ideathmean);

	/*
	 * If we have an array of island populations, make sure that
	 * each has more than two islanders.
	 * While here, revert to a single island population size if our
	 * sizes are, in fact, uniform.
	 * Calculate our total population size as well.
	 */
	if (NULL != islandpops) {
		g_assert(0 == islandpop);
		for (i = 0; i < islands; i++) {
			if (islandpops[i] > 1)
				continue;
			gtk_label_set_text(err, "Error: need at "
				"least two islanders per island.");
			gtk_widget_show_all(GTK_WIDGET(err));
			goto cleanup;
		}
		for (totalpop = i = 0; i < islands; i++)
			totalpop += islandpops[i];

		/* Check for uniformity. */
		for (i = 1; i < islands; i++)
			if (islandpops[i] != islandpops[0])
				break;
		if (i == islands && 0 == ideathmean) {
			g_debug("Reverting to uniform island "
				"populations: all islands have "
				"the same (with no death process): "
				"%zu", islandpops[0]);
			islandpop = islandpops[0];
			g_free(islandpops);
			islandpops = NULL;
		}
	}

	/*
	 * Handle uniform island sizes.
	 * (This isn't part of the conditional above because we might
	 * set ourselves to be uniform whilst processing.)
	 */
	if (NULL == islandpops && islandpop < 2) {
		gtk_label_set_text(err, "Error: need at "
			"least two islanders per island.");
		gtk_widget_show_all(GTK_WIDGET(err));
		goto cleanup;
	} else if (NULL == islandpops)
		totalpop = islands * islandpop;

	for (mutants = 0; mutants < MUTANTS__MAX; mutants++)
		if (gtk_toggle_button_get_active
			(GTK_TOGGLE_BUTTON(b->wins.mutants[mutants])))
			break;

	if ( ! entry2double(b->wins.xmin, &xmin, err))
		goto cleanup;
	if ( ! entry2double(b->wins.xmax, &xmax, err))
		goto cleanup;
	if ( ! entryworder(b->wins.xmin, b->wins.xmax, xmin, xmax, err))
		goto cleanup;
	ymin = xmin;
	ymax = xmax;
	if (MUTANTS_GAUSSIAN == mutants) {
		if ( ! entry2double(b->wins.ymin, &ymin, err))
			goto cleanup;
		if ( ! entry2double(b->wins.ymax, &ymax, err))
			goto cleanup;
		if ( ! entryorder(b->wins.ymin,
				b->wins.xmin, ymin, xmin, err))
			goto cleanup;
		if ( ! entryorder(b->wins.xmax,
				b->wins.ymax, xmax, ymax, err))
			goto cleanup;
		if ( ! entryworder(b->wins.ymin,
				b->wins.ymax, ymin, ymax, err))
			goto cleanup;
		if ( ! entry2double(b->wins.mutantsigma, &sigma, err))
			goto cleanup;
	}
	if ( ! entry2double(b->wins.alpha, &alpha, err))
		goto cleanup;
	if ( ! entry2double(b->wins.delta, &delta, err))
		goto cleanup;
	if ( ! entry2double(b->wins.migrate[input], &m, err))
		goto cleanup;
	if ( ! entry2size(b->wins.incumbents, &slices, err, 1))
		goto cleanup;
	if ( ! entry2func(b->wins.func, &exp, err))
		goto cleanup;

	func = gtk_entry_get_text(b->wins.func);
	name = gtk_entry_get_text(b->wins.name);

	if ('\0' == *name)
		name = "unnamed";

	gtk_widget_hide(GTK_WIDGET(err));

	/*
	 * All parameters check out!
	 * Allocate the simulation now.
	 */
	if (button == b->wins.buttonrange) {
		if (0 == b->rangeid) {
			g_debug("Starting rangefinder");
			b->rangeid = g_idle_add
				((GSourceFunc)on_rangefind_idle, b);
		} else
			g_debug("Re-using rangefinder");

		hnode_free(b->range.exp);
		if (NULL != islandpops) {
			b->range.n = 0;
			for (i = 0; i < islands; i++)
				b->range.n = islandpops[i] > b->range.n ?
					islandpops[i] : b->range.n;
		} else
			b->range.n = islandpop;

		b->range.exp = exp;
		b->range.alpha = alpha;
		b->range.delta = delta;
		b->range.slices = slices;
		b->range.slicex = b->range.slicey = 0;
		b->range.piaggr = 0.0;
		b->range.picount = 0;
		b->range.pimin = DBL_MAX;
		b->range.pimax = -DBL_MAX;
		b->range.xmin = b->range.ymin = xmin;
		b->range.xmax = b->range.ymax = xmax;
		if (MUTANTS_GAUSSIAN == mutants) {
			b->range.ymin = ymin;
			b->range.ymax = ymax;
		}
		exp = NULL;
		file = g_strdup_printf
			("%s, X=[%g, %g), Y=[%g, %g), n=%zu",
			func, b->range.xmin, b->range.xmax,
			b->range.ymin, b->range.ymax,
			b->range.n);
		gtk_label_set_text(b->wins.rangefunc, file);
		g_free(file);

		file = g_strdup_printf
			("%g(1 + %gπ)",
			 b->range.alpha, b->range.delta);
		gtk_label_set_markup(b->wins.rangeparms, file);
		g_free(file);

		gtk_widget_hide(GTK_WIDGET(b->wins.rangeerrorbox));
		gtk_widget_set_visible
			(GTK_WIDGET(b->wins.rangefind), TRUE);
		goto cleanup;
	}

	sim = g_malloc0(sizeof(struct sim));
	sim->dims = slices;
	sim->islands = islands;
	sim->mutants = mutants;
	sim->mutantsigma = sigma;
	sim->func = g_strdup(func);
	sim->name = g_strdup(name);
	sim->fitpoly = gtk_adjustment_get_value(b->wins.fitpoly);
	sim->weighted = gtk_toggle_button_get_active(b->wins.weighted);
	sim->smoothing = gtk_adjustment_get_value(b->wins.smoothing);
	sim->ideathmean = ideathmean;
	sim->ideathcoef = idcoef;
	sim->kml = kml;
	sim->migrant = migrants;
	sim->maptop = maptop;
	sim->mapindex = mapindex;
	if ((sim->mapindexfix = mapindexfix) > sim->islands) {
		g_debug("Clamping map index to %zu", sim->islands - 1);
		sim->mapindexfix = sim->islands - 1;
	}
	g_mutex_init(&sim->hot.mux);
	g_cond_init(&sim->hot.cond);

	/* Source for the fraction of mutants. */
	sim->bufs.fractions = kdata_array_alloc(NULL, slices);
	g_assert(NULL != sim->bufs.fractions);
	sim->bufs.ifractions = kdata_array_alloc(NULL, islands);
	g_assert(NULL != sim->bufs.ifractions);
	sim->bufs.mutants = kdata_array_alloc(NULL, slices);
	g_assert(NULL != sim->bufs.mutants);
	sim->bufs.incumbents = kdata_array_alloc(NULL, slices);
	g_assert(NULL != sim->bufs.incumbents);
	sim->bufs.islands = kdata_array_alloc(NULL, islands);
	g_assert(NULL != sim->bufs.islands);
	sim->bufs.meanmins = kdata_array_alloc(NULL, slices);
	g_assert(NULL != sim->bufs.meanmins);
	sim->bufs.mextinctmaxs = kdata_array_alloc(NULL, slices);
	g_assert(NULL != sim->bufs.mextinctmaxs);
	sim->bufs.iextinctmins = kdata_array_alloc(NULL, slices);
	g_assert(NULL != sim->bufs.iextinctmins);
	sim->bufs.fitpoly = kdata_array_alloc(NULL, slices);
	g_assert(NULL != sim->bufs.fitpoly);
	sim->bufs.fitpolybuf = kdata_buffer_alloc(slices);
	g_assert(NULL != sim->bufs.fitpolybuf);
	sim->bufs.fitpolymins = kdata_array_alloc(NULL, slices);
	g_assert(NULL != sim->bufs.fitpolymins);
	sim->bufs.meanminqbuf = kdata_array_alloc(NULL, CQUEUESZ);
	g_assert(NULL != sim->bufs.meanminqbuf);
	sim->bufs.fitminqbuf = kdata_array_alloc(NULL, CQUEUESZ);
	g_assert(NULL != sim->bufs.fitminqbuf);

	for (i = 0; i < slices; i++) {
		strat = xmin + (xmax - xmin) * (i / (double)slices);
		kdata_array_set(sim->bufs.fractions, i, strat, 0);
		kdata_array_set(sim->bufs.mutants, i, strat, 0);
		kdata_array_set(sim->bufs.incumbents, i, strat, 0);
		kdata_array_set(sim->bufs.meanmins, i, strat, 0);
		kdata_array_set(sim->bufs.mextinctmaxs, i, strat, 0);
		kdata_array_set(sim->bufs.iextinctmins, i, strat, 0);
		kdata_array_set(sim->bufs.fitpoly, i, strat, 0);
		kdata_array_set(sim->bufs.fitpolymins, i, strat, 0);
	}

	sim->bufs.times = simbuf_alloc
		(kdata_array_alloc(NULL, stop + 1), stop + 1);
	sim->bufs.islandmeans = simbuf_alloc
		(kdata_mean_alloc(sim->bufs.islands), islands);
	sim->bufs.islandstddevs = simbuf_alloc
		(kdata_stddev_alloc(sim->bufs.islands), islands);
	sim->bufs.imeans = simbuf_alloc
		(kdata_mean_alloc(sim->bufs.ifractions), islands);
	sim->bufs.istddevs = simbuf_alloc
		(kdata_stddev_alloc(sim->bufs.ifractions), islands);
	sim->bufs.means = simbuf_alloc
		(kdata_mean_alloc(sim->bufs.fractions), slices);
	sim->bufs.stddevs = simbuf_alloc
		(kdata_stddev_alloc(sim->bufs.fractions), slices);
	sim->bufs.mextinct = simbuf_alloc
		(kdata_mean_alloc(sim->bufs.mutants), slices);
	sim->bufs.iextinct = simbuf_alloc
		(kdata_mean_alloc(sim->bufs.incumbents), slices);

	/*
	 * Conditionally allocate our fitness polynomial structures.
	 * These are per-simulation as they're only run by one thread at
	 * any one time during the simulation.
	 */
	if (sim->fitpoly) {
		sim->work.X = gsl_matrix_alloc
			(sim->dims, sim->fitpoly + 1);
		sim->work.y = gsl_vector_alloc(sim->dims);
		sim->work.w = gsl_vector_alloc(sim->dims);
		sim->work.c = gsl_vector_alloc(sim->fitpoly + 1);
		sim->work.cov = gsl_matrix_alloc
			(sim->fitpoly + 1, sim->fitpoly + 1);
		sim->work.work = gsl_multifit_linear_alloc
			(sim->dims, sim->fitpoly + 1);
		sim->work.coeffs = g_malloc0_n
			(sim->fitpoly + 1, sizeof(double));
	}

	sim->nprocs = gtk_adjustment_get_value(b->wins.nthreads);
	sim->totalpop = totalpop;
	sim->stop = stop;
	sim->alpha = alpha;
	sim->colour = b->nextcolour++;
	sim->delta = delta;
	sim->m = m;
	sim->ms = ms;
	sim->pop = islandpop;
	sim->pops = islandpops;
	sim->input = input;
	sim->exp = exp;
	sim->xmin = xmin;
	sim->xmax = xmax;
	sim->ymin = ymin;
	sim->ymax = ymax;
	b->sims = g_list_append(b->sims, sim);
	sim_ref(sim, NULL);
	sim->threads = g_malloc0_n(sim->nprocs, sizeof(struct simthr));
	g_debug("%p: Simulation created", sim);

	/* Create the simulation threads. */
	for (i = 0; i < sim->nprocs; i++) {
		sim->threads[i].rank = i;
		sim->threads[i].sim = sim;
		sim->threads[i].thread = g_thread_new
			(NULL, simulation, &sim->threads[i]);
	}

	/* Create the simulation window. */
	cur = g_malloc0(sizeof(struct curwin));
	cur->view = VIEW_MEAN;
	window_init(b, cur, g_list_append(NULL, sim));

	/* Initialise the name of our simulation. */
	donamefill(&b->wins);
	return;
cleanup:
	if (NULL != ms)
		for (i = 0; i < islands; i++)
			g_free(ms[i]);
	g_free(islandpops);
	g_free(ms);
	hnode_free(exp);
	kml_free(kml);
}

ext: c

max_stars_repo_name: kristapsdz/bmigrate

max_issues_repo_name: kristapsdz/bmigrate

max_forks_repo_name: kristapsdz/bmigrate


text:
#ifndef PRIMAL_DUAL_LOOPLESS_H
#define PRIMAL_DUAL_LOOPLESS_H

#include "problem_data.h"

#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <sstream>
#include <omp.h>

//This class implements the loopless variance reduced type methods with arbitrary sampling

/*
The optimization problem to solve is:

\sum_{i=1}^n \lambda_i\phi_i(A_i^{\top} w)+ g(w)
Assumption 1: For each i, \phi_i is 1-smooth
By default, lambda_i=1/n for all i.
*
* The dual problem is
*         g*(-sum_{i=1}^n \lambda_i \alpha_i A_i)+\sum_{i=1}^n \lambda_i \phi*_i(\alpha_i)
*/

template<typename L, typename D>
class Primal_Dual_LOOPLESS: public problem_data<L, D>
{

private:

  // involved variables

  std::vector<D> proba_vector;

  std::vector<D>  tilde_proba_vector;

  std::vector<D> group_C;

  std::vector<D> index_group_C;

  std::vector<D> maxp_group_C;

  std::vector<D> isolated_I;

  std::vector<D> sump_group_C;
  vector<D> theta_S;   //theta_S in the paper

  std::vector<D> Li;

  std::vector<D> gk;     // the vector g^k in the paper

  // auxiliary variables
  std::vector<D> deltagk;
  std::vector<D> index_to_do_prox;
  std::vector<D> whether_or_not_to_do_prox;

  L nb_of_iters_per_loop;

  D primal_value;

  D dual_value;

  D epsilon;

  L max_nb_loops;

  D max_p;

  D maxv;

  D lambda_1;
  D lambda_2;

  D running_time;

  L nb_loops;

  L print_every_N;

  D delta;

  vector<L> batch_i;
  vector<L> my_batch;
  vector<D> batch_deltaalphai;

  L batch_size;

  L nb_groups;    //number of groups in group sampling

protected:

  std::vector<D> primal_x;     // the x variable

  string uniform;

  D Lf;

  D L2;

  D L1;

  D sumLi;

  D p; // the probability of changing w to x as in the paper
  L tau;  //number of threads on each node/computer

  D scaler;

  std::vector<D> dual_alpha;  //  dual_alpha[i]=\phi_i'(A_i't)

  std::vector<D> lambda_f;    // lambda_f[i]=lambda_i
  std::vector<D> baralpha;   // baralpha=sum_{i=1}^{nsamples} lambda_f[i]* dual_alpha[i]*A_i

  L current_nb_iters;

  L nb_iters;

  std::vector<D> last_update_i;

  ofstream samp;

public:

  gsl_rng * rng;
  virtual inline D gradient_of_phi_i(D, L){return D(NULL);}
  virtual inline D gradient_of_gstar_j(D, L){return D(NULL);}
  virtual inline D value_of_phi_i(D, L) {return D(NULL);}
  virtual inline D value_of_g_j(D, L){return D(NULL);}
  virtual inline D value_of_phistar_i(D,L) {return D(NULL);}
  virtual inline D value_of_gstar(vector<D> &){return D(NULL);}
  virtual inline D value_of_gstar_minus(vector<D> &, D){return D(NULL);}
  virtual inline D feasible_dual(vector<D> &){return D(NULL);}
  virtual inline D compute_delta_alpha(D,D,L){return D(NULL);}
  virtual inline void set_auxiliary_v(){}
  virtual inline D get_lambda1(){return D(NULL);}
  virtual inline D get_lambda2(){return D(NULL);}
  virtual inline D compute_current_xj_value(D , L, L , L){return D(NULL);}
  virtual inline D compute_current_xj_value_without_update(D , L, L , L){return D(NULL);}
  virtual inline void set_stepsize(){}
  virtual inline void update_x(D,L){}
  virtual inline void update_baralpha(){}

  Primal_Dual_LOOPLESS()
  : problem_data<L,D>()
  {

  }

  Primal_Dual_LOOPLESS(const char* matrix_file, const char* vector_file)
  : problem_data<L,D>(matrix_file, vector_file)
  {

  }

  Primal_Dual_LOOPLESS(const char* matrix_file)
  : problem_data<L,D>(matrix_file)
  {

  }

  void set_rng()
  {
    gsl_rng_env_setup();
    const gsl_rng_type * T;
    T = gsl_rng_default;
    rng = gsl_rng_alloc(T);
    gsl_rng_set(rng,time(NULL));
    //gsl_rng_set(rng, 27432042);

  }

  void set_print_every_N(L i){print_every_N=i;}

  void set_L1(L p){
    if(p==0){ //batch sampling mode  (stochastic process in the paper)
      L1=(1-1./tau)*Lf+L2;
    }else{ //group sampling mode
      L1=Lf+L2;
    }
    //  alpha_l=1/(4.*L1+2*L2);
  }

  void set_Li_Lf()
  {
    sumLi=0;
    maxv=0;
    D minv=std::numeric_limits<double>::max();
    Li.resize(this->nsamples,0);
    for(L i=0;i<this->nsamples;i++)
    {
      D vi=0;
      for (L k = this->ptr[i]; k < this->ptr[i + 1];k++)
      {
        vi+=lambda_f[i]*this->nsamples*this->A[k]*this->A[k];
      }
      Li[i]=vi;
      if(maxv<vi) maxv=vi;
      if(minv>vi) minv=vi;
    }
    for(L i=0;i<this->nsamples;i++){
      if(Li[i]<maxv*1e-15)
      {
        for (L k = this->ptr[i]; k < this->ptr[i + 1];k++)
        { this->A[k]=0;
          L j=this->row_idx[k];
          for (L kj = this->ptr_t[j]; kj < this->ptr_t[j + 1];	kj++)
          {
            L i2=this->col_idx[kj];
            if(i2==i) this->A_t[kj]=0;
          }
        }
        Li[i]=maxv*1e-15;
      }
      sumLi+=Li[i];
    }
    Lf=compute_lambda_max(10);
    cout<<"  max of v: "<<maxv<<"  min of v: "<<minv<<" sumof Li="<<sumLi<<" Lf= "<<Lf<<endl;
  }

  //compute the maximal eigenvalue of sum_{i=1}^m \lambda_f[i]*A_i*A_i' by power iteration
  D compute_lambda_max(L K){

    std::vector<D> bk(this->nfeatures);
    for (L j=0;j<this->nfeatures;j++)
    {
      bk[j]=1;
    }
    std::vector<D> yk(this->nsamples);
    D normk;
    D tmp;
    for(L kk=0;kk<K;kk++){
      for (L i=0;i<this->nsamples;i++){
        tmp=0;
        for (L k = this->ptr[i]; k < this->ptr[i + 1];	k++)
        {
          L j=this->row_idx[k];
          tmp+=this->A[k]*bk[j];
        }
        yk[i]=tmp;
      }
      normk=0;
      for (L j=0;j<this->nfeatures;j++){
        bk[j]=0;
        for (L k = this->ptr_t[j]; k < this->ptr_t[j + 1];	k++)
        {
          L i=this->col_idx[k];
          bk[j]+=this->A_t[k]*yk[i]*lambda_f[i];
        }
        normk+=bk[j]*bk[j];
      }
      normk=sqrt(normk);
      for (L j=0;j<this->nfeatures;j++)
      {bk[j]=bk[j]/normk; }
    }
    cout<<endl;
    D res=0;
    normk=0;
    for (L i=0;i<this->nsamples;i++){
      tmp=0;
      for (L k = this->ptr[i]; k < this->ptr[i + 1];	k++)
      {
        L j=this->row_idx[k];
        tmp+=this->A[k]*bk[j];
      }
      yk[i]=tmp;
      normk+=yk[i]*yk[i];
    }
    std::vector<D> bk2(this->nfeatures);
    for (L j=0;j<this->nfeatures;j++){
      bk2[j]=0;
      for (L k = this->ptr_t[j]; k < this->ptr_t[j + 1];	k++)
      {
        L i=this->col_idx[k];
        bk2[j]+=this->A_t[k]*yk[i]*lambda_f[i];
      }
    }
    for (L j=0;j<this->nfeatures;j++)
    res+=bk2[j]*bk[j];
    return res;
  }

  void set_optimal_probability()
  {
    tilde_proba_vector.clear();
    tilde_proba_vector.resize(this->nsamples,0);
    proba_vector.clear();
    proba_vector.resize(this->nsamples,0);
    D sum=0;
    for(L i=0; i<this->nsamples;i++)
    {
      tilde_proba_vector[i]=Li[i];
      sum+=Li[i];
    }
    max_p=0;
    for(L i=0; i<this->nsamples;i++)
    {
      tilde_proba_vector[i]=tilde_proba_vector[i]/sum;
      proba_vector[i]=1-pow(1-tilde_proba_vector[i],tau);
      if(max_p<tilde_proba_vector[i])
      {
        max_p=tilde_proba_vector[i];
      }
    }
    cout<<"sum="<<sum<<"; proba"<<tilde_proba_vector[0]<<endl;
  }

  void set_group_sampling_probability(){
    proba_vector.clear();
    proba_vector.resize(this->nsamples,0);
    std::vector<D> q(this->nsamples);
    D sumq=0;
    cout<<"start setting group sampling probablity"<<endl;
    for(L i=0; i<this->nsamples;i++)
    {
      q[i]=Li[i];
      sumq+=q[i];
    }
    D maxq=0;
    L nb=0;
    D tmp=0;
    for(L i=0; i<this->nsamples;i++)
    {
      q[i]=q[i]/sumq*tau;
      if(q[i]>maxq) maxq=q[i];
      if(q[i]>1) {
        nb++; //count the number of elements larger than 1
        tmp+=q[i]-1;
      }
    }
    cout<<"maxq="<<maxq<<endl;
    if(maxq<=1){
      for(L i=0; i<this->nsamples;i++)
      {
        proba_vector[i]=q[i];
      }
    }else{
      for(L i=0; i<this->nsamples;i++)
      {
        if(q[i]>1) q[i]=1;
        else {
          D deltaq=min(1-q[i],tmp);
          tmp=tmp-deltaq;
          q[i]+=deltaq;
        }
        proba_vector[i]=q[i];
      }
    }
  }

  void set_L2(L p_mod, L u){
    if(p_mod==0){ //batch sampling mode
      D tmp=0;
      D st;
      for(L i=0;i<this->nsamples;i++){
        st=Li[i]/tilde_proba_vector[i];
        if (st>tmp) tmp=st;
      }
      L2=tmp/this->nsamples/tau;
      cout<<"L2="<<L2<<"; "<<tilde_proba_vector[0]<<endl;
    }
    else{ //group sampling mode
      D tmp=0;
      D st;
      for(L i=0;i<this->nsamples;i++){
        st=Li[i]/proba_vector[i];
        if(isolated_I[i]==1) st=st-Li[i];
        if(st>tmp) tmp=st;
      }
      L2=tmp/this->nsamples;
    }
    cout<<"set_L2="<<L2<<endl;
  }

  void set_uniform_probability(L p_mod)
  {
    if(p_mod==0) {
      //batch sampling (stochastic process as defined in the paper)
      proba_vector.clear();
      tilde_proba_vector.clear();
      tilde_proba_vector.resize(this->nsamples,1./this->nsamples);
      proba_vector.resize(this->nsamples,1-pow(1-1.0/this->nsamples,tau));
      max_p=1.0/this->nsamples;
      cout<<"uniform="<<tilde_proba_vector[0]<<endl;
    }
    else{ //group sampling
      D pi=(tau*this->c+0.0)/this->nsamples;
      proba_vector.clear();
      proba_vector.resize(this->nsamples,pi);
    }
  }

  void sort_p(){
    vector<pair<D,L> >a;
    for (L i = 0 ;i < this->nsamples ; i++) {
      a.push_back(make_pair(proba_vector[i],i)); // k = value, i = original index
    }
    sort(a.begin(),a.end());
    group_C.clear();
    group_C.resize(this->nsamples);
    isolated_I.clear();
    isolated_I.resize(this->nsamples,0);
    D tmp=0;
    D maxpi=0;
    D previous_i=0;
    index_group_C.clear();
    index_group_C.push_back(0);
    maxp_group_C.clear();
    sump_group_C.clear();

    for (L i = 0 ;i < this->nsamples ; i++){
      group_C[i]=a[this->nsamples-1-i].second;
      D pi=a[this->nsamples-1-i].first;
      if(tmp+pi<=1){
        if(pi>maxpi) maxpi=pi;
        tmp+=pi;
      }
      else{
        index_group_C.push_back(i);
        maxp_group_C.push_back(maxpi);
        sump_group_C.push_back(tmp);
        tmp=pi;
        maxpi=pi;
        if(i-previous_i==1) isolated_I[group_C[previous_i]]=1;
        previous_i=i;
      }
    }
    index_group_C.push_back(this->nsamples);
    maxp_group_C.push_back(maxpi);
    sump_group_C.push_back(tmp);
    nb_groups=sump_group_C.size();
    if(previous_i==this->nsamples-1) isolated_I[group_C[previous_i]]=1;
    cout<<"size of group="<<nb_groups<<endl;
  }

  inline L sampling(L n)
  {
    L i=(floor)(gsl_rng_uniform(rng)*n);
    D y=gsl_rng_uniform(rng);
    while(y*max_p>tilde_proba_vector[i])
    {
      i=(floor)(gsl_rng_uniform(rng)*n);
      y=gsl_rng_uniform(rng);
    }
    return i;
  }

  void set_initial_dual(){
    for(L i=0;i<this->nsamples;i++)
    {
      D res=0;
      for (L k = this->ptr[i]; k < this->ptr[i + 1];k++)
      {
        L j=this->row_idx[k];
        res+=this->A[k]*primal_x[j];
      }
      dual_alpha[i]=gradient_of_phi_i(res,i);
    }
  }

  void update_primal()
  {
    L nb_indices=0;
    //#pragma omp parallel for
    for(L i_t=0;i_t<batch_size;i_t++)
    {
      L i=batch_i[i_t];
      D deltaalphai=batch_deltaalphai[i_t];
      for (L k = this->ptr[i]; k < this->ptr[i + 1];k++)
      {
        L j=this->row_idx[k];
        if (whether_or_not_to_do_prox[j]==0) {
          //#pragma omp critical
          {
          whether_or_not_to_do_prox[j]=1;
          index_to_do_prox[nb_indices]=j;
          nb_indices++;
          deltagk[j]=0;
          }
        }
        deltagk[j]+=lambda_f[i]*deltaalphai*this->A[k]*theta_S[i];
      }
    }
    #pragma omp parallel for
    for(L i_d=0;i_d<nb_indices;i_d++){
      L j=index_to_do_prox[i_d];
      D xj=compute_current_xj_value(baralpha[j]+deltagk[j], j,last_update_i[j],current_nb_iters);
      //update_x(xj,j);
      last_update_i[j]=current_nb_iters;
      whether_or_not_to_do_prox[j]=0;
      index_to_do_prox[i_d]=-1;
    }

  }

  void update_baralpha(L  i, D deltaalphai)
  {
    for (L k = this->ptr[i]; k < this->ptr[i + 1];	k++)
    {
      L j=this->row_idx[k];
      baralpha[j]+=lambda_f[i]*deltaalphai*this->A[k];
    }
  }

  D compute_AiTxk(L i) //This function compute A_i^\top x^k
  {
    D res=0;
    D xj;
    for (L k = this->ptr[i]; k < this->ptr[i + 1];k++)
    {
      L j=this->row_idx[k];
      if(current_nb_iters!=last_update_i[j])
      {
        xj=compute_current_xj_value(baralpha[j], j,last_update_i[j],current_nb_iters);
        //update_x(xj,j);
      }
      res+=this->A[k]*primal_x[j];
      last_update_i[j]=current_nb_iters;
    }
    return res;
  }

  D compute_AiTxk_without_update(L i) //This function compute A_i^\top x^k  without updating x^k
  {
    D res=0;
    D xj;
    //#pragma omp parallel for reduction(+:res)
    for (L k = this->ptr[i]; k < this->ptr[i + 1];k++)
    {
      L j=this->row_idx[k];
      if(current_nb_iters!=last_update_i[j])
      {
        xj=compute_current_xj_value_without_update(baralpha[j], j,last_update_i[j],current_nb_iters);
      }
      else {
        xj=primal_x[j];
      }
      res+=this->A[k]*xj;
    }
    return res;
  }

  void compute_dual_value()
  {
    D res=0;
    //D alphaAx= 0;
    std::vector<D> alpha_tmp(this->nsamples,0);
    std::vector<D> baralpha_tmp(this->nfeatures,0);
    //#pragma omp parallel for schedule(dynamic)
    for(L i=0;i<this->nsamples;i++)
    {
      D aitx=compute_AiTxk_without_update(i);
      alpha_tmp[i]=gradient_of_phi_i(aitx,i);
      //alphaAx+= lambda_f[i]*alpha_tmp[i]*aitx;
      for (L k = this->ptr[i]; k < this->ptr[i + 1];	k++)
      {
        L j=this->row_idx[k];
        baralpha_tmp[j]+=lambda_f[i]*alpha_tmp[i]*this->A[k];
      }
    }

    D scal=feasible_dual(baralpha_tmp);
    //#pragma omp parallel for reduction(-:res)
    for(L i=0;i<this->nsamples;i++)
    {
      res-=value_of_phistar_i(scal*alpha_tmp[i],i)*lambda_f[i];
    }
    res-=value_of_gstar_minus(baralpha_tmp,scal); // should have used value_of_gstar(-baralpha)
    dual_value=res;
  }

  void compute_primal_value()
  {
    D res=0;
    //#pragma omp parallel for reduction(+:res)
    for(L i=0;i<this->nsamples;i++)
    {
      D aitx=compute_AiTxk_without_update(i);
      res+=lambda_f[i]*value_of_phi_i(aitx,i);
    }

    D res2=0;
    //#pragma omp parallel for reduction(+:res2)
    for(L j=0; j<this->nfeatures; j++)
    {
      D xj=compute_current_xj_value_without_update(baralpha[j], j,last_update_i[j],current_nb_iters);
      D gj=value_of_g_j(xj,j);
      res2+=gj;
    }
    primal_value= res+res2;
  }

  void compute_initial_dual_value()
  {

    D res=0;
    D scal=feasible_dual(baralpha);
    #pragma omp parallel for reduction(-:res)
    for(L i=0;i<this->nsamples;i++)
    {
      res-=value_of_phistar_i(scal*dual_alpha[i],i)*lambda_f[i];
    }
    //cout<< "d_1(x)= "<< res<< " d_2(x)= "<< value_of_gstar_minus(baralpha,scal)<< endl;
    res-=value_of_gstar_minus(baralpha,scal);  // should have used value_of_gstar(-baralpha)
    dual_value=res;
  }

  void compute_initial_primal_value()
  {
    D res=0;
    #pragma omp parallel for reduction(+:res)
    for(L j=0; j<this->nfeatures; j++)
    {
      res+=value_of_g_j(primal_x[j],j);
    }
    D res2=0;
    #pragma omp parallel for reduction(+:res2)
    for(L i=0;i<this->nsamples;i++)
    {
      D aitx=0;
      for (L k = this->ptr[i]; k < this->ptr[i + 1];	k++)
      {
        L j=this->row_idx[k];
        aitx+=this->A[k]*primal_x[j];
      }
      res2+=lambda_f[i]*value_of_phi_i(aitx,i);
    }
    //cout<< "p_1(x)= "<< res<< " p_2(x)= "<< res2<< endl;
    primal_value= res+res2;
  }

  void set_p(D scal_p){

    if(scal_p==4.6){
      p=sqrt(this->mu/L1*tau/this->nsamples);
    }
    else if (scal_p==5.6){
      p=this->mu/L1;
    }
    else{
      p=scal_p*tau/(0.0+this->nsamples);
    }
    cout<<"changing probablity: "<<p<<endl;
  }

  void initialize(vector<D> & x0, vector<D> & w0, vector<D> & L_phi, D val_mu, D val_epsilon, L max_nb, L nb_tau, L nb_c, L u, L p_mod, D scal_p)
  {
    cout<<"start initializing"<<" u="<<u<<endl;
    this->distributed_set_up(nb_c);
    if(nb_tau>this->noverc) perror("tau should be less than n over c");
    tau=nb_tau;
    nb_of_iters_per_loop=floor(this->nsamples/(this->c*(tau+0.0)));
    batch_i.clear();
    batch_deltaalphai.clear();
    batch_i.resize(this->nsamples,0);
    batch_deltaalphai.resize(this->nsamples,0);

    deltagk.clear();
    deltagk.resize(this->nfeatures,0);
    index_to_do_prox.clear();
    index_to_do_prox.resize(this->nfeatures,-1);
    whether_or_not_to_do_prox.clear();
    whether_or_not_to_do_prox.resize(this->nfeatures,0);
    gk.clear();
    gk.resize(this->nfeatures,0);

    /**setup parameters**/
    epsilon=val_epsilon;
    max_nb_loops=max_nb;
    this->mu=val_mu;
    cout<<"mu="<<this->mu<<endl;
    lambda_f.clear();
    lambda_f.resize(this->nsamples,0); //L_phi is the Lipschitz constant of the function \phi_i.
    for(int i=0;i< this->nsamples; i++){
      lambda_f[i]= L_phi[i]/this->nsamples;
    }

    set_rng();
    set_auxiliary_v();

    lambda_1=get_lambda1();
    lambda_2=get_lambda2();

    set_Li_Lf();

    running_time=0;
    /**setup probability**/
    if(p_mod==0)   // batch sampling (stochastic process as defined in the paper)
    {
      if(u==0)   // uniform sampling
      {
        set_uniform_probability(p_mod);
        this->uniform="uniform";
      }
      else{
        set_optimal_probability();
        this->uniform="nonuniform";
      }
    }
    else{  //group sampling
      if(u==0){
        set_uniform_probability(p_mod);
        this->uniform="uniform";
      }else{
        this->set_group_sampling_probability();
        this->uniform="nonuniform";
      }
      sort_p();
      cout<<"sort p"<<endl;
    }

    set_L2(p_mod,u);
    scaler=scal_p;

    set_L1(p_mod);
    cout<<"set L1="<<L1<<endl;
    set_p(scal_p);
    set_thetaS(p_mod);
    cout<<"setthetaS"<<endl;

    primal_x=w0;
    dual_alpha=x0;
    set_stepsize();
    cout<<"set stepsize"<<endl;

    set_initial_dual();

    baralpha.clear();
    baralpha.resize(this->nfeatures,0);
    last_update_i.clear();
    last_update_i.resize(this->nfeatures,0);
    for(L i=0;i<this->nsamples;i++)
    {
      update_baralpha(i, dual_alpha[i]);
    }

    current_nb_iters=0;
    D start_time_initial=std::clock();
    compute_initial_primal_value();
    compute_initial_dual_value();
    D end_time_initial=std::clock();
    cout<<"computation of initial primal dual value time="<<(end_time_initial-start_time_initial) / (double) CLOCKS_PER_SEC<<endl;

    cout<<"primal value="<<primal_value<<endl;
    cout<<"dual_value="<<dual_value<<endl;

    cout<<"Initialization is finished!"<<endl;
  }

  void result_record()
  {
    cout<<primal_value<<endl;
    cout<<dual_value<<endl;
    ofstream result;
    result.open("results/x0.dat");
    for (L j=0;j<this->nfeatures;j++){
      result <<setprecision(20)<< primal_x[j] << " ";
      //cout<< primal_x[j] << " ";
    }
    result << endl;
    ofstream resulta;
    resulta.open("results/alpha0.dat");
    for (L i=0;i<this->nsamples;i++)
    resulta << setprecision(20)<<dual_alpha[i] << " ";
    resulta << endl;
    result.close();
  }

  void read_w0(){
    ifstream recordw0("results/x0.dat");
    for(L j=0;j<this->nfeatures;j++)
    {
      recordw0>>primal_x[j];
    }
    recordw0.close();

    ifstream recordalpha0("results/alpha0.dat");
    for(L i=0;i<this->nsamples;i++)
    {
      recordalpha0>>dual_alpha[i];
    }
    recordalpha0.close();
  }

  void compute_and_record_result()
  {
    if(nb_loops%print_every_N==0){
      D start = std::clock();
      compute_primal_value();
      D endt=std::clock();
      cout<<"computation of primal value time="<<(endt-start) / (double) CLOCKS_PER_SEC<<endl;
      compute_dual_value();
      delta=primal_value-dual_value;
      cout<<setprecision(16)<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<";;;  "<<running_time<<"; primal dual gap="<<delta<<" primal value="<<primal_value<<"; dual value="<<dual_value<<endl;
      samp<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<" "<<delta<<" "<<running_time<<" "<<primal_value<<endl;
    }

  }

  void batch_sampling()
  {
    batch_size=0;
    L i;
    for(L it_t=0;it_t<tau;it_t++)
    {
      i=sampling(this->nsamples);
      batch_i[batch_size]=i;
      batch_size++;
    }
  }

  void set_thetaS(L p_mod){
    theta_S.clear();
    theta_S.resize(this->nsamples,0);
    if(p_mod==0){
      for(L i=0;i<this->nsamples;i++)
      theta_S[i]=1.0/(tau*tilde_proba_vector[i]);
    }
    else{
      for(L i=0;i<this->nsamples;i++)
      theta_S[i]=1.0/proba_vector[i];
    }
  }

  void group_sampling(){
    batch_size=0;
    std::vector<L> sampled_groups;
    for(L i=0;i<nb_groups;i++){
      D y=gsl_rng_uniform(rng);
      if(y<sump_group_C[i])
      {
        sampled_groups.push_back(i);
        batch_size++;
      }
    }
    for(L t=0;t<batch_size;t++){
      L group_i=sampled_groups[t];
      L s1=index_group_C[group_i];
      L s2=index_group_C[group_i+1];
      L i=s1+(floor)(gsl_rng_uniform(rng)*(s2-s1));
      D y=gsl_rng_uniform(rng);
      D maxpi=maxp_group_C[group_i];
      while(y*maxpi>proba_vector[group_C[i]])
      {
        i=s1+(floor)(gsl_rng_uniform(rng)*(s2-s1));
        y=gsl_rng_uniform(rng);
      }
      batch_i[t]=group_C[i];
    }
  }

  void loopless(vector<D> & x0, vector<D> & w0, string filename, vector<D> & L_phi, D val_mu, D val_epsilon, L max_nb, L nb_tau, L nb_c, L u, L p_mod, D scal_p)
  {
    initialize(x0,  w0,  L_phi, val_mu, val_epsilon,  max_nb,  nb_tau, nb_c,  u, p_mod, scal_p);
    string sampname="results/L_"+filename+uniform;
    if(p_mod==0) sampname=sampname+"_batch";
    else sampname=sampname+"_group";
    string scal_str;
    stringstream  scal_convert;
    scal_convert<<scal_p;
    scal_str=scal_convert.str();
    sampname+=scal_str;
    cout<<"running Loopless"<<" ; "<<sampname<<endl;
    samp.open(sampname.c_str());

    delta=primal_value-dual_value;
    nb_loops=0;
    nb_iters=0;

    cout<<setprecision(9)<<"initial: "<<" delta: "<<delta<<" primal: "<<primal_value<<"; dual: "<<dual_value<<"epsilon="<<epsilon<<endl;
    samp<<((0.0+nb_iters)*this->c*tau/(this->nsamples))<<" "<<delta<<running_time<<" "<<primal_value<<endl;
    srand48(27432042);
    //srand(time(NULL));
    printf("There are %d threads\n",omp_get_num_threads());

    D start;
    while(delta>=epsilon && nb_loops<max_nb_loops)
    {
      compute_and_record_result();

      if(delta<epsilon) break;
      start = std::clock();

      nb_loops++;

      //cout<<"before the loop time elapsed="<<duration<<endl;
      start = std::clock();

      //cout<<"nb_of_iters_per_loop="<<nb_of_iters_per_loop<<endl;
      for(L it=0;it<nb_of_iters_per_loop;it++)
      {
        if(p_mod==0)
        batch_sampling();
        else
        group_sampling();

        //start = std::clock();
        start=omp_get_wtime();
        #pragma omp parallel
        {
          #pragma omp for
          for(L i_t=0;i_t<batch_size;i_t++)
          {
          L i=batch_i[i_t];
          D aitg=compute_AiTxk(i);
          D deltaalphai=(gradient_of_phi_i(aitg,i)-dual_alpha[i]);
          batch_deltaalphai[i_t]=deltaalphai;
         }
      }
        //cout<<omp_get_wtime()-start<<endl;
        current_nb_iters++;
        update_primal();
        update_baralpha();
        nb_iters++;
        D end_time = omp_get_wtime();
        //running_time+= ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
        running_time+=end_time-start;

      }

      //cout<<"after the loop time elapsed="<<duration<<endl;
    }
    if(delta>=epsilon) {compute_primal_value();
      compute_dual_value();
      delta=primal_value-dual_value;
      cout<<setprecision(9)<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<";  "<<running_time<<"; primal dual gap="<<delta<<" primal value="<<primal_value<<"; dual value="<<dual_value<<endl;
      samp<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<" "<<delta<<" "<<running_time<<" "<<primal_value<<endl;
    }
    cout<<setprecision(9)<<"nb_loops: "<<nb_loops<<"nb_iters"<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<";  delta: "<<delta<<" primal: "<<primal_value<<"; dual: "<<dual_value<<endl;
    samp.close();
  }

  void loopless2(vector<D> & x0, vector<D> & w0, string filename, vector<D> & L_phi, D val_mu, D val_epsilon, L max_nb, L nb_tau, L nb_c, L u, L p_mod, D scal_p)
  {
    initialize(x0,  w0,  L_phi, val_mu, val_epsilon,  max_nb,  nb_tau, nb_c,  u, p_mod, scal_p);
    string sampname="results/L_"+filename+uniform;
    if(p_mod==0) sampname=sampname+"_batch";
    else sampname=sampname+"_group";
    string scal_str;
    stringstream  scal_convert;
    scal_convert<<scal_p;
    scal_str=scal_convert.str();
    sampname+=scal_str;
    cout<<"running Loopless"<<" ; "<<sampname<<endl;
    samp.open(sampname.c_str());

    delta=primal_value-dual_value;
    nb_loops=0;
    nb_iters=0;

    cout<<setprecision(9)<<"initial: "<<" delta: "<<delta<<" primal: "<<primal_value<<"; dual: "<<dual_value<<endl;
    samp<<((0.0+nb_iters)*this->c*tau/(this->nsamples))<<" "<<delta<<running_time<<" "<<primal_value<<endl;
    //srand48(27432042);
    srand(time(NULL));

    D start;
    while(nb_loops<max_nb_loops)
    {
      compute_and_record_result();

      start = std::clock();

      nb_loops++;

      //cout<<"before the loop time elapsed="<<duration<<endl;
      start = std::clock();
      //cout<<"nb_of_iters_per_loop="<<nb_of_iters_per_loop<<endl;
      for(L it=0;it<nb_of_iters_per_loop;it++)
      {
        if(p_mod==0)
        batch_sampling();
        else
        group_sampling();

        start = std::clock();
        for(L i_t=0;i_t<batch_size;i_t++)
        {
          L i=batch_i[i_t];
          D aitg=compute_AiTxk(i);
          D deltaalphai=(gradient_of_phi_i(aitg,i)-dual_alpha[i]);
          batch_deltaalphai[i_t]=deltaalphai;
        }
        current_nb_iters++;
        update_primal();
        update_baralpha();
        nb_iters++;
        running_time+= ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
      }
      //cout<<"after the loop time elapsed="<<duration<<endl;
    }
    cout<<setprecision(9)<<"nb_loops: "<<floor(((0.0+nb_iters)*this->c*tau/(this->nsamples)))<<";  delta: "<<delta<<" primal: "<<primal_value<<"; dual: "<<dual_value<<endl;
    samp.close();
  }

};

#endif /* MIN_SMOOTH_CONVEX_H */


text:
#pragma once

#include <assert.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

#include "defines.h"
#include "matrix_vector_ops.h"

#define DEFAULT_DIFF_ERROR_TOLERANCE 5e-5
#define DEFAULT_DIFF_STEPSIZE 1e-2

#define MAX_DIFF_CORRECTION_ATTEMPTS 100

#define RUNGE_KUTTA_ORDER 5

typedef int (*DynamicsFunction) (gsl_vector *dy,
                                 double t,
                                 const gsl_vector *y,
                                 const gsl_vector *u);

// Runge-Kutte-Dormand-Prince embedded method modified for controlled systems
//TODO implement adaptive step size
int rungeKutteStep(DynamicsFunction dyn,
                   double t,
                   gsl_vector *y_next,
                   gsl_vector *rk_e_next,
                   const VehicleState *vehicle,
                   const Controller *controller,
                   double h);

int rungeKutteAdaptiveStep(DynamicsFunction dyn,
                           double t,
                           gsl_vector *y_next,
                           gsl_vector *rk_e_next,
                           const VehicleState *vehicle,
                           const Controller *controller,
                           double *stepSize,
                           double tolerance);

//TODO implement implicit RK method as well for stiff equations

max_stars_repo_name: umd-agrc/SimpleControlSim

max_issues_repo_name: umd-agrc/SimpleControlSim

max_forks_repo_name: umd-agrc/SimpleControlSim


text:
/* Implementation for next reaction method */

#include <math.h>
#include <stdint.h>
#include <time.h>
#include <stdio.h>
#include <stdint.h>
#include <float.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "ssa.h"
#include "nrm.h"
#include "priorityq.h"
#include "depgraph.h"

const gsl_rng_type *SSARNGT;
gsl_rng *SSARNG;

void update_reaction(REACTION *r, SYSTEM *s, double t)
{
	double newp = ssa_h(s->R[r->number], s->x, s->n) * s->k[r->number];
	double oldp = s->propensities[r->number];

	if (newp > 0.0) {
		if (oldp > 0.0) {
			pq_update(s->pq, r, (oldp / newp) * (r->tau - t) + t);
		} else {
			r->tau = gsl_ran_gamma(SSARNG, s->steps[r->number], 1 / newp) + t;
			pq_insert(s->pq, r);
		}
	} else if (oldp > 0.0) {
		pq_delete(s->pq, r);
	}

	s->propensities[r->number] = newp;
}

double ssa_nrmstep(SYSTEM *s)
{
	REACTION *r = pq_min(s->pq);
        ssa_doreaction(s->R[r->number], s->P[r->number], s->x, s->n);
	double t = r->tau;

	s->propensities[r->number] = ssa_h(s->R[r->number], s->x, s->n) * s->k[r->number];
	if (s->propensities[r->number] > 0.0)
		pq_update(s->pq, r, gsl_ran_gamma(SSARNG, s->steps[r->number], 1 / s->propensities[r->number]) + t);
	else
		pq_delete(s->pq, r);

	LLIST *head = r->affects;
	while (head != NULL) {
		REACTION *update = (REACTION *) head->data;
		update_reaction(update, s, t);
		head = head->next;
	}

	return t;
}

REACTION **setup_reactions(SYSTEM *s)
{
	REACTION **reactions = malloc(sizeof(REACTION*) * s->m);
	char **adjmat = adjacencymatrix(s->R, s->P, s->n, s->m);

	INDEX i, j;
	for (i = 0; i < s->m; i++) {
		REACTION *reaction = reaction_make(i);
		reactions[i] = reaction;
	}

	for (i = 0; i < s->m; i++) {
		REACTION *r = reactions[i];
		for (j = 0; j < s->m; j++)
			if (i != j && adjmat[i][j] != 0)
				r->affects = llist_push(r->affects, reactions[j]);
	}

        for (i = 0; i < s->m; i++) {
		REACTION *r = reactions[i];
		for (j = 0; j < s->m; j++) {
			if (s->creates[i][j] != 0)
				r->creates = llist_push(r->creates, reactions[j]);
                        if (s->destroys[i][j] != 0)
                                r->destroys = llist_push(r->destroys, reactions[j]);
                }
	}

	for (i = 0; i < s->m; i++)
		free(adjmat[i]);

	free(adjmat);
	return reactions;
}

void ssa_nrm(INDEX **R, INDEX **P, INDEX n, INDEX m, double *k, COUNT *x, COUNT *steps, INDEX **creates, INDEX **destroys, double T)
{
	gsl_rng_env_setup();
	SSARNGT = gsl_rng_default;
	SSARNG = gsl_rng_alloc(SSARNGT);
	gsl_rng_set(SSARNG, time(NULL));

        SYSTEM *s = malloc(sizeof(SYSTEM));
	s->m = m;
        s->n = n;
        s->R = R;
        s->P = P;
        s->creates = creates;
        s->destroys = destroys;

	REACTION **reactions = setup_reactions(s);
	double *propensities = malloc(sizeof(double) * m);
	PQ *pq = pq_make();

        s->pq = pq;
        s->x = x;
        s->steps = steps;
        s->k = k;
        s->propensities = propensities;

	INDEX i;
	for (i = 0; i < m; i++) {
		REACTION *r = reactions[i];
		propensities[i] = ssa_h(R[i], x, n) * k[i];
		if (propensities[i] > 0.0) {
			r->tau = gsl_ran_gamma(SSARNG, steps[i], 1 / propensities[i]);
			pq_insert(pq, r);
		}
	}

	ssa_printstate(0.0, x, n);

	while(!pq_isempty(pq) && pq_min(pq)->tau < T)
		ssa_printstate(ssa_nrmstep(s), x, n);

	pq_free(pq);
	gsl_rng_free(SSARNG);
}

max_stars_repo_name: lgrozinger/ssapy

max_issues_repo_name: lgrozinger/ssapy

max_forks_repo_name: lgrozinger/ssapy