text: #include <math.h>
  #include <stdio.h>
  #include <string>
  #include <random>
  #include <time.h>
  #include <unordered_map>
  #include "mex.h"
  #include "matrix.h"
  
  #define maxR 4096
  #define maxN 524288
  #define maxe 8192
  #define maxK 512
  #define accr 1e-12
  
  void swap_double(double * a, int p, int q){
  	double temp = a[p];
  	a[p] = a[q];
  	a[q] = temp;
  }
  void swap_int(int * a, int p, int q){
  	int temp = a[p];
  	a[p] = a[q];
  	a[q] = temp;
  }
  
  void MYsort(double * a, double * b, int len, char type){
  	// sort a according to rule type, permutate b according to a 
  	if(len <= 1) return;
   	std::default_random_engine generator;
  	std::uniform_int_distribution<int> distribution(0, len-1);
  	int pos = distribution(generator), p = 0;
  	double pivot = a[pos];
      swap_double(a, pos, len-1);
      swap_double(b, pos, len-1);
  	if(type == 'd'){
  		for(int i = 0; i < len-1; i++){
  			if(a[i] > pivot){
  				swap_double(a, i, p);
  				swap_double(b, i, p);
  			    p++;
  			}
  		}
  		swap_double(a, p, len-1);
  		swap_double(b, p, len-1);
  	}
  	if(type == 'a'){
  		for(int i = 0; i < len-1; i++){
  			if(a[i] < pivot){
  				swap_double(a, i, p);
  				swap_double(b, i, p);
  			    p++;
  			}
  		}
  		swap_double(a, p, len-1);
  		swap_double(b, p, len-1);
  	}
  	MYsort(a, b, p, type);
  	MYsort(a+p+1, b+p+1, len-p-1, type);
  }
  
  void projection_edge_weight(double * y, double * a, double * W, const double * para){
      y[0] = (W[0]*a[0] - W[1]*a[1])/(W[0] + W[1]);
      if(y[0] < 0) y[0] = fmax(-para[0], y[0]);
      else y[0] = fmin(para[0], y[0]);
      y[1] = - y[0];
  }
  
  void projection_para_edge_weight(double * y, double * a, double * W, const double * para, int len){
      int n = len/2, i, twoi;
      for(i = 0; i < n;i++){
          twoi = 2*i;
          y[twoi] = (W[twoi]*a[twoi] - W[twoi+1]*a[twoi+1])/(W[twoi] + W[twoi+1]);
          if(y[twoi] < 0) y[twoi] = fmax(-para[i], y[twoi]);
          else y[twoi] = fmin(para[i], y[twoi]);
          y[twoi+1] = - y[twoi];
      }
  }
  
  
  void projection_homo_weight(double * y, double * a, double * W, const double * para, int len){
  	double * invW = new double[len];
      double suma = 0, suminvW = 0, lambda = 0, ypossum = 0, invWpossum = 0,
              ynegsum = 0, invWnegsum = 0, dlambda_thre = 0, tempy = 0;
      int *maskpos = new int [len];
      int *maskneg = new int [len];
      int maskposlen = 0, maskneglen = 0;
      int i;
      for(i = 0; i < len; i++){
          suma += a[i];
          invW[i] = 1/W[i];
          suminvW += invW[i];
          //mexPrintf("%f, %f", a[i], W[i]);
      }
      lambda = suma/suminvW;
      for(i = 0; i < len; i++){
          y[i] = -lambda * invW[i] + a[i];
          if(y[i] > accr){
              maskpos[maskposlen++] = i;
              ypossum += y[i];
              invWpossum += invW[i];
          }
          if(y[i] < -accr){
              maskneg[maskneglen++] = i;
              ynegsum += y[i];
              invWnegsum += invW[i];            
          }
      }
      while(ypossum > para[0] + accr && maskposlen >= 1){
          dlambda_thre = (ypossum - para[0])/invWpossum;
          for(i = maskposlen-1; i >=0; i--){
              //dlambda = min(y[maskpos[i]]*W[maskpos[i]],dlambda_thre);
              tempy = fmax(0, y[maskpos[i]]-dlambda_thre * invW[maskpos[i]]);
              ypossum -= y[maskpos[i]] - tempy;
              y[maskpos[i]] = tempy;
              if(y[maskpos[i]] < accr){
                  invWpossum -= invW[maskpos[i]];
                  swap_int(maskpos, i, maskposlen-1);
                  maskposlen--;
              }
          }
      }
      while(ynegsum < -para[0] - accr && maskneglen >= 1){
          dlambda_thre = (-ynegsum - para[0])/invWnegsum;
          for(i = maskneglen-1; i >=0; i--){
              //dlambda = min(y[maskpos[i]]*W[maskpos[i]],dlambda_thre);
              tempy = fmin(0, y[maskneg[i]] + dlambda_thre * invW[maskneg[i]]);
              ynegsum += tempy - y[maskneg[i]];
              y[maskneg[i]] = tempy;
              if(y[maskneg[i]] > -accr){
                  invWnegsum -= invW[maskneg[i]];
                  swap_int(maskneg, i, maskneglen-1);
                  maskneglen--;
              }
          }
      }
      delete[] invW;
      delete[] maskpos;
      delete[] maskneg;
  }
  
  void concave_card_subroutine_weight(double * x, double lambdamin, double lambdamax, double * a, double * W, const double * para, int i, int len){
  	int j;
  	double minval = INFINITY, tempsum = 0, lambda = 0, sumW =0;
      double * index = new double [len] , * val = new double [len], * invindex,
              * newa, * newW;
  	int minindex;
  //     printf("%d", len);
  //     for(j = 0; j < len; j++){
  //         printf("a[%d]:%f\n", j, a[j]);
  //         printf("W[%d]:%f\n", j, W[j]);
  //         printf("para[%d]:%f\n", j, para[j]);
  //     }
      for(j = 0; j < len; j++){
          index[j] = j;
      }
  	if(i%2 == 0){
  		for(j = 0; j < len; j++){
              lambda -= para[j] + a[j];
              sumW += W[j];
          }
          lambda = lambda/sumW;
  		for(j = 0; j < len; j++){
  			val[j] = lambda*W[j] + a[j];
  		}
          MYsort(val, index, len, 'a');
          for(j = 0; j < len; j++){
          	tempsum = tempsum + val[j] + para[j];
  			if(tempsum < minval){
  				minval = tempsum;
  				minindex = j;
  			}
          }
  		if(minval > 0 || minindex == len - 1){
  			for(j = 0; j < len; j++){
  				x[j] = lambda;
  			}
              delete[] index;
              delete[] val;
  			return;
  		}
  	}else{
  		lambda = (lambdamin + lambdamax)/2;
  		for(j = 0; j < len; j++){
  			val[j] = lambda*W[j] + a[j];
  		}
          MYsort(val, index, len, 'a');
          for(j = 0; j < len; j++){
          	tempsum = tempsum + val[j] + para[j];
  			if(tempsum < minval){
  				minval = tempsum;
  				minindex = j;
  			}
          }       
  		if(minval > 0){
  			concave_card_subroutine_weight(x, lambdamin, lambda, a, W, para, i+1, len);
              delete[] index;
              delete[] val;
  			return;
  		}
  		if(minindex == len-1){
  			concave_card_subroutine_weight(x, lambda, lambdamax, a, W, para, i+1, len);
              delete[] index;
              delete[] val;
  			return;
  		}		
  	}
      invindex = new double [len];
      newa = new double [len];
      newW = new double [len];
      for(j = 0; j < len; j++){
          invindex[(int)index[j]] = j;
          newa[j] = a[(int)index[j]];
          newW[j] = W[(int)index[j]];
      }
  	concave_card_subroutine_weight(x, lambda, lambdamax, newa, newW, para, i+1, minindex+1);
  	concave_card_subroutine_weight(x + minindex+1, lambdamin, lambda, newa + minindex + 1, newW + minindex + 1, para + minindex + 1, i+1, len - minindex - 1);
      double * newx = new double [len];
      for(j = 0; j < len; j++){
          newx[j] = x[(int)invindex[j]];
      }
      for(j = 0; j < len; j++){
          x[j] = newx[j];
      }    
      delete[] newx;
      delete[] index;
      delete[] invindex;
      delete[] newW;
      delete[] newa;
      delete[] val;
  
  }
  
  
  void projection_concave_card_weight(double * y, double * a, double * W, const double * para, int len){
  
      double * invW = new double [len];
      double * nega = new double [len];
      double lambdamax  = -INFINITY, lambdamin = INFINITY;
      int i;
  	for(i = 0; i < len; i++){
  		invW[i] = 1/W[i];
          nega[i] = -a[i];
          lambdamax = fmax(lambdamax, (-para[len-1]+a[i])*W[i]);
          lambdamin = fmin(lambdamin, (-para[0]+a[i])*W[i]);
  	}    
      
  	double * x = new double[len];
      // nega , invW may change
  	concave_card_subroutine_weight(x, lambdamin, lambdamax, nega, invW, para, 0, len);
  	for(i = 0; i < len; i++){
  		y[i] = a[i] - x[i]/W[i];
  	}
  	delete[] invW;
  	delete[] nega;
  	delete[] x;
  }
  
  double gap_eff(int ** incidence_list, int * incidence_list_size,
          double ** parameter_list, int * parameter_list_size,
          char ** submodular_type, int * submodular_type_size,
          double * bias_vec, double * y_sum, double * x,
              int N, int R){
      double gap = 0, dv = 0, tempmax, tempmin;
      int i,j;
      for(i = 0; i < N; i++){
          gap -= y_sum[i]*x[i];
      }
      for(i = 0; i < R; i++){
          if(submodular_type[i][0] == 'h'){
              tempmax = -INFINITY;
              tempmin = INFINITY;
              for(j = 0; j < incidence_list_size[i]; j++){
                  tempmax = fmax(tempmax, x[incidence_list[i][j]]);
                  tempmin = fmin(tempmin, x[incidence_list[i][j]]);
              }
              gap += (tempmax - tempmin)*parameter_list[i][0];
          }
          if(submodular_type[i][0] == 'e'){
              gap += fabs(x[incidence_list[i][0]] - x[incidence_list[i][1]])*parameter_list[i][0];
          }
          if(submodular_type[i][0] == 'p'){
              for(j = 0; j < incidence_list_size[i]/2; j++){
                  gap += fabs(x[incidence_list[i][2*j]] - x[incidence_list[i][2*j+1]])*parameter_list[i][j];
              }
          }
          if(submodular_type[i][0] == 'c'){
              double * tempx = new double [incidence_list_size[i]];
              double * index = new double [incidence_list_size[i]]();
              for(j = 0; j < incidence_list_size[i]; j++){
                  tempx[j] = x[incidence_list[i][j]];
              }
              MYsort(tempx, index, incidence_list_size[i], 'd');
              for(j = 0; j < incidence_list_size[i]; j++){
                  gap += tempx[j]*parameter_list[i][j];
              }
              delete[] tempx;
              delete[] index;
          }
      }
      
      return gap;
  }
  
  double gap_d_eff(int ** incidence_list, int * incidence_list_size,
          double ** parameter_list, int * parameter_list_size,
          char ** submodular_type, int * submodular_type_size,
          double * bias_vec, double * x,
              int N, int R){
      double * xsort = new double [N], * index = new double [N], * invindex = new double [N];
      double * level_cost = new double [N](), level_cost_min = 0;
      double pos_x_sum = 0, indexmax, indexmin;
      int i,j;
      for(i = 0; i < N; i++){
          xsort[i] = x[i];
          index[i] = i;
          if(x[i] > 0) pos_x_sum += x[i];
      }
      MYsort(xsort, index, N, 'd');
      for(i = 0; i < N; i++){
          invindex[(int)index[i]] = i;
      }
      for(i = 0; i < R; i++){
          if(submodular_type[i][0] == 'h'){
              indexmax = 0;
              indexmin = N-1;
              for(j = 0; j < incidence_list_size[i]; j++){
                  indexmax = fmax(indexmax, invindex[incidence_list[i][j]]);
                  indexmin = fmin(indexmin, invindex[incidence_list[i][j]]);
              }
              level_cost[(int)indexmin] += parameter_list[i][0];
              level_cost[(int)indexmax] -= parameter_list[i][0];
          }
          if(submodular_type[i][0] == 'e'){
              level_cost[(int)fmin(invindex[incidence_list[i][0]], invindex[incidence_list[i][1]])] += parameter_list[i][0];
              level_cost[(int)fmax(invindex[incidence_list[i][0]], invindex[incidence_list[i][1]])] -= parameter_list[i][0];
          }
          if(submodular_type[i][0] == 'p'){
              for(j = 0; j < incidence_list_size[i]/2; j++){
                  level_cost[(int)fmin(invindex[incidence_list[i][2*j]], invindex[incidence_list[i][2*j+1]])] += parameter_list[i][j];
                  level_cost[(int)fmax(invindex[incidence_list[i][2*j]], invindex[incidence_list[i][2*j+1]])] -= parameter_list[i][j];
              }
          }
          if(submodular_type[i][0] == 'c'){
              double * tempindex = new double [incidence_list_size[i]];
              double * nouseindex = new double [incidence_list_size[i]];
              for(j = 0; j < incidence_list_size[i]; j++){
                  tempindex[j] = invindex[incidence_list[i][j]];
                  nouseindex[j] = j;
              }
              MYsort(tempindex, nouseindex, incidence_list_size[i], 'a');
              for(j = 0; j < incidence_list_size[i]; j++){
                  level_cost[(int)tempindex[j]] += parameter_list[i][j];
              }
              delete[] tempindex;
              delete[] nouseindex;
          }
      }
      for(i = 0; i < N; i++){
          level_cost[i] += bias_vec[(int)index[i]];
      }
      level_cost_min = fmin(level_cost[0], level_cost_min);
      for(i = 1; i < N; i++){
          level_cost[i] += level_cost[i-1];
          if(level_cost[i]< level_cost_min) level_cost_min = level_cost[i];
                  
      }
      delete[] index;
      delete[] invindex;
      delete[] xsort;
      delete[] level_cost;
      //mexPrintf("%f\n", pos_x_sum);
      return level_cost_min + pos_x_sum;
  }
  void degree_comp_full(double * degree_vec, int ** incidence_list, int * incidence_list_size, int R){
      int i, j;
      for(i = 0; i < R; i++){
          for(j = 0; j < incidence_list_size[i]; j++){
              degree_vec[incidence_list[i][j]]++;
          }
      }    
  }
  
  void main_func(double * y_sum, double * record, double * record_discrete,
          int ** incidence_list, int * incidence_list_size,
          double ** parameter_list, int * parameter_list_size,
          char ** submodular_type, int * submodular_type_size,
          double * bias_vec, int N, int R, unsigned int T,
          unsigned int record_dis){
      double ** y = new double * [R];
      double * y_ave = new double [N], * temp_y_ave = new double [N];
      double a[maxe];
      double W[maxe];
      double newy[maxe];
      double * degree_vec = new double [N]();
      double * x = new double [N]();
      int i, j, t;
      double gap;
   	std::default_random_engine generator(time(0));
  	std::uniform_int_distribution<int> distribution(0, R-1);
      
      degree_comp_full(degree_vec, incidence_list, incidence_list_size, R);
      for(i = 0; i < R; i++){
          y[i] = new double [incidence_list_size[i]]();
      }
      for(i = 0; i < N; i++){
          y_ave[i] = bias_vec[i]/degree_vec[i];
      }
  
      for(t = 0; t < T; t++){
          for(i = 0; i < N; i++){
              temp_y_ave[i] = y_ave[i];
          }
          for(j = 0; j < R; j++){
              for(i = 0; i < incidence_list_size[j]; i++){
                  a[i] = y[j][i] -temp_y_ave[incidence_list[j][i]];
                  W[i] = degree_vec[incidence_list[j][i]];
              }
              if(submodular_type[j][0] == 'h'){
                  projection_homo_weight(newy, a, W, parameter_list[j], incidence_list_size[j]);
              }
              if(submodular_type[j][0] == 'e'){
                  projection_edge_weight(newy, a, W, parameter_list[j]);
              }
              if(submodular_type[j][0] == 'p'){
                  projection_para_edge_weight(newy, a, W, parameter_list[j], incidence_list_size[j]);
              }
              if(submodular_type[j][0] == 'c'){
                  projection_concave_card_weight(newy, a, W, parameter_list[j], incidence_list_size[j]);
              }
              for(i = 0; i < incidence_list_size[j]; i++){
                  y_ave[incidence_list[j][i]] += (newy[i] - y[j][i])/W[i];
                  y[j][i] = newy[i];
              }
          }
          if((t+1)%record_dis == 0){
              for(i = 0; i < N; i++){
                  x[i] = -degree_vec[i]* y_ave[i];
                  y_sum[i] = - x[i] - bias_vec[i];
              }
              gap = gap_eff(incidence_list, incidence_list_size,
                      parameter_list, parameter_list_size,
                      submodular_type, submodular_type_size,
                      bias_vec, y_sum, x,
                      N, R);
              mexPrintf("%d, gap: %f\n", t, gap);
              record[(t+1)/record_dis] = gap;
              gap = gap_d_eff(incidence_list, incidence_list_size,
                  parameter_list, parameter_list_size,
                  submodular_type, submodular_type_size,
                  bias_vec, x, N, R);
              record_discrete[(t+1)/record_dis] = gap;
          }
      }
      delete[] y;
      delete[] y_ave;
      delete[] temp_y_ave;
      delete[] degree_vec;
      delete[] x;
  }
  
  void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
  	if (nlhs != 3 || nrhs != 8) {
      mexWarnMsgTxt("Check Parameters");
      return;
      }
      // read incidence_list
      int R = *(mxGetPr(prhs[5]));
      const mxArray * incidence_list_org = prhs[0];
      mxArray * incidence_Element;
      //const int * R_org = mxGetDimensions(incidence_list_org);
      int * * incidence_list = new int * [R];
      int * incidence_list_size = new int [R];
      double * templist;
      int j, k;
      for(j = 0; j < R;j++){
         incidence_Element = mxGetCell(incidence_list_org, j);
         incidence_list_size[j] = (int)mxGetN(incidence_Element);
         incidence_list[j] = new int [incidence_list_size[j]];
         templist = mxGetPr(incidence_Element);
         for(k = 0; k < incidence_list_size[j]; k++){
             incidence_list[j][k] = templist[k] - 1;
         }
      }
      // read parameter_list
      const mxArray * parameter_list_org = prhs[1];
      mxArray * parameter_Element;
      double * * parameter_list = new double * [R];
      int * parameter_list_size = new int [R];
      for(j = 0; j < R;j++){
         parameter_Element = mxGetCell(parameter_list_org, j);
         parameter_list_size[j] = (int)mxGetN(parameter_Element);
         parameter_list[j] = new double[parameter_list_size[j]];
         templist = mxGetPr(parameter_Element);
         for(k = 0; k < parameter_list_size[j]; k++){
             parameter_list[j][k] = templist[k];
             //mexPrintf("%d, %d, %f\n", j, k, parameter_list[j][k]);
         }
      }
      // read submodular_type
      const mxArray * submodular_type_org = prhs[2];
      mxArray * submodular_type_Element;
      char ** submodular_type = new char * [R];
      int * submodular_type_size = new int [R];
      char * temp_type;
      for(j = 0; j < R;j++){
         submodular_type_Element = mxGetCell(submodular_type_org, j);
         temp_type = (char *)mxGetPr(submodular_type_Element);
         submodular_type_size[j] = (int)mxGetN(submodular_type_Element);
         submodular_type[j] = new char[submodular_type_size[j]];
         for(k=0; k < submodular_type_size[j]; k++){
             submodular_type[j][k] = (char)temp_type[2*k];
             //mexPrintf("%d, %d, %c\n", j, k, temp_type[2*k]);
         }
      }
      
      double * bias_vec = mxGetPr(prhs[3]);
      int N = *(mxGetPr(prhs[4]));
      unsigned int T = *(mxGetPr(prhs[6]));
      unsigned int record_dis = *(mxGetPr(prhs[7]));
      int record_num = T/record_dis + 1;
      
      plhs[0] = mxCreateDoubleMatrix(1, N, mxREAL);
      double * y_sum = mxGetPr(plhs[0]);
      plhs[1] = mxCreateDoubleMatrix(record_num, 1, mxREAL);
      double * record = mxGetPr(plhs[1]);
      plhs[2] = mxCreateDoubleMatrix(record_num, 1, mxREAL);
      double * record_discrete = mxGetPr(plhs[2]);    
      main_func(y_sum, record, record_discrete,
              incidence_list, incidence_list_size,
              parameter_list, parameter_list_size,
              submodular_type, submodular_type_size,
              bias_vec, N, R, T, record_dis);
      delete[] incidence_list;
      delete[] incidence_list_size;
      delete[] parameter_list;
      delete[] parameter_list_size;
      delete[] submodular_type;
      delete[] submodular_type_size;
  }
  
  
  
  
  
  
  
  