text:
#include <math.h>
#include "2d_array.h"
#include "const.h"
#include "misc.h"
#include "utilities.h"
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_fit.h>

/******************************************************************************
MODULE:  rmse_from_square_root_mean

PURPOSE:  simulate matlab calculate rmse from square root mean

RETURN VALUE:
Type = void
Value           Description
-----           -----------

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
5/28/2019   Su Ye            Original Development

******************************************************************************/
void rmse_from_square_root_mean
(
    float **array,      /* I: input array */
    float fit_cft,      /* I: input fit_cft value */
    int dim1_index,     /* I: dimension 1 index in input array */
    int dim2_len,       /* I: dimension 2 length */
    float *rmse         /* O: output rmse */
)
{
    int i;
    float sum = 0.0;

    for (i = 0; i < dim2_len; i++)
    {
         sum += (array[dim1_index][i] - fit_cft) *
                (array[dim1_index][i] - fit_cft);
    }
    *rmse = sqrt(sum / dim2_len);
}

/******************************************************************************
MODULE:  dofit

PURPOSE: Declare data type and allocate memory and do multiple linear robust
         fit used for auto_robust_fit

RETURN VALUE: None

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
5/28/2019   Su Ye            Original Development

NOTES:
******************************************************************************/
void dofit(const gsl_multifit_robust_type *T,
      const gsl_matrix *X, const gsl_vector *y,
      gsl_vector *c, gsl_matrix *cov)
{
  gsl_multifit_robust_workspace * work
    = gsl_multifit_robust_alloc (T, X->size1, X->size2);
  gsl_multifit_robust (X, y, c, cov, work);
  gsl_multifit_robust_free (work);
  work = NULL;
}

void dofit_linear(const gsl_matrix *X, const gsl_vector *y,
      gsl_vector *c, gsl_matrix *cov)
{
    double chisq;
    gsl_multifit_linear_workspace * work
        = gsl_multifit_linear_alloc ( X->size1, X->size2);
    gsl_multifit_linear(X, y, c, cov, &chisq,work);
    gsl_multifit_linear_free (work);
    // work = NULL; // SY 03242019
}

/******************************************************************************
MODULE:  auto_robust_fit

PURPOSE:  Robust fit for one band

RETURN VALUE: None

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
5/28/2019   Su Ye            Original Development

NOTES:
******************************************************************************/
void auto_robust_fit
(
    float **clrx,
    float **clry,
    int nums,
    int start,
    int band_index,
    float *coefs
)
{
    int i, j;
    const int p = 2; /* linear fit */
    gsl_matrix *x, *cov;
    gsl_vector *y, *c;

    /******************************************************************/
    /*                                                                */
    /* Defines the inputs/outputs for robust fitting                  */
    /*                                                                */
    /******************************************************************/

    x = gsl_matrix_alloc (nums, p);
    y = gsl_vector_alloc (nums);

    c = gsl_vector_alloc (p);
    cov = gsl_matrix_alloc (p, p);

    /******************************************************************/
    /*                                                                */
    /* construct design matrix x for linear fit                       */
    /*                                                                */
    /******************************************************************/

    for (i = 0; i < nums; ++i)
    {
        for (j = 0; j < p; j++)
        {
            if (j == 0)
            {
                gsl_matrix_set (x, i, j, 1.0);
            }
            else
            {
                gsl_matrix_set (x, i, j, (double)clrx[i][j-1]);
            }
        }
        gsl_vector_set(y,i,(double)clry[band_index][i+start]);
    }

    /******************************************************************/
    /*                                                                */
    /* perform robust fit                                             */
    /*                                                                */
    /******************************************************************/

    dofit(gsl_multifit_robust_bisquare, x, y, c, cov);

    for (j = 0; j < (int)c->size; j++)
    {
        coefs[j] = gsl_vector_get(c, j);
    }

    /******************************************************************/
    /*                                                                */
    /* Free the memories                                              */
    /*                                                                */
    /******************************************************************/

    gsl_matrix_free (x);
    x = NULL;
    gsl_vector_free (y);
    y = NULL;
    gsl_vector_free (c);
    c = NULL;
    gsl_matrix_free (cov);
    cov = NULL;
}

/******************************************************************************
MODULE:  auto_mask

PURPOSE:  Multitemporal cloud, cloud shadow, & snow masks (global version)

RETURN VALUE:
Type = int
ERROR error out due to memory allocation
SUCCESS no error encounted

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
05282019    Su Ye

NOTES:
******************************************************************************/
int auto_mask
(
    int *clrx,
    float **clry,
    int start,
    int end,
    float years,
    float t_b1,
    float t_b2,
    float n_t,
    int *bl_ids
)
{
    char FUNC_NAME[] = "auto_mask";
    int i;
    float **x;
    float pred_b2, pred_b4;
    int nums;
    float coefs[ROBUST_COEFFS];
    float coefs2[ROBUST_COEFFS];

    nums = end - start + 1;
    /* Allocate memory */
    x = (float **)allocate_2d_array(nums, ROBUST_COEFFS - 1, sizeof(float));
    if (x == NULL)
    {
        RETURN_ERROR("ERROR allocating x memory", FUNC_NAME, ERROR);
    }

    for (i = 0; i < nums; i++)
    {
        x[i][0] = (float)clrx[i+start];
    }

    /******************************************************************/
    /*                                                                */
    /* Do robust fitting for band 2 */
    /*                                                                */
    /******************************************************************/

    auto_robust_fit(x, clry, nums, start, 1, coefs);

    /******************************************************************/
    /*                                                                */
    /* Do robust fitting for band 4 */
    /*                                                                */
    /******************************************************************/

    auto_robust_fit(x, clry, nums, start, 3, coefs2);

    /******************************************************************/
    /*                                                                */
    /* predict band 2 and band 4 refs, bl_ids value of 0 is clear and   */
    /* 1 otherwise                                                    */
    /*                                                                */
    /******************************************************************/

    for (i = 0; i < nums; i++)
    {
        pred_b2 = coefs[0] + coefs[1] * (float)clrx[i+start];
        pred_b4 = coefs2[0] + coefs2[1] * (float)clrx[i+start];
        if (((clry[1][i+start]-pred_b2) > (n_t * t_b1)) ||
            ((clry[3][i+start]-pred_b4) < -(n_t * t_b2)))
        {
            bl_ids[i] = 1;
        }
        else
        {
            bl_ids[i] = 0;
        }
    }

    /* Free allocated memory */
    if (free_2d_array ((void **) x) != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: x\n", FUNC_NAME, ERROR);
    }

    return (SUCCESS);
}

/******************************************************************************
MODULE:  adjust_median_variogram

PURPOSE:  calculate absolute variogram for auto_mask

RETURN VALUE:
Type = void
Value           Description
-----           -----------

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
28/05/2019   Su Ye         Original Development

NOTES:
******************************************************************************/
int adjust_median_variogram
(
    int *clrx,                  /* I: dates                                          */
    float **array,              /* I: input array                                    */
    int dim1_len,               /* I: dimension 1 length in input array              */
    int dim2_start,             /* I: dimension 2 start index                        */
    int dim2_end,               /* I: dimension 2 end index                          */
    float *date_vario,          /* O: outputted median variogran for dates           */
    float *max_neighdate_diff,  /* O: maximum difference for two neighbor times       */
    float *output_array        /* O: output array                                   */
)
{
    int i, j;           /* loop indecies                                     */
    float *var;         /* pointer for allocation variable memory            */
    int dim2_len = dim2_end - dim2_start + 1; /* perhaps should get defined  */
    char FUNC_NAME[] = "adjust_median_variogram"; /* for error messages             */
    int max_freq;
    int m;

//    for (i = 0; i < dim2_len; i++)
//    {
//        for (j = 0; j < TOTAL_IMAGE_BANDS; j++)
//        {
//           printf("%f\n", (float)array[j][i]);
//        }

//    }

    if (dim2_len == 1)
    {
        for (i = 0; i < dim1_len; i++)
        {
            output_array[i] = array[i][dim2_start];
            *date_vario = clrx[dim2_start];
            return (SUCCESS);
        }
    }

    var = malloc((dim2_len-1) * sizeof(float));
    if (var == NULL)
    {
        RETURN_ERROR ("Allocating var memory", FUNC_NAME, ERROR);
    }

    for (j = dim2_start; j < dim2_end; j++)
    {
        var[j] = abs(clrx[j+1] - clrx[j]);
    }
    quick_sort_float(var, dim2_start, dim2_end-1);
    m = (dim2_len-1) / 2;
    if ((dim2_len-1) % 2 == 0)
    {

        *date_vario = (var[m-1] + var[m]) / 2.0;
    }
    else
        *date_vario = var[m];

    *max_neighdate_diff = var[dim2_len-2];

    for (i = 0; i < dim1_len; i++)
    {
        for (j = dim2_start; j < dim2_end; j++)
        {
            var[j] = abs(array[i][j+1] - array[i][j]);
            //printf("%d var for band %d: %f\n", j, i+1, (float)var[j]);

        }
        quick_sort_float(var, dim2_start, dim2_end-1);
//            for (j = 0; j < dim2_end; j++)
//            {
//               printf("%f\n", var[j]);
//            }
        m = (dim2_len-1) / 2;
        if ((dim2_len-1) % 2 == 0)
        {
            //printf("%f\n", var[m-1]);
           //printf("%f\n", var[m]);
            output_array[i] = (var[m-1] + var[m]) / 2.0;
        }
        else
            output_array[i] = var[m];

    }

    free(var);

    return (SUCCESS);
}

/******************************************************************************
 *
MODULE:  single_median_variogram
PURPOSE:  calculate absolute variogram for a single band

RETURN VALUE:
Type = void
Value           Description
-----           -----------

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
22/06/2019   Su Ye         Original Development

NOTES:
******************************************************************************/
int single_median_variogram
(
    short int *array,              /* I: input array                                    */
    int i_start,
    int i_end,
    short int *variogram,
    short int *mediam_value          /* O: outputted mediam value                         */
)
{
    int j, m;           /* loop indecies                                     */
    char FUNC_NAME[] = "single_median_variogram"; /* for error messages             */
    short int *var;         /* pointer for allocation variable memory            */
    int obs_num;
    int var_count = 0;
    short int *array_cpy;

    obs_num = i_end - i_start + 1;

//    for (j = 0; j < obs_num; j++)
//    {
//        printf("%i\n", array[j]);

//    }

    if (obs_num == 1)
    {
        *variogram = 0;
        *mediam_value = array[0];
        return (SUCCESS);

    }

    var = malloc((obs_num-1) * sizeof(short int));
    if (var == NULL)
    {
        RETURN_ERROR ("Allocating var memory", FUNC_NAME, ERROR);
    }

    array_cpy = malloc(obs_num * sizeof(short int));
    if (array_cpy == NULL)
    {
        RETURN_ERROR ("Allocating array_cpy memory", FUNC_NAME, ERROR);
    }

    for (j = i_start; j < i_end; j++)
    {
        var[j] = (short int)abs(array[j+1] - array[j]);
        //printf("%d var for band %d: %f\n", j, i+1, (float)var[j]);

    }

    for (j = i_start; j < i_end+1; j++)
    {
        array_cpy[j] = array[j];
        //printf("%d var for band %d: %f\n", j, i+1, (float)var[j]);

    }
    quick_sort_shortint(var, 0, obs_num-2);
//            for (j = 0; j < obs_num-1; j++)
//            {
//               printf("%d\n", var[j]);
//            }

    /* compute variogram */
    m = (obs_num-1) / 2;
    if ((obs_num-1) % 2 == 0)
    {
        *variogram = (short int) ((var[m-1] + var[m]) / 2.0);
    }
    else
        *variogram = var[m];

    /* compute mediam value */
    quick_sort_shortint(array_cpy, i_start, i_end);
//    for (j = 0; j < obs_num; j++)
//    {
//       printf("%d\n", array_cpy[j]);
//    }

    m = obs_num / 2;
    if (obs_num % 2 == 0)
    {
        *mediam_value = (short int)((array_cpy[m-1] + array_cpy[m]) / 2.0);
    }
    else
        *mediam_value = array_cpy[m];

    free(var);
    free(array_cpy);

    return (SUCCESS);
}

/******************************************************************************
 *
MODULE:  single_median_variogram
PURPOSE:  calculate absolute variogram for a single band

RETURN VALUE:
Type = void
Value           Description
-----           -----------

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
22/06/2019   Su Ye         Original Development

NOTES:
******************************************************************************/
int single_median_quantile
(
    short int *array,              /* I: input array                                    */
    int i_start,
    int i_end,
    short int *quantile,          /* O: outputted quantile value                         */
    short int *mediam_value
)
{
    int j, m;           /* loop indecies                                     */
    char FUNC_NAME[] = "single_median_quantile"; /* for error messages             */
    short int *var;         /* pointer for allocation variable memory            */
    int obs_num;
    int var_count = 0;
    short int *array_cpy;

    obs_num = i_end - i_start + 1;

//    for (j = 0; j < obs_num; j++)
//    {
//        printf("%i\n", array[j]);

//    }

    if (obs_num == 1)
    {
        *mediam_value = array[0];
        *quantile = array[0];
        return (SUCCESS);

    }

    var = malloc((obs_num-1) * sizeof(short int));
    if (var == NULL)
    {
        RETURN_ERROR ("Allocating var memory", FUNC_NAME, ERROR);
    }

    array_cpy = malloc(obs_num * sizeof(short int));
    if (array_cpy == NULL)
    {
        RETURN_ERROR ("Allocating array_cpy memory", FUNC_NAME, ERROR);
    }

    for (j = i_start; j < i_end+1; j++)
    {
        array_cpy[j] = array[j];
        //printf("%d var for band %d: %f\n", j, i+1, (float)var[j]);

    }

    /* compute mediam value */
    quick_sort_shortint(array_cpy, i_start, i_end);
//    for (j = 0; j < obs_num; j++)
//    {
//       printf("%d\n", array_cpy[j]);
//    }

    m = obs_num / 2;
    if (obs_num % 2 == 0)
    {
        *mediam_value = (short int)((array_cpy[m-1] + array_cpy[m]) / 2.0);
    }
    else
        *mediam_value = array_cpy[m];

    *quantile = array_cpy[obs_num / 4];

    free(var);
    free(array_cpy);

    return (SUCCESS);
}
/******************************************************************************
 *
MODULE:  single_mean_rmse
PURPOSE:  calculate absolute variogram for a single band

RETURN VALUE:
Type = void
Value           Description
-----           -----------

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
22/06/2019   Su Ye         Original Development

NOTES:
******************************************************************************/
int single_mean_rmse
(
    short int *array,              /* I: input array                                    */
    int i_start,
    int i_end,
    float *rmse,
    float *mean          /* O: outputted mediam value                         */
)
{
    int j, m;           /* loop indecies                                     */
    char FUNC_NAME[] = "single_mean_rse"; /* for error messages             */
    short int *var;         /* pointer for allocation variable memory            */
    int obs_num;
    int sum = 0;
    int rmse_sum = 0;

    obs_num = i_end - i_start + 1;

//    for (j = 0; j < obs_num; j++)
//    {
//        printf("%i\n", array[j]);

//    }

    if (obs_num == 1)
    {
        *rmse = 0;
        *mean = array[0];
        return (SUCCESS);

    }

    var = malloc((obs_num-1) * sizeof(short int));
    if (var == NULL)
    {
        RETURN_ERROR ("Allocating var memory", FUNC_NAME, ERROR);
    }

    for (j = i_start; j < i_end + 1; j++)
    {
        sum = sum + (int)array[j];
    }

    *mean = (float)((float)sum / (float)obs_num);

    for (j = i_start; j < i_end + 1; j++)
    {
        rmse_sum = rmse_sum + (array[j] - *mean) * (array[j] - *mean);
    }

    *rmse = (float)sqrt((double)rmse_sum/(double)obs_num);
    return (SUCCESS);
}

/******************************************************************************
MODULE:  linear_fit_centerdate

PURPOSE:  Robust fit for one band

RETURN VALUE: None

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
5/31/2019   Su Ye            Original Development

NOTES:
******************************************************************************/
void linear_fit_centerdate
(
    int *clrx,
    float **clry,
    int nums,
    int start,
    int center_date,
    int i_col,
    short int **composites,
    int bweighted,
    float *C0,
    float *C1
)
{
    char FUNC_NAME[] = "linear_fit_centerdate";
    double c0;
    double c1;
    double cov00;
    double cov01;
    double cov11;
    double sumsq;
    int i;
    double* x;
    double* y_b1;
    double* y_b2;
    double* y_b3;
    double* y_b4;
    double* w;
    float coefs[ROBUST_COEFFS];
    float** x_t;

    x = (double*)malloc(nums * sizeof(double));
    if (x == NULL)
    {
        RETURN_ERROR("ERROR allocating x memory", FUNC_NAME, FAILURE);
    }

    y_b1 = (double*)malloc(nums * sizeof(double));
    if (y_b1 == NULL)
    {
        RETURN_ERROR("ERROR allocating y_b1 memory", FUNC_NAME, FAILURE);
    }

    y_b2 = (double*)malloc(nums * sizeof(double));
    if (y_b2 == NULL)
    {
        RETURN_ERROR("ERROR allocating y_b2 memory", FUNC_NAME, FAILURE);
    }

    y_b3 = (double*)malloc(nums * sizeof(double));
    if (y_b3 == NULL)
    {
        RETURN_ERROR("ERROR allocating y_b3 memory", FUNC_NAME, FAILURE);
    }

    y_b4 = (double*)malloc(nums * sizeof(double));
    if (y_b4 == NULL)
    {
        RETURN_ERROR("ERROR allocating y_b4 memory", FUNC_NAME, FAILURE);
    }

    w = (double*)malloc(nums * sizeof(double));
    if (w == NULL)
    {
        RETURN_ERROR("ERROR allocating y_b4 memory", FUNC_NAME, FAILURE);
    }

    x_t = (float **)allocate_2d_array(nums, ROBUST_COEFFS - 1, sizeof(float));
    if (x_t == NULL)
    {
        RETURN_ERROR("ERROR allocating x_t memory", FUNC_NAME, ERROR);
    }

    for (i = 0; i < nums; i++)
    {
        x_t[i][0] = (float)clrx[i+start];
    }

    /******************************************************************/
    /*                                                                */
    /* Defines the inputs/outputs for robust fitting                  */
    /*                                                                */
    /******************************************************************/

    /******************************************************************/
    /*                                                                */
    /* construct design matrix x for linear fit                       */
    /*                                                                */
    /******************************************************************/

    for (i = 0; i < nums; ++i)
    {
        x[i] = (double)clrx[i];
        y_b1[i] = (double)clry[0][i];
        y_b2[i] = (double)clry[1][i];
        y_b3[i] = (double)clry[2][i];
        y_b4[i] = (double)clry[3][i];

    }

    /******************************************************************/
    /*                                                                */
    /* perform robust fit                                             */
    /*                                                                */
    /******************************************************************/

    if (bweighted == TRUE)
    {
        /* create hot weights */
        for (i = 0; i < nums; ++i)
        {
             w[i] = (double)1/(abs(y_b1[i] - 0.5 * y_b3[i]));
        }

        gsl_fit_wlinear (x, 1, w, 1, y_b1, 1, nums, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
        composites[0][i_col] = (short int)(c0 + c1 * center_date);
        C0[0] = (float)c0;
        C1[0] = (float)c1;

        gsl_fit_wlinear (x, 1, w, 1, y_b2, 1, nums, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
        composites[1][i_col] = (short int)(c0 + c1 * center_date);
        C0[1] = (float)c0;
        C1[1] = (float)c1;

        gsl_fit_wlinear (x, 1, w, 1, y_b3, 1, nums, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
        composites[2][i_col] = (short int)(c0 + c1 * center_date);
        C0[2] = (float)c0;
        C1[2] = (float)c1;

        gsl_fit_wlinear (x, 1, w, 1, y_b4, 1, nums, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
        composites[3][i_col] = (short int)(c0 + c1 * center_date);
        C0[3] = (float)c0;
        C1[3] = (float)c1;
    }
    else
    {
        /* OLS*/
//        gsl_fit_linear (x, 1, y_b1, 1, nums, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
//        composites[0][i_col] = (short int)(c0 + c1 * center_date);
//        C0[0] = (float)c0;
//        C1[0] = (float)c1;

//        gsl_fit_linear (x, 1, y_b2, 1, nums, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
//        composites[1][i_col] = (short int)(c0 + c1 * center_date);
//        C0[1] = (float)c0;
//        C1[1] = (float)c1;

//        gsl_fit_linear (x, 1, y_b3, 1, nums, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
//        composites[2][i_col] = (short int)(c0 + c1 * center_date);
//        C0[2] = (float)c0;
//        C1[2] = (float)c1;

//        gsl_fit_linear (x, 1, y_b4, 1, nums, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
//        composites[3][i_col] = (short int)(c0 + c1 * center_date);
//        C0[3] = (float)c0;
//        C1[3] = (float)c1;
        /* robust regression */
        for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        {
            auto_robust_fit(x_t, clry, nums, start, i, coefs);
            composites[i][i_col] = (short int)(coefs[0] + coefs[1] * center_date);
            C0[i] = (float)coefs[0];
            C1[i] = (float)coefs[1];

        }
    }

    /******************************************************************/
    /*                                                                */
    /* Free the memories                                              */
    /*                                                                */
    /******************************************************************/

    free (x);
    free (y_b1);
    free (y_b2);
    free (y_b3);
    free (y_b4);
    free (w);
    /* Free allocated memory */
    if (free_2d_array ((void **) x_t) != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: x_t\n", FUNC_NAME, ERROR);
    }

}

/******************************************************************************
MODULE:  median_filter

PURPOSE:  calculate mediam filtering for three buf scanline

RETURN VALUE: None

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
6/12/2019   Su Ye            Original Development

NOTES:
******************************************************************************/
int median_filter
(
    short int **buf1,
    short int **buf2,
    short int **buf3,
    short int **buf,
    int *valid_datecount_scanline,       /* I: the number of valid dates               */
    int num_scenes,
    int n_cols
)
{
    int i, b, j;
    short int window[9];

    for(i = 0; i < n_cols; i++)
    {
        for(j = 0; j < valid_datecount_scanline[i]; j++)
        {
            if(i==0)
            {
                for(b = 0; b < TOTAL_IMAGE_BANDS; b++)
                {
                    window[0] = buf1[b][i * num_scenes + j];
                    window[1] = buf1[b][i * num_scenes + j];
                    window[2] = buf1[b][(i+1) * num_scenes + j];
                    window[3] = buf2[b][i * num_scenes + j];
                    window[4] = buf2[b][i * num_scenes + j];
                    window[5] = buf2[b][(i+1) * num_scenes + j];
                    window[6] = buf3[b][i * num_scenes + j];
                    window[7] = buf3[b][i * num_scenes + j];
                    window[8] = buf3[b][(i+1) * num_scenes + j];

                    quick_sort_shortint(window, 0, 8);

                    buf[b][i * num_scenes + j] = window[4];

                }
            }
            else if (i == n_cols - 1)
            {
                for(b = 0; b < TOTAL_IMAGE_BANDS; b++)
                {
                    window[0] = buf1[b][(i-1) * num_scenes + j];
                    window[1] = buf1[b][i * num_scenes + j];
                    window[2] = buf1[b][i * num_scenes + j];
                    window[3] = buf2[b][(i-1) * num_scenes + j];
                    window[4] = buf2[b][i * num_scenes + j];
                    window[5] = buf2[b][i * num_scenes + j];
                    window[6] = buf3[b][(i-1) * num_scenes + j];
                    window[7] = buf3[b][i * num_scenes + j];
                    window[8] = buf3[b][i * num_scenes + j];

                    quick_sort_shortint(window, 0, 8);

                    buf[b][i * num_scenes + j] = window[4];

                }
            }
            else
            {
                for(b = 0; b < TOTAL_IMAGE_BANDS; b++)
                {
                    window[0] = buf1[b][(i-1) * num_scenes + j];
                    window[1] = buf1[b][i * num_scenes + j];
                    window[2] = buf1[b][(i+1) * num_scenes + j];
                    window[3] = buf2[b][(i-1) * num_scenes + j];
                    window[4] = buf2[b][i * num_scenes + j];
                    window[5] = buf2[b][(i+1) * num_scenes + j];
                    window[6] = buf3[b][(i-1) * num_scenes + j];
                    window[7] = buf3[b][i * num_scenes + j];
                    window[8] = buf3[b][(i+1) * num_scenes + j];

                    quick_sort_shortint(window, 0, 8);

                    buf[b][i * num_scenes + j] = window[4];

                }
            }

        }

    }

}

max_stars_repo_name: agroimpacts/imager

max_issues_repo_name: agroimpacts/imager

max_forks_repo_name: agroimpacts/imager


text:
#ifndef _read_h_included_
#define _read_h_included_

#include <assert.h>
#include <stdlib.h>
#include <algorithm>
#include <cctype>
#include <fstream>
#include <functional>
#include <iostream>
#include <map>
#include <numeric>
#include <set>
#include <sstream>
#include <vector>
#include <sstream>
#include <string>
#include <boost/config.hpp>
#include <boost/foreach.hpp>
#include <boost/random.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <jsc/bioinfo/gene_anno.hpp>
#include <jsc/util/log.hpp>

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

#include "splicing_graph.h"
#include "accessible_read_starts.h"

using namespace std;
using namespace boost;
using namespace boost::numeric;
using namespace jsc::bioinfo;
using namespace jsc::util;

using boost::shared_ptr;

// returns the index (in the isoform) of the first matching exon in the isoform
long is_connected_exons_compatible_with_isoform(
		vector<unsigned long> const & exon_indices,
		vector<unsigned long> const & iso_exon_indices) {
	vector<unsigned long>::const_iterator itr, itr2;
	itr = exon_indices.begin();
	itr2 = iso_exon_indices.begin();
	bool is_compatible = false;
	bool found_first_match = false;
	long first_match_idx = 0;
	while (itr != exon_indices.end()) {
		if (itr2 == iso_exon_indices.end()) {
			is_compatible = false;
			break;
		}
		if (*itr != *itr2) {
			if (!found_first_match) {
				itr2++;
				first_match_idx++;
			} else {
				is_compatible = false;
				break;
			}
		} else {
			if (!found_first_match) {
				found_first_match = true;
				is_compatible = true;
				itr++;
				itr2++;
			} else {
				itr++;
				itr2++;
			}
		}
	}
	return (is_compatible ? first_match_idx : -1);
}

class ReadType {
	public:
		static const int LONG_READ = 0;
		static const int MEDIUM_READ = 1;
		static const int SHORT_READ = 2;
		static const int SHORT_CONTIG = 3;
		static const int CONTIG = 4;
		static const int LONG_READ_PE = 5;
		static const int MEDIUM_READ_PE = 6;
		static const int SHORT_READ_PE = 7;
		static const string TYPES[];
};
const string ReadType::TYPES[] = {"LONG_READ", "MEDIUM_READ", "SHORT_READ", "SHORT_CONTIG", "CONTIG", "LONG_READ_PE", "MEDIUM_READ_PE", "SHORT_READ_PE"};

class Read {
	public:
		string read_type;
		vector<unsigned long> exon_indices;
		unsigned long start_at_first_exon;
		unsigned long end_at_last_exon;

		Read(shared_ptr<AccessibleReadStarts> const & ars
				= shared_ptr<AccessibleReadStarts>())
			: ars(ars) {
			read_length = 0;
			expected_read_length = 0;
			read_type = "General read";
		}

		unsigned long get_num_isoforms() const {
			return ars->get_num_isoforms();
		}

		unsigned long get_read_length() const {
			return read_length;
		}

		shared_ptr<AccessibleReadStarts> get_ARS() const {
			return ars;
		}

		virtual unsigned long get_read_span() const {
			return read_length;
		}

		unsigned long get_expected_read_length() const {
			return expected_read_length;
		}

		virtual unsigned long get_expected_read_span() const {
			return expected_read_length;
		}

		virtual bool is_compatible_with_iso(
				unsigned long const & iso_idx) const = 0;

		virtual void build(
				interval_list<long> const & read_il,
				SetExonp const & exonps,
				interval_list<long> const & read_il2 = interval_list<long>()) = 0;

		/// Generates a read from an isoform.
		//
		/// The start position of the read in the isoform is either
		/// determined by \a fixed_read_start or chosen randomly
		/// according to the actual type of the read.
		//
		/// \param[in] iso_idx							index of the isoform from which
		///																	the read will be generated.
		///																	exons in the corresponding gene.
		/// \param[in] fixed_read_start			When >= 0, specifies the start
		///																	position (in the isoform) of the
		///																	read to be generated.
		///	\returns	The actual start position (in the isoform)
		///						of the generated read.
		virtual unsigned long generate_read(
				unsigned long const & iso_idx,
				long fixed_read_start = -1) = 0;

		virtual double prob_generated_by_iso(
				unsigned long const & iso_idx) const = 0;

		virtual string to_string() const = 0;

		virtual ~Read() {}

	protected:
		shared_ptr<AccessibleReadStarts> ars;
		unsigned long read_length;
		unsigned long expected_read_length;
};

ostream& operator << (ostream& os, Read const & r) {
	os << r.to_string();
	return os;
};

typedef shared_ptr<Read> Read_ptr;

class Read_stats {
	public:
		static const unsigned long medium_single_expected_read_length = 250;
		static const double medium_cost_per_bp = 7E-5;
		static const unsigned long short_single_expected_read_length = 30;
		static const unsigned long short_paired_expected_insert_size = 300;
		static const double short_cost_per_bp = 7E-6;
};

class Read_single : public Read {
	public:
		Read_single(shared_ptr<AccessibleReadStarts> ars
				= shared_ptr<AccessibleReadStarts>())
			: Read(ars)	{
			rng = NULL;
			read_type = "General single read";
		}

		virtual bool is_compatible_with_iso(
				unsigned long const & iso_idx) const {
			return (is_connected_exons_compatible_with_isoform(
						exon_indices, ars->iso_exon_indices[iso_idx]) >= 0);
		}

		virtual void build(
				interval_list<long> const & read_il,
				SetExonp const & exonps,
				interval_list<long> const & read_il2 = interval_list<long>()) {
			L_(debug2) << "Read::build";
			unsigned long matching_length = 0;
			set<unsigned long> index_set;
			exon_indices.clear();
			start_at_first_exon = end_at_last_exon = 0;
			bool found_start = false;
			unsigned long il_idx = 0;
			SetExonp::const_iterator itr = exonps.begin();
			long cur_start, cur_end, cur_exon_start;
			cur_exon_start = 0;
			L_(debug2) << "Read: " << read_il;
			L_(debug2) << "Exons: " << exonps;
			while (il_idx < read_il.get_num_intervals()) {
				cur_start = read_il.get_starts()[il_idx];
				cur_end = read_il.get_ends()[il_idx];
				L_(debug2) << "Current interval: [" << cur_start << ","
					<< cur_end << ")";
				// find all the exons that overlap with/contain the interval
				while (itr != exonps.end() && (*itr)->start < cur_end) {
					L_(debug2) << "Current exon: [" << (*itr)->start << ","
						<< (*itr)->end << ")";
					cur_exon_start = max(cur_exon_start, (*itr)->start);
					if (cur_start >= cur_exon_start && cur_start < (*itr)->end) {
						if (!found_start) {
							L_(debug2) << "Found start";
							found_start = true;
							start_at_first_exon = cur_start - (*itr)->start;
						} else if (cur_start > cur_exon_start) {
							break;
						}

						L_(debug2) << "Add exon #" << distance(exonps.begin(), itr);
						index_set.insert(distance(exonps.begin(), itr));
						cur_exon_start = min((*itr)->end, cur_end);
						matching_length += cur_exon_start - cur_start;

						if (cur_end < (*itr)->end) {
							cur_start = cur_end;
							end_at_last_exon = cur_end - (*itr)->start;
							++il_idx;
							break;
						} else if (cur_end == (*itr)->end) {
							cur_start = cur_end;
							end_at_last_exon = cur_end - (*itr)->start;
							++il_idx;
						} else {
							cur_start = (*itr)->end;
							end_at_last_exon = (*itr)->end - (*itr)->start;
						}
					} else if (cur_exon_start > (*itr)->start
							&& cur_exon_start < (*itr)->end) {
						break;
					}
					itr++;
				}
				if (cur_start == cur_end) {
					continue;
				} else {
					break;
				}
			}

			BOOST_FOREACH (unsigned long i, index_set) {
				exon_indices.push_back(i);
			}
			read_length = matching_length;
		}

		virtual unsigned long generate_read(
				unsigned long const & iso_idx,
				long fixed_read_start = -1) {
			assert(rng != NULL || fixed_read_start >= 0);
			vector<unsigned long> const & iso_exon_indices =
				ars->iso_exon_indices[iso_idx];
			vector<unsigned long> const & iso_exon_total_lengths =
				ars->iso_exon_total_lengths[iso_idx];

			exon_indices.clear();
			unsigned long read_start;
			if (fixed_read_start < 0) {	// randomly generate a read start
				read_start = (unsigned long)(
						gsl_ran_flat(rng, 0,
							ars->get_iso_ARS_total_length(iso_idx)));
				read_start = ars->ARStart2IsoStart(iso_idx, read_start);
			} else {
				read_start = fixed_read_start;
			}

			unsigned long read_end = read_start + get_expected_read_length();
			read_length = read_end - read_start;
			L_(debug2) << "Generating read [" << read_start
				<< ", " << read_end << ") length=" << read_length << "bp";
			vector<unsigned long>::const_iterator itr =
				upper_bound(iso_exon_total_lengths.begin(),
						iso_exon_total_lengths.end(),
						read_start);
			assert(itr != iso_exon_total_lengths.end());
			unsigned long start_exon_idx = distance(
					iso_exon_total_lengths.begin(), itr);
			itr = lower_bound(iso_exon_total_lengths.begin(),
					iso_exon_total_lengths.end(),
					read_end);
			assert(itr != iso_exon_total_lengths.end());
			unsigned long end_exon_idx = distance(
					iso_exon_total_lengths.begin(), itr);
			L_(debug2) << "Start exon: "
				<< iso_exon_indices[start_exon_idx]
				<< ", end exon: "
				<< iso_exon_indices[end_exon_idx];
			for (unsigned long i = start_exon_idx;
					i <= end_exon_idx; ++i) {
				exon_indices.push_back(iso_exon_indices[i]);
			}
			start_at_first_exon = read_start -
				(start_exon_idx == 0 ?
				 0 : iso_exon_total_lengths[start_exon_idx - 1]);
			end_at_last_exon = read_end -
				(end_exon_idx == 0 ?
				 0 : iso_exon_total_lengths[end_exon_idx - 1]);

			return read_start;
		}

		virtual double prob_generated_by_iso(
				unsigned long const & iso_idx) const {
			double num_diff_reads = ars->get_iso_ARS_total_length(iso_idx);
			if (num_diff_reads <= 0) {
				L_(warning) << read_type << ": supported read length > isoform length" << endl;
				return 0;
			}
			double prob = (double)1.0 / num_diff_reads;
			return prob;
		}

		virtual string to_string() const {
			std::ostringstream os;
			os << "[" << read_type << "] "
				<< "(Exons) ";
			BOOST_FOREACH(unsigned long const & idx, exon_indices) {
				os << idx << ",";
			}
			os << " (Start at first exon) " << start_at_first_exon << " ";
			os << "(End at last exon) " << end_at_last_exon;
			return os.str();
		}
	protected:
		gsl_rng * rng;
};

class Read_medium_single : public Read_single {
	public:
		Read_medium_single(shared_ptr<AccessibleReadStarts> const & ars
				= shared_ptr<AccessibleReadStarts>(),
				gsl_rng * r = NULL,
				unsigned long expected_length =
				Read_stats::medium_single_expected_read_length)
			: Read_single(ars) {
			rng = r;
			expected_read_length = expected_length;
			read_type = "Medium single read";
		}
};

class Read_short_single : public Read_single {
	public:
		Read_short_single(shared_ptr<AccessibleReadStarts> const & ars
				= shared_ptr<AccessibleReadStarts>(),
				gsl_rng * r = NULL,
				unsigned long expected_length =
				Read_stats::short_single_expected_read_length)
			: Read_single(ars) {
			rng = r;
			expected_read_length = expected_length;
			read_type = "Short single read";
		}
};

class Read_paired: public Read {
	public:
		Read_single end1;
		Read_single end2;
		Read_single span;

		Read_paired(shared_ptr<AccessibleReadStarts> ars
				= shared_ptr<AccessibleReadStarts>(),
				double tolerance = 0)
			: Read(ars), tolerance(tolerance) {
			insert_size = expected_insert_size = 0;
			read_type = "General paired-end read";
			end1 = Read_single(ars);
			end2 = Read_single(ars);
			span = Read_single(ars);
		}

		unsigned long get_expected_insert_size() const {
			return expected_insert_size;
		}

		virtual unsigned long get_read_span() const {
			return read_length + insert_size;
		}

		virtual unsigned long get_expected_read_span() const {
			return (end1.get_expected_read_length()
				+ end2.get_expected_read_length()
				+ expected_insert_size);
		}

		virtual bool is_compatible_with_iso(
				unsigned long const & iso_idx) const {
			assert(end1.exon_indices.size() > 0
					&& end2.exon_indices.size() > 0);
			if (end1.is_compatible_with_iso(iso_idx)
					&& end2.is_compatible_with_iso(iso_idx)) {
				if (expected_insert_size == 0) {
					return true;
				}
				vector<unsigned long> const & iso_exon_indices =
					ars->iso_exon_indices[iso_idx];
				vector<unsigned long> const & exon_total_lengths =
					ars->iso_exon_total_lengths[iso_idx];

				vector<unsigned long>::const_iterator itr;
				unsigned long end1_last_exon_idx = end1.exon_indices
					[end1.exon_indices.size() - 1];
				itr = lower_bound(iso_exon_indices.begin(),
						iso_exon_indices.end(),
						end1_last_exon_idx);
				unsigned long exon_idx = distance(
						iso_exon_indices.begin(), itr);
				unsigned long end1_end =
					(exon_idx > 0 ? exon_total_lengths[exon_idx - 1] : 0)
					+ end1.end_at_last_exon;
				unsigned long end2_first_exon_idx = end2.exon_indices[0];
				itr = lower_bound(iso_exon_indices.begin(),
						iso_exon_indices.end(),
						end2_first_exon_idx);
				exon_idx = distance(iso_exon_indices.begin(), itr);
				unsigned long end2_start =
					(exon_idx > 0 ? exon_total_lengths[exon_idx - 1] : 0)
					+ end2.start_at_first_exon;
				double gap = end2_start - end1_end;
				double diff = gap / (double)(expected_insert_size);
				if (diff <= (1.0 + tolerance + epsilon)
						&& diff >= (1.0 - tolerance - epsilon)) {
					return true;
				}
			}
			return false;
		}

		virtual double prob_generated_by_iso(
				unsigned long const & iso_idx) const {
			double num_diff_reads = ars->get_iso_ARS_total_length(iso_idx);
			if (num_diff_reads <= 0) {
				L_(warning) << read_type << ": supported paired-end read length (including insert) > isoform length" << endl;
				return 0;
			}
			double prob = (double)1.0 / num_diff_reads;
			return prob;
		}

		virtual string to_string() const {
			std::ostringstream os;
			os << "[" << read_type << "]"
				<< " End #1: " << end1.to_string()
				<< "; End #2: " << end2.to_string();
			return os.str();
		}

		virtual void build(
				interval_list<long> const & read_il,
				SetExonp const & exonps,
				interval_list<long> const & read_il2 = interval_list<long>()) {
			L_(debug2) << "Read::build";
			end1.build(read_il, exonps);
			end2.build(read_il2, exonps);
			read_length = end1.get_read_length() + end2.get_read_length();
			/// \todo compute insert size and construct span, exon_indices, etc.
		}

		virtual unsigned long generate_read(
				unsigned long const & iso_idx,
				long fixed_read_start = -1) {
			assert(rng != NULL || fixed_read_start >= 0);

			unsigned long read_start;
			if (fixed_read_start < 0) {	// randomly generate a read start
				read_start = (unsigned long)(
						gsl_ran_flat(rng, 0,
							ars->get_iso_ARS_total_length(iso_idx)));
				read_start = ars->ARStart2IsoStart(iso_idx, read_start);
			} else {
				read_start = fixed_read_start;
			}

			unsigned long end1_start =
				end1.generate_read(iso_idx, read_start);
			unsigned long end2_start =
				end1_start + end1.get_read_length()
				+ expected_insert_size;
			end2.generate_read(iso_idx, end2_start);
			read_length = end1.get_read_length() + end2.get_read_length();

			span.generate_read(iso_idx, read_start);
			exon_indices = span.exon_indices;
			start_at_first_exon = span.start_at_first_exon;
			end_at_last_exon = span.end_at_last_exon;

			return end1_start;
		}

	protected:
		unsigned long insert_size;
		unsigned long expected_insert_size;
		double tolerance;
		static const double epsilon = 1E-5;
		gsl_rng * rng;
};

void generate_reads(unsigned long const & num_reads,
		Isoforms const & iso,
		shared_ptr<Read> readp,
		gsl_rng * rng,
		ublas::matrix<double> & m_delta_G,
		bool known_isoforms_only = false) {
	gsl_ran_discrete_t * discrete_g = gsl_ran_discrete_preproc(
			iso.num_known_isoforms, iso.known_iso_probs);

	for (unsigned long i = 0; i < num_reads; ++i) {
		unsigned long iso_idx = gsl_ran_discrete(rng, discrete_g);
		L_(debug2) << "Selected known isoform " << iso_idx;
		if (!known_isoforms_only) {
			iso_idx = iso.known_iso_indices[iso_idx];
		}
		readp->generate_read(iso_idx);
		L_(debug2) << *readp;
		if (known_isoforms_only) {
			for (unsigned long j = 0;
					j < iso.num_known_isoforms; ++j) {
				if (readp->is_compatible_with_iso(j)) {
					m_delta_G(i,j) = readp->prob_generated_by_iso(j);
				}
			}
		} else {
			for (unsigned long j = 0;
					j < iso.num_possible_isoforms; ++j) {
				if (readp->is_compatible_with_iso(j)) {
					m_delta_G(i,j) = readp->prob_generated_by_iso(j);
				}
			}
		}
	}

	gsl_ran_discrete_free(discrete_g);
}

class Read_short_paired: public Read_paired {
	public:
		Read_short_paired(shared_ptr<AccessibleReadStarts> ars
				= shared_ptr<AccessibleReadStarts>(),
				gsl_rng * r = NULL,
				unsigned long expected_length =
				2 * Read_stats::short_single_expected_read_length,
				unsigned long expected_insert =
				Read_stats::short_paired_expected_insert_size,
				double tolerance = 0)
			: Read_paired(ars, tolerance) {
			rng = r;
			expected_read_length = expected_length;
			insert_size = expected_insert_size = expected_insert;
			read_type = "Short paired-end read";
			shared_ptr<AccessibleReadStarts> end_ars(
					new AccessibleReadStarts(
						ars->iso_exon_indices,
						ars->iso_exon_total_lengths,
						ars->exon_lengths,
						expected_length / 2));
			end1 = Read_short_single(end_ars, r, expected_length / 2);
			end2 = Read_short_single(end_ars, r, expected_length / 2);
			span = Read_medium_single(ars, r, expected_length + expected_insert);
		}
};

void compute_iso_probs_EM_step(
		unsigned long num_possible_isoforms,
		vector<ublas::matrix<double> > const & m_delta_Gs,
		vector<double> const & old_theta,
		vector<double> & new_theta) {
	for (unsigned k = 0; k < num_possible_isoforms; ++k) {
		double sum_zeta = 0;
		double num_total_reads = 0;
		BOOST_FOREACH (ublas::matrix<double> const & m_delta_G, m_delta_Gs) {
			unsigned long num_reads = m_delta_G.size1();
			num_total_reads += (double)num_reads;
			for (unsigned long i = 0; i < num_reads; ++i) {
				double sum_local_probs = 0;
				for (unsigned k2 = 0; k2 < num_possible_isoforms; ++k2) {
					sum_local_probs += old_theta[k2] * m_delta_G(i, k2);
				}
				if (sum_local_probs > 0) {
					double local_prob = old_theta[k] * m_delta_G(i, k);
					if (local_prob > 0) {
						sum_zeta += local_prob / sum_local_probs;
					}
				}
			}
		}
		new_theta[k] = sum_zeta / num_total_reads;
	}
}

double reads_log_likelihood(
		unsigned long num_possible_isoforms,
		vector<ublas::matrix<double> > const & m_delta_Gs,
		vector<double> const & theta) {
	double ll = 0;
	BOOST_FOREACH (ublas::matrix<double> const & m_delta_G, m_delta_Gs) {
		unsigned long num_reads = m_delta_G.size1();
		for (unsigned long i = 0; i < num_reads; ++i) {
			double sum_local_probs = 0;
			for (unsigned k = 0; k < num_possible_isoforms; ++k) {
				sum_local_probs += theta[k] * m_delta_G(i, k);
			}
			ll += log(sum_local_probs);
		}
	}
	return ll;
}

void compute_iso_probs_EM(
		unsigned long num_possible_isoforms,
		vector<ublas::matrix<double> > const & m_delta_Gs,
		vector<double> & theta) {
	theta.resize(num_possible_isoforms,
			(double)1.0 / (double)num_possible_isoforms);
	double ll, old_ll;
	vector<double> old_theta;
	unsigned long count = 0;
	do {
		old_theta = theta;
		old_ll = reads_log_likelihood(
				num_possible_isoforms, m_delta_Gs, old_theta);
		compute_iso_probs_EM_step(num_possible_isoforms,
				m_delta_Gs, old_theta, theta);
		ll = reads_log_likelihood(
				num_possible_isoforms, m_delta_Gs, theta);
		if (++count == 100) {
			L_(debug) << "log likelihood: " << old_ll << " -> " << ll;
			count = 0;
		}
	} while (abs(1.0 - old_ll / ll) > 1E-6);
}

void err_known_iso_probs(Isoforms const & iso,
		vector<double> const & theta, double & abs_err,
		double & sse) {
	abs_err = sse = 0;
	for (unsigned long i = 0; i < iso.num_known_isoforms; ++i) {
		abs_err += abs(iso.known_iso_probs[i]
				- theta[iso.known_iso_indices[i]]);
		sse += (iso.known_iso_probs[i] - theta[iso.known_iso_indices[i]]) *
			(iso.known_iso_probs[i] - theta[iso.known_iso_indices[i]]);
	}
}

void err_possible_iso_probs(Isoforms const & iso,
		vector<double> const & theta, double & abs_err,
		double & sse) {
	abs_err = sse = 0;
	for (unsigned long i = 0; i < iso.num_possible_isoforms; ++i) {
		abs_err += theta[i];
		sse += theta[i] * theta[i];
	}
	for (unsigned long i = 0; i < iso.num_known_isoforms; ++i) {
		abs_err -= theta[iso.known_iso_indices[i]];
		abs_err += abs(iso.known_iso_probs[i]
				- theta[iso.known_iso_indices[i]]);
		sse -= theta[iso.known_iso_indices[i]] * theta[iso.known_iso_indices[i]];
		sse += (iso.known_iso_probs[i] - theta[iso.known_iso_indices[i]]) *
			(iso.known_iso_probs[i] - theta[iso.known_iso_indices[i]]);
	}
}

#endif

ext: h

max_stars_repo_name: gersteinlab/LESSeq

max_issues_repo_name: gersteinlab/LESSeq

max_forks_repo_name: gersteinlab/LESSeq


text:
/*
 *  mvdens.h
 *  likely
 *
 *  Created by Karim Benabed on 10/03/08.
 *  Copyright 2008 __MyCompanyName__. All rights reserved.
 *
 */

#ifndef __MVDENS_H
#define __MVDENS_H
#ifdef __PLANCK__

#include "HL2_likely/tools/errorlist.h"
#include "HL2_likely/tools/io.h"
#include "HL2_likely/tools/maths_base.h"

#else

#include "errorlist.h"
#include "io.h"
#include "maths_base.h"

#endif

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_int.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_gamma.h>

/*#ifdef HAS_LAPACK
#ifdef HAS_MKL
#include "mkl_blas.h"
#include "mkl_lapack.h"
#else
#include "blas.h"
#include "lapack.h"
#endif
#endif
*/

#define __MVDENS_PARANOID_DEBUG__

/* errors */
#define mv_base       -700
#define mv_allocate   -1 + mv_base
#define mv_serialize  -2 + mv_base
#define mv_outOfBound -3 + mv_base
#define mv_badComm    -4 + mv_base
#define mv_negWeight  -5 + mv_base
#define mv_cholesky   -6 + mv_base
#define mv_negative   -7 + mv_base
#define mv_undef      -8 + mv_base
#define mv_file       -9 + mv_base
#define mv_io        -10 + mv_base
#define mv_tooManySteps -11 + mv_base
#define mv_dimension -12 + mv_base
#define mv_type      -13 + mv_base
#define mv_negHatCl  -14 + mv_base

#define PI         3.141592653589793
#define LOGSQRT2PI 0.918938533204673
#ifndef DYNMAX
#define DYNMAX 500.0
#endif
#define mv_INF 1e99

#define mv_EXP  2
#define mv_NORM 1
#define mv_UNORM 0

#define MVT_DF 3

typedef struct {
  size_t ndim;
  void *buf;
  int own_buf;
  double *mean;
  double *std;
  double *x_tmp;
  gsl_vector_view  mean_view_container;
  gsl_vector *mean_view;
  gsl_vector_view  x_tmp_view_container;
  gsl_vector *x_tmp_view;
  gsl_matrix_view  std_view_container;
  gsl_matrix  *std_view;
  size_t band_limit;
  int chol;
  int df;
  double detL;   /* Determinant of L, where L^t L = Covariance matrix */
} mvdens;

typedef struct {
  size_t ncomp, ndim;
  int init_cwght;
  double *wght, *cwght;
  gsl_vector_view wght_view_container, cwght_view_container;
  gsl_vector *wght_view, *cwght_view;
  mvdens **comp;
  void *buf_comp;
} mix_mvdens;

// init and stuff
mvdens * mvdens_alloc(size_t ndim, error **err);
mvdens * mvdens_init(mvdens *g, size_t nndim, double* mn, double* std, error **err);
void mvdens_empty(mvdens *g);
void mvdens_free(mvdens ** g);
void mvdens_free_void(void **g);
void mvdens_print(FILE* where, mvdens* what);
mvdens *mvdens_dwnp(FILE* where, error **err);
mvdens *mvdens_read_and_chol(FILE *where, error **err);
void mvdens_dump(FILE* where, mvdens* what);
void mvdens_chdump(const char *name, mvdens* what, error **err);
void mvdens_set_band_limit(mvdens *self, size_t value);
void mvdens_from_meanvar(mvdens *m, const double *pmean, const double *pvar, double correct);

//computations
double determinant(const double *L, size_t ndim);
void mvdens_cholesky_decomp(mvdens* self, error **err);
double* mvdens_ran(double* dest, mvdens * g, gsl_rng * r,error **err);
double scalar_product(mvdens *g, const double *x, error **err);
double mvdens_log_pdf(mvdens *g, const double * x, error ** err);
double mvdens_log_pdf_void(void *g, const double *x, error **err);
double mvdens_inverse(mvdens *m, error **err);

/* mix_mvdens functions */
mix_mvdens *mix_mvdens_alloc(size_t ncomp, size_t size,error **err);
void mix_mvdens_copy(mix_mvdens *target, const mix_mvdens *source, error **err);
void mix_mvdens_free_void(void **m);
void mix_mvdens_free(mix_mvdens **m);
void mix_mvdens_print(FILE* where,mix_mvdens* what);
void mix_mvdens_dump(FILE* where,mix_mvdens* what);
mix_mvdens * mix_mvdens_dwnp(FILE* where,error **err);
void mix_mvdens_set_band_limit(mix_mvdens *self, size_t value);
double effective_number_of_components(const mix_mvdens *self, error **err);

size_t mix_mvdens_size(size_t ncomp,size_t ndim);
void* serialize_mix_mvdens(mix_mvdens* self, size_t *sz,error** err);
mix_mvdens* deserialize_mix_mvdens(void* serialized, size_t sz,error **err);

void mix_mvdens_cholesky_decomp(mix_mvdens* self, error **err);
double* mix_mvdens_ran(double* dest,size_t *index, mix_mvdens *m, gsl_rng * r, error **err);
size_t ranindx(double *cw, size_t nE, gsl_rng * r, error ** err);
double mix_mvdens_log_pdf(mix_mvdens *m, const double *x, error **err);
double mix_mvdens_log_pdf_void(void *m, const double *x, error **err);

#define MALLOC_IF_NEEDED(ret,dest,size,err) { \
  if (dest==NULL) { \
    ret= (double*) malloc_err(size,err); \
  } else { \
    ret=dest; \
  } \
}

#ifdef HAS_HDF5
#include "hdf5.h"
#define hdf5_base mv_io
void mvdens_hdfdump_infile(mvdens* mv,hid_t loc_id,error **err);
void mvdens_hdfdump(mvdens* mv,char *fname,error **err);
void mix_mvdens_hdfdump(mix_mvdens* mmv,char *fname,error **err);
mvdens* mvdens_hdfdwnp_infile(mvdens* mv,hid_t loc_id,error **err);
mvdens* mvdens_hdfdwnp(char *fname,error **err);
mix_mvdens* mix_mvdens_hdfdwnp(char *fname,error **err);
#endif

#endif


text:
/*
  Copyright [2017-2021] [IBM Corporation]
  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at
  http://www.apache.org/licenses/LICENSE-2.0
  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License.
*/

#ifndef _MCAS_POOL_SESSION_H_
#define _MCAS_POOL_SESSION_H_

#include <gsl/pointers>
#include <memory>

class Pool_instance;

struct pool_session {
  pool_session(gsl::not_null<std::shared_ptr<Pool_instance>> ph) : pool(ph), canary(0x45450101) {
  }

  ~pool_session() {
  }

  bool check() const { return canary == 0x45450101; }
  gsl::not_null<std::shared_ptr<Pool_instance>> pool;
  const unsigned canary;
};

#endif

max_stars_repo_name: IBM/artemis

max_issues_repo_name: IBM/artemis

max_forks_repo_name: IBM/artemis


text:
#ifndef __SAMPLE_H__
#define __SAMPLE_H__

#include <gsl/gsl_rng.h>
#include <pthread.h>
#include <stdint.h>

#include "parse_mmaps.h"

#define NUM_USER_LOCKS 2000

struct sample_thread_info;
struct index_update;
struct sample_threads;

/* Thread utility functions */

void initialize_threads(struct sample_threads* sample_threads, int num_threads,
			const struct mmap_info* mmap_info);
void destroy_threads(struct sample_threads* sample_threads);

/* Resampling functions */

/* Set all topic and POV assignments to -1, not doing any index
   updates. Used before initializing topics and POVs.*/
void resample_null(struct sample_threads* sample_threads);
/* Initialize topic and POV assignments by sampling them uniformly at
   random and performing index updates. */
void resample_uniform(struct sample_threads* sample_threads);

/* Re-sample topics and POVs. One iteration of Gibbs sampling. */
void resample(struct sample_threads* sample_threads);
/* Re-sample assignments, but (approximately) in the opposite order as
   resample(). Used during model selection. */
void resample_reverse(struct sample_threads* sample_threads);

/* Rather than re-sampling randomly, always choose the maximum
   probability assignment. Useful for finding a high probability
   assignment of topics and POVs after random sampling. */
void resample_maximize(struct sample_threads* sample_threads);
/* Restore a set of assignments, updating indexes as necessary */
void resample_restore(struct sample_threads* sample_threads,
		      const struct revision_assignment* revision_assignments);
/* Set all topic and POV assignments to -1, subtracting the
   assignments from indexes. Primarily useful for verifying that
   indexes were correct. */
void resample_zero(struct sample_threads* sample_threads);

/* Parallel probability computations */

/* Compute the log probability of transitioning to the specified
   assignments from the current assignments (don't actually transition
   to them). */
double transition_probability(struct sample_threads* sample_threads,
			      const struct revision_assignment* revision_assignments);
/* Compute the likelihood of the current model and topic/POV
   assignments */
double parallel_log_likelihood(struct sample_threads* sample_threads);

/* Other miscellaneous parallel utility routines */

/* Set all of the user topic/POV distributions to alpha. Used during
   initialization. */
void parallel_initialize_user_topics(struct sample_threads* sample_threads);

/* Structs to hold synchronization and thread information. */

struct sample_thread_info {
  /* Each thread gets its own random number generator, initialized
     with a different seed. */
  gsl_rng* rand_gen;

  /* A num_topics * num_povs_per_topic array which threads use to
     sample new topic and POV assignments. */
  double* sampling_array;

  /* Each thread gets its own copy of mmap_info, allowing us to
     replace topic_index_mmap with a private copy which is manually
     synchronized. */
  struct mmap_info mmap_info;

  // Sample only pages numbered (page# % mod_n) == sample_pages
  int sample_pages;
  int mod_n;

  // Information about the current thread.
  pthread_t thread;

  /* An array of locks to synchronize access and updates to user
     topic/POV distributions. For a given user's topic and POV
     distribution, the lock to use is userid % NUM_USER_LOCKS. */
  pthread_rwlock_t* user_locks;

  /* Memory accounting for the topic_index_mmap copy, if any. One
     thread updates the original, and so does not have a specially
     allocated topic index, in which case allocated_topic_dist will be
     NULL. */
  char* allocated_topic_dist;

  /* The index of the first position in the topic index update queue
     that we have not read. */
  int64_t last_queue_position;

  /* Pointers to the global topic index update queue, and the global
     position of the first empty slot in the queue. Access to the
     queue location is synchronized with queue_lock. */
  struct index_update* index_update_queue;
  int64_t* queue_location;
  pthread_mutex_t* queue_lock;

  /* Generally +1 or -1, indicating whether sampling is forward or
     backward. */
  int increment;

  /* Helper functions to specify sampling behavior, replacing small or
     large parts of the sampling routine. */
  void (*revision_callback) (struct sample_thread_info*, int64_t);
  void (*index_update_function) (const struct mmap_info* mmap_info,
				 const struct index_update* patch);
  void (*sample_function) (const double* sampling_array,
			   double probability_sum,
			   int num_topics, int pov_per_topic,
			   int* chosen_topic, int* chosen_pov,
			   gsl_rng* rand_gen);

  /* When computing transition probabilities, these assignments
     specify which assignments the transition is to. */
  const struct revision_assignment* reference_assignments;

  double output;
};

struct sample_threads {
  struct sample_thread_info* thread_info;
  pthread_mutex_t queue_lock;
  pthread_rwlock_t user_locks[NUM_USER_LOCKS];
  int num_threads;

  struct index_update* index_update_queue;
  // One past the last value in the queue
  int64_t queue_location;
};

#endif

max_stars_repo_name: allenlavoie/topic-pov

max_issues_repo_name: allenlavoie/topic-pov

max_forks_repo_name: allenlavoie/topic-pov


text:
#pragma once
#include <autograph/Core/Types.h>
#include <gsl/span>
#include <gsl/string_span>

namespace ag
{
	// until a standard solution appears
	using gsl::span;
	using gsl::string_span;
}

max_stars_repo_name: ennis/autograph-pipelines

max_issues_repo_name: ennis/autograph-pipelines

max_forks_repo_name: ennis/autograph-pipelines


text:
#pragma once

#include <gsl\gsl>
#include "Texture.h"
#include "Rectangle.h"

namespace Library
{
	class Texture2D final : public Texture
	{
		RTTI_DECLARATIONS(Texture2D, Texture)

	public:
		Texture2D(const winrt::com_ptr<ID3D11ShaderResourceView>& shaderResourceView, std::uint32_t width, std::uint32_t height);
		Texture2D(const Texture2D&) = default;
		Texture2D& operator=(const Texture2D&) = default;
		Texture2D(Texture2D&&) = default;
		Texture2D& operator=(Texture2D&&) = default;
		~Texture2D() = default;

		static std::shared_ptr<Texture2D> CreateTexture2D(gsl::not_null<ID3D11Device*> device, const D3D11_TEXTURE2D_DESC& textureDesc);
		static std::shared_ptr<Texture2D> CreateTexture2D(gsl::not_null<ID3D11Device*> device, std::uint32_t width, std::uint32_t height, std::uint32_t mipLevels = 1, std::uint32_t arraySize = 1, DXGI_FORMAT format = DXGI_FORMAT_R8G8B8A8_UNORM, DXGI_SAMPLE_DESC sampleDesc = { 1, 0 }, std::uint32_t bindFlags = D3D11_BIND_UNORDERED_ACCESS | D3D11_BIND_SHADER_RESOURCE, std::uint32_t cpuAccessFlags = 0);

		std::uint32_t Width() const;
		std::uint32_t Height() const;
		Rectangle Bounds() const;

	private:
		std::uint32_t mWidth;
		std::uint32_t mHeight;
	};
}


text:
/*
Copyright (c) 2014, Kai Klindworth
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#ifndef BLAS_WRAPPER_H
#define BLAS_WRAPPER_H

#include <cblas.h>

namespace blas {

inline void gemm(double* Y, const double* A, bool transposeA, int rowsA, int colsA, const double* B, bool transposeB, int rowsB, int colsB, double alpha = 1.0, double beta = 0.0)
{
	cblas_dgemm(CblasRowMajor, transposeA ? CblasTrans : CblasNoTrans, transposeB ? CblasTrans : CblasNoTrans, transposeA ? colsA : rowsA, transposeB ? rowsB : colsB, transposeB ? colsB : rowsB, alpha, A, colsA, B, colsB, beta, Y, transposeB ? rowsB : colsB);
}

inline void gemm(float* Y, const float* A, bool transposeA, int rowsA, int colsA, const float* B, bool transposeB, int rowsB, int colsB, float alpha = 1.0f, float beta = 0.0f)
{
	cblas_sgemm(CblasRowMajor, transposeA ? CblasTrans : CblasNoTrans, transposeB ? CblasTrans : CblasNoTrans, transposeA ? colsA : rowsA, transposeB ? rowsB : colsB, transposeB ? colsB : rowsB, alpha, A, colsA, B, colsB, beta, Y, transposeB ? rowsB : colsB);
}

inline void gemv(double* Y, const double* A, bool transposeA, int rowsA, int colsA, const double* x, double alpha = 1.0, double beta = 0.0)
{
	cblas_dgemv(CblasRowMajor, transposeA ? CblasTrans : CblasNoTrans, rowsA, colsA, alpha, A, colsA, x, 1, beta, Y, 1);
}

inline void gemv(float* Y, const float* A, bool transposeA, int rowsA, int colsA, const float* x, float alpha = 1.0, float beta = 0.0)
{
	cblas_sgemv(CblasRowMajor, transposeA ? CblasTrans : CblasNoTrans, rowsA, colsA, alpha, A, colsA, x, 1, beta, Y, 1);
}

inline void ger(double* Y, const double * A, int rowsA, const double* X, int rowsX, double alpha = 1.0)
{
	cblas_dger(CblasRowMajor, rowsA, rowsX, alpha, A, 1, X, 1, Y, rowsX);
}

inline void ger(float* Y, const float * A, int rowsA, const float* X, int rowsX, float alpha = 1.0)
{
	cblas_sger(CblasRowMajor, rowsA, rowsX, alpha, A, 1, X, 1, Y, rowsX);
}

inline double norm2(double* X, int n, int stride = 1)
{
	return cblas_dnrm2(n, X, stride);
}

inline float norm2(float* X, int n, int stride = 1)
{
	return cblas_snrm2(n, X, stride);
}

inline void scale(float alpha, float* X, int n, int stride = 1)
{
	return cblas_sscal(n, alpha, X, stride);
}

inline void scale(double alpha, double* X, int n, int stride = 1)
{
	return cblas_dscal(n, alpha, X, stride);
}
}

#endif


text:
/* specfunc/mathieu_angfunc.c
 *
 * Copyright (C) 2002 Lowell Johnson
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  L. Johnson */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_mathieu.h>

int gsl_sf_mathieu_ce_e(int order, double qq, double zz, gsl_sf_result *result)
{
  int even_odd, ii, status;
  double coeff[GSL_SF_MATHIEU_COEFF], norm, fn, factor;
  gsl_sf_result aa;

  norm = 0.0;
  even_odd = 0;
  if (order % 2 != 0)
      even_odd = 1;

  /* Handle the trivial case where q = 0. */
  if (qq == 0.0)
  {
      norm = 1.0;
      if (order == 0)
          norm = sqrt(2.0);

      fn = cos(order*zz)/norm;

      result->val = fn;
      result->err = 2.0*GSL_DBL_EPSILON;
      factor = fabs(fn);
      if (factor > 1.0)
          result->err *= factor;

      return GSL_SUCCESS;
  }

  /* Use symmetry characteristics of the functions to handle cases with
     negative order. */
  if (order < 0)
      order *= -1;

  /* Compute the characteristic value. */
  status = gsl_sf_mathieu_a_e(order, qq, &aa);
  if (status != GSL_SUCCESS)
  {
      return status;
  }

  /* Compute the series coefficients. */
  status = gsl_sf_mathieu_a_coeff(order, qq, aa.val, coeff);
  if (status != GSL_SUCCESS)
  {
      return status;
  }

  if (even_odd == 0)
  {
      fn = 0.0;
      norm = coeff[0]*coeff[0];
      for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
      {
          fn += coeff[ii]*cos(2.0*ii*zz);
          norm += coeff[ii]*coeff[ii];
      }
  }
  else
  {
      fn = 0.0;
      for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
      {
          fn += coeff[ii]*cos((2.0*ii + 1.0)*zz);
          norm += coeff[ii]*coeff[ii];
      }
  }

  norm = sqrt(norm);
  fn /= norm;

  result->val = fn;
  result->err = 2.0*GSL_DBL_EPSILON;
  factor = fabs(fn);
  if (factor > 1.0)
      result->err *= factor;

  return GSL_SUCCESS;
}

int gsl_sf_mathieu_se_e(int order, double qq, double zz, gsl_sf_result *result)
{
  int even_odd, ii, status;
  double coeff[GSL_SF_MATHIEU_COEFF], norm, fn, factor;
  gsl_sf_result aa;

  norm = 0.0;
  even_odd = 0;
  if (order % 2 != 0)
      even_odd = 1;

  /* Handle the trivial cases where order = 0 and/or q = 0. */
  if (order == 0)
  {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
  }

  if (qq == 0.0)
  {
      norm = 1.0;
      fn = sin(order*zz);

      result->val = fn;
      result->err = 2.0*GSL_DBL_EPSILON;
      factor = fabs(fn);
      if (factor > 1.0)
          result->err *= factor;

      return GSL_SUCCESS;
  }

  /* Use symmetry characteristics of the functions to handle cases with
     negative order. */
  if (order < 0)
      order *= -1;

  /* Compute the characteristic value. */
  status = gsl_sf_mathieu_b_e(order, qq, &aa);
  if (status != GSL_SUCCESS)
  {
      return status;
  }

  /* Compute the series coefficients. */
  status = gsl_sf_mathieu_b_coeff(order, qq, aa.val, coeff);
  if (status != GSL_SUCCESS)
  {
      return status;
  }

  if (even_odd == 0)
  {
      fn = 0.0;
      for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
      {
          norm += coeff[ii]*coeff[ii];
          fn += coeff[ii]*sin(2.0*(ii + 1)*zz);
      }
  }
  else
  {
      fn = 0.0;
      for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
      {
          norm += coeff[ii]*coeff[ii];
          fn += coeff[ii]*sin((2.0*ii + 1)*zz);
      }
  }
  norm = sqrt(norm);
  fn /= norm;

  result->val = fn;
  result->err = 2.0*GSL_DBL_EPSILON;
  factor = fabs(fn);
  if (factor > 1.0)
      result->err *= factor;

  return GSL_SUCCESS;
}

int gsl_sf_mathieu_ce_array(int nmin, int nmax, double qq, double zz,
                            gsl_sf_mathieu_workspace *work,
                            double result_array[])
{
  int even_odd, order, ii, jj, status;
  double coeff[GSL_SF_MATHIEU_COEFF], *aa = work->aa, norm;

  /* Initialize the result array to zeroes. */
  for (ii=0; ii<nmax-nmin+1; ii++)
      result_array[ii] = 0.0;

  /* Ensure that the workspace is large enough to accomodate. */
  if (work->size < (unsigned int)nmax)
  {
      GSL_ERROR("Work space not large enough", GSL_EINVAL);
  }

  if (nmin < 0 || nmax < nmin)
  {
      GSL_ERROR("domain error", GSL_EDOM);
  }

  /* Compute all of the eigenvalues up to nmax. */
  gsl_sf_mathieu_a_array(0, nmax, qq, work, aa);

  for (ii=0, order=nmin; order<=nmax; ii++, order++)
  {
      norm = 0.0;
      even_odd = 0;
      if (order % 2 != 0)
          even_odd = 1;

      /* Handle the trivial case where q = 0. */
      if (qq == 0.0)
      {
          norm = 1.0;
          if (order == 0)
              norm = sqrt(2.0);

          result_array[ii] = cos(order*zz)/norm;

          continue;
      }

      /* Compute the series coefficients. */
      status = gsl_sf_mathieu_a_coeff(order, qq, aa[order], coeff);
      if (status != GSL_SUCCESS)
          return status;

      if (even_odd == 0)
      {
          norm = coeff[0]*coeff[0];
          for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
          {
              result_array[ii] += coeff[jj]*cos(2.0*jj*zz);
              norm += coeff[jj]*coeff[jj];
          }
      }
      else
      {
          for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
          {
              result_array[ii] += coeff[jj]*cos((2.0*jj + 1.0)*zz);
              norm += coeff[jj]*coeff[jj];
          }
      }

      norm = sqrt(norm);
      result_array[ii] /= norm;
  }

  return GSL_SUCCESS;
}

int gsl_sf_mathieu_se_array(int nmin, int nmax, double qq, double zz,
                            gsl_sf_mathieu_workspace *work,
                            double result_array[])
{
  int even_odd, order, ii, jj, status;
  double coeff[GSL_SF_MATHIEU_COEFF], *bb = work->bb, norm;

  /* Initialize the result array to zeroes. */
  for (ii=0; ii<nmax-nmin+1; ii++)
      result_array[ii] = 0.0;

  /* Ensure that the workspace is large enough to accomodate. */
  if (work->size < (unsigned int)nmax)
  {
      GSL_ERROR("Work space not large enough", GSL_EINVAL);
  }

  if (nmin < 0 || nmax < nmin)
  {
      GSL_ERROR("domain error", GSL_EDOM);
  }

  /* Compute all of the eigenvalues up to nmax. */
  gsl_sf_mathieu_b_array(0, nmax, qq, work, bb);

  for (ii=0, order=nmin; order<=nmax; ii++, order++)
  {
      norm = 0.0;
      even_odd = 0;
      if (order % 2 != 0)
          even_odd = 1;

      /* Handle the trivial cases where order = 0 and/or q = 0. */
      if (order == 0)
      {
          norm = 1.0;
          result_array[ii] = 0.0;
          continue;
      }

      if (qq == 0.0)
      {
          norm = 1.0;
          result_array[ii] = sin(order*zz);

          continue;
      }

      /* Compute the series coefficients. */
      status = gsl_sf_mathieu_b_coeff(order, qq, bb[order], coeff);
      if (status != GSL_SUCCESS)
      {
          return status;
      }

      if (even_odd == 0)
      {
          for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
          {
              result_array[ii] += coeff[jj]*sin(2.0*(jj + 1)*zz);
              norm += coeff[jj]*coeff[jj];
          }
      }
      else
      {
          for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
          {
              result_array[ii] += coeff[jj]*sin((2.0*jj + 1.0)*zz);
              norm += coeff[jj]*coeff[jj];
          }
      }

      norm = sqrt(norm);
      result_array[ii] /= norm;
  }

  return GSL_SUCCESS;
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_mathieu_ce(int order, double qq, double zz)
{
	EVAL_RESULT(gsl_sf_mathieu_ce_e(order, qq, zz, &result));
}

double gsl_sf_mathieu_se(int order, double qq, double zz)
{
	EVAL_RESULT(gsl_sf_mathieu_se_e(order, qq, zz, &result));
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
** ECM - eigenvector centrality mapping
**
** G.Lohmann, MPI-KYB, updated Oct 2018
*/

#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

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

#define SQR(x) ((x) * (x))
#define ReLU(x) ((x) > 0 ? (x) : 0)

extern void   VMatrixECM(gsl_matrix_float *,float *,int,int,VImage);
extern void   VMatrixProjection(gsl_matrix_float *X,float *ev,int,int seed);
extern VImage VoxelMap(VAttrList list,VImage mask,size_t nvox);
extern size_t NumVoxels(VImage);
extern gsl_matrix_float *ReadDataECM(VAttrList,VImage,VImage,VShort,VShort,int,size_t);
extern VImage WriteOutput(VAttrList,VImage,float *,size_t);
extern int Irreducible(gsl_matrix_float *X,int type);

void NormVec(float *x,size_t nvox)
{
 size_t i;
 float sum=0;
 sum = 0;
 for (i=0; i<nvox; i++) sum += (x[i]*x[i]);
 sum = sqrt(sum);
 for (i=0; i<nvox; i++) x[i] /= sum;
}

/* needed for ECM-RLC, negative */
void SignFlip(gsl_matrix_float *X)
{
  float u=0;
  size_t nt = X->size2;
  size_t i,j,n=nt/2;
  for (i=0; i<X->size1; i++) {
    for (j=n; j<nt; j++) {
      u = gsl_matrix_float_get(X,i,j);
      gsl_matrix_float_set(X,i,j,-u);
    }
  }
}

/* fast, memory efficient power iteration, used in ECM-RLC, ECM-Add1 */
void PowerIteration(gsl_matrix_float *X,float *ev,int type,int maxiter)
{
  int iter;
  float vsum=0,d=0;
  size_t i;
  size_t nvox = X->size1;
  size_t nt = X->size2;
  float scale = 1.0/(float)nt;

  fprintf(stderr," power iteration...\n");
  gsl_vector_float *y = gsl_vector_float_calloc(nvox);
  gsl_vector_float *z = gsl_vector_float_calloc(nt);

  Irreducible(X,type);  /* check if corr matrix is irreducible */

  NormVec(ev,nvox);
  for (i=0; i<nvox; i++) y->data[i] = ev[i];

  /* power iterations */
  for (iter=0; iter<maxiter; iter++) {

    /* R = (X^t X)/n */
    gsl_blas_sgemv (CblasTrans,1.0,X,y,0.0,z);
    if (type == 7) SignFlip(X);
    gsl_blas_sgemv (CblasNoTrans,1.0,X,z,0.0,y);
    if (type == 7) SignFlip(X);
    gsl_vector_float_scale(y,scale);

    /* ECM-Add1,  (R+E)v = Rv + Ev */
    if (type == 1) {
      vsum = 0;
      for (i=0; i<nvox; i++) vsum += ev[i];
      for (i=0; i<nvox; i++) y->data[i] += vsum;
    }

    /* normalize */
    NormVec(y->data,nvox);

    /* check convergence */
    d = 0;
    for (i=0; i<nvox; i++) {
      d += SQR(ev[i] - y->data[i]);
      ev[i] = y->data[i];
    }
    fprintf(stderr," %5d   %.6f\n",(int)iter,d);

    if (iter > 2 && d < 1.0e-6) break;
  }
  gsl_vector_float_free(y);
  gsl_vector_float_free(z);
}

float *VECM(gsl_matrix_float *X,int type,VBoolean project,int seed,int maxiter,VImage map)
{
  size_t i,j;
  size_t nvox = X->size1;
  size_t nt = X->size2;
  if (type == 0 || type == 7) nt /= 2;

  /* ini eigenvector */
  float u=0;
  float *eigvec = (float *) VCalloc(nvox,sizeof(float));
  for (i=0; i<nvox; i++) eigvec[i] = 1.0/(float)nt;

  /* no matrix projection */
  if (!project) {

    switch (type) {

    case 0:   /*  ECM-RLC positive */
      for (i=0; i<X->size1; i++) {
	for (j=0; j<nt; j++) {
	  u = gsl_matrix_float_get(X,i,j);
	  gsl_matrix_float_set(X,i,j+nt,fabs(u));
	}
      }
      PowerIteration(X,eigvec,type,maxiter);
      break;

    case 1:   /* ECM-Add1 */
      PowerIteration(X,eigvec,type,maxiter);
      break;

    case 6:   /* no filter, risky */
      PowerIteration(X,eigvec,1,maxiter);  /* initialize using ecm-add1 */
      PowerIteration(X,eigvec,0,maxiter);
      break;

      /* ECM-RLC negative */
    case 7:
      for (i=0; i<X->size1; i++) {
	for (j=0; j<nt; j++) {
	  u = gsl_matrix_float_get(X,i,j);
	  gsl_matrix_float_set(X,i,j,fabs(u));
	  gsl_matrix_float_set(X,i,j+nt,u);
	}
      }
      PowerIteration(X,eigvec,type,maxiter);
      break;

    default:  /* dense matrix representation */
      VMatrixECM(X,eigvec,type,maxiter,map);
      break;

    }
  }

  /* Matrix projection, ECM-project */
  else {
    VMatrixProjection(X,eigvec,type,seed);
  }

  /* rescale for better visualization */
  float nx = (float)nvox;
  nx = sqrt(nx);
  for (i=0; i<nvox; i++) eigvec[i] *= nx;

  return eigvec;
}

VDictEntry TYPDict[] = {
  { "rlc",     0, 0,0,0,0 },
  { "add",     1, 0,0,0,0 },
  { "pos",     2, 0,0,0,0 },
  { "abs",     3, 0,0,0,0 },
  { "neg",     4, 0,0,0,0 },
  { "gauss",   5, 0,0,0,0 },
  { NULL, 0,0,0,0,0 }
};

int main (int argc,char *argv[])
{
  static VString  mask_filename = "";
  static VShort   first  = 0;
  static VShort   length = 0;
  static VShort   type   = 0;
  static VShort   maxiter= 20;
  static VBoolean project = FALSE;
  static VInteger  seed   = 99402622;
  static VShort   nproc  = 0;
  static VOptionDescRec  options[] = {
    {"mask",VStringRepn,1,(VPointer) &mask_filename,VRequiredOpt,NULL,"Region of interest mask"},
    {"first",VShortRepn,1,(VPointer) &first,VOptionalOpt,NULL,"First timestep to use"},
    {"length",VShortRepn,1,(VPointer) &length,VOptionalOpt,NULL,"Length of time series to use, '0' to use full length"},
    {"metric",VShortRepn,1,(VPointer) &type,VOptionalOpt,TYPDict,"Type of correlation metric"},
    {"iterations",VShortRepn,1,(VPointer) &maxiter,VOptionalOpt,NULL,"Maximum number of power iterations"},
    {"project",VBooleanRepn,1,(VPointer) &project,VOptionalOpt,NULL,"Whether to do matrix projection"},
    {"seed",VIntegerRepn,1,(VPointer) &seed,VOptionalOpt,NULL,"Seed for random number generation (only for ECM-project)"},
    {"j",VShortRepn,1,(VPointer) &nproc,VOptionalOpt,NULL,"Number of processors to use, '0' to use all"}
  };
  FILE *in_file=NULL,*out_file=NULL;
  VString in_filename=NULL;
  char *prg_name=GetLipsiaName("vecm");
  fprintf(stderr, "%s\n", prg_name);

  /* parse command line */
  VParseFilterCmdZ (VNumber (options),options,argc,argv,&in_file,&out_file,&in_filename);
  if (type < 0 || type > 8) VError(" unknown metric");
  if (type < 6) fprintf(stderr," Correlation metric %d:  %s\n",type,TYPDict[type].keyword);
  if (project && (type < 2 || type == 7))
    VError(" matrix projection not needed for metric %d (%s)\n",type,TYPDict[type].keyword);

  /* omp-stuff */
#ifdef _OPENMP
  int num_procs=omp_get_num_procs();
  if (nproc > 0 && nproc < num_procs) num_procs = nproc;
  if (type > 3 || project == TRUE) fprintf(stderr,"using %d cores\n",(int)num_procs);
  omp_set_num_threads(num_procs);
#endif /* _OPENMP */

  /* read functional data */
  VAttrList list = VReadAttrListZ(in_file,in_filename,0L,TRUE,FALSE);
  if (list == NULL) VError(" error reading input file %s",in_file);
  VAttrList geolist = VGetGeoInfo(list);

  /* read mask */
  VAttrList listm = VReadAttrList(mask_filename,0L,TRUE,FALSE);
  VImage mask = VReadImage(listm);
  if (mask == NULL) VError(" no mask found");
  size_t nvox = NumVoxels(mask);

  /* voxel map */
  VImage map = VoxelMap(list,mask,nvox);

  /* read data into X */
  gsl_matrix_float *X = ReadDataECM(list,mask,map,first,length,type,nvox);

  /* main process */
  float *eigvec = VECM(X,(int)type,project,(VLong)seed,(int)maxiter,map);

  /* create output image */
  VImage dest = WriteOutput(list,map,eigvec,nvox);
  VSetAttr(VImageAttrList(dest),"name",NULL,VStringRepn,"ECM");
  VAttrList out_list = VCreateAttrList();
  VAppendAttr(out_list,"image",NULL,VImageRepn,dest);

  /* update geoinfo, 4D to 3D */
  if (geolist != NULL) {
    double *D = VGetGeoDim(geolist,NULL);
    D[0] = 3;  /* 3D */
    D[4] = 1;  /* just one timestep */
    VSetGeoDim(geolist,D);
    VSetGeoInfo(geolist,out_list);
  }

  /* write to disk */
  VHistory(VNumber(options),options,prg_name,&list,&out_list);
  if (! VWriteFile (out_file, out_list)) exit (1);
  fprintf (stderr, "%s: done.\n", argv[0]);
  return 0;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#ifndef __GSL_VERSION_H__
#define __GSL_VERSION_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_types.h>

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

#define GSL_VERSION "2.4"
#define GSL_MAJOR_VERSION 2
#define GSL_MINOR_VERSION 4

GSL_VAR const char * gsl_version;

__END_DECLS

#endif /* __GSL_VERSION_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
/* movstat/qnacc.c
 *
 * Moving window Q_n accumulator
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

typedef double qnacc_type_t;
typedef qnacc_type_t ringbuf_type_t;

#include "ringbuf.c"

typedef struct
{
  qnacc_type_t *window; /* linear array for current window */
  qnacc_type_t *work;   /* workspace, length 3*n */
  int *work_int;        /* integer workspace, length 5*n */
  ringbuf *rbuf;        /* ring buffer storing current window */
} qnacc_state_t;

static size_t
qnacc_size(const size_t n)
{
  size_t size = 0;

  size += sizeof(qnacc_state_t);
  size += n * sizeof(qnacc_type_t);     /* window */
  size += 3 * n * sizeof(qnacc_type_t); /* work */
  size += 5 * n * sizeof(int);          /* work_int */
  size += ringbuf_size(n);

  return size;
}

static int
qnacc_init(const size_t n, void * vstate)
{
  qnacc_state_t * state = (qnacc_state_t *) vstate;

  state->window = (qnacc_type_t *) ((unsigned char *) vstate + sizeof(qnacc_state_t));
  state->work = (qnacc_type_t *) ((unsigned char *) state->window + n * sizeof(qnacc_type_t));
  state->work_int = (int *) ((unsigned char *) state->work + 3 * n * sizeof(qnacc_type_t));
  state->rbuf = (ringbuf *) ((unsigned char *) state->work_int + 5 * n * sizeof(int));

  ringbuf_init(n, state->rbuf);

  return GSL_SUCCESS;
}

static int
qnacc_insert(const qnacc_type_t x, void * vstate)
{
  qnacc_state_t * state = (qnacc_state_t *) vstate;

  /* add new element to ring buffer */
  ringbuf_insert(x, state->rbuf);

  return GSL_SUCCESS;
}

static int
qnacc_delete(void * vstate)
{
  qnacc_state_t * state = (qnacc_state_t *) vstate;

  if (!ringbuf_is_empty(state->rbuf))
    ringbuf_pop_back(state->rbuf);

  return GSL_SUCCESS;
}

/* FIXME XXX: this is inefficient - could be improved by maintaining a sorted ring buffer */
static int
qnacc_get(void * params, qnacc_type_t * result, const void * vstate)
{
  const qnacc_state_t * state = (const qnacc_state_t *) vstate;
  size_t n = ringbuf_copy(state->window, state->rbuf);

  (void) params;

  gsl_sort(state->window, 1, n);
  *result = gsl_stats_Qn_from_sorted_data(state->window, 1, n, state->work, state->work_int);

  return GSL_SUCCESS;
}

static const gsl_movstat_accum qn_accum_type =
{
  qnacc_size,
  qnacc_init,
  qnacc_insert,
  qnacc_delete,
  qnacc_get
};

const gsl_movstat_accum *gsl_movstat_accum_Qn = &qn_accum_type;

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#pragma once

#include <string>
#include <utility>
#include <vector>
#include <unordered_map>
#include <gsl/span>

#include <lmdb.h>
#include <functional>
#include "Transaction.h"
#include "LmdbExceptions.h"

namespace DeepestScatter
{
    class Dataset
    {
    public:
        struct Settings;
        class Example;

        explicit Dataset(std::shared_ptr<Settings> settings);
        ~Dataset();

        template<class T>
        size_t getRecordsCount();

        template<class T>
        T getRecord(int32_t recordId);

        template<class T>
        void dropTable();

        template<class T>
        void append(const T& example);

        template<class T>
        void batchAppend(const gsl::span<T>& examples, int32_t startId);

        struct Settings
        {
            Settings(std::string path) :
                path(std::move(path))
            {}

            std::string path;
        };

    private:
        using TableName = std::string;

        template<class T>
        void tryAppend(const T& example, int entryId);

        template<class T>
        void tryBatchAppend(const gsl::span<T>& examples, int nextExampleId);

        void increaseSizeIfNeededWhile(const std::function<void(void)>& action);
        void increaseSize();

        MDB_dbi getTable(const TableName& name);
        MDB_dbi openTable(const TableName& name);
        void closeTable(const TableName& name);

        std::unordered_map<TableName, MDB_dbi> openedTables;
        std::unordered_map<TableName, int32_t> nextIds;

        MDB_env* mdbEnv = nullptr;
    };

    template <class T>
    size_t Dataset::getRecordsCount()
    {
        const TableName tableName = T::descriptor()->name();

        MDB_dbi dbi = getTable(tableName);

        return Transaction::withTransaction<size_t>(mdbEnv, nullptr, 0, [&](Transaction& transaction)
        {
            MDB_stat stats;
            LmdbExceptions::checkError(mdb_stat(transaction, dbi, &stats));
            return stats.ms_entries;
        });
    }

    template <class T>
    T Dataset::getRecord(int32_t recordId)
    {
        const TableName tableName = T::descriptor()->name();
        MDB_dbi dbi = getTable(tableName);

        return Transaction::withTransaction<T>(mdbEnv, nullptr, 0, [&](Transaction& transaction)
        {
            MDB_val mdbKey
            {
                sizeof(int32_t),
                &recordId
            };

            MDB_val mdbVal;
            LmdbExceptions::checkError(mdb_get(transaction, dbi, &mdbKey, &mdbVal));

            T record{};
            record.ParseFromArray(mdbVal.mv_data, mdbVal.mv_size);

            return record;
        });
    }

    template <class T>
    void Dataset::dropTable()
    {
        const TableName tableName = T::descriptor()->name();
        MDB_dbi dbi = getTable(tableName);

        if (getRecordsCount<T>() == 0)
        {
            return;
        }

        while (true)
        {
            std::string tableConfirmation;

            std::cout << "YOU ARE GOING TO DELETE " << getRecordsCount<T>() << " RECORDS!!!" << std::endl;
            std::cout << "Type table name: \"" << T::descriptor()->name() << "\" to drop it" << std::endl;
            std::cin >> tableConfirmation;

            if (tableConfirmation == T::descriptor()->name())
            {
                break;
            }

            std::cout << "Mismatched table name. Try again.";
        }

        Transaction::withTransaction<void>(mdbEnv, nullptr, 0, [&](Transaction& transaction)
        {
            mdb_drop(transaction, dbi, 0);
            nextIds[tableName] = 0;
        });
    }

    template <class T>
    void Dataset::append(const T& example)
    {
        const TableName tableName = T::descriptor()->name();
        // Returns zero if not initialized;
        int entryId = nextIds[tableName];

        increaseSizeIfNeededWhile([&]()
        {
            return tryAppend(example, entryId);
        });

        nextIds[tableName] = entryId + 1;
    }

    template <class T>
    void Dataset::batchAppend(const gsl::span<T>& examples, int32_t startId)
    {
        const TableName tableName = T::descriptor()->name();

        increaseSizeIfNeededWhile([&]()
        {
            return tryBatchAppend(examples, startId);
        });

        nextIds[tableName] = startId + examples.size();
    }

    template <class T>
    void Dataset::tryAppend(const T& example, int entryId)
    {
        const TableName tableName = T::descriptor()->name();
        MDB_dbi dbi = getTable(tableName);

        Transaction::withTransaction<void>(mdbEnv, nullptr, 0, [&](Transaction& transaction)
        {

            MDB_val mdbKey
            {
                sizeof(int32_t),
                &entryId
            };

            const size_t length = example.ByteSizeLong();
            std::vector<uint8_t> serialized(length);
            example.SerializeWithCachedSizesToArray(&serialized[0]);

            MDB_val mdbVal
            {
                length * sizeof(uint8_t),
                static_cast<void*>(&serialized[0])
            };

            LmdbExceptions::checkError(mdb_put(transaction, dbi, &mdbKey, &mdbVal, 0));
        });
    }

    template <class T>
    void Dataset::tryBatchAppend(const gsl::span<T>& examples, int nextExampleId)
    {
        const TableName tableName = T::descriptor()->name();
        MDB_dbi dbi = getTable(tableName);

        Transaction::withTransaction<void>(mdbEnv, nullptr, 0, [&](Transaction& transaction)
        {
            for (const auto& example : examples)
            {
                MDB_val mdbKey
                {
                    sizeof(int32_t),
                    &nextExampleId
                };

                const size_t length = example.ByteSizeLong();
                std::vector<uint8_t> serialized(length);
                example.SerializeWithCachedSizesToArray(&serialized[0]);

                MDB_val mdbVal
                {
                    length * sizeof(uint8_t),
                    static_cast<void*>(&serialized[0])
                };

                LmdbExceptions::checkError(mdb_put(transaction, dbi, &mdbKey, &mdbVal, 0));
                nextExampleId++;
            }
        });
    }
}

max_stars_repo_name: marsermd/DeepestScatter

max_issues_repo_name: marsermd/DeepestScatter

max_forks_repo_name: marsermd/DeepestScatter


text:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>

int
func_f (const gsl_vector * x, void *params, gsl_vector * f)
{
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);

  gsl_vector_set(f, 0, 100.0 * (x2 - x1*x1));
  gsl_vector_set(f, 1, 1.0 - x1);

  return GSL_SUCCESS;
}

int
func_df (const gsl_vector * x, void *params, gsl_matrix * J)
{
  double x1 = gsl_vector_get(x, 0);

  gsl_matrix_set(J, 0, 0, -200.0*x1);
  gsl_matrix_set(J, 0, 1, 100.0);
  gsl_matrix_set(J, 1, 0, -1.0);
  gsl_matrix_set(J, 1, 1, 0.0);

  return GSL_SUCCESS;
}

int
func_fvv (const gsl_vector * x, const gsl_vector * v,
          void *params, gsl_vector * fvv)
{
  double v1 = gsl_vector_get(v, 0);

  gsl_vector_set(fvv, 0, -200.0 * v1 * v1);
  gsl_vector_set(fvv, 1, 0.0);

  return GSL_SUCCESS;
}

void
callback(const size_t iter, void *params,
         const gsl_multifit_nlinear_workspace *w)
{
  gsl_vector * x = gsl_multifit_nlinear_position(w);

  /* print out current location */
  printf("%f %f\n",
         gsl_vector_get(x, 0),
         gsl_vector_get(x, 1));
}

void
solve_system(gsl_vector *x0, gsl_multifit_nlinear_fdf *fdf,
             gsl_multifit_nlinear_parameters *params)
{
  const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
  const size_t max_iter = 200;
  const double xtol = 1.0e-8;
  const double gtol = 1.0e-8;
  const double ftol = 1.0e-8;
  const size_t n = fdf->n;
  const size_t p = fdf->p;
  gsl_multifit_nlinear_workspace *work =
    gsl_multifit_nlinear_alloc(T, params, n, p);
  gsl_vector * f = gsl_multifit_nlinear_residual(work);
  gsl_vector * x = gsl_multifit_nlinear_position(work);
  int info;
  double chisq0, chisq, rcond;

  /* initialize solver */
  gsl_multifit_nlinear_init(x0, fdf, work);

  /* store initial cost */
  gsl_blas_ddot(f, f, &chisq0);

  /* iterate until convergence */
  gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol,
                              callback, NULL, &info, work);

  /* store final cost */
  gsl_blas_ddot(f, f, &chisq);

  /* store cond(J(x)) */
  gsl_multifit_nlinear_rcond(&rcond, work);

  /* print summary */

  fprintf(stderr, "NITER         = %zu\n", gsl_multifit_nlinear_niter(work));
  fprintf(stderr, "NFEV          = %zu\n", fdf->nevalf);
  fprintf(stderr, "NJEV          = %zu\n", fdf->nevaldf);
  fprintf(stderr, "NAEV          = %zu\n", fdf->nevalfvv);
  fprintf(stderr, "initial cost  = %.12e\n", chisq0);
  fprintf(stderr, "final cost    = %.12e\n", chisq);
  fprintf(stderr, "final x       = (%.12e, %.12e)\n",
          gsl_vector_get(x, 0), gsl_vector_get(x, 1));
  fprintf(stderr, "final cond(J) = %.12e\n", 1.0 / rcond);

  printf("\n\n");

  gsl_multifit_nlinear_free(work);
}

int
main (void)
{
  const size_t n = 2;
  const size_t p = 2;
  gsl_vector *f = gsl_vector_alloc(n);
  gsl_vector *x = gsl_vector_alloc(p);
  gsl_multifit_nlinear_fdf fdf;
  gsl_multifit_nlinear_parameters fdf_params =
    gsl_multifit_nlinear_default_parameters();

  /* print map of Phi(x1, x2) */
  {
    double x1, x2, chisq;
    double *f1 = gsl_vector_ptr(f, 0);
    double *f2 = gsl_vector_ptr(f, 1);

    for (x1 = -1.2; x1 < 1.3; x1 += 0.1)
      {
        for (x2 = -0.5; x2 < 2.1; x2 += 0.1)
          {
            gsl_vector_set(x, 0, x1);
            gsl_vector_set(x, 1, x2);
            func_f(x, NULL, f);

            chisq = (*f1) * (*f1) + (*f2) * (*f2);
            printf("%f %f %f\n", x1, x2, chisq);
          }
        printf("\n");
      }
    printf("\n\n");
  }

  /* define function to be minimized */
  fdf.f = func_f;
  fdf.df = func_df;
  fdf.fvv = func_fvv;
  fdf.n = n;
  fdf.p = p;
  fdf.params = NULL;

  /* starting point */
  gsl_vector_set(x, 0, -0.5);
  gsl_vector_set(x, 1, 1.75);

  fprintf(stderr, "=== Solving system without acceleration ===\n");
  fdf_params.trs = gsl_multifit_nlinear_trs_lm;
  solve_system(x, &fdf, &fdf_params);

  fprintf(stderr, "=== Solving system with acceleration ===\n");
  fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
  solve_system(x, &fdf, &fdf_params);

  gsl_vector_free(f);
  gsl_vector_free(x);

  return 0;
}