text:
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

// y' = 2y + x
// https://www.wolframalpha.com/input/?i=y%27+%3D+2y+%2B+x%2C+y%280%29+%3D+1%2C+y%281%29
// https://www.wolframalpha.com/input/?i=1%2F4+*+%285+*+e%5E2+-+3%29

// 저항값: 1Kohm
// 인덕턴스: 300uH
// DC 전원: 5V

// ---DC---R---L---GND---
// 직렬 회로의 특성상 전류값이 일정함 i(t)
// V = i(t) * R + L * [di(t) / dt]
// 5 = y * R + L * y'
// 미분 방정식의 표준형은 무조건 y'의 계수가 1이어야 합니다!!!
// 5 / L = (R / L) * y + y'
// 5 / 0.0003 = (1000 / 0.0003) * y + y'

// y' = 5 / 0.0003 - (1000 / 0.0003) * y   <<<=== 표준형
int odefunc (double x, const double y[], double f[], void *params)
{
    // f[0] = x+2*y[0];
	//f[0] = pow(M_E, x);
	//f[0] = -(1000 / 0.0003) * x - (5 / 0.0003);
	f[0] = -(1000 / 0.0003) * y[0] + (5 / 0.0003);

    return GSL_SUCCESS;
}

int * jac;

int main(void)
{
    int dim = 1;
	// 미방을 풀기 위한 데이터 타입
    gsl_odeiv2_system sys = {odefunc, NULL, dim, NULL};

	// 드라이버는 미방을 쉽게 풀 수 있도록 시간에 따른 변화, 제어, 스텝에 대한 래퍼
	// gsl_odeiv2_step_rkf45는 룬게쿠타 4-5차를 의미함
	// 다음 인자는 미분방정식의 시작지점
	// 절대적인 에러 바운드와 상대적인 에러 바운드
    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rkf45, 1e-6, 1e-6, 0.0);
    int i;
    double x0 = 0.0,  xf = 10.0; /* start and end of integration interval */
    double x = x0;
    double y[1] = { 0.5 };  /* initial value */
	double tmp = 0.0;
	int cnt = 0;

    for (i = 0; tmp <= xf; tmp += 0.0000001)
    {
        double xi = x0 + tmp * (xf-x0) / xf;
        int status = gsl_odeiv2_driver_apply (d, &x, xi, y);

        if (status != GSL_SUCCESS)
        {
            printf ("error, return value=%d\n", status);
            break;
        }

		if (cnt++ < 101)
		{
        	printf ("%.8e %.8e\n", x, y[0]);
		}
    }

    gsl_odeiv2_driver_free (d);
    return 0;
}


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

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

struct f
{
  char *name;
  double (*f) (gsl_complex z);
  double x;
  double y;
  double fx;
  double fy;
};

struct fz
{
  char *name;
  gsl_complex (*f) (gsl_complex z);
  double x;
  double y;
  double fx;
  double fy;
};

#define FN(x) "gsl_complex_" #x, gsl_complex_ ## x
#define ARG(x,y) x, y
#define RES(x,y) x, y

struct f list[] =
{
#include "results1.h"
  {"", 0, 0, 0, 0, 0}
};

struct fz listz[] =
{
#include "results.h"
  {"", 0, 0, 0, 0, 0}
};

int
main (void)
{
  size_t i = 0;

  gsl_ieee_env_setup();

  for (i = 0 ; i < 10; i++)
    {
      double r = (i - 5.0) * 0.3 ;
      double t = 2.0 * M_PI * i / 5 ;
      double x = r * cos(t), y = r * sin(t) ;
      gsl_complex z = gsl_complex_polar (r, t) ;
      gsl_test_rel (GSL_REAL(z), x, 10 * GSL_DBL_EPSILON, "gsl_complex_polar real part at (r=%g,t=%g)", r, t);

      gsl_test_rel (GSL_IMAG(z), y, 10 * GSL_DBL_EPSILON, "gsl_complex_polar imag part at (r=%g,t=%g)", r, t);
    }

    i = 0;

  while (list[i].f)
    {
      struct f t = list[i];
      gsl_complex z = gsl_complex_rect (t.x, t.y);
      double f = (t.f) (z);
      gsl_test_rel (f, t.fx, 10 * GSL_DBL_EPSILON, "%s at (%g,%g)", t.name, t.x, t.y);
      i++;
    }

  i = 0;

  while (listz[i].f)
    {
      struct fz t = listz[i];
      gsl_complex z = gsl_complex_rect (t.x, t.y);
      gsl_complex fz = (t.f) (z);
      double fx = GSL_REAL (fz), fy = GSL_IMAG (fz);

#ifdef DEBUG
      printf("x = "); gsl_ieee_fprintf_double (stdout, &t.x); printf("\n");
      printf("y = "); gsl_ieee_fprintf_double (stdout, &t.y); printf("\n");
      printf("fx = "); gsl_ieee_fprintf_double (stdout, &fx); printf("\n");
      printf("ex = "); gsl_ieee_fprintf_double (stdout, &t.fx); printf("\n");
      printf("fy = "); gsl_ieee_fprintf_double (stdout, &fy); printf("\n");
      printf("ey = "); gsl_ieee_fprintf_double (stdout, &t.fy); printf("\n");
#endif

      gsl_test_rel (fx, t.fx, 10 * GSL_DBL_EPSILON, "%s real part at (%g,%g)", t.name, t.x, t.y);
      gsl_test_rel (fy, t.fy, 10 * GSL_DBL_EPSILON, "%s imag part at (%g,%g)", t.name, t.x, t.y);
      i++;
    }

  exit (gsl_test_summary ());
}

ext: c


text:
/* [response orientation nms filterBank] = steerableDetector3D(image, filterOrder, sigma);
 *
 * (c) Francois Aguet, 30/08/2012 (last modified 09/02/2012).
 *
 * Compilation:
 * Mac/Linux: mex -I/usr/local/include -I../../mex/include /usr/local/lib/libgsl.a /usr/local/lib/libgslcblas.a steerableDetector3D.cpp
 * Windows: mex COMPFLAGS="$COMPFLAGS /TP /MT" -I"..\..\..\extern\mex\include\gsl-1.14" -I"..\..\mex\include" "..\..\..\extern\mex\lib\gsl.lib" "..\..\..\extern\mex\lib\cblas.lib" -output steerableDetector3D steerableDetector3D.cpp
 */

#include <cstring>
#include <gsl_include/gsl_poly.h>
#include <gsl_include/gsl_math.h>
#include <gsl_include/gsl_eigen.h>
#include <gsl_include/gsl_cblas.h>
#include <algorithm>
#include <vector>
#include <fstream>
#include <string>
#include <math.h>
#include "stackutil.h"
#include "../../../v3d_main/jba/newmat11/newmatap.h"
#include "../../../v3d_main/jba/newmat11/newmatio.h"

#include "convolver3D.h"

double zhi_abs(double num);
int swapthree(double& dummya, double& dummyb, double& dummyc);

class Filter {

public:
    Filter(const double voxels[], const int nx, const int ny, const int nz, const int M, const double sigma);
    Filter(const double voxels[], const int nx, const int ny, const int nz, const int M, const double sigma, const double zxRatio);
    ~Filter();
    void init();

    double* getResponse();
    double** getOrientation();
    double* getNMS();

private:
    const double* voxels_;
    double* response_;
    double** orientation_;
    double* nms_;
    int nx_, ny_, nz_;
    int M_;
    double sigma_, sigmaZ_;
    int N_;
    double alpha_, sign_, c_;

    double *gxx_, *gxy_, *gxz_, *gyy_, *gyz_, *gzz_;

    void calculateTemplates();
    double interpResponse(const double x, const double y, const double z);
    static int mirror(const int x, const int nx);
    void normalize(double v[], const int k);
    void cross(const double v1[], const double v2[], double r[]);

    void computeCurveNMS();
    void computeSurfaceNMS();
    void run();
};

ext: h


text:
#ifndef __EMULATOR_INC_
#define __EMULATOR_INC_
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sf.h>

/**
 * the fn ptr to the covariance function, this is the most called function in libEmu
 * you can change this when you setup the optstruct.
 */
double (*covariance_fn)(gsl_vector*, gsl_vector*, gsl_vector*, int, int) ;

void print_matrix(gsl_matrix* m, int nx, int ny);

double covariance_fn_gaussian(gsl_vector *xm, gsl_vector* xn, gsl_vector* thetas, int nthetas, int nparams);

void derivative_l_gauss(gsl_matrix *dCdTheta, gsl_matrix* xmodel,
												double thetaLength, int index, int nmodel_points, int nparams);

// same as above but without clamping on small values
double covariance_fn_gaussian_exact(gsl_vector *xm, gsl_vector* xn, gsl_vector* thetas, int nthetas, int nparams);

double covariance_fn_matern_three(gsl_vector *xm, gsl_vector* xn, gsl_vector* thetas, int nthetas, int nparams);

void derivative_l_matern_three(gsl_matrix *dCdTheta, gsl_matrix* xmodel, double thetaLength,
															  int index, int nmodel_points, int nparams);

double covariance_fn_matern_five(gsl_vector *xm, gsl_vector* xn, gsl_vector* thetas, int nthetas, int nparams);

void derivative_l_matern_five(gsl_matrix *dCdTheta, gsl_matrix* xmodel, double thetaLength,
															int index, int nmodel_points, int nparams);

void makeKVector(gsl_vector* kvector, gsl_matrix *xmodel, gsl_vector *xnew, gsl_vector* thetas, int nmodel_points, int nthetas, int nparams);

void makeKVector_fnptr(gsl_vector* kvector, gsl_matrix *xmodel, gsl_vector *xnew, gsl_vector* thetas, int nmodel_points, int nthetas, int nparams,
											 double (*covariance_fn_ptr)(gsl_vector*, gsl_vector*, gsl_vector*, int, int));

void makeCovMatrix(gsl_matrix *cov_matrix, gsl_matrix *xmodel, gsl_vector* thetas, int nmodel_points, int nthetas, int nparams);

void makeCovMatrix_fnptr(gsl_matrix *cov_matrix, gsl_matrix *xmodel, gsl_vector* thetas, int nmodel_points, int nthetas, int nparams,
												 double (*covariance_fn_ptr)(gsl_vector*, gsl_vector*, gsl_vector*, int, int));

double makeEmulatedMean(gsl_matrix *inverse_cov_matrix, gsl_vector *training_vector, gsl_vector *kplus_vector, gsl_vector* h_vector, gsl_matrix* h_matrix, gsl_vector* beta_vector,  int nmodel_points);

double makeEmulatedVariance(gsl_matrix *inverse_cov_matrix, gsl_vector *kplus_vector, gsl_vector *h_vector, gsl_matrix *h_matrix, double kappa, int nmodel_points, int nregression_fns);

void initialise_new_x(gsl_matrix* new_x, int nparams, int nemulate_points, double emulate_min, double emulate_max);

/* deprecated */
/* double covariance_fn_gaussian_nondiag(gsl_vector* xm, gsl_vector*xn, gsl_vector*thetas, int nthetas, int nparams); */
/* double covariance_fn_matern(gsl_vector *xm, gsl_vector* xn, gsl_vector* thetas, int nthetas, int nparams); */

#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator


text:
/*
 *  vbHmmPc.c
 *  Model-specific core functions for VB-HMM-PC.
 *
 *  Created by OKAMOTO Kenji, SAKO Yasushi and RIKEN
 *  Copyright 2011-2015
 *  Cellular Informatics Laboratory, Advance Science Institute, RIKEN, Japan.
 *  All rights reserved.
 *
 *  Ver. 1.1.0
 *  Last modified on 2015.09.17
 */

#include <math.h>
#include <string.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_pow_int.h>
#include <string.h>
#include "vbHmmPc.h"
#include "rand.h"

#ifdef _OPENMP
#include "omp.h"
#endif

#define  MAX(a,b)  ((a)>(b)?(a):(b))
#define  MIN(a,b)  ((a)<(b)?(a):(b))

//static int isGlobalAnalysis = 0;

void setFunctions_pc(){
    commonFunctions funcs;
    funcs.newModelParameters    = newModelParameters_pc;
    funcs.freeModelParameters   = freeModelParameters_pc;
    funcs.newModelStats         = newModelStats_pc;
    funcs.freeModelStats        = freeModelStats_pc;
    funcs.initializeVbHmm       = initializeVbHmm_pc;
    funcs.pTilde_z1             = pTilde_z1_pc;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_pc;
    funcs.pTilde_xn_zn          = pTilde_xn_zn_pc;
    funcs.calcStatsVars         = calcStatsVars_pc;
    funcs.maximization          = maximization_pc;
    funcs.varLowerBound         = varLowerBound_pc;
    funcs.reorderParameters     = reorderParameters_pc;
    funcs.outputResults         = outputResults_pc;
    setFunctions( funcs );
}

//void setGFunctions_pc(){
//    gCommonFunctions funcs;
//    funcs.newModelParameters    = newModelParameters_pc;
//    funcs.freeModelParameters   = freeModelParameters_pc;
//    funcs.newModelStats         = newModelStats_pc;
//    funcs.freeModelStats        = freeModelStats_pc;
//    funcs.newModelStatsG        = newModelStatsG_pc;
//    funcs.freeModelStatsG       = freeModelStatsG_pc;
//    funcs.initializeVbHmmG      = initializeVbHmmG_pc;
//    funcs.pTilde_z1             = pTilde_z1_pc;
//    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_pc;
//    funcs.pTilde_xn_zn          = pTilde_xn_zn_pc;
//    funcs.calcStatsVarsG        = calcStatsVarsG_pc;
//    funcs.maximizationG         = maximizationG_pc;
//    funcs.varLowerBoundG        = varLowerBoundG_pc;
//    funcs.reorderParametersG    = reorderParametersG_pc;
//    funcs.outputResultsG        = outputResultsG_pc;
//    setGFunctions( funcs );
//    isGlobalAnalysis = 1;
//}

void outputResults_pc( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    outputPcResults( xn, gv, iv, logFP );
}

//void outputResultsG_pc( xns, gv, ivs, logFP )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//FILE *logFP;
//{
//    outputPcResultsG( xns, gv, ivs, logFP );
//}

void *newModelParameters_pc( xn, sNo )
xnDataSet *xn;
int sNo;
{
    int i;
    pcParameters *p = (void*)malloc( sizeof(pcParameters) );

    p->uPiArr = (double*)malloc( sNo * sizeof(double) );
    p->sumUPi = 0.0;
    p->uAMat = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->uAMat[i] = (double*)malloc( sNo * sizeof(double) );
    }
    p->sumUAArr = (double*)malloc( sNo * sizeof(double) );

    p->avgPi = (double *)malloc( sNo * sizeof(double) );
    p->avgLnPi = (double *)malloc( sNo * sizeof(double) );
    p->avgA = (double **)malloc( sNo * sizeof(double*) );
    p->avgLnA = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->avgA[i] = (double *)malloc( sNo * sizeof(double) );
        p->avgLnA[i] = (double *)malloc( sNo * sizeof(double) );
    }

    p->avgI = (double *)malloc( sNo * sizeof(double) );
    p->avgLnI = (double *)malloc( sNo * sizeof(double) );

    p->aIArr = (double*)malloc( sNo * sizeof(double) );
    p->bIArr = (double*)malloc( sNo * sizeof(double) );

    return p;
}

void freeModelParameters_pc( p, xn, sNo )
void **p;
xnDataSet *xn;
int sNo;
{
    pcParameters *gp = *p;
    int i;

    free( gp->uPiArr );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->uAMat[i] );
    }
    free( gp->uAMat );
    free( gp->sumUAArr );

    free( gp->avgPi );
    free( gp->avgLnPi );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->avgA[i] );
        free( gp->avgLnA[i] );
    }
    free( gp->avgA );
    free( gp->avgLnA );

    free( gp->avgI );
    free( gp->avgLnI );

    free( gp->aIArr );
    free( gp->bIArr );

    free( *p );
    *p = NULL;
}

void *newModelStats_pc( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
//    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        pcStats *s = (pcStats*)malloc( sizeof(pcStats) );

        int i;
        s->Ni = (double *)malloc( sNo * sizeof(double) );
        s->Nij = (double **)malloc( sNo * sizeof(double*) );
        for( i = 0 ; i < sNo ; i++ )
        {   s->Nij[i] = (double *)malloc( sNo * sizeof(double) );   }
        s->Ci = (double *)malloc( sNo * sizeof(double) );
        s->Mi = (double *)malloc( sNo * sizeof(double) );

        return s;

//    } else {
//
//        return NULL;
//
//    }
}

void freeModelStats_pc( s, xn, gv, iv )
void **s;
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
//    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        pcStats *gs = *s;
        int i;
        free( gs->Ni );
        for( i = 0 ; i < sNo ; i++ )
        {   free( gs->Nij[i] );   }
        free( gs->Nij );
        free( gs->Ci );
        free( gs->Mi );

        free( gs );
        *s = NULL;
//    }
}

//void *newModelStatsG_pc( xns, gv, ivs)
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//    int sNo = gv->sNo;
//    pcGlobalStats *gs = (pcGlobalStats*)malloc( sizeof(pcGlobalStats) );
//
//    return gs;
//}

//void freeModelStatsG_pc( gs, xns, gv, ivs )
//void **gs;
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//    int sNo = gv->sNo;
//    pcGlobalStats *ggs = *gs;
//
//    free( *gs );
//    *gs = NULL;
//}

void initializeVbHmm_pc( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    pcData *d = xn->data;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    pcParameters *p = gv->params;

    int i, j;
    size_t totalC = 0;
    for( i = 0 ; i < dLen ; i++ ){
        totalC += d->counts[i];
    }
    double meanI = (double)totalC / (double)dLen;

    // hyperparameter for p( pi(k) )
    p->sumUPi = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        p->uPiArr[i] = 1.0;
        p->sumUPi += p->uPiArr[i];
    }

    // hyperparameter for p( A(i,j) )
    for( i = 0 ; i < sNo ; i++ ){
        p->sumUAArr[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            if( j == i ){
                p->uAMat[i][j] = 100.0;
            } else {
                p->uAMat[i][j] = 1.0;
            }
            p->sumUAArr[i] += p->uAMat[i][j];
        }
    }

    // hyperparameter for p( I(k) )
    for( i = 0 ; i < sNo ; i++ ){
        p->aIArr[i] = 1.0;
        p->bIArr[i] = 1.0 / meanI;
    }

    initialize_indVars_pc( xn, gv, iv );

    calcStatsVars_pc( xn, gv, iv );
    maximization_pc( xn, gv, iv );
}

//void initializeVbHmmG_pc( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void initialize_indVars_pc( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;

    int i;
    size_t n;
    double sumPar;
    for( n = 0 ; n < dLen ; n++ ){
        sumPar = 0.0;
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] = enoise(1.0) + 1.0;
            sumPar += gmMat[n][i];
        }
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] /= sumPar;
        }
    }
}

xnDataSet *newXnDataSet_pc( filename )
const char *filename;
{
    xnDataSet *xn = (xnDataSet*)malloc( sizeof(xnDataSet) );
    xn->name = (char*)malloc( strlen(filename) + 2 );
    strncpy( xn->name, filename, strlen(filename)+1 );
    xn->data = (pcData*)malloc( sizeof(pcData) );
    pcData *d = (pcData*)xn->data;
    d->binSize = 0.0;
    d->counts = NULL;
    return xn;
}

void freeXnDataSet_pc( xn )
xnDataSet **xn;
{
    pcData *d = (pcData*)(*xn)->data;
    free( d->counts );
    free( (*xn)->data );
    free( (*xn)->name );
    free( *xn );
    *xn = NULL;
}

double pTilde_z1_pc( i, params )
int i;
void *params;
{
    pcParameters *p = (pcParameters*)params;
    return exp( p->avgLnPi[i] );
}

double pTilde_zn_zn1_pc( i, j, params )
int i, j;
void *params;
{
    pcParameters *p = (pcParameters*)params;
    return exp( p->avgLnA[i][j] );
}

double pTilde_xn_zn_pc( xnWv, n, i, params )
xnDataSet *xnWv;
size_t n;
int i;
void *params;
{
    pcParameters *p = (pcParameters*)params;
    pcData *xn = (pcData*)xnWv->data;
    return exp( p->avgLnI[i] * (double)xn->counts[n] - p->avgI[i] ) / gsl_sf_fact( xn->counts[n] );
}

void calcStatsVars_pc( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    pcData *d = (pcData*)xn->data;
    pcStats *s = (pcStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *Ni = s->Ni, **Nij = s->Nij;
    double *Mi = s->Mi, *Ci = s->Ci;
    size_t n;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        Ni[i] = 1e-10;
        Ci[i] = 1e-10;
        Mi[i] = 1e-10;
        for( j = 0 ; j < sNo ; j++ ){
            Nij[i][j] = 1e-10;
        }
    }
    for( n = 0 ; n < dLen ; n++ ){
        for( i = 0 ; i < sNo ; i++ ){
            Ni[i]  += gmMat[n][i];
            Ci[i]  += gmMat[n][i] * (double)d->counts[n];
            for( j = 0 ; j < sNo ; j++ ){
                Mi[i]     += xiMat[n][i][j];
                Nij[i][j] += xiMat[n][i][j];
            }
        }
    }
}

//void calcStatsVarsG_pc( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void maximization_pc( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    pcParameters *p = (pcParameters*)gv->params;
    pcStats *s = (pcStats*)iv->stats;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgI = p->avgI, *avgLnI = p->avgLnI;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi, *aIArr = p->aIArr, *bIArr = p->bIArr;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *Ni = s->Ni, **Nij = s->Nij;
    double *Mi = s->Mi, *Ci = s->Ci;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        avgPi[i] = ( uPiArr[i] + gmMat[0][i] ) / ( sumUPi + 1.0 );
        avgLnPi[i] = gsl_sf_psi( uPiArr[i] + gmMat[0][i] ) - gsl_sf_psi( sumUPi + 1.0 );

        for( j = 0 ; j < sNo ; j++ ){
            avgA[i][j] = ( uAMat[i][j] + Nij[i][j] ) / ( sumUAArr[i] + Mi[i] );
            avgLnA[i][j] = gsl_sf_psi( uAMat[i][j] + Nij[i][j] ) - gsl_sf_psi( sumUAArr[i] + Mi[i] );
        }

        avgI[i] = (Ci[i] + aIArr[i]) / (Ni[i] + bIArr[i]);
        avgLnI[i] = gsl_sf_psi( Ci[i] + aIArr[i] ) - log( Ni[i] + bIArr[i] );
#ifdef  INTENSITY_CAP
        size_t n, totalC = 0;
        pcData *pc = xnWv->data;
        for( n = 0 ; n < xnWv->N ; n++ ){
            totalC += pc->counts[n];
        }
        double meanI = (double)totalC / (double)xnWv->N;
        avgI[i] = MIN( avgI[i], maxIntensityRatio * meanI );
        avgLnI[i] = MIN( avgLnI[i], log(maxIntensityRatio * meanI) );
#endif
    }
}

//void maximizationG_pc( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

double varLowerBound_pc( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    pcParameters *p = (pcParameters*)gv->params;
    pcStats *s = (pcStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, *cn = iv->cn;
    double *avgLnPi = p->avgLnPi, **avgLnA = p->avgLnA;
    double *avgI = p->avgI, *avgLnI = p->avgLnI;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi, *aIArr = p->aIArr, *bIArr = p->bIArr;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *Ni = s->Ni, **Nij = s->Nij;
    double *Mi = s->Mi, *Ci = s->Ci;
    size_t n;
    int i, j;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double lnpA = 0.0;
    double lnpI = 0.0;
    double lnqPi = gsl_sf_lngamma(sumUPi + 1.0);
    double lnqA = 0.0;
    double lnqI = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        lnpPi += (uPiArr[i]-1.0) * avgLnPi[i] - gsl_sf_lngamma(uPiArr[i]);

        lnpI += aIArr[i] * log(bIArr[i]) - gsl_sf_lngamma(aIArr[i]);
        lnpI += (aIArr[i] - 1.0) * avgLnI[i] - bIArr[i] * avgI[i];

        lnqPi += (uPiArr[i]+gmMat[0][i]-1.0) * (gsl_sf_psi(uPiArr[i]+gmMat[0][i]) - gsl_sf_psi(sumUPi+1.0));
        lnqPi -= gsl_sf_lngamma(uPiArr[i] + gmMat[0][i]);

        lnqI += (Ci[i] + aIArr[i]) * log(Ni[i] + bIArr[i]) - gsl_sf_lngamma(Ci[i] + aIArr[i]);
        lnqI += (Ci[i] + aIArr[i] - 1.0) * avgLnI[i] - (Ni[i] + bIArr[i]) * avgI[i];

        lnpA += gsl_sf_lngamma(sumUAArr[i]);
        lnqA += gsl_sf_lngamma(sumUAArr[i] + Mi[i]);
        for( j = 0 ; j < sNo ; j++ ){
            lnpA += (uAMat[i][j]-1.0) * avgLnA[i][j] - gsl_sf_lngamma(uAMat[i][j]);

            lnqA += (uAMat[i][j] + Nij[i][j] - 1.0) * (gsl_sf_psi(uAMat[i][j]+Nij[i][j]) - gsl_sf_psi( sumUAArr[i]+Mi[i]) );
            lnqA -= gsl_sf_lngamma( uAMat[i][j] + Nij[i][j] );
        }
    }

    double lnpX = 0.0;
    for( n = 0 ; n < dLen ; n++ ){
        lnpX += log( cn[n] );
    }

    double val;
    val  = lnpPi + lnpA + lnpI;
    val -= lnqPi + lnqA + lnqI;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

//double varLowerBoundG_pc( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void reorderParameters_pc( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    pcParameters *p = (pcParameters*)gv->params;
    pcStats *s = (pcStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgI = p->avgI, *avgLnI = p->avgLnI;
    double *Ni = s->Ni, *Ci = s->Ci;
    size_t n;
    int i, j;

    int *index = (int*)malloc( sNo * sizeof(int) );
    double *store = (double*)malloc( sNo * sizeof(double) );
    double **s2D = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   s2D[i] = (double*)malloc( sNo * sizeof(double) );   }

    // index indicates order of avgI values (0=biggest avgI -- sNo=smallest avgI).
    for( i = 0 ; i < sNo ; i++ ){
        index[i] = sNo - 1;
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                if( avgI[i] < avgI[j] ){
                    index[i]--;
                } else if( avgI[i] == avgI[j] ){
                    if( i < j )
                    {   index[i]--;   }
                }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgI[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgI[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnI[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnI[i] = store[i];   }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgA[i][j] = s2D[i][j];   }
    }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgLnA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgLnA[i][j] = s2D[i][j];   }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = Ni[i];   }
    for( i = 0 ; i < sNo ; i++ ){   Ni[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = Ci[i];   }
    for( i = 0 ; i < sNo ; i++ ){   Ci[i] = store[i];   }

    for( n = 0 ; n < dLen ; n++ ){
        for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = gmMat[n][i];   }
        for( i = 0 ; i < sNo ; i++ ){   gmMat[n][i] = store[i];   }
    }

    for( n = 0 ; n < dLen ; n++ ){
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = xiMat[n][i][j];   }
        }
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   xiMat[n][i][j] = s2D[i][j];   }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   free( s2D[i] );   }
    free( s2D );
    free( store );
    free( index );
}

//void reorderParametersG_pc( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void outputPcResults( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    pcParameters *p = (pcParameters*)gv->params;
    int sNo = gv->sNo;

    int i, j;
    fprintf(logFP, "  results: K = %d \n", sNo);

    fprintf(logFP, "   intensity: ( %g", p->avgI[0]);
    for( i = 1 ; i < sNo ; i++ )
    {   fprintf(logFP, ", %g", p->avgI[i]);   }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   pi: ( %g", p->avgPi[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgPi[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   A_matrix: [");
    for( i = 0 ; i < sNo ; i++ ){
        fprintf(logFP, " ( %g", p->avgA[i][0]);
        for( j = 1 ; j < sNo ; j++ )
        {   fprintf(logFP, ", %g", p->avgA[i][j]);   }
        fprintf(logFP, ")");
    }
    fprintf(logFP, " ] \n\n");

    char fn[256];
    FILE *fp;
    size_t n;

    sprintf( fn, "%s.param%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        fprintf(fp, "I, pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", A%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g, %g", p->avgI[i], p->avgPi[i]);
            for( j = 0 ; j < sNo ; j++ )
            {   fprintf(fp, ", %g", p->avgA[j][i]);   }
            fprintf(fp, "\n");
        }
        fclose(fp);
    }

    sprintf( fn, "%s.Lq%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < gv->iteration ; n++ ){
            fprintf( fp, "%24.20e\n", gv->LqArr[n] );
        }
        fclose(fp);
    }

    sprintf( fn, "%s.maxS%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < xn->N ; n++ ){
            fprintf( fp, "%d\n", iv->stateTraj[n] );
        }
        fclose(fp);
    }

}

//void outputPcResultsG( xns, gv, ivs, logFP )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//FILE *logFP;
//{
//}

//

max_stars_repo_name: okamoto-kenji/varBayes-HMM

max_issues_repo_name: okamoto-kenji/varBayes-HMM

max_forks_repo_name: okamoto-kenji/varBayes-HMM


text:
static char help[] = "Newton's method for arctan x = 0.  Run with -snes_fd or -snes_mf.\n\n";

#include <petsc.h>

extern PetscErrorCode FormFunction(SNES, Vec, Vec, void*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    SNES      snes;          // nonlinear solver
    Vec       x, r;          // solution, residual vectors
    PetscReal x0 = 2.0;

    PetscInitialize(&argc,&argv,(char*)0,help);

    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","options to atan","");CHKERRQ(ierr);
    ierr = PetscOptionsReal("-x0","initial value","atan.c",x0,&x0,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd();CHKERRQ(ierr);

    ierr = VecCreate(PETSC_COMM_WORLD,&x); CHKERRQ(ierr);
    ierr = VecSetSizes(x,PETSC_DECIDE,1); CHKERRQ(ierr);
    ierr = VecSetFromOptions(x); CHKERRQ(ierr);
    ierr = VecSet(x,x0); CHKERRQ(ierr);
    ierr = VecDuplicate(x,&r); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetFunction(snes,r,FormFunction,NULL); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
    ierr = SNESSolve(snes,NULL,x); CHKERRQ(ierr);
    ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

    VecDestroy(&x);  VecDestroy(&r);  SNESDestroy(&snes);
    PetscFinalize();
    return 0;
}

PetscErrorCode FormFunction(SNES snes, Vec x, Vec F, void *ctx) {
    PetscErrorCode ierr;
    const PetscReal  *ax;
    PetscReal        *aF;

    ierr = VecGetArrayRead(x,&ax);CHKERRQ(ierr);
    ierr = VecGetArray(F,&aF);CHKERRQ(ierr);
    aF[0] = atan(ax[0]);
    ierr = VecRestoreArrayRead(x,&ax);CHKERRQ(ierr);
    ierr = VecRestoreArray(F,&aF);CHKERRQ(ierr);
    return 0;
}


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

#ifndef __GSL_WAVELET2D_H__
#define __GSL_WAVELET2D_H__
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_wavelet.h>
#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

GSL_EXPORT int gsl_wavelet2d_transform (const gsl_wavelet * w,
                                        double *data,
                                        size_t tda, size_t size1, size_t size2,
                                        gsl_wavelet_direction dir,
                                        gsl_wavelet_workspace * work);

GSL_EXPORT int gsl_wavelet2d_transform_forward (const gsl_wavelet * w,
                                                double *data,
                                                size_t tda, size_t size1, size_t size2,
                                                gsl_wavelet_workspace * work);

GSL_EXPORT int gsl_wavelet2d_transform_inverse (const gsl_wavelet * w,
                                                double *data,
                                                size_t tda, size_t size1, size_t size2,
                                                gsl_wavelet_workspace * work);

GSL_EXPORT int gsl_wavelet2d_nstransform (const gsl_wavelet * w,
                                          double *data,
                                          size_t tda, size_t size1, size_t size2,
                                          gsl_wavelet_direction dir,
                                          gsl_wavelet_workspace * work);

GSL_EXPORT int gsl_wavelet2d_nstransform_forward (const gsl_wavelet * w,
                                                  double *data,
                                                  size_t tda, size_t size1, size_t size2,
                                                  gsl_wavelet_workspace * work);

GSL_EXPORT int gsl_wavelet2d_nstransform_inverse (const gsl_wavelet * w,
                                                  double *data,
                                                  size_t tda, size_t size1, size_t size2,
                                                  gsl_wavelet_workspace * work);

GSL_EXPORT
int
gsl_wavelet2d_transform_matrix (const gsl_wavelet * w,
                                gsl_matrix * a,
                                gsl_wavelet_direction dir,
                                gsl_wavelet_workspace * work);

GSL_EXPORT
int
gsl_wavelet2d_transform_matrix_forward (const gsl_wavelet * w,
                                        gsl_matrix * a,
                                        gsl_wavelet_workspace * work);

GSL_EXPORT
int
gsl_wavelet2d_transform_matrix_inverse (const gsl_wavelet * w,
                                        gsl_matrix * a,
                                        gsl_wavelet_workspace * work);

GSL_EXPORT
int
gsl_wavelet2d_nstransform_matrix (const gsl_wavelet * w,
                                  gsl_matrix * a,
                                  gsl_wavelet_direction dir,
                                  gsl_wavelet_workspace * work);

GSL_EXPORT
int
gsl_wavelet2d_nstransform_matrix_forward (const gsl_wavelet * w,
                                          gsl_matrix * a,
                                          gsl_wavelet_workspace * work);

GSL_EXPORT
int
gsl_wavelet2d_nstransform_matrix_inverse (const gsl_wavelet * w,
                                          gsl_matrix * a,
                                          gsl_wavelet_workspace * work);

__END_DECLS

#endif /* __GSL_WAVELET2D_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program                         */
/*          GCG --- Generic Column Generation                                */
/*                  a Dantzig-Wolfe decomposition based extension            */
/*                  of the branch-cut-and-price framework                    */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/* Copyright (C) 2010-2018 Operations Research, RWTH Aachen University       */
/*                         Zuse Institute Berlin (ZIB)                       */
/*                                                                           */
/* This program is free software; you can redistribute it and/or             */
/* modify it under the terms of the GNU Lesser 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 Lesser General Public License for more details.                       */
/*                                                                           */
/* You should have received a copy of the GNU Lesser General Public License  */
/* along with this program; if not, write to the Free Software               */
/* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.*/
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/**@file   sepa_basis.c
 * @brief  basis separator
 * @author Jonas Witt
 */

/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/

/*#define SCIP_DEBUG*/

#include <assert.h>
#include <stdio.h>
#include <string.h>

#include "scip/scip.h"
#include "scip/scipdefplugins.h"
#include "sepa_basis.h"
#include "sepa_master.h"
#include "gcg.h"
#include "relax_gcg.h"
#include "pricer_gcg.h"
#include "pub_gcgvar.h"

#ifdef GSL
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#endif

#define SEPA_NAME              "basis"
#define SEPA_DESC              "separator calculates a basis of the orig problem to generate cuts, which cut off the master lp sol"
#define SEPA_PRIORITY                100
#define SEPA_FREQ                     0
#define SEPA_MAXBOUNDDIST           1.0
#define SEPA_USESSUBSCIP           FALSE /**< does the separator use a secondary SCIP instance? */
#define SEPA_DELAY                FALSE /**< should separation method be delayed, if other separators found cuts? */

#define STARTMAXCUTS 50       /**< maximal cuts used at the beginning */
#define MAXCUTSINC   20       /**< increase of allowed number of cuts */

/*
 * Data structures
 */

/** separator data */
struct SCIP_SepaData
{
   SCIP_ROW**            mastercuts;         /**< cuts in the master problem */
   SCIP_ROW**            origcuts;           /**< cuts in the original problem */
   int                   norigcuts;          /**< number of cuts in the original problem */
   int                   nmastercuts;        /**< number of cuts in the master problem */
   int                   maxcuts;            /**< maximal number of allowed cuts */
   SCIP_ROW**            newcuts;            /**< new cuts to tighten original problem */
   int                   nnewcuts;           /**< number of new cuts */
   int                   maxnewcuts;         /**< maximal number of allowed new cuts */
   SCIP_ROW*             objrow;             /**< row with obj coefficients */
   SCIP_Bool             enable;             /**< parameter returns if basis separator is enabled */
   SCIP_Bool             enableobj;          /**< parameter returns if objective constraint is enabled */
   SCIP_Bool             enableobjround;     /**< parameter returns if rhs/lhs of objective constraint is rounded, when obj is int */
   SCIP_Bool             enableppcuts;       /**< parameter returns if cuts generated during pricing are added to newconss array */
   SCIP_Bool             enableppobjconss;   /**< parameter returns if objective constraint for each redcost of pp is enabled */
   SCIP_Bool             enableppobjcg;      /**< parameter returns if objective constraint for each redcost of pp is enabled during pricing */
   int                   separationsetting;  /**< parameter returns which parameter setting is used for separation */
   SCIP_Bool             chgobj;             /**< parameter returns if basis is searched with different objective */
   SCIP_Bool             chgobjallways;      /**< parameter returns if obj is not only changed in first iteration */
   SCIP_Bool             genobjconvex;       /**< parameter returns if objconvex is generated dynamically */
   SCIP_Bool             enableposslack;     /**< parameter returns if positive slack should influence the probing objective function */
   SCIP_Bool             forcecuts;          /**< parameter returns if cuts are forced to enter the LP */
   int                   posslackexp;        /**< parameter returns exponent of usage of positive slack */
   SCIP_Bool             posslackexpgen;     /**< parameter returns if exponent should be automatically generated */
   SCIP_Real             posslackexpgenfactor; /**< parameter returns factor for automatically generated exponent */
   int                   iterations;         /**< parameter returns number of new rows adding iterations (rows just cut off probing lp sol) */
   int                   mincuts;            /**< parameter returns number of minimum cuts needed to return *result = SCIP_Separated */
   SCIP_Real             objconvex;          /**< parameter return convex combination factor */
};

/*
 * Local methods
 */

/** allocates enough memory to hold more cuts */
static
SCIP_RETCODE ensureSizeCuts(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data data structure */
   int                   size                /**< new size of cut arrays */
   )
{
   assert(scip != NULL);
   assert(sepadata != NULL);
   assert(sepadata->mastercuts != NULL);
   assert(sepadata->origcuts != NULL);
   assert(sepadata->norigcuts <= sepadata->maxcuts);
   assert(sepadata->norigcuts >= 0);
   assert(sepadata->nmastercuts <= sepadata->maxcuts);
   assert(sepadata->nmastercuts >= 0);

   if( sepadata->maxcuts < size )
   {
      while ( sepadata->maxcuts < size )
      {
         sepadata->maxcuts += MAXCUTSINC;
      }
      SCIP_CALL( SCIPreallocMemoryArray(scip, &(sepadata->mastercuts), sepadata->maxcuts) );
      SCIP_CALL( SCIPreallocMemoryArray(scip, &(sepadata->origcuts), sepadata->maxcuts) );
   }
   assert(sepadata->maxcuts >= size);

   return SCIP_OKAY;
}

/** allocates enough memory to hold more new cuts */
static
SCIP_RETCODE ensureSizeNewCuts(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_SEPADATA*        sepadata,           /**< separator data data structure */
   int                   size                /**< new size of cut arrays */
   )
{
   assert(scip != NULL);
   assert(sepadata != NULL);
   assert(sepadata->newcuts != NULL);
   assert(sepadata->nnewcuts <= sepadata->maxnewcuts);
   assert(sepadata->nnewcuts >= 0);

   if( sepadata->maxnewcuts < size )
   {
      while ( sepadata->maxnewcuts < size )
      {
         sepadata->maxnewcuts += MAXCUTSINC;
      }
      SCIP_CALL( SCIPreallocMemoryArray(scip, &(sepadata->newcuts), sepadata->maxnewcuts) );
   }
   assert(sepadata->maxnewcuts >= size);

   return SCIP_OKAY;
}

/** returns the result of the exponentiation for given exponent and basis (basis^exponent) */
static
SCIP_Real exponentiate(
   SCIP_Real            basis,               /**< basis for exponentiation */
   int                  exponent            /**< exponent for exponentiation */
   )
{
   SCIP_Real result;
   int i;

   assert(exponent >= 0);

   result = 1.0;
   for( i = 0; i < exponent; ++i )
   {
      result *= basis;
   }

   return result;
}

/**< Initialize probing objective coefficient for each variable with original objective. */
static
SCIP_RETCODE initProbingObjWithOrigObj(
   SCIP*                origscip,           /**< orig scip problem */
   SCIP_Bool            enableobj,          /**< returns if objective row was added to the lp */
   SCIP_Real            objfactor           /**< factor, the objective is multiplied with */
)
{
   SCIP_VAR** origvars;
   int norigvars;
   SCIP_VAR* origvar;

   SCIP_Real newobj;
   int i;

   assert(SCIPinProbing(origscip));

   origvars = SCIPgetVars(origscip);
   norigvars = SCIPgetNVars(origscip);

   /** loop over original variables */
   for( i = 0; i < norigvars; ++i )
   {
      /* get variable information */
      origvar = origvars[i];
      newobj = 0.0;

      /* if objective row is enabled consider also the original objective value */
      if( enableobj )
         newobj = objfactor * SCIPvarGetObj(origvar);

      SCIP_CALL( SCIPchgVarObjProbing(origscip, origvar, newobj) );
   }
   return SCIP_OKAY;
}

/**< Change probing objective coefficient for each variable by adding original objective
 *   to the probing objective.
 */
static
SCIP_RETCODE chgProbingObjAddingOrigObj(
   SCIP*                origscip,           /**< orig scip problem */
   SCIP_Real            objfactor,          /**< factor the additional part of the objective is multiplied with */
   SCIP_Real            objdivisor          /**< factor the additional part of the objective is divided with */
)
{
   SCIP_VAR** origvars;
   int norigvars;
   SCIP_VAR* origvar;

   SCIP_Real newobj;
   int i;

   assert(SCIPinProbing(origscip));

   origvars = SCIPgetVars(origscip);
   norigvars = SCIPgetNVars(origscip);

   /** loop over original variables */
   for( i = 0; i < norigvars; ++i )
   {
      /* get variable information */
      origvar = origvars[i];

      newobj = SCIPgetVarObjProbing(origscip, origvar) + (objfactor * SCIPvarGetObj(origvar))/ objdivisor ;

      SCIP_CALL( SCIPchgVarObjProbing(origscip, origvar, newobj) );
   }
   return SCIP_OKAY;
}

/**< Initialize probing objective coefficient for each variable depending on the current origsol.
 *
 *   If variable is at upper bound set objective to -1, if variable is at lower bound set obj to 1,
 *   else set obj to 0.
 *   Additionally, add original objective to the probing objective if this is enabled.
 */
static
SCIP_RETCODE initProbingObjUsingVarBounds(
   SCIP*                origscip,           /**< orig scip problem */
   SCIP_SEPADATA*       sepadata,           /**< separator specific data */
   SCIP_SOL*            origsol,            /**< orig solution */
   SCIP_Bool            enableobj,          /**< returns if objective row was added to the lp */
   SCIP_Real            objfactor           /**< factor the objective is multiplied with */
)
{
   SCIP_Bool enableposslack;
   int posslackexp;

   SCIP_VAR** origvars;
   int norigvars;
   SCIP_VAR* origvar;

   SCIP_Real lb;
   SCIP_Real ub;
   SCIP_Real solval;
   SCIP_Real newobj;
   SCIP_Real distance;

   int i;

   origvars = SCIPgetVars(origscip);
   norigvars = SCIPgetNVars(origscip);

   enableposslack = sepadata->enableposslack;
   posslackexp = sepadata->posslackexp;

   /** loop over original variables */
   for( i = 0; i < norigvars; ++i )
   {
      /* get variable information */
      origvar = origvars[i];
      lb = SCIPvarGetLbLocal(origvar);
      ub = SCIPvarGetUbLocal(origvar);
      solval = SCIPgetSolVal(origscip, origsol, origvar);

      assert(SCIPisFeasLE(origscip, solval, ub));
      assert(SCIPisFeasGE(origscip, solval, lb));

      /* if solution value of variable is at ub or lb initialize objective value of the variable
       * such that the difference to this bound is minimized
       */
      if( SCIPisFeasEQ(origscip, lb, ub) )
      {
         newobj = 0.0;
      }
      else if( SCIPisLT(origscip, ub, SCIPinfinity(origscip)) && SCIPisFeasLE(origscip, ub, solval) )
      {
         newobj = -1.0;
      }
      else if( SCIPisGT(origscip, lb, -SCIPinfinity(origscip)) && SCIPisFeasGE(origscip, lb, solval) )
      {
         newobj = 1.0;
      }
      else if( enableposslack )
      {
         /* compute distance from solution to variable bound */
         distance = MIN(solval - lb, ub - solval);

         assert(SCIPisFeasPositive(origscip, distance));

         /* check if distance is lower than 1 and compute factor */
         if( SCIPisLT(origscip, distance, 1.0) )
         {
            newobj = exponentiate(MAX(0.0, 1.0 - distance), posslackexp);

            /* check if algebraic sign has to be changed */
            if( SCIPisLT(origscip, distance, solval - lb) )
               newobj = -newobj;
         }
         else
         {
            newobj = 0.0;
         }
      }
      else
      {
         newobj = 0.0;
      }

      newobj = newobj * (int) SCIPgetObjsense(origscip);

      /* if objective row is enabled consider also the original objective value */
      if( enableobj )
         newobj = newobj + SCIPvarGetObj(origvar);

      SCIP_CALL( SCIPchgVarObjProbing(origscip, origvar, objfactor*newobj) );
   }

   return SCIP_OKAY;
}

/**< Change probing objective depending on the current origsol.
 *
 * Loop over all constraints lhs <= sum a_i*x_i <= rhs. If lhs == sum a_i*x_i^* add a_i to objective
 * of variable i and if rhs == sum a_i*x_i^* add -a_i to objective of variable i.
 */
static
SCIP_RETCODE chgProbingObjUsingRows(
   SCIP*                origscip,           /**< orig scip problem */
   SCIP_SEPADATA*       sepadata,           /**< separator data */
   SCIP_SOL*            origsol,            /**< orig solution */
   SCIP_Real            objfactor,          /**< factor the objective is multiplied with */
   SCIP_Real            objdivisor          /**< factor the objective is divided with */
)
{
   SCIP_Bool enableposslack;
   int posslackexp;

   SCIP_ROW** rows;
   int nrows;
   SCIP_ROW* row;
   SCIP_Real* vals;
   SCIP_VAR** vars;
   SCIP_COL** cols;
   int nvars;

   SCIP_Real lhs;
   SCIP_Real rhs;
   SCIP_Real* solvals;
   SCIP_Real activity;
   SCIP_Real factor;
   SCIP_Real objadd;
   SCIP_Real obj;
   SCIP_Real norm;
   SCIP_Real distance;

   int i;
   int j;

   rows = SCIPgetLPRows(origscip);
   nrows = SCIPgetNLPRows(origscip);

   enableposslack = sepadata->enableposslack;
   posslackexp = sepadata->posslackexp;

   assert(SCIPinProbing(origscip));

   SCIP_CALL( SCIPallocBufferArray(origscip, &solvals, SCIPgetNVars(origscip)) );
   SCIP_CALL( SCIPallocBufferArray(origscip, &vars, SCIPgetNVars(origscip)) );

   /** loop over constraint and check activity */
   for( i = 0; i < nrows; ++i )
   {
      row = rows[i];
      lhs = SCIProwGetLhs(row);
      rhs = SCIProwGetRhs(row);

      nvars = SCIProwGetNNonz(row);
      if( nvars == 0 || (sepadata->objrow != NULL && strcmp(SCIProwGetName(row),SCIProwGetName(sepadata->objrow)) == 0 ) )
         continue;

      /* get values, variables and solution values */
      vals = SCIProwGetVals(row);
      cols = SCIProwGetCols(row);
      for( j = 0; j < nvars; ++j )
      {
         vars[j] = SCIPcolGetVar(cols[j]);
      }

      activity = SCIPgetRowSolActivity(origscip, row, origsol);

      if( SCIPisFeasEQ(origscip, rhs, lhs) )
      {
         continue;
      }
      if( SCIPisLT(origscip, rhs, SCIPinfinity(origscip)) && SCIPisFeasLE(origscip, rhs, activity) )
      {
         factor = -1.0;
      }
      else if( SCIPisGT(origscip, lhs, -SCIPinfinity(origscip)) && SCIPisFeasGE(origscip, lhs, activity) )
      {
         factor = 1.0;
      }
      else if( enableposslack )
      {
         assert(!(SCIPisInfinity(origscip, rhs) && SCIPisInfinity(origscip, lhs)));
         assert(!(SCIPisInfinity(origscip, activity) && SCIPisInfinity(origscip, -activity)));

         /* compute distance from solution to row */
         if( SCIPisInfinity(origscip, rhs) && SCIPisGT(origscip, lhs, -SCIPinfinity(origscip)) )
            distance = activity - lhs;
         else if( SCIPisInfinity(origscip, lhs) && SCIPisLT(origscip, rhs, SCIPinfinity(origscip)) )
            distance = rhs - activity;
         else
            distance = MIN(activity - lhs, rhs - activity);

         assert(SCIPisFeasPositive(origscip, distance) || !SCIPisCutEfficacious(origscip, origsol, row));

         /* check if distance is lower than 1 and compute factor */
         if( SCIPisLT(origscip, distance, 1.0) )
         {
            factor = exponentiate(MAX(0.0, 1.0 - distance), posslackexp);

            /* check if algebraic sign has to be changed */
            if( SCIPisLT(origscip, distance, activity - lhs) )
               factor = -1.0*factor;
         }
         else
         {
            continue;
         }
      }
      else
      {
         continue;
      }

      norm = SCIProwGetNorm(row);

      /** loop over variables of the constraint and change objective */
      for( j = 0; j < nvars; ++j )
      {
         obj = SCIPgetVarObjProbing(origscip, vars[j]);
         objadd = (factor * vals[j]) / norm;

         objadd = objadd * (int) SCIPgetObjsense(origscip);

         SCIP_CALL( SCIPchgVarObjProbing(origscip, vars[j], obj + (objfactor * objadd) / objdivisor) );
      }
   }

   SCIPfreeBufferArray(origscip, &solvals);
   SCIPfreeBufferArray(origscip, &vars);

   return SCIP_OKAY;
}

#ifdef GSL
/**< Get matrix (including nrows and ncols) of rows that are satisfied with equality by sol */
static
SCIP_RETCODE getEqualityMatrixGsl(
   SCIP*                scip,               /**< SCIP data structure */
   SCIP_SOL*            sol,                /**< solution */
   gsl_matrix**         matrix,             /**< pointer to store equality matrix */
   int*                 nrows,              /**< pointer to store number of rows */
   int*                 ncols,              /**< pointer to store number of columns */
   int*                 prerank             /**< pointer to store preprocessed rank */
)
{
   int* var2col;
   int* delvars;

   int nvar2col;
   int ndelvars;

   SCIP_ROW** lprows;
   int nlprows;

   SCIP_COL** lpcols;
   int nlpcols;

   int i;
   int j;

   *ncols = SCIPgetNLPCols(scip);
   nlprows = SCIPgetNLPRows(scip);
   lprows = SCIPgetLPRows(scip);
   nlpcols = SCIPgetNLPCols(scip);
   lpcols = SCIPgetLPCols(scip);

   *nrows = 0;

   ndelvars = 0;
   nvar2col = 0;

   SCIP_CALL( SCIPallocBufferArray(scip, &var2col, nlpcols) );
   SCIP_CALL( SCIPallocBufferArray(scip, &delvars, nlpcols) );

   /* loop over lp cols and check if it is at one of its bounds */
   for( i = 0; i < nlpcols; ++i )
   {
      SCIP_COL* lpcol;
      SCIP_VAR* lpvar;

      lpcol = lpcols[i];

      lpvar = SCIPcolGetVar(lpcol);

      if( SCIPisEQ(scip, SCIPgetSolVal(scip, sol, lpvar), SCIPcolGetUb(lpcol) )
         || SCIPisEQ(scip, SCIPgetSolVal(scip, sol, lpvar), SCIPcolGetLb(lpcol)) )
      {
         int ind;

         ind = SCIPcolGetIndex(lpcol);

         delvars[ndelvars] = ind;

         ++ndelvars;

         var2col[ind] = -1;
      }
      else
      {
         int ind;

         ind = SCIPcolGetIndex(lpcol);

         var2col[ind] = nvar2col;

         ++nvar2col;
      }
   }

   SCIPsortInt(delvars, ndelvars);

   *matrix = gsl_matrix_calloc(nlprows, nvar2col);

   *ncols = nvar2col;

   /* loop over lp rows and check if solution feasibility is equal to zero */
   for( i = 0; i < nlprows; ++i )
   {
      SCIP_ROW* lprow;

      lprow = lprows[i];

      /* if solution feasiblity is equal to zero, add row to matrix */
      if( SCIPisEQ(scip, SCIPgetRowSolFeasibility(scip, lprow, sol), 0.0) )
      {
         SCIP_COL** cols;
         SCIP_Real* vals;
         int nnonz;

         cols = SCIProwGetCols(lprow);
         vals = SCIProwGetVals(lprow);
         nnonz = SCIProwGetNNonz(lprow);

         /* get nonzero coefficients of row */
         for( j = 0; j < nnonz; ++j )
         {
            int ind;
            int pos;

            ind = SCIPcolGetIndex(cols[j]);
            assert(ind >= 0 && ind < nlpcols);

            if( !SCIPsortedvecFindInt(delvars, ind, ndelvars, &pos) )
            {
               gsl_matrix_set(*matrix, *nrows, var2col[ind], vals[j]);
            }
         }
         ++(*nrows);
      }
   }
   *nrows = nlprows;
   *prerank = ndelvars;

   SCIPfreeBufferArray(scip, &delvars);
   SCIPfreeBufferArray(scip, &var2col);

   return SCIP_OKAY;
}

/**< get the rank of a given matrix */
static
SCIP_RETCODE getRank(
   SCIP*                scip,
   gsl_matrix*          matrix,
   int                  nrows,
   int                  ncols,
   int*                 rank
)
{
   gsl_matrix* matrixq;
   gsl_matrix* matrixr;

   gsl_vector* tau;
   gsl_vector* norm;
   gsl_permutation* permutation;

   int ranktmp;
   int signum;

   int i;

   matrixq = gsl_matrix_alloc(nrows, nrows);
   matrixr = gsl_matrix_alloc(nrows, ncols);

   norm = gsl_vector_alloc(ncols);
   tau = gsl_vector_alloc(MIN(nrows, ncols));

   permutation = gsl_permutation_alloc(ncols);

   gsl_linalg_QRPT_decomp(matrix, tau, permutation, &signum, norm);

   gsl_linalg_QR_unpack(matrix, tau, matrixq, matrixr);

   ranktmp = 0;

   for( i = 0; i < MIN(nrows, ncols); ++i )
   {
      SCIP_Real val;

      val = gsl_matrix_get(matrixr, i, i);

      if( SCIPisZero(scip, val) )
      {
         break;
      }
      ++(ranktmp);
   }

   *rank = ranktmp;

   gsl_matrix_free(matrixq);
   gsl_matrix_free(matrixr);

   gsl_vector_free(tau);
   gsl_vector_free(norm);

   gsl_permutation_free(permutation);

   return SCIP_OKAY;
}

/**< Get rank (number of linear independent rows) of rows that are satisfied
 *   with equality by solution sol */
static
SCIP_RETCODE getEqualityRankGsl(
   SCIP*                scip,               /**< SCIP data structure */
   SCIP_SOL*            sol,                /**< solution */
   int*                 equalityrank        /**< pointer to store rank of rows with equality */
   )
{
   gsl_matrix* matrix;
   int nrows;
   int ncols;

   int prerank;
   int rowrank;

   SCIP_CALL( getEqualityMatrixGsl(scip, sol, &matrix, &nrows, &ncols, &prerank) );

   SCIP_CALL( getRank(scip, matrix, nrows, ncols, &rowrank) );

   gsl_matrix_free(matrix);

   *equalityrank = rowrank + prerank;

   return SCIP_OKAY;
}
#endif

/** add cuts which are due to the latest objective function of the pricing problems
 *  (reduced cost non-negative) */
static
SCIP_RETCODE addPPObjConss(
   SCIP*                scip,               /**< SCIP data structure */
   SCIP_SEPA*           sepa,               /**< separator basis */
   int                  ppnumber,           /**< number of pricing problem */
   SCIP_Real            dualsolconv,        /**< dual solution corresponding to convexity constraint */
   SCIP_Bool            newcuts,            /**< add cut to newcuts in sepadata? (otherwise add it just to the cutpool) */
   SCIP_Bool            probing             /**< add cut to probing LP? */
)
{
   SCIP_SEPADATA* sepadata;

   SCIP* pricingscip;

   SCIP_VAR** pricingvars;
   SCIP_VAR* var;

   int npricingvars;
   int nvars;

   char name[SCIP_MAXSTRLEN];

   int j;
   int k;

   SCIP_OBJSENSE objsense;

   SCIP_Real lhs;
   SCIP_Real rhs;

   sepadata = SCIPsepaGetData(sepa);

   nvars = 0;
   pricingscip = GCGgetPricingprob(scip, ppnumber);
   pricingvars = SCIPgetOrigVars(pricingscip);
   npricingvars = SCIPgetNOrigVars(pricingscip);

   if( !GCGisPricingprobRelevant(scip, ppnumber) || pricingscip == NULL )
      return SCIP_OKAY;

   objsense = SCIPgetObjsense(pricingscip);

   if( objsense == SCIP_OBJSENSE_MINIMIZE )
   {
      lhs = dualsolconv;
      rhs = SCIPinfinity(scip);
   }
   else
   {
      rhs = dualsolconv;
      lhs = -SCIPinfinity(scip);
   }

   for( k = 0; k < GCGgetNIdenticalBlocks(scip, ppnumber); ++k )
   {
      SCIP_ROW* origcut;

      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "newconstraint_%d_%d_%d", SCIPsepaGetNCalls(sepa), ppnumber, k);

      SCIP_CALL( SCIPcreateEmptyRowUnspec(scip, &origcut, name, lhs, rhs, FALSE, FALSE, TRUE) );

      nvars = 0;

      for( j = 0; j < npricingvars ; ++j )
      {
         assert(GCGvarIsPricing(pricingvars[j]));

         if( !SCIPisEQ(scip, SCIPvarGetObj(pricingvars[j]), 0.0) )
         {
            var = GCGpricingVarGetOrigvars(pricingvars[j])[k];
            assert(var != NULL);
            SCIP_CALL( SCIPaddVarToRow(scip, origcut, var, SCIPvarGetObj(pricingvars[j])) );
            ++nvars;
         }
      }

      if( nvars > 0 )
      {
         if( newcuts )
         {
            SCIP_CALL( ensureSizeNewCuts(scip, sepadata, sepadata->nnewcuts + 1) );

            sepadata->newcuts[sepadata->nnewcuts] = origcut;
            SCIP_CALL( SCIPcaptureRow(scip, sepadata->newcuts[sepadata->nnewcuts]) );
            ++(sepadata->nnewcuts);

            SCIPdebugMessage("cut added to new cuts in relaxdata\n");
         }
         else
         {
            SCIP_CALL( SCIPaddPoolCut(scip, origcut) );
            SCIPdebugMessage("cut added to orig cut pool\n");
         }

         if( probing )
         {
            SCIP_CALL( SCIPaddRowProbing(scip, origcut) );
            SCIPdebugMessage("cut added to probing\n");
         }

      }
      SCIP_CALL( SCIPreleaseRow(scip, &origcut) );
   }

   return SCIP_OKAY;
}
/*
 * Callback methods of separator
 */

/** copy method for separator plugins (called when SCIP copies plugins) */
#if 0
static
SCIP_DECL_SEPACOPY(sepaCopyBasis)
{  /*lint --e{715}*/
   SCIPerrorMessage("method of basis separator not implemented yet\n");
   SCIPABORT(); /*lint --e{527}*/

   return SCIP_OKAY;
}
#else
#define sepaCopyBasis NULL
#endif

/** destructor of separator to free user data (called when SCIP is exiting) */
static
SCIP_DECL_SEPAFREE(sepaFreeBasis)
{  /*lint --e{715}*/
   SCIP_SEPADATA* sepadata;

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);

   SCIPfreeMemory(scip, &sepadata);

   return SCIP_OKAY;
}

/** initialization method of separator (called after problem was transformed) */
static
SCIP_DECL_SEPAINIT(sepaInitBasis)
{  /*lint --e{715}*/
   SCIP*   origscip;
   SCIP_SEPADATA* sepadata;

   SCIP_VAR** origvars;
   int norigvars;

   char name[SCIP_MAXSTRLEN];

   SCIP_Real obj;
   int i;

   SCIP_Bool enable;
   SCIP_Bool enableobj;

   assert(scip != NULL);

   origscip = GCGmasterGetOrigprob(scip);
   assert(origscip != NULL);

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);

   origvars = SCIPgetVars(origscip);
   norigvars = SCIPgetNVars(origscip);

   SCIPdebugMessage("sepaInitBasis\n");

   enable = sepadata->enable;
   enableobj = sepadata->enableobj;

   sepadata->maxcuts = STARTMAXCUTS;
   sepadata->norigcuts = 0;
   sepadata->maxnewcuts = 0;
   sepadata->nnewcuts = 0;
   sepadata->objrow = NULL;

   /* if separator is disabled do nothing */
   if( !enable )
   {
      return SCIP_OKAY;
   }

   SCIP_CALL( SCIPallocMemoryArray(scip, &(sepadata->origcuts), STARTMAXCUTS) ); /*lint !e506*/
   SCIP_CALL( SCIPallocMemoryArray(scip, &(sepadata->mastercuts), STARTMAXCUTS) ); /*lint !e506*/
   SCIP_CALL( SCIPallocMemoryArray(scip, &(sepadata->newcuts), STARTMAXCUTS) ); /*lint !e506*/

   /* if objective row is enabled create row with objective coefficients */
   if( enableobj )
   {
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "objrow");
      SCIP_CALL( SCIPcreateEmptyRowUnspec(origscip, &(sepadata->objrow), name, -SCIPinfinity(origscip), SCIPinfinity(origscip), TRUE, FALSE, TRUE) );

      for( i = 0; i < norigvars; ++i )
      {
         obj = SCIPvarGetObj(origvars[i]);
         SCIP_CALL( SCIPaddVarToRow(origscip, sepadata->objrow, origvars[i], obj) );
      }
   }

   return SCIP_OKAY;
}

/** deinitialization method of separator (called before transformed problem is freed) */
static
SCIP_DECL_SEPAEXIT(sepaExitBasis)
{  /*lint --e{715}*/
   SCIP* origscip;
   SCIP_SEPADATA* sepadata;
   SCIP_Bool enableobj;

   int i;

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);
   enableobj = sepadata->enableobj;
   assert(sepadata->nmastercuts == sepadata->norigcuts);

   origscip = GCGmasterGetOrigprob(scip);
   assert(origscip != NULL);

   for( i = 0; i < sepadata->norigcuts; i++ )
   {
      SCIP_CALL( SCIPreleaseRow(origscip, &(sepadata->origcuts[i])) );
   }

   for( i = 0; i < sepadata->nnewcuts; ++i )
   {
      if( sepadata->newcuts[i] != NULL )
         SCIP_CALL( SCIPreleaseRow(origscip, &(sepadata->newcuts[i])) );
   }

   if( enableobj )
      SCIP_CALL( SCIPreleaseRow(origscip, &(sepadata->objrow)) );

   SCIPfreeMemoryArrayNull(scip, &(sepadata->origcuts));
   SCIPfreeMemoryArrayNull(scip, &(sepadata->mastercuts));
   SCIPfreeMemoryArrayNull(scip, &(sepadata->newcuts));

   return SCIP_OKAY;
}

/** solving process initialization method of separator (called when branch and bound process is about to begin) */
static
SCIP_DECL_SEPAINITSOL(sepaInitsolBasis)
{  /*lint --e{715}*/
   SCIP_SEPADATA* sepadata;

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);

   sepadata->nmastercuts = 0;

   return SCIP_OKAY;
}

/** solving process deinitialization method of separator (called before branch and bound process data is freed) */
static
SCIP_DECL_SEPAEXITSOL(sepaExitsolBasis)
{  /*lint --e{715}*/
   SCIP_SEPADATA* sepadata;
   int i;

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);
   assert(sepadata->nmastercuts == sepadata->norigcuts);

   assert(GCGmasterGetOrigprob(scip) != NULL);

   for( i = 0; i < sepadata->nmastercuts; i++ )
   {
      SCIP_CALL( SCIPreleaseRow(scip, &(sepadata->mastercuts[i])) );
   }

   return SCIP_OKAY;
}

/**< Initialize objective due to generation of convex combination */
static
SCIP_RETCODE initGenconv(
   SCIP*                origscip,           /**< original SCIP data structure */
   SCIP_SEPADATA*       sepadata,           /**< separator data structure */
   SCIP_SOL*            origsol,            /**< current original solution */
   int                  nbasis,             /**< rank of constraint matrix */
   SCIP_Real*           convex              /**< pointer to store convex combination coefficient */
)
{  /*lint --e{715}*/
#ifdef GSL
   int rank;

   SCIP_CALL( getEqualityRankGsl(origscip, origsol, &rank) );

   *convex = 1.0* rank/nbasis;

   SCIPdebugMessage("use generic coefficient %d/%d = %f\n", rank, nbasis, *convex);

#else
   SCIPwarningMessage(origscip, "Gnu Scientific Library is not enabled! \n"
      "either set sepa/basis/genobjconvex = FALSE sepa/basis/posslackexpgen = FALSE \n"
      "or compile with GSL=true and include Gnu Scientific Library\n");
   *convex = sepadata->objconvex;
#endif

   return SCIP_OKAY;
}

/**< Initialize objective due to generation of convex combination */
static
SCIP_RETCODE initConvObj(
   SCIP*                origscip,           /**< original SCIP data structure */
   SCIP_SEPADATA*       sepadata,           /**< separator data structure */
   SCIP_SOL*            origsol,            /**< current original solution */
   SCIP_Real            convex,             /**< convex coefficient to initialize objective */
   SCIP_Bool            genericconv         /**< was convex coefficient calculated generically? */
)
{
   SCIP_Real objnormnull;
   SCIP_Real objnormcurrent;

   objnormnull = 1.0;
   objnormcurrent = 1.0;

   if( SCIPisEQ(origscip, convex, 0.0) )
   {
      SCIP_CALL( initProbingObjWithOrigObj(origscip, TRUE, 1.0) );
   }
   else if( SCIPisLT(origscip, convex, 1.0) )
   {
      SCIP_CALL( initProbingObjWithOrigObj(origscip, TRUE, 1.0) );
      objnormnull = SCIPgetObjNorm(origscip);

      SCIP_CALL( initProbingObjUsingVarBounds(origscip, sepadata, origsol, FALSE, convex) );
      SCIP_CALL( chgProbingObjUsingRows(origscip, sepadata, origsol, convex, 1.0) );

      objnormcurrent = SCIPgetObjNorm(origscip)/(convex);

      if( SCIPisEQ(origscip, objnormcurrent, 0.0) )
         SCIP_CALL( initProbingObjWithOrigObj(origscip, TRUE, 1.0) );
      else if( SCIPisGT(origscip, objnormnull, 0.0) )
         SCIP_CALL( chgProbingObjAddingOrigObj(origscip, (1.0 - convex) * objnormcurrent, objnormnull) );
   }
   else if( SCIPisEQ(origscip, convex, 1.0) )
   {
      SCIP_CALL( initProbingObjUsingVarBounds(origscip, sepadata, origsol, !genericconv && sepadata->enableobj, 1.0) );
      SCIP_CALL( chgProbingObjUsingRows(origscip, sepadata, origsol, 1.0, 1.0) );
   }

   return SCIP_OKAY;
}

/** LP solution separation method of separator */
static
SCIP_DECL_SEPAEXECLP(sepaExeclpBasis)
{  /*lint --e{715}*/

   SCIP* origscip;
   SCIP_SEPADATA* sepadata;

   SCIP_SEPA** sepas;
   int nsepas;

   SCIP_ROW** cuts;
   SCIP_ROW* mastercut;
   SCIP_ROW* origcut;
   SCIP_COL** cols;
   SCIP_VAR** roworigvars;
   SCIP_VAR** mastervars;
   SCIP_Real* mastervals;
   int ncols;
   int ncuts;
   SCIP_Real* vals;
   int nmastervars;

   SCIP_RESULT resultdummy;
   SCIP_OBJSENSE objsense;
   SCIP_SOL* origsol;
   SCIP_Bool lperror;
   SCIP_Bool delayed;
   SCIP_Bool cutoff;
   SCIP_Bool infeasible;
   SCIP_Real obj;

   SCIP_Bool enable;
   SCIP_Bool enableobj;
   SCIP_Bool enableobjround;
   SCIP_Bool enableppobjconss;

   char name[SCIP_MAXSTRLEN];

   int i;
   int j;
   int iteration;
   int nbasis;
   int nlprowsstart;
   int nlprows;
   SCIP_ROW** lprows;

   assert(scip != NULL);
   assert(result != NULL);

   origscip = GCGmasterGetOrigprob(scip);
   assert(origscip != NULL);

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);

   SCIPdebugMessage("calling sepaExeclpBasis\n");

   *result = SCIP_DIDNOTFIND;

   enable = sepadata->enable;
   enableobj = sepadata->enableobj;
   enableobjround = sepadata->enableobjround;
   enableppobjconss = sepadata->enableppobjconss;

   /* if separator is disabled do nothing */
   if( !enable )
   {
      SCIPdebugMessage("separator is not enabled\n");
      *result = SCIP_DIDNOTRUN;
      return SCIP_OKAY;
   }

   /* ensure master LP is solved to optimality */
   if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
   {
      SCIPdebugMessage("master LP not solved to optimality, do no separation!\n");
      *result = SCIP_DIDNOTRUN;
      return SCIP_OKAY;
   }

   if( GCGgetNRelPricingprobs(origscip) < GCGgetNPricingprobs(origscip) )
   {
      SCIPdebugMessage("aggregated pricing problems, do no separation!\n");
      *result = SCIP_DIDNOTRUN;
      return SCIP_OKAY;
   }

   if( GCGrelaxIsOrigSolFeasible(origscip) )
   {
      SCIPdebugMessage("Current solution is feasible, no separation necessary!\n");
      *result = SCIP_DIDNOTRUN;
      return SCIP_OKAY;
   }

   /* get current original solution */
   origsol = GCGrelaxGetCurrentOrigSol(origscip);

   /* get obj and objsense */
   objsense = SCIPgetObjsense(origscip);
   obj = SCIPgetSolOrigObj(origscip, origsol);

   /** get number of linearly independent rows needed for basis */
   nbasis = SCIPgetNLPCols(origscip);

   *result = SCIP_DIDNOTFIND;

   /* init iteration counter */
   iteration = 0;

   /* set parameter setting for separation */
   SCIP_CALL( SCIPsetSeparating(origscip, (SCIP_PARAMSETTING) sepadata->separationsetting, TRUE) );

   /* start diving */
   SCIP_CALL( SCIPstartProbing(origscip) );

   SCIP_CALL( SCIPnewProbingNode(origscip) );

   SCIP_CALL( SCIPconstructLP(origscip, &cutoff) );

   /** add origcuts to probing lp */
   for( i = 0; i < GCGsepaGetNCuts(scip); ++i )
   {
      if( SCIProwGetLPPos(GCGsepaGetOrigcuts(scip)[i]) == -1 )
         SCIP_CALL( SCIPaddRowProbing(origscip, GCGsepaGetOrigcuts(scip)[i]) );
   }

   /** add new cuts which did not cut off master sol to probing lp */
   for( i = 0; i < sepadata->nnewcuts; ++i )
   {
      if( SCIProwGetLPPos(sepadata->newcuts[i]) == -1 )
         SCIP_CALL( SCIPaddRowProbing(origscip, sepadata->newcuts[i]) );
   }

   /* store number of lp rows in the beginning */
   nlprowsstart = SCIPgetNLPRows(origscip);

   /* while the counter is smaller than the number of allowed iterations,
    * try to separate origsol via probing lp sol */
   /* TODO: while z*(T) = 0 like Range suggests? But then we have to adjust which cuts are added */
   while( iteration < sepadata->iterations )
   {
      SCIPdebugMessage("iteration %d of at most %d iterations\n", iteration + 1, sepadata->iterations);

      SCIP_CALL( SCIPapplyCutsProbing(origscip, &cutoff) );

      /* add new constraints if this is enabled  */
      if( enableppobjconss && iteration == 0 )
      {
         SCIP_Real* dualsolconv;

         SCIPdebugMessage("add reduced cost cut for relevant pricing problems\n");

         SCIP_CALL( SCIPallocMemoryArray(scip, &dualsolconv, GCGgetNPricingprobs(origscip)) );
         SCIP_CALL( GCGsetPricingObjs(scip, dualsolconv) );

         for( i = 0; i < GCGgetNPricingprobs(origscip); ++i )
         {
            SCIP_CALL( addPPObjConss(origscip, sepa, i, dualsolconv[i], FALSE, TRUE) );
         }

         SCIPfreeMemoryArray(scip, &dualsolconv);
      }

      /* init objective */
      if( sepadata->chgobj && (iteration == 0 || sepadata->chgobjallways) )
      {
         SCIPdebugMessage("initialize objective function\n");
         if( sepadata->genobjconvex )
         {
            SCIP_Real genconvex;

            SCIP_CALL( initGenconv(origscip, sepadata, origsol, nbasis, &genconvex) );

            SCIP_CALL( initConvObj(origscip, sepadata, origsol, genconvex, TRUE) );
         }
         else
         {
            SCIPdebugMessage("use given coefficient %g\n", sepadata->objconvex);

            if( sepadata->enableposslack && sepadata->posslackexpgen )
            {
               SCIP_Real genconvex;
               SCIP_Real factor;

               factor = sepadata->posslackexpgenfactor;

               SCIP_CALL( initGenconv(origscip, sepadata, origsol, nbasis, &genconvex) );

               sepadata->posslackexp = (int) (SCIPceil(origscip, factor/(1.0 - genconvex)) + 0.5);

               SCIPdebugMessage("exponent = %d\n", sepadata->posslackexp);

            }
            SCIP_CALL( initConvObj(origscip, sepadata, origsol, sepadata->objconvex, FALSE) );
         }
      }

      /* update rhs/lhs of objective constraint and add it to probing LP, if it exists (only in first iteration) */
      if( enableobj && iteration == 0 )
      {
         SCIPdebugMessage("initialize original objective cut\n");

         /* round rhs/lhs of objective constraint, if it exists, obj is integral and this is enabled */
         if( SCIPisObjIntegral(origscip) && enableobjround )
         {
            if( objsense == SCIP_OBJSENSE_MAXIMIZE )
            {
               SCIPdebugMessage("round rhs down\n");
               obj = SCIPfloor(origscip, obj);
            }
            else
            {
               SCIPdebugMessage("round lhs up\n");
               obj = SCIPceil(origscip, obj);
            }
         }

         /* update rhs/lhs of objective constraint */
         if( objsense == SCIP_OBJSENSE_MAXIMIZE )
         {
            SCIP_CALL( SCIPchgRowRhs(origscip, sepadata->objrow, obj) );
            SCIP_CALL( SCIPchgRowLhs(origscip, sepadata->objrow, -1.0*SCIPinfinity(origscip)) );
         }
         else
         {
            SCIP_CALL( SCIPchgRowLhs(origscip, sepadata->objrow, obj) );
            SCIP_CALL( SCIPchgRowRhs(origscip, sepadata->objrow, SCIPinfinity(origscip)) );
         }
         SCIPdebugMessage("add original objective cut to probing LP\n");

         /** add row to probing lp */
         SCIP_CALL( SCIPaddRowProbing(origscip, sepadata->objrow) );
      }

      SCIPdebugMessage("solve probing LP\n");

      /* solve probing lp */
      SCIP_CALL( SCIPsolveProbingLP(origscip, -1, &lperror, &cutoff) );

      assert(!lperror);

      /* get separators of origscip */
      sepas = SCIPgetSepas(origscip);
      nsepas = SCIPgetNSepas(origscip);

      SCIPdebugMessage("set parameters of separators\n");

      /* loop over sepas and enable/disable sepa */
      for( i = 0; i < nsepas; ++i )
      {
         const char* sepaname;
         char paramname[SCIP_MAXSTRLEN];

         sepaname = SCIPsepaGetName(sepas[i]);

         (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "separating/%s/freq", sepaname);

         /* disable intobj, closecuts, rapidlearning and cgmip separator*/
         if( strcmp(sepaname, "intobj") == 0 || strcmp(sepaname, "closecuts") == 0
            || strcmp(sepaname, "rapidlearning") == 0
            || (strcmp(sepaname, "cgmip") == 0))
         {
            SCIP_CALL( SCIPsetIntParam(origscip, paramname, -1) );
         }
         else
         {
            SCIP_CALL( SCIPsetIntParam(origscip, paramname, 0) );
         }
      }

      SCIPdebugMessage("separate current LP solution\n");

      /** separate current probing lp sol of origscip */
      SCIP_CALL( SCIPseparateSol(origscip, NULL, TRUE, FALSE, TRUE, &delayed, &cutoff) );

      if( delayed && !cutoff )
      {
         SCIPdebugMessage("call delayed separators\n");

         SCIP_CALL( SCIPseparateSol(origscip, NULL, TRUE, TRUE, TRUE, &delayed, &cutoff) );
      }

      /* if cut off is detected set result pointer and return SCIP_OKAY */
      if( cutoff )
      {
         *result = SCIP_CUTOFF;
         SCIP_CALL( SCIPendProbing(origscip) );

         /* disable separating again */
         SCIP_CALL( SCIPsetSeparating(origscip, SCIP_PARAMSETTING_OFF, TRUE) );

         return SCIP_OKAY;
      }

      SCIPdebugMessage("separate current LP solution and current original solution in cutpool\n");

      /* separate cuts in cutpool */
      SCIP_CALL( SCIPseparateSolCutpool(origscip, SCIPgetGlobalCutpool(origscip), NULL, &resultdummy) );
      SCIP_CALL( SCIPseparateSolCutpool(origscip, SCIPgetDelayedGlobalCutpool(origscip), NULL, &resultdummy) );

      /* separate cuts in cutpool */
      SCIP_CALL( SCIPseparateSolCutpool(origscip, SCIPgetGlobalCutpool(origscip), origsol, &resultdummy) );
      SCIP_CALL( SCIPseparateSolCutpool(origscip, SCIPgetDelayedGlobalCutpool(origscip), origsol, &resultdummy) );

      assert(sepadata->norigcuts == sepadata->nmastercuts);

      SCIPdebugMessage("%d cuts are in the original sepastore!\n", SCIPgetNCuts(origscip));

      /* get separated cuts */
      cuts = SCIPgetCuts(origscip);
      ncuts = SCIPgetNCuts(origscip);

      SCIP_CALL( ensureSizeCuts(scip, sepadata, sepadata->norigcuts + ncuts) );

      mastervars = SCIPgetVars(scip);
      nmastervars = SCIPgetNVars(scip);
      SCIP_CALL( SCIPallocBufferArray(scip, &mastervals, nmastervars) );

      /** loop over cuts and transform cut to master problem (and safe cuts) if it seperates origsol */
      for( i = 0; i < ncuts; i++ )
      {
         SCIP_Bool colvarused;

         colvarused = FALSE;
         origcut = cuts[i];

         /* get columns and vals of the cut */
         ncols = SCIProwGetNNonz(origcut);
         cols = SCIProwGetCols(origcut);
         vals = SCIProwGetVals(origcut);

         /* get the variables corresponding to the columns in the cut */
         SCIP_CALL( SCIPallocBufferArray(scip, &roworigvars, ncols) );
         for( j = 0; j < ncols; j++ )
         {
            roworigvars[j] = SCIPcolGetVar(cols[j]);
            assert(roworigvars[j] != NULL);
            if( !GCGvarIsOriginal(roworigvars[j]) )
            {
               colvarused = TRUE;
               break;
            }
         }

         if( colvarused )
         {
            SCIPwarningMessage(origscip, "colvar used in original cut %s\n", SCIProwGetName(origcut));
            SCIPfreeBufferArray(scip, &roworigvars);
            continue;
         }

         if( !SCIPisCutEfficacious(origscip, origsol, origcut) )
         {
            if( !SCIProwIsLocal(origcut) )
               SCIP_CALL( SCIPaddPoolCut(origscip, origcut) );

            SCIPfreeBufferArray(scip, &roworigvars);

            continue;
         }

         /* add the cut to the original cut storage */
         sepadata->origcuts[sepadata->norigcuts] = origcut;
         SCIP_CALL( SCIPcaptureRow(origscip, sepadata->origcuts[sepadata->norigcuts]) );
         sepadata->norigcuts++;

         /* create new cut in the master problem */
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "mc_basis_%s", SCIProwGetName(origcut));
         SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &mastercut, sepa, name,
               ( SCIPisInfinity(scip, -SCIProwGetLhs(origcut)) ?
                  SCIProwGetLhs(origcut) : SCIProwGetLhs(origcut) - SCIProwGetConstant(origcut)),
               ( SCIPisInfinity(scip, SCIProwGetRhs(origcut)) ?
                  SCIProwGetRhs(origcut) : SCIProwGetRhs(origcut) - SCIProwGetConstant(origcut)),
                  SCIProwIsLocal(origcut), TRUE, FALSE) );

         /* transform the original variables to master variables and add them to the cut */
         GCGtransformOrigvalsToMastervals(origscip, roworigvars, vals, ncols, mastervars, mastervals, nmastervars);
         SCIP_CALL( SCIPaddVarsToRow(scip, mastercut, nmastervars, mastervars, mastervals) );

         /* add the cut to the master problem and to the master cut storage */
         SCIP_CALL( SCIPaddRow(scip, mastercut, sepadata->forcecuts, &infeasible) );
         sepadata->mastercuts[sepadata->nmastercuts] = mastercut;
         SCIP_CALL( SCIPcaptureRow(scip, sepadata->mastercuts[sepadata->nmastercuts]) );
         sepadata->nmastercuts++;
         SCIP_CALL( GCGsepaAddMastercuts(scip, origcut, mastercut) );

         SCIP_CALL( SCIPreleaseRow(scip, &mastercut) );
         SCIPfreeBufferArray(scip, &roworigvars);
      }

      if( SCIPgetNCuts(scip) >= sepadata->mincuts )
      {
         *result = SCIP_SEPARATED;

         iteration = sepadata->iterations;
      }
      else if( SCIPgetNCuts(origscip) == 0 )
      {
         iteration = sepadata->iterations;
      }
      else
      {
         ++iteration;
      }

      SCIPdebugMessage("%d cuts are in the master sepastore!\n", SCIPgetNCuts(scip));

      SCIPfreeBufferArray(scip, &mastervals);

      assert(sepadata->norigcuts == sepadata->nmastercuts );
   }

   SCIP_CALL( SCIPclearCuts(origscip) );

   lprows = SCIPgetLPRows(origscip);
   nlprows = SCIPgetNLPRows(origscip);

   assert(nlprowsstart <= nlprows);

   SCIP_CALL( ensureSizeNewCuts(scip, sepadata, sepadata->nnewcuts + nlprows - nlprowsstart) );

   for( i = nlprowsstart; i < nlprows; ++i )
   {
      if( SCIProwGetOrigintype(lprows[i]) == SCIP_ROWORIGINTYPE_SEPA )
	   {
         sepadata->newcuts[sepadata->nnewcuts] = lprows[i];
	      SCIP_CALL( SCIPcaptureRow(origscip, sepadata->newcuts[sepadata->nnewcuts]) );
	      ++(sepadata->nnewcuts);
	   }
   }

   /* end diving */
   SCIP_CALL( SCIPendProbing(origscip) );

   if( SCIPgetNCuts(scip) > 0 )
   {
      *result = SCIP_SEPARATED;
   }

   /* disable separating again */
   SCIP_CALL( SCIPsetSeparating(origscip, SCIP_PARAMSETTING_OFF, TRUE) );

   SCIPdebugMessage("exiting sepaExeclpBasis\n");

   return SCIP_OKAY;
}

/** arbitrary primal solution separation method of separator */
#if 0
static
SCIP_DECL_SEPAEXECSOL(sepaExecsolBasis)
{  /*lint --e{715}*/
   SCIPerrorMessage("method of basis separator not implemented yet\n");
   SCIPABORT(); /*lint --e{527}*/

   return SCIP_OKAY;
}
#else
#define sepaExecsolBasis NULL
#endif

/*
 * separator specific interface methods
 */

/** creates the basis separator and includes it in SCIP */
SCIP_RETCODE SCIPincludeSepaBasis(
   SCIP*                 scip                /**< SCIP data structure */
   )
{
   SCIP_SEPADATA* sepadata;

   /* create master separator data */
   SCIP_CALL( SCIPallocMemory(scip, &sepadata) );

   sepadata->mastercuts = NULL;
   sepadata->origcuts = NULL;
   sepadata->norigcuts = 0;
   sepadata->nmastercuts = 0;
   sepadata->maxcuts = 0;
   sepadata->newcuts = NULL;
   sepadata->nnewcuts = 0;
   sepadata->maxnewcuts = 0;
   sepadata->objrow = NULL;

   /* include separator */
   SCIP_CALL( SCIPincludeSepa(scip, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
         SEPA_USESSUBSCIP, SEPA_DELAY,
         sepaCopyBasis, sepaFreeBasis, sepaInitBasis, sepaExitBasis, sepaInitsolBasis, sepaExitsolBasis, sepaExeclpBasis, sepaExecsolBasis,
         sepadata) );

   /* add basis separator parameters */
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/enable", "is basis separator enabled?",
         &(sepadata->enable), FALSE, TRUE, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/enableobj", "is objective constraint of separator enabled?",
         &(sepadata->enableobj), FALSE, FALSE, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/enableobjround", "round obj rhs/lhs of obj constraint if obj is int?",
         &(sepadata->enableobjround), FALSE, FALSE, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/enableppcuts", "add cuts generated during pricing to newconss array?",
         &(sepadata->enableppcuts), FALSE, FALSE, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/enableppobjconss", "is objective constraint for redcost of each pp of "
      "separator enabled?", &(sepadata->enableppobjconss), FALSE, FALSE, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/enableppobjcg", "is objective constraint for redcost of each pp during "
      "pricing of separator enabled?", &(sepadata->enableppobjcg), FALSE, FALSE, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/genobjconvex", "generated obj convex dynamically",
         &(sepadata->genobjconvex), FALSE, FALSE, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/enableposslack", "should positive slack influence the probing objective "
      "function?", &(sepadata->enableposslack), FALSE, FALSE, NULL, NULL) );
   SCIP_CALL( SCIPaddIntParam(GCGmasterGetOrigprob(scip), "sepa/basis/posslackexp", "exponent of positive slack usage",
         &(sepadata->posslackexp), FALSE, 1, 1, INT_MAX, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/posslackexpgen", "automatically generated exponent?",
            &(sepadata->posslackexpgen), FALSE, FALSE, NULL, NULL) );
   SCIP_CALL( SCIPaddRealParam(GCGmasterGetOrigprob(scip), "sepa/basis/posslackexpgenfactor", "factor for automatically generated exponent",
            &(sepadata->posslackexpgenfactor), FALSE, 0.1, SCIPepsilon(GCGmasterGetOrigprob(scip)),
            SCIPinfinity(GCGmasterGetOrigprob(scip)), NULL, NULL) );
   SCIP_CALL( SCIPaddRealParam(GCGmasterGetOrigprob(scip), "sepa/basis/objconvex", "convex combination factor",
         &(sepadata->objconvex), FALSE, 1.0, 0.0, 1.0, NULL, NULL) );
   SCIP_CALL( SCIPaddIntParam(GCGmasterGetOrigprob(scip), "sepa/basis/paramsetting", "parameter returns which parameter setting is used for "
      "separation (default = 0, aggressive = 1, fast = 2", &(sepadata->separationsetting), FALSE, 1, 0, 2, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/chgobj", "parameter returns if basis is searched with different objective",
      &(sepadata->chgobj), FALSE, TRUE, NULL, NULL) );
   SCIP_CALL( SCIPaddIntParam(GCGmasterGetOrigprob(scip), "sepa/basis/iterations", "parameter returns if number new rows adding"
      "iterations (rows just cut off probing lp sol)", &(sepadata->iterations), FALSE, 100, 1, 10000000 , NULL, NULL) );
   SCIP_CALL( SCIPaddIntParam(GCGmasterGetOrigprob(scip), "sepa/basis/mincuts", "parameter returns number of minimum cuts needed to "
      "return *result = SCIP_Separated", &(sepadata->mincuts), FALSE, 1, 1, 100, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/chgobjallways", "parameter returns if obj is changed not only in the "
      "first iteration", &(sepadata->chgobjallways), FALSE, FALSE, NULL, NULL) );
   SCIP_CALL( SCIPaddBoolParam(GCGmasterGetOrigprob(scip), "sepa/basis/forcecuts", "parameter returns if cuts are forced to enter the LP ",
      &(sepadata->forcecuts), FALSE, FALSE, NULL, NULL) );

   return SCIP_OKAY;
}

/** returns the array of original cuts saved in the separator data */
SCIP_ROW** GCGsepaBasisGetOrigcuts(
   SCIP*                 scip                /**< SCIP data structure */
   )
{
   SCIP_SEPA* sepa;
   SCIP_SEPADATA* sepadata;

   assert(scip != NULL);

   sepa = SCIPfindSepa(scip, SEPA_NAME);
   assert(sepa != NULL);

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);

   return sepadata->origcuts;
}

/** returns the number of original cuts saved in the separator data */
int GCGsepaBasisGetNOrigcuts(
   SCIP*                 scip                /**< SCIP data structure */
   )
{
   SCIP_SEPA* sepa;
   SCIP_SEPADATA* sepadata;

   assert(scip != NULL);

   sepa = SCIPfindSepa(scip, SEPA_NAME);
   assert(sepa != NULL);

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);

   return sepadata->norigcuts;
}

/** returns the array of master cuts saved in the separator data */
SCIP_ROW** GCGsepaBasisGetMastercuts(
   SCIP*                 scip                /**< SCIP data structure */
   )
{
   SCIP_SEPA* sepa;
   SCIP_SEPADATA* sepadata;

   assert(scip != NULL);

   sepa = SCIPfindSepa(scip, SEPA_NAME);
   assert(sepa != NULL);

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);

   return sepadata->mastercuts;
}

/** returns the number of master cuts saved in the separator data */
int GCGsepaBasisGetNMastercuts(
   SCIP*                 scip                /**< SCIP data structure */
   )
{
   SCIP_SEPA* sepa;
   SCIP_SEPADATA* sepadata;

   assert(scip != NULL);

   sepa = SCIPfindSepa(scip, SEPA_NAME);
   assert(sepa != NULL);

   sepadata = SCIPsepaGetData(sepa);
   assert(sepadata != NULL);

   return sepadata->nmastercuts;
}

/** transforms cut in pricing variables to cut in original variables and adds it to newcuts array */
SCIP_RETCODE GCGsepaBasisAddPricingCut(
   SCIP*                scip,
   int                  ppnumber,
   SCIP_ROW*            cut
   )
{
   SCIP* origscip;
   SCIP_SEPA* sepa;
   SCIP_SEPADATA* sepadata;

   SCIP* pricingprob;
   SCIP_Real* vals;
   SCIP_COL** cols;
   SCIP_VAR** pricingvars;
   int nvars;

   int i;
   int j;
   int k;

   char name[SCIP_MAXSTRLEN];

   assert(GCGisMaster(scip));

   sepa = SCIPfindSepa(scip, SEPA_NAME);

   if( sepa == NULL )
   {
      SCIPerrorMessage("sepa basis not found\n");
      return SCIP_OKAY;
   }

   sepadata = SCIPsepaGetData(sepa);
   origscip = GCGmasterGetOrigprob(scip);
   pricingprob = GCGgetPricingprob(origscip, ppnumber);

   if( !sepadata->enableppcuts )
   {
      return SCIP_OKAY;
   }

   assert(!SCIProwIsLocal(cut));

   nvars = SCIProwGetNNonz(cut);
   cols = SCIProwGetCols(cut);
   vals = SCIProwGetVals(cut);

   if( nvars == 0 )
   {
      return SCIP_OKAY;
   }

   SCIP_CALL( SCIPallocMemoryArray(scip, &pricingvars, nvars) );

   for( i = 0; i < nvars; ++i )
   {
      pricingvars[i] = SCIPcolGetVar(cols[i]);
      assert(pricingvars[i] != NULL);
   }

   for( k = 0; k < GCGgetNIdenticalBlocks(origscip, ppnumber); ++k )
   {
      SCIP_ROW* origcut;

      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "ppcut_%d_%d_%d", SCIPsepaGetNCalls(sepa), ppnumber, k);

      SCIP_CALL( SCIPcreateEmptyRowUnspec(origscip, &origcut, name,
         ( SCIPisInfinity(pricingprob, -SCIProwGetLhs(cut)) ?
            -SCIPinfinity(origscip) : SCIProwGetLhs(cut) - SCIProwGetConstant(cut)),
         ( SCIPisInfinity(pricingprob, SCIProwGetRhs(cut)) ?
            SCIPinfinity(origscip) : SCIProwGetRhs(cut) - SCIProwGetConstant(cut)),
             FALSE, FALSE, TRUE) );

      for( j = 0; j < nvars ; ++j )
      {
         SCIP_VAR* var;

         if( !GCGvarIsPricing(pricingvars[j]) )
         {
            nvars = 0;
            break;
         }
         assert(GCGvarIsPricing(pricingvars[j]));

         var = GCGpricingVarGetOrigvars(pricingvars[j])[k];
         assert(var != NULL);

         SCIP_CALL( SCIPaddVarToRow(origscip, origcut, var, vals[j]) );
      }

      if( nvars > 0 )
      {
         SCIP_CALL( ensureSizeNewCuts(scip, sepadata, sepadata->nnewcuts + 1) );

         sepadata->newcuts[sepadata->nnewcuts] = origcut;
         SCIP_CALL( SCIPcaptureRow(scip, sepadata->newcuts[sepadata->nnewcuts]) );
         ++(sepadata->nnewcuts);

         SCIPdebugMessage("cut added to orig cut pool\n");
      }
      SCIP_CALL( SCIPreleaseRow(origscip, &origcut) );
   }

   SCIPfreeMemoryArray(scip, &pricingvars);

   return SCIP_OKAY;
}

/** add cuts which are due to the latest objective function of the pricing problems
 *  (reduced cost non-negative) */
SCIP_RETCODE SCIPsepaBasisAddPPObjConss(
   SCIP*                scip,               /**< SCIP data structure */
   int                  ppnumber,           /**< number of pricing problem */
   SCIP_Real            dualsolconv,        /**< dual solution corresponding to convexity constraint */
   SCIP_Bool            newcuts             /**< add cut to newcuts in sepadata? (otherwise add it just to the cutpool) */
)
{
   SCIP_SEPA* sepa;

   assert(GCGisMaster(scip));

   sepa = SCIPfindSepa(scip, SEPA_NAME);

   if( sepa == NULL )
   {
      SCIPerrorMessage("sepa basis not found\n");
      return SCIP_OKAY;
   }

   SCIP_CALL( addPPObjConss(GCGmasterGetOrigprob(scip), sepa, ppnumber, dualsolconv, newcuts, FALSE) );

   return SCIP_OKAY;
}


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"
#include "error_cblas_l2.h"

void
cblas_zgbmv (const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA,
             const int M, const int N, const int KL, const int KU,
             const void *alpha, const void *A, const int lda, const void *X,
             const int incX, const void *beta, void *Y, const int incY)
{
#define BASE double
#include "source_gbmv_c.h"
#undef BASE
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE