/*
 *  call to numeric_deriv from R -
 *  .Call("numeric_deriv", expr, theta, rho)
 *  Returns: ans
 */
SEXP
numeric_deriv(SEXP expr, SEXP theta, SEXP rho, SEXP dir)
{
    SEXP ans, gradient, pars;
    double eps = sqrt(DOUBLE_EPS), *rDir;
    int start, i, j, k, lengthTheta = 0;

    if(!isString(theta))
	error(_("'theta' should be of type character"));
    if (isNull(rho)) {
	error(_("use of NULL environment is defunct"));
	rho = R_BaseEnv;
    } else
	if(!isEnvironment(rho))
	    error(_("'rho' should be an environment"));
    PROTECT(dir = coerceVector(dir, REALSXP));
    if(TYPEOF(dir) != REALSXP || LENGTH(dir) != LENGTH(theta))
	error(_("'dir' is not a numeric vector of the correct length"));
    rDir = REAL(dir);

    PROTECT(pars = allocVector(VECSXP, LENGTH(theta)));

    PROTECT(ans = duplicate(eval(expr, rho)));

    if(!isReal(ans)) {
	SEXP temp = coerceVector(ans, REALSXP);
	UNPROTECT(1);
	PROTECT(ans = temp);
    }
    for(i = 0; i < LENGTH(ans); i++) {
	if (!R_FINITE(REAL(ans)[i]))
	    error(_("Missing value or an infinity produced when evaluating the model"));
    }
    const void *vmax = vmaxget();
    for(i = 0; i < LENGTH(theta); i++) {
	const char *name = translateChar(STRING_ELT(theta, i));
	SEXP s_name = install(name);
	SEXP temp = findVar(s_name, rho);
	if(isInteger(temp))
	    error(_("variable '%s' is integer, not numeric"), name);
	if(!isReal(temp))
	    error(_("variable '%s' is not numeric"), name);
	if (MAYBE_SHARED(temp)) /* We'll be modifying the variable, so need to make sure it's unique PR#15849 */
	    defineVar(s_name, temp = duplicate(temp), rho);
	MARK_NOT_MUTABLE(temp);
	SET_VECTOR_ELT(pars, i, temp);
	lengthTheta += LENGTH(VECTOR_ELT(pars, i));
    }
    vmaxset(vmax);
    PROTECT(gradient = allocMatrix(REALSXP, LENGTH(ans), lengthTheta));

    for(i = 0, start = 0; i < LENGTH(theta); i++) {
	for(j = 0; j < LENGTH(VECTOR_ELT(pars, i)); j++, start += LENGTH(ans)) {
	    SEXP ans_del;
	    double origPar, xx, delta;

	    origPar = REAL(VECTOR_ELT(pars, i))[j];
	    xx = fabs(origPar);
	    delta = (xx == 0) ? eps : xx*eps;
	    REAL(VECTOR_ELT(pars, i))[j] += rDir[i] * delta;
	    PROTECT(ans_del = eval(expr, rho));
	    if(!isReal(ans_del)) ans_del = coerceVector(ans_del, REALSXP);
	    UNPROTECT(1);
	    for(k = 0; k < LENGTH(ans); k++) {
		if (!R_FINITE(REAL(ans_del)[k]))
		    error(_("Missing value or an infinity produced when evaluating the model"));
		REAL(gradient)[start + k] =
		    rDir[i] * (REAL(ans_del)[k] - REAL(ans)[k])/delta;
	    }
	    REAL(VECTOR_ELT(pars, i))[j] = origPar;
	}
    }
    setAttrib(ans, install("gradient"), gradient);
    UNPROTECT(4);
    return ans;
}

SEXP lapack_qr(SEXP Xin, SEXP tl)
{
    SEXP ans, Givens, Gcpy, nms, pivot, qraux, X;
    int i, n, nGivens = 0, p, trsz, *Xdims, rank;
    double rcond = 0., tol = asReal(tl), *work;

    if (!(isReal(Xin) & isMatrix(Xin)))
	error(_("X must be a real (numeric) matrix"));
    if (tol < 0.) error(_("tol, given as %g, must be non-negative"), tol);
    if (tol > 1.) error(_("tol, given as %g, must be <= 1"), tol);
    ans = PROTECT(allocVector(VECSXP,5));
    SET_VECTOR_ELT(ans, 0, X = duplicate(Xin));
    Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
    n = Xdims[0]; p = Xdims[1];
    SET_VECTOR_ELT(ans, 2, qraux = allocVector(REALSXP, (n < p) ? n : p));
    SET_VECTOR_ELT(ans, 3, pivot = allocVector(INTSXP, p));
    for (i = 0; i < p; i++) INTEGER(pivot)[i] = i + 1;
    trsz = (n < p) ? n : p;	/* size of triangular part of decomposition */
    rank = trsz;
    Givens = PROTECT(allocVector(VECSXP, rank - 1));
    setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 5));
    SET_STRING_ELT(nms, 0, mkChar("qr"));
    SET_STRING_ELT(nms, 1, mkChar("rank"));
    SET_STRING_ELT(nms, 2, mkChar("qraux"));
    SET_STRING_ELT(nms, 3, mkChar("pivot"));
    SET_STRING_ELT(nms, 4, mkChar("Givens"));
    if (n > 0 && p > 0) {
	int  info, *iwork, lwork;
	double *xpt = REAL(X), tmp;

	lwork = -1;
	F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), &tmp, &lwork, &info);
	if (info)
	    error(_("First call to dgeqrf returned error code %d"), info);
	lwork = (int) tmp;
	work = (double *) R_alloc((lwork < 3*trsz) ? 3*trsz : lwork,
				  sizeof(double));
	F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), work, &lwork, &info);
	if (info)
	    error(_("Second call to dgeqrf returned error code %d"), info);
	iwork = (int *) R_alloc(trsz, sizeof(int));
	F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
			 work, iwork, &info);
	if (info)
	    error(_("Lapack routine dtrcon returned error code %d"), info);
	while (rcond < tol) {	/* check diagonal elements */
	    double minabs = (xpt[0] < 0.) ? -xpt[0]: xpt[0];
	    int jmin = 0;
	    for (i = 1; i < rank; i++) {
		double el = xpt[i*(n+1)];
		el = (el < 0.) ? -el: el;
		if (el < minabs) {
		    jmin = i;
		    minabs = el;
		}
	    }
	    if (jmin < (rank - 1)) {
		SET_VECTOR_ELT(Givens, nGivens, getGivens(xpt, n, jmin, rank));
		nGivens++;
	    }
	    rank--;
	    F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
			     work, iwork, &info);
	    if (info)
		error(_("Lapack routine dtrcon returned error code %d"), info);
	}
    }
    SET_VECTOR_ELT(ans, 4, Gcpy = allocVector(VECSXP, nGivens));
    for (i = 0; i < nGivens; i++)
	SET_VECTOR_ELT(Gcpy, i, VECTOR_ELT(Givens, i));
    SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
    setAttrib(ans, install("useLAPACK"), ScalarLogical(1));
    setAttrib(ans, install("rcond"), ScalarReal(rcond));
    UNPROTECT(2);
    return ans;
}

SEXP DropDims(SEXP x)
{
    SEXP dims, dimnames, newnames = R_NilValue;
    int i, n, ndims;

    PROTECT(x);
    dims = getDimAttrib(x);
    dimnames = getDimNamesAttrib(x);

    /* Check that dropping will actually do something. */
    /* (1) Check that there is a "dim" attribute. */

    if (dims == R_NilValue) {
	UNPROTECT(1);
	return x;
    }
    ndims = LENGTH(dims);

    /* (2) Check whether there are redundant extents */
    n = 0;
    for (i = 0; i < ndims; i++)
	if (INTEGER(dims)[i] != 1) n++;
    if (n == ndims) {
	UNPROTECT(1);
	return x;
    }

    if (n <= 1) {
	/* We have reduced to a vector result.
	   If that has length one, it is ambiguous which dimnames to use,
	   so use it if there is only one (as from R 2.7.0).
	 */
	if (dimnames != R_NilValue) {
	    if(XLENGTH(x) != 1) {
		for (i = 0; i < LENGTH(dims); i++) {
		    if (INTEGER(dims)[i] != 1) {
			newnames = VECTOR_ELT(dimnames, i);
			break;
		    }
		}
	    } else { /* drop all dims: keep names if unambiguous */
		int cnt;
		for(i = 0, cnt = 0; i < LENGTH(dims); i++)
		    if(VECTOR_ELT(dimnames, i) != R_NilValue) cnt++;
		if(cnt == 1)
		    for (i = 0; i < LENGTH(dims); i++) {
			newnames = VECTOR_ELT(dimnames, i);
			if(newnames != R_NilValue) break;
		    }
	    }
	}
	PROTECT(newnames);
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
	setAttrib(x, R_DimSymbol, R_NilValue);
	setAttrib(x, R_NamesSymbol, newnames);
	/* FIXME: the following is desirable, but pointless as long as
	   subset.c & others have a contrary version that leaves the
	   S4 class in, incorrectly, in the case of vectors.  JMC
	   3/3/09 */
/* 	if(IS_S4_OBJECT(x)) {/\* no longer valid subclass of array or
 	matrix *\/ */
/* 	    setAttrib(x, R_ClassSymbol, R_NilValue); */
/* 	    UNSET_S4_OBJECT(x); */
/* 	} */
	UNPROTECT(1);
    } else {
	/* We have a lower dimensional array. */
	SEXP newdims, dnn, newnamesnames = R_NilValue;
	dnn = getNamesAttrib(dimnames);
	PROTECT(newdims = allocVector(INTSXP, n));
	for (i = 0, n = 0; i < ndims; i++)
	    if (INTEGER(dims)[i] != 1)
		INTEGER(newdims)[n++] = INTEGER(dims)[i];
	if (!isNull(dimnames)) {
	    int havenames = 0;
	    for (i = 0; i < ndims; i++)
		if (INTEGER(dims)[i] != 1 &&
		    VECTOR_ELT(dimnames, i) != R_NilValue)
		    havenames = 1;
	    if (havenames) {
		PROTECT(newnames = allocVector(VECSXP, n));
		PROTECT(newnamesnames = allocVector(STRSXP, n));
		for (i = 0, n = 0; i < ndims; i++) {
		    if (INTEGER(dims)[i] != 1) {
			if(!isNull(dnn))
			    SET_STRING_ELT(newnamesnames, n,
					   STRING_ELT(dnn, i));
			SET_VECTOR_ELT(newnames, n++, VECTOR_ELT(dimnames, i));
		    }
		}
	    }
	    else dimnames = R_NilValue;
	}
	PROTECT(dimnames);
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
	setAttrib(x, R_DimSymbol, newdims);
	if (dimnames != R_NilValue)
	{
	    if(!isNull(dnn))
		setAttrib(newnames, R_NamesSymbol, newnamesnames);
//.........这里部分代码省略.........

SEXP pollSocket(SEXP sockets_, SEXP events_, SEXP timeout_) {
    SEXP result;
    
    if(TYPEOF(timeout_) != INTSXP) {
        error("poll timeout must be an integer.");
    }

    if(TYPEOF(sockets_) != VECSXP || LENGTH(sockets_) == 0) {
        error("A non-empy list of sockets is required as first argument.");
    }

    int nsock = LENGTH(sockets_);
    PROTECT(result = allocVector(VECSXP, nsock));

    if (TYPEOF(events_) != VECSXP) {
        error("event list must be a list of strings or a list of vectors of strings.");
    }
    if(LENGTH(events_) != nsock) {
        error("event list must be the same length as socket list.");
    }

    zmq_pollitem_t *pitems = (zmq_pollitem_t*)R_alloc(nsock, sizeof(zmq_pollitem_t));
    if (pitems == NULL) {
        error("failed to allocate memory for zmq_pollitem_t array.");
    }

    try {
        for (int i = 0; i < nsock; i++) {
            zmq::socket_t* socket = reinterpret_cast<zmq::socket_t*>(checkExternalPointer(VECTOR_ELT(sockets_, i), "zmq::socket_t*"));
            pitems[i].socket = (void*)*socket;
            pitems[i].events = rzmq_build_event_bitmask(VECTOR_ELT(events_, i));
        }

        int rc = zmq::poll(pitems, nsock, *INTEGER(timeout_));

        if(rc >= 0) {
            for (int i = 0; i < nsock; i++) {
                SEXP events, names;

                // Pre count number of polled events so we can
                // allocate appropriately sized lists.
                short eventcount = 0;
                if (pitems[i].events & ZMQ_POLLIN) eventcount++;
                if (pitems[i].events & ZMQ_POLLOUT) eventcount++;
                if (pitems[i].events & ZMQ_POLLERR) eventcount++;

                PROTECT(events = allocVector(VECSXP, eventcount));
                PROTECT(names = allocVector(VECSXP, eventcount));

                eventcount = 0;
                if (pitems[i].events & ZMQ_POLLIN) {
                    SET_VECTOR_ELT(events, eventcount, ScalarLogical(pitems[i].revents & ZMQ_POLLIN));
                    SET_VECTOR_ELT(names, eventcount, mkChar("read"));
                    eventcount++;
                }

                if (pitems[i].events & ZMQ_POLLOUT) {
                    SET_VECTOR_ELT(names, eventcount, mkChar("write"));

                    SET_VECTOR_ELT(events, eventcount, ScalarLogical(pitems[i].revents & ZMQ_POLLOUT));
                    eventcount++;
                }

                if (pitems[i].events & ZMQ_POLLERR) {
                    SET_VECTOR_ELT(names, eventcount, mkChar("error"));
                    SET_VECTOR_ELT(events, eventcount, ScalarLogical(pitems[i].revents & ZMQ_POLLERR));
                }
                setAttrib(events, R_NamesSymbol, names);
                SET_VECTOR_ELT(result, i, events);
            }
        } else {
            error("polling zmq sockets failed.");
        }
    } catch(std::exception& e) {
        error(e.what());
    }
    // Release the result list (1), and per socket
    // events lists with associated names (2*nsock).
    UNPROTECT(1 + 2*nsock);
    return result;
}

//.........这里部分代码省略.........
              }
            }
            gap=0;
          }
        }
        if((int)gap > (int)maxgap) {  /* check that we don't have excessive trailing gap */
          for(ii = i-1; ii > i-gap-1; ii--) {
            int_result[ii] = NA_INTEGER; 
          }
        }
      } else {
        /* nr-2 is first position to fill fromLast=TRUE */
        int_result[nr-1] = int_x[nr-1];
        for(i=nr-2; i>=0; i--) {
          int_result[i] = int_x[i];
          if(int_result[i] == NA_INTEGER) {
            if(limit > gap)
              int_result[i] = int_result[i+1];
            gap++;
          } else {
            if((int)gap > (int)maxgap) {
              for(ii = i+1; ii < i+gap+1; ii++) {
                int_result[ii] = NA_INTEGER; 
              }
            }
            gap=0;
          }
        }
        if((int)gap > (int)maxgap) {  /* check that we don't have leading trailing gap */
          for(ii = i+1; ii < i+gap+1; ii++) {
            int_result[ii] = NA_INTEGER; 
          }
        }
      }
      break;
    case REALSXP:
      real_x = REAL(x);
      real_result = REAL(result);
      if(!LOGICAL(fromLast)[0]) {   /* fromLast=FALSE */
        for(i=0; i < (_first+1); i++) {
          real_result[i] = real_x[i];
        }
        for(i=_first+1; i<nr; i++) {
          real_result[i] = real_x[i];
          if( ISNA(real_result[i]) || ISNAN(real_result[i])) {
            if(limit > gap)
              real_result[i] = real_result[i-1];
            gap++;
          } else {
            if((int)gap > (int)maxgap) {
              for(ii = i-1; ii > i-gap-1; ii--) {
                real_result[ii] = NA_REAL; 
              }
            }
            gap=0;
          }
        }
        if((int)gap > (int)maxgap) {  /* check that we don't have excessive trailing gap */
          for(ii = i-1; ii > i-gap-1; ii--) {
            real_result[ii] = NA_REAL; 
          }
        }
      } else {                      /* fromLast=TRUE */
        real_result[nr-1] = real_x[nr-1];
        for(i=nr-2; i>=0; i--) {
          real_result[i] = real_x[i];
          if(ISNA(real_result[i]) || ISNAN(real_result[i])) {
            if(limit > gap)
              real_result[i] = real_result[i+1];
            gap++;
          } else {
            if((int)gap > (int)maxgap) {
              for(ii = i+1; ii < i+gap+1; ii++) {
                real_result[ii] = NA_REAL; 
              }
            }
            gap=0;
          }
        }
        if((int)gap > (int)maxgap) {  /* check that we don't have leading trailing gap */
          for(ii = i+1; ii < i+gap+1; ii++) {
            real_result[ii] = NA_REAL; 
          }
        }
      }
      break;
    default:
      error("unsupported type");
      break;
  }
  if(isXts(x)) {
    setAttrib(result, R_DimSymbol, getAttrib(x, R_DimSymbol));
    setAttrib(result, R_DimNamesSymbol, getAttrib(x, R_DimNamesSymbol));
    setAttrib(result, xts_IndexSymbol, getAttrib(x, xts_IndexSymbol));
    copy_xtsCoreAttributes(x, result);
    copy_xtsAttributes(x, result);
  }
  UNPROTECT(P);
  return(result);
}

//.........这里部分代码省略.........
	for (i=0; i<n_dimensions; ++i){
		REAL(xDimSteps)[i] = dim_steps[i];
	}
	list_index++;
	SET_VECTOR_ELT(rtnList, list_index, xDimSteps);
	SET_STRING_ELT(listNames, list_index, mkChar("dimSteps"));


	/* Loop over the dimensions to grab the remaining info ... */
	for( i=0; i < n_dimensions; ++i ){
	//
	/* get (and print) the dimension names for all dimensions*
	... remember that since miget_dimension_name calls strdup which, in turn,
	... calls malloc to get memory for the new string -- we need to call "mifree" on
	... our pointer to release that memory.  */
		result = miget_dimension_name(dimensions[i], &dim_name);
		
		// do we have a time dimension?
		if ( !strcmp(dim_name, "time") ) { 
			time_dim_exists = TRUE;
			n_frames = ( time_dim_exists ) ? dim_sizes[0] : 0;
		}
		
		// store the dimension name and units
		SET_STRING_ELT(xDimNames, i, mkChar(dim_name));
		mifree_name(dim_name);
		
		result = miget_dimension_units(dimensions[i], &dim_units);
		SET_STRING_ELT(xDimUnits, i, mkChar(dim_units));
		mifree_name(dim_units);
		
	}
	/* add number of frames to return list */
	list_index++;
	SET_VECTOR_ELT(rtnList, list_index, ScalarInteger(n_frames));
	SET_STRING_ELT(listNames, list_index, mkChar("nFrames"));
	
	// add dim names to return list
	list_index++;
	SET_VECTOR_ELT(rtnList, list_index, xDimNames);
	SET_STRING_ELT(listNames, list_index, mkChar("dimNames"));
	// add dim units
	list_index++;
	SET_VECTOR_ELT(rtnList, list_index, xDimUnits);
	SET_STRING_ELT(listNames, list_index, mkChar("dimUnits"));


	/* get the dimension OFFSETS values for the TIME dimension */
	if ( time_dim_exists ) {

		PROTECT( xTimeOffsets=allocVector(REALSXP,n_frames) );
		n_protects++;
		result = miget_dimension_offsets(dimensions[0], n_frames, 0, time_offsets);
		if ( result == MI_ERROR ) { error("Error returned from miget_dimension_offsets.\n"); }
		/* add to R vector ... */
		for (i=0; i < n_frames; ++i) {
			REAL(xTimeOffsets)[i] = time_offsets[i];
//			if (R_DEBUG_mincIO) Rprintf("Time offset[%d] =  %g\n", i, time_offsets[i]);
		}
		list_index++;
		SET_VECTOR_ELT(rtnList, list_index, xTimeOffsets);
		SET_STRING_ELT(listNames, list_index, mkChar("timeOffsets"));

		/* get the dimension WIDTH values for the TIME dimension */
		PROTECT( xTimeWidths=allocVector(REALSXP,n_frames) );
		n_protects++;
	
		result = miget_dimension_widths(dimensions[0], MI_ORDER_FILE, n_frames, 0, time_widths);
		if ( result == MI_ERROR ) { error("Error returned from miget_dimension_widths.\n"); }
		/* add to R vector ... */
		for (i=0; i<n_frames; ++i) {
			REAL(xTimeWidths)[i] = time_widths[i];
//			if (R_DEBUG_mincIO) Rprintf("Time width[%d] =  %g\n", i, time_widths[i]);
		}
		list_index++;
		SET_VECTOR_ELT(rtnList, list_index, xTimeWidths);
		SET_STRING_ELT(listNames, list_index, mkChar("timeWidths"));
	}



	// free heap memory
	free(dimensions);


	/* close volume */
	miclose_volume(minc_volume);


	/* attach the list component names to the list */
	setAttrib(rtnList, R_NamesSymbol, listNames);


	/* remove gc collection protection */
	UNPROTECT(n_protects);

   /* return */
	if ( R_DEBUG_mincIO ) Rprintf("get_volume_info: returning ...\n");
   return(rtnList);
}

SEXP do_rgb(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP c, r, g, b, a, nam;
    int OP, i, l_max, nr, ng, nb, na;
    Rboolean max_1 = FALSE;
    double mV = 0.0; /* -Wall */

    checkArity(op, args);
    OP = PRIMVAL(op);
    if(OP) {/* op == 1:  rgb256() :*/
	PROTECT(r = coerceVector(CAR(args), INTSXP)); args = CDR(args);
	PROTECT(g = coerceVector(CAR(args), INTSXP)); args = CDR(args);
	PROTECT(b = coerceVector(CAR(args), INTSXP)); args = CDR(args);
	PROTECT(a = coerceVector(CAR(args), INTSXP)); args = CDR(args);
    }
    else {
	PROTECT(r = coerceVector(CAR(args), REALSXP)); args = CDR(args);
	PROTECT(g = coerceVector(CAR(args), REALSXP)); args = CDR(args);
	PROTECT(b = coerceVector(CAR(args), REALSXP)); args = CDR(args);
	PROTECT(a = coerceVector(CAR(args), REALSXP)); args = CDR(args);
	mV = asReal(CAR(args));			       args = CDR(args);
	max_1 = (mV == 1.);
    }

    nr = LENGTH(r); ng = LENGTH(g); nb = LENGTH(b); na = LENGTH(a);
    if (nr <= 0 || ng <= 0 || nb <= 0 || na <= 0) {
	UNPROTECT(4);
	return(allocVector(STRSXP, 0));
    }
    l_max = nr;
    if (l_max < ng) l_max = ng;
    if (l_max < nb) l_max = nb;
    if (l_max < na) l_max = na;

    PROTECT(nam = coerceVector(CAR(args), STRSXP)); args = CDR(args);
    if (length(nam) != 0 && length(nam) != l_max)
	errorcall(call, _("invalid names vector"));
    PROTECT(c = allocVector(STRSXP, l_max));

#define _R_set_c_RGBA(_R,_G,_B,_A)				\
    for (i = 0; i < l_max; i++)				\
	SET_STRING_ELT(c, i, mkChar(RGBA2rgb(_R,_G,_B,_A)))

    if(OP) { /* OP == 1:  rgb256() :*/
	_R_set_c_RGBA(CheckColor(INTEGER(r)[i%nr]),
		      CheckColor(INTEGER(g)[i%ng]),
		      CheckColor(INTEGER(b)[i%nb]),
		      CheckAlpha(INTEGER(a)[i%na]));
    }
    else if(max_1) {
	_R_set_c_RGBA(ScaleColor(REAL(r)[i%nr]),
		      ScaleColor(REAL(g)[i%ng]),
		      ScaleColor(REAL(b)[i%nb]),
		      ScaleAlpha(REAL(a)[i%na]));
    }
    else { /* maxColorVal not in {1, 255} */
	_R_set_c_RGBA(ScaleColor(REAL(r)[i%nr] / mV),
		      ScaleColor(REAL(g)[i%ng] / mV),
		      ScaleColor(REAL(b)[i%nb] / mV),
		      ScaleAlpha(REAL(a)[i%na] / mV));
    }
    if (length(nam) != 0)
	setAttrib(c, R_NamesSymbol, nam);
    UNPROTECT(6);
    return c;
}

/* generate a graph with given node ordering and arc probability. */
SEXP ordered_graph(SEXP nodes, SEXP num, SEXP prob) {

int i = 0, j = 0, k = 0, nnodes = LENGTH(nodes), *a = NULL, *n = INTEGER(num);
double *p = REAL(prob);
SEXP list, res, args, argnames, amat, arcs, cached, debug2, null, temp;

  /* a fake debug argument (set to FALSE) for cache_structure(). */
  PROTECT(debug2 = allocVector(LGLSXP, 1));
  LOGICAL(debug2)[0] = FALSE;

  /* the list of optional arguments. */
  PROTECT(argnames = allocVector(STRSXP, 1));
  SET_STRING_ELT(argnames, 0, mkChar("prob"));

  PROTECT(args = allocVector(VECSXP, 1));
  setAttrib(args, R_NamesSymbol, argnames);
  SET_VECTOR_ELT(args, 0, prob);

  /* allocate and initialize the adjacency matrix. */
  PROTECT(amat = allocMatrix(INTSXP, nnodes, nnodes));
  a = INTEGER(amat);
  memset(a, '\0', nnodes * nnodes * sizeof(int));

  GetRNGstate();

#define ORDERED_AMAT(prob) \
      for (i = 0; i < nnodes; i++) \
        for (j = i + 1; j < nnodes; j++) \
          if (unif_rand() < prob) \
            a[CMC(i, j, nnodes)] = 1; \
          else \
            a[CMC(i, j, nnodes)] = 0; \

  /* return a list if more than one bn is generated. */
  if (*n > 1) {

    PROTECT(list = allocVector(VECSXP, *n));
    PROTECT(null = allocVector(NILSXP, 1));

    /* generate the "bn" structure, with dummy NULLs for the "arcs" and
     * "nodes" elements (which will be initialized later on). */
    PROTECT(res = bn_base_structure(nodes, args, null, null, 0, "none", "ordered"));

    for (k = 0; k < *n; k++) {

      /* sample each arc in the upper-triangular portion of the adjacency matrix
       * (so that node ordering is conserved) with the specified probability. */
      ORDERED_AMAT(*p);

      /* generate the arc set and the cached information form the adjacency
       * matrix. */
      PROTECT(arcs = amat2arcs(amat, nodes));
      PROTECT(cached = cache_structure(nodes, amat, debug2));
      SET_VECTOR_ELT(res, 1, cached);
      SET_VECTOR_ELT(res, 2, arcs);

      /* save the structure in the list. */
      PROTECT(temp = duplicate(res));
      SET_VECTOR_ELT(list, k, temp);

      UNPROTECT(3);

    }/*FOR*/

    PutRNGstate();

    UNPROTECT(7);
    return list;

  }/*THEN*/
  else {

    /* sample each arc in the upper-triangular portion of the adjacency matrix
     * (so that node ordering is conserved) with the specified probability. */
    ORDERED_AMAT(*p);

    /* generate the arc set and the cached information form the adjacency
     * matrix. */
    PROTECT(arcs = amat2arcs(amat, nodes));
    PROTECT(cached = cache_structure(nodes, amat, debug2));

    /* generate the "bn" structure. */
    PROTECT(res = bn_base_structure(nodes, args, arcs, cached, 0, "none", "ordered"));

    PutRNGstate();

    UNPROTECT(7);
    return res;

  }/*ELSE*/

}/*ORDERED_GRAPH*/

SEXP Scatt( const SEXP data_sxp,		// data matrix  
			const SEXP clusters_sxp,	// vector with information about clusters (object num. -> cluster)
			const SEXP clust_num_sxp, // number of clusters (table with one element)
			const SEXP choosen_metric_sxp	// number representing choosen metric
		  )
{

	// additional variables especially needed in loops and in functions as parameters
	int i, j, obj_num, dim_num, clust_num, protect_num;
	// define distance between choosen objects
	double dist;

	protect_num = 0;
	
	// declaration of intracluster distances (vectors)
	SEXP mean_sxp, variance_sxp;
	double *mean, *variance;
	
	// matrix (and pointer to the table) with cluster centers
	SEXP cluster_center_sxp, cluster_variance_sxp, cluster_size_sxp;
	double *cluster_center, *cluster_variance;
	int *cluster_size;

	// table - which object belongs to which cluster
	int *cluster_tab = INTEGER(clusters_sxp);

	// get information about data matrix
	SEXP dim = NILSXP;
	PROTECT( dim = getAttrib(data_sxp, R_DimSymbol) );
	protect_num++;

	obj_num = INTEGER(dim)[0];
	dim_num = INTEGER(dim)[1];

	// and number of clusters
	clust_num = INTEGER(clust_num_sxp)[0];

	// compute mean 
	PROTECT( mean_sxp = clv_mean(data_sxp, obj_num, dim_num) );
	protect_num++;

	// and variance
	PROTECT( variance_sxp = clv_variance(data_sxp, obj_num, dim_num, mean_sxp) );
	protect_num++;
	variance = REAL(variance_sxp);

	// vector with information about size of each cluster
	PROTECT( cluster_size_sxp = clv_clustersSize(clusters_sxp, clust_num) );
	protect_num++;
	cluster_size = INTEGER(cluster_size_sxp);
	
	PROTECT( cluster_center_sxp = clv_clusterCenters(data_sxp, obj_num, dim_num, clust_num, cluster_tab, cluster_size) );
	protect_num++;

	PROTECT( cluster_variance_sxp = clv_clusterVariance(data_sxp, obj_num, dim_num, clust_num, cluster_tab, cluster_size, cluster_center_sxp) );
	protect_num++;
	cluster_variance = REAL(cluster_variance_sxp);

	double sum_cls_var_norm = 0, tmp;
	int pos;

	// compute "stdev" value ( sum[ forall k in {1, ... ,cluster num.} ] ||sigma(C_k)||)
	for(i=0; i<clust_num; i++)
	{
		sum_cls_var_norm += clv_norm(cluster_variance, i, dim_num, clust_num); 
	}
	
	// compute norm of variance of dataset (||sigma(X)||)
	double var_norm = clv_norm(variance, 0, dim_num, 1); ;
	
	SEXP Scatt, stdev;
	PROTECT( Scatt = allocVector(REALSXP, 1) );
	protect_num++;
	PROTECT( stdev = allocVector(REALSXP, 1) );
	protect_num++;
	
	REAL(Scatt)[0] = sum_cls_var_norm/(clust_num*var_norm);
	REAL(stdev)[0] = sqrt(sum_cls_var_norm)/clust_num;
	
	// time to gather all particular indicies into one result list 
	int list_elem_num = 3;
	SEXP result_list;
	PROTECT( result_list = allocVector(VECSXP, list_elem_num) );
	protect_num++;
	
	SEXP names;
	PROTECT( names = allocVector(STRSXP, list_elem_num) );
	protect_num++;
	
	pos = 0;
	SET_STRING_ELT( names, pos++, mkChar("Scatt") );
	SET_STRING_ELT( names, pos++, mkChar("stdev") );
	SET_STRING_ELT( names, pos++, mkChar("cluster.center") );

	setAttrib( result_list, R_NamesSymbol, names );	

	pos = 0;
	SET_VECTOR_ELT( result_list, pos++, Scatt );
	SET_VECTOR_ELT( result_list, pos++, stdev );
	SET_VECTOR_ELT( result_list, pos++, cluster_center_sxp );
//.........这里部分代码省略.........

/* an Ide-Cozman alternative for 2-nodes graphs. */
static SEXP ic_2nodes(SEXP nodes, SEXP num, SEXP burn_in, SEXP max_in_degree,
    SEXP max_out_degree, SEXP max_degree, SEXP connected, SEXP debug) {

int i = 0, *n = INTEGER(num), *a = NULL;
int *debuglevel = LOGICAL(debug);
double u = 0;
SEXP list, resA, resB, arcsA, arcsB, cachedA, cachedB;
SEXP amatA, amatB, args, argnames, false;

  /* the list of optional arguments. */
  PROTECT(argnames = allocVector(STRSXP, 4));
  SET_STRING_ELT(argnames, 0, mkChar("burn.in"));
  SET_STRING_ELT(argnames, 1, mkChar("max.in.degree"));
  SET_STRING_ELT(argnames, 2, mkChar("max.out.degree"));
  SET_STRING_ELT(argnames, 3, mkChar("max.degree"));

  PROTECT(args = allocVector(VECSXP, 4));
  setAttrib(args, R_NamesSymbol, argnames);
  SET_VECTOR_ELT(args, 0, burn_in);
  SET_VECTOR_ELT(args, 1, max_in_degree);
  SET_VECTOR_ELT(args, 2, max_out_degree);
  SET_VECTOR_ELT(args, 3, max_degree);

  /* allocate a FALSE variable. */
  PROTECT(false = allocVector(LGLSXP, 1));
  LOGICAL(false)[0] = FALSE;

  /* allocate and initialize the tow adjacency matrices. */
  PROTECT(amatA = allocMatrix(INTSXP, 2, 2));
  a = INTEGER(amatA);
  memset(a, '\0', sizeof(int) * 4);
  a[2] = 1;
  PROTECT(amatB = allocMatrix(INTSXP, 2, 2));
  a = INTEGER(amatB);
  memset(a, '\0', sizeof(int) * 4);
  a[1] = 1;
  /* generates the arc sets. */
  PROTECT(arcsA = amat2arcs(amatA, nodes));
  PROTECT(arcsB = amat2arcs(amatB, nodes));
  /* generate the cached node information. */
  PROTECT(cachedA = cache_structure(nodes, amatA, false));
  PROTECT(cachedB = cache_structure(nodes, amatB, false));
  /* generate the two "bn" structures. */
  PROTECT(resA = bn_base_structure(nodes, args, arcsA, cachedA, 0, "none", "empty"));
  PROTECT(resB = bn_base_structure(nodes, args, arcsB, cachedB, 0, "none", "empty"));

  if (*debuglevel > 0)
    Rprintf("* no burn-in required.\n");

  GetRNGstate();

  /* return a list if more than one bn is generated. */
  if (*n > 1) {

    PROTECT(list = allocVector(VECSXP, *n));
    for (i = 0; i < *n; i++) {

      if (*debuglevel > 0)
        Rprintf("* current model (%d):\n", i + 1);

      /* sample which graph to return. */
      u = unif_rand();

      if (u <= 0.5) {

        /* pick the graph with A -> B. */
        SET_VECTOR_ELT(list, i, resA);

        /* print the model string to allow a sane debugging experience. */
        if (*debuglevel > 0)
          print_modelstring(resA);

      }/*THEN*/
      else {

        /* pick the graph with B -> A. */
        SET_VECTOR_ELT(list, i, resB);

        /* print the model string to allow a sane debugging experience. */
        if (*debuglevel > 0)
          print_modelstring(resB);

      }/*ELSE*/

    }/*FOR*/

    PutRNGstate();

    UNPROTECT(12);
    return list;

  }/*THEN*/
  else {

    if (*debuglevel > 0)
      Rprintf("* current model (1):\n");

    /* sample which graph to return. */
    u = unif_rand();
//.........这里部分代码省略.........

static SEXP c_ide_cozman(SEXP nodes, SEXP num, SEXP burn_in, SEXP max_in_degree,
    SEXP max_out_degree, SEXP max_degree, SEXP connected, SEXP debug) {

int i = 0, k = 0, nnodes = LENGTH(nodes), *n = INTEGER(num);
int changed = 0, *work = NULL, *arc = NULL, *a = NULL, *burn = INTEGER(burn_in);
int *degree = NULL, *in_degree = NULL, *out_degree = NULL;
int *debuglevel = LOGICAL(debug), *cozman = LOGICAL(connected);
double *max_in = REAL(max_in_degree), *max_out = REAL(max_out_degree),
  *max = REAL(max_degree);
SEXP list, res, args, argnames, amat, arcs, cached, debug2, null, temp;
char *label = (*cozman > 0) ? "ic-dag" : "melancon";

  /* a fake debug argument (set to FALSE) for cache_structure(). */
  PROTECT(debug2 = allocVector(LGLSXP, 1));
  LOGICAL(debug2)[0] = FALSE;

  /* the list of optional arguments. */
  PROTECT(argnames = allocVector(STRSXP, 4));
  SET_STRING_ELT(argnames, 0, mkChar("burn.in"));
  SET_STRING_ELT(argnames, 1, mkChar("max.in.degree"));
  SET_STRING_ELT(argnames, 2, mkChar("max.out.degree"));
  SET_STRING_ELT(argnames, 3, mkChar("max.degree"));

  PROTECT(args = allocVector(VECSXP, 4));
  setAttrib(args, R_NamesSymbol, argnames);
  SET_VECTOR_ELT(args, 0, burn_in);
  SET_VECTOR_ELT(args, 1, max_in_degree);
  SET_VECTOR_ELT(args, 2, max_out_degree);
  SET_VECTOR_ELT(args, 3, max_degree);

  /* allocate and initialize the adjacency matrix. */
  PROTECT(amat = allocMatrix(INTSXP, nnodes, nnodes));
  a = INTEGER(amat);
  memset(a, '\0', nnodes * nnodes * sizeof(int));

  /* initialize a simple ordered tree with n nodes, where all nodes
   * have just one parent, except the first one that does not have
   * any parent. */
  for (i = 1; i < nnodes; i++)
    a[CMC(i - 1, i, nnodes)] = 1;

  /* allocate the arrays needed by SampleNoReplace. */
  arc = alloc1dcont(2);
  work = alloc1dcont(nnodes);

  /* allocate and initialize the degree arrays. */
  degree = alloc1dcont(nnodes);
  in_degree = alloc1dcont(nnodes);
  out_degree = alloc1dcont(nnodes);

  for (i = 0; i < nnodes; i++) {

    in_degree[i] = out_degree[i] = 1;
    degree[i] = 2;

  }/*FOR*/
  in_degree[0] = out_degree[nnodes - 1] = 0;
  degree[0] = degree[nnodes - 1] = 1;

  GetRNGstate();

  /* wait for the markov chain monte carlo simulation to reach stationarity. */
  for (k = 0; k < *burn; k++) {

    if (*debuglevel > 0)
      Rprintf("* current model (%d):\n", k + 1);

    changed = ic_logic(a, nodes, &nnodes, arc, work, degree, max, in_degree, max_in,
                out_degree, max_out, cozman, debuglevel);

    /* print the model string to allow a sane debugging experience; note that this
     * has a huge impact on performance, so use it with care. */
    if ((*debuglevel > 0) && (changed)) {

      PROTECT(null = allocVector(NILSXP, 1));
      PROTECT(res = bn_base_structure(nodes, args, null, null, 0, "none", label));
      PROTECT(arcs = amat2arcs(amat, nodes));
      PROTECT(cached = cache_structure(nodes, amat, debug2));
      SET_VECTOR_ELT(res, 1, cached);
      SET_VECTOR_ELT(res, 2, arcs);
      print_modelstring(res);
      UNPROTECT(4);

    }/*THEN*/

  }/*FOR*/

#define UPDATE_NODE_CACHE(cur) \
          if (*debuglevel > 0) \
            Rprintf("  > updating cached information about node %s.\n", NODE(cur)); \
          memset(work, '\0', nnodes * sizeof(int)); \
          PROTECT(temp = c_cache_partial_structure(cur, nodes, amat, work, debug2)); \
          SET_VECTOR_ELT(cached, cur, temp); \
          UNPROTECT(1);

  /* return a list if more than one bn is generated. */
  if (*n > 1) {

    if (*debuglevel > 0)
      Rprintf("* end of the burn-in iterations.\n");
//.........这里部分代码省略.........

/* generate an empty graph. */
SEXP empty_graph(SEXP nodes, SEXP num) {

int i = 0, nnodes = LENGTH(nodes), *n = INTEGER(num);
SEXP list, res, args, arcs, cached;
SEXP dimnames, colnames, elnames, base, base2;

  /* an empty list of optional arguments. */
  PROTECT(args = allocVector(VECSXP, 0));

  /* names for the arc set columns. */
  PROTECT(dimnames = allocVector(VECSXP, 2));
  PROTECT(colnames = allocVector(STRSXP, 2));
  SET_STRING_ELT(colnames, 0, mkChar("from"));
  SET_STRING_ELT(colnames, 1, mkChar("to"));
  SET_VECTOR_ELT(dimnames, 1, colnames);

  /* names for the cached information. */
  PROTECT(elnames = allocVector(STRSXP, 4));
  SET_STRING_ELT(elnames, 0, mkChar("mb"));
  SET_STRING_ELT(elnames, 1, mkChar("nbr"));
  SET_STRING_ELT(elnames, 2, mkChar("parents"));
  SET_STRING_ELT(elnames, 3, mkChar("children"));

  /* allocate and initialize the arc set. */
  PROTECT(arcs = allocMatrix(STRSXP, 0, 2));
  setAttrib(arcs, R_DimNamesSymbol, dimnames);

  /* allocate and initialize nodes' cached information. */
  PROTECT(base2 = allocVector(STRSXP, 0));
  PROTECT(base = allocVector(VECSXP, 4));
  setAttrib(base, R_NamesSymbol, elnames);

  PROTECT(cached = allocVector(VECSXP, nnodes));
  setAttrib(cached, R_NamesSymbol, nodes);

  for (i = 0; i < 4; i++)
    SET_VECTOR_ELT(base, i, base2);
  for (i = 0; i < nnodes; i++)
    SET_VECTOR_ELT(cached, i, base);

  /* generate the "bn" structure. */
  PROTECT(res = bn_base_structure(nodes, args, arcs, cached, 0, "none", "empty"));

  /* return a list if more than one bn is generated. */
  if (*n > 1) {

    PROTECT(list = allocVector(VECSXP, *n));
    for (i = 0; i < *n; i++)
      SET_VECTOR_ELT(list, i, res);

    UNPROTECT(10);
    return list;

  }/*THEN*/
  else {

    UNPROTECT(9);
    return res;

  }/*ELSE*/

}/*EMPTY_GRAPH*/

//.........这里部分代码省略.........
				imat[j][i] +=  (wtave/d2)*
				    ((cmat[i][j] - temp*cmat2[i][j]) -
				    temp2*(a[j]-temp*a2[j]));
    		            }
    		        }

		    for (i=0; i<nvar; i++) { /*in anticipation */
			a2[i] =0;
			for (j=0; j<nvar; j++) cmat2[i][j] =0;
		        }
	            }
		}
	    }   /* end  of accumulation loop  */

	/* am I done?
	**   update the betas and test for convergence
	*/
	*flag = cholesky2(imat, nvar, toler);

	if (fabs(1-(loglik[1]/newlk))<= eps && halving==0) { /* all done */
	    loglik[1] = newlk;
	    chinv2(imat, nvar);     /* invert the information matrix */
	    for (i=0; i<nvar; i++) {
		beta[i] = newbeta[i]*scale[i];
		u[i] /= scale[i];
		imat[i][i] *= scale[i]*scale[i];
		for (j=0; j<i; j++) {
		    imat[j][i] *= scale[i]*scale[j];
		    imat[i][j] = imat[j][i];
		    }
	    }
	    goto finish;
	}

	if (*iter== maxiter) break;  /*skip the step halving calc*/

	if (newlk < loglik[1])   {    /*it is not converging ! */
		halving =1;
		for (i=0; i<nvar; i++)
		    newbeta[i] = (newbeta[i] + beta[i]) /2; /*half of old increment */
		}
	else {
	    halving=0;
	    loglik[1] = newlk;
	    chsolve2(imat,nvar,u);
	    j=0;
	    for (i=0; i<nvar; i++) {
		beta[i] = newbeta[i];
		newbeta[i] = newbeta[i] +  u[i];
	        }
	    }
	}   /* return for another iteration */

    /*
    ** We end up here only if we ran out of iterations 
    */
    loglik[1] = newlk;
    chinv2(imat, nvar);
    for (i=0; i<nvar; i++) {
	beta[i] = newbeta[i]*scale[i];
	u[i] /= scale[i];
	imat[i][i] *= scale[i]*scale[i];
	for (j=0; j<i; j++) {
	    imat[j][i] *= scale[i]*scale[j];
	    imat[i][j] = imat[j][i];
	    }
	}
    *flag = 1000;


finish:
    /*
    ** create the output list
    */
    PROTECT(rlist= allocVector(VECSXP, 8));
    SET_VECTOR_ELT(rlist, 0, beta2);
    SET_VECTOR_ELT(rlist, 1, means2);
    SET_VECTOR_ELT(rlist, 2, u2);
    SET_VECTOR_ELT(rlist, 3, imat2);
    SET_VECTOR_ELT(rlist, 4, loglik2);
    SET_VECTOR_ELT(rlist, 5, sctest2);
    SET_VECTOR_ELT(rlist, 6, iter2);
    SET_VECTOR_ELT(rlist, 7, flag2);
    

    /* add names to the objects */
    PROTECT(rlistnames = allocVector(STRSXP, 8));
    SET_STRING_ELT(rlistnames, 0, mkChar("coef"));
    SET_STRING_ELT(rlistnames, 1, mkChar("means"));
    SET_STRING_ELT(rlistnames, 2, mkChar("u"));
    SET_STRING_ELT(rlistnames, 3, mkChar("imat"));
    SET_STRING_ELT(rlistnames, 4, mkChar("loglik"));
    SET_STRING_ELT(rlistnames, 5, mkChar("sctest"));
    SET_STRING_ELT(rlistnames, 6, mkChar("iter"));
    SET_STRING_ELT(rlistnames, 7, mkChar("flag"));
    setAttrib(rlist, R_NamesSymbol, rlistnames);

    unprotect(nprotect+2);
    return(rlist);
    }

//.........这里部分代码省略.........
	    }
	}
    }

    PROTECT(args);

    drop = 1;
    ExtractDropArg(args, &drop);
    x = CAR(args);

    /* This was intended for compatibility with S, */
    /* but in fact S does not do this. */
    /* FIXME: replace the test by isNull ... ? */

    if (x == R_NilValue) {
	UNPROTECT(1);
	return x;
    }
    subs = CDR(args);
    nsubs = length(subs); /* Will be short */
    type = TYPEOF(x);

    /* Here coerce pair-based objects into generic vectors. */
    /* All subsetting takes place on the generic vector form. */

    ax = x;
    if (isVector(x))
	PROTECT(ax);
    else if (isPairList(x)) {
	SEXP dim = getAttrib(x, R_DimSymbol);
	int ndim = length(dim);
	if (ndim > 1) {
	    PROTECT(ax = allocArray(VECSXP, dim));
	    setAttrib(ax, R_DimNamesSymbol, getAttrib(x, R_DimNamesSymbol));
	    setAttrib(ax, R_NamesSymbol, getAttrib(x, R_DimNamesSymbol));
	}
	else {
	    PROTECT(ax = allocVector(VECSXP, length(x)));
	    setAttrib(ax, R_NamesSymbol, getAttrib(x, R_NamesSymbol));
	}
	for(px = x, i = 0 ; px != R_NilValue ; px = CDR(px))
	    SET_VECTOR_ELT(ax, i++, CAR(px));
    }
    else errorcall(call, R_MSG_ob_nonsub, type2char(TYPEOF(x)));

    /* This is the actual subsetting code. */
    /* The separation of arrays and matrices is purely an optimization. */

    if(nsubs < 2) {
	SEXP dim = getAttrib(x, R_DimSymbol);
	int ndim = length(dim);
	PROTECT(ans = VectorSubset(ax, (nsubs == 1 ? CAR(subs) : R_MissingArg),
				   call));
	/* one-dimensional arrays went through here, and they should
	   have their dimensions dropped only if the result has
	   length one and drop == TRUE
	*/
	if(ndim == 1) {
	    SEXP attr, attrib, nattrib;
	    int len = length(ans);

	    if(!drop || len > 1) {
		// must grab these before the dim is set.
		SEXP nm = PROTECT(getAttrib(ans, R_NamesSymbol));
		PROTECT(attr = allocVector(INTSXP, 1));
		INTEGER(attr)[0] = length(ans);

//.........这里部分代码省略.........
        {
            _glfwInputError(GLFW_API_UNAVAILABLE,
                            "WGL: OpenGL ES requested but WGL_ARB_create_context_es2_profile is unavailable");
            return GLFW_FALSE;
        }
    }

    if (_glfw.wgl.ARB_create_context)
    {
        int index = 0, mask = 0, flags = 0;

        if (ctxconfig->client == GLFW_OPENGL_API)
        {
            if (ctxconfig->forward)
                flags |= WGL_CONTEXT_FORWARD_COMPATIBLE_BIT_ARB;

            if (ctxconfig->profile == GLFW_OPENGL_CORE_PROFILE)
                mask |= WGL_CONTEXT_CORE_PROFILE_BIT_ARB;
            else if (ctxconfig->profile == GLFW_OPENGL_COMPAT_PROFILE)
                mask |= WGL_CONTEXT_COMPATIBILITY_PROFILE_BIT_ARB;
        }
        else
            mask |= WGL_CONTEXT_ES2_PROFILE_BIT_EXT;

        if (ctxconfig->debug)
            flags |= WGL_CONTEXT_DEBUG_BIT_ARB;

        if (ctxconfig->robustness)
        {
            if (_glfw.wgl.ARB_create_context_robustness)
            {
                if (ctxconfig->robustness == GLFW_NO_RESET_NOTIFICATION)
                {
                    setAttrib(WGL_CONTEXT_RESET_NOTIFICATION_STRATEGY_ARB,
                              WGL_NO_RESET_NOTIFICATION_ARB);
                }
                else if (ctxconfig->robustness == GLFW_LOSE_CONTEXT_ON_RESET)
                {
                    setAttrib(WGL_CONTEXT_RESET_NOTIFICATION_STRATEGY_ARB,
                              WGL_LOSE_CONTEXT_ON_RESET_ARB);
                }

                flags |= WGL_CONTEXT_ROBUST_ACCESS_BIT_ARB;
            }
        }

        if (ctxconfig->release)
        {
            if (_glfw.wgl.ARB_context_flush_control)
            {
                if (ctxconfig->release == GLFW_RELEASE_BEHAVIOR_NONE)
                {
                    setAttrib(WGL_CONTEXT_RELEASE_BEHAVIOR_ARB,
                              WGL_CONTEXT_RELEASE_BEHAVIOR_NONE_ARB);
                }
                else if (ctxconfig->release == GLFW_RELEASE_BEHAVIOR_FLUSH)
                {
                    setAttrib(WGL_CONTEXT_RELEASE_BEHAVIOR_ARB,
                              WGL_CONTEXT_RELEASE_BEHAVIOR_FLUSH_ARB);
                }
            }
        }

        if (ctxconfig->noerror)
        {
            if (_glfw.wgl.ARB_create_context_no_error)