int
main()
{
#if N & 1
	long double	value = 0;
#else
	double		value = 0;
#endif
#if N < 5
	int		exp = 0;
#endif

#if N == 1
	return ldexpl(value, exp) != 0;
#endif
#if N == 2
	return ldexp(value, exp) != 0;
#endif
#if N == 3
	return frexpl(value, &exp) != 0;
#endif
#if N == 4
	return frexp(value, &exp) != 0;
#endif
#if N == 5
	return isnan(value);
#endif
#if N == 6
	return isnan(value);
#endif
#if N == 7
	return copysign(1.0, value) < 0;
#endif
#if N == 8
	return signbit(value);
#endif
}

/* Tests along the real and imaginary axes. */
void
test_axes(void)
{
	static const long double nums[] = {
		-2, -1, -0.5, 0.5, 1, 2
	};
	long double complex z;
	int i;

	for (i = 0; i < sizeof(nums) / sizeof(nums[0]); i++) {
		/* Real axis */
		z = CMPLXL(nums[i], 0.0);
		if (fabs(nums[i]) <= 1) {
			testall_tol(cacosh, z, CMPLXL(0.0, acos(nums[i])), 1);
			testall_tol(cacos, z, CMPLXL(acosl(nums[i]), -0.0), 1);
			testall_tol(casin, z, CMPLXL(asinl(nums[i]), 0.0), 1);
			testall_tol(catanh, z, CMPLXL(atanh(nums[i]), 0.0), 1);
		} else {
			testall_tol(cacosh, z,
				    CMPLXL(acosh(fabs(nums[i])),
					   (nums[i] < 0) ? pi : 0), 1);
			testall_tol(cacos, z,
				    CMPLXL((nums[i] < 0) ? pi : 0,
					   -acosh(fabs(nums[i]))), 1);
			testall_tol(casin, z,
				    CMPLXL(copysign(pi / 2, nums[i]),
					   acosh(fabs(nums[i]))), 1);
			testall_tol(catanh, z,
				    CMPLXL(atanh(1 / nums[i]), pi / 2), 1);
		}
		testall_tol(casinh, z, CMPLXL(asinh(nums[i]), 0.0), 1);
		testall_tol(catan, z, CMPLXL(atan(nums[i]), 0), 1);

		/* TODO: Test the imaginary axis. */
	}
}

double CGaussian::Deviance(const CDataset& kData, const Bag& kBag,
                           const double* kFuncEstimate) {
  double loss = 0.0;
  double weight = 0.0;

  unsigned long num_rows_in_set = kData.get_size_of_set();
#pragma omp parallel for schedule(static, get_array_chunk_size()) \
    reduction(+ : loss, weight) num_threads(get_num_threads())
  for (unsigned long i = 0; i < num_rows_in_set; i++) {
    const double tmp =
        (kData.y_ptr()[i] - kData.offset_ptr()[i] - kFuncEstimate[i]);
    loss += kData.weight_ptr()[i] * tmp * tmp;
    weight += kData.weight_ptr()[i];
  }

  // TODO: Check if weights are all zero for validation set
  if ((weight == 0.0) && (loss == 0.0)) {
    return nan("");
  } else if (weight == 0.0) {
    return copysign(HUGE_VAL, loss);
  }

  return loss / weight;
}

static Handle powerOf(TaskData *mdTaskData, Handle args)
{
    double x = real_arg1(args), y = real_arg2(args);
    /* Some of the special cases are defined and don't seem to match
       the C pow function (at least as implemented in MS C). */
    /* Maybe handle all this in ML? */
    if (isnan(x))
    {
        if (y == 0.0) return real_result(mdTaskData, 1.0);
        else return real_result(mdTaskData, notANumber);
    }
    else if (isnan(y)) return real_result(mdTaskData, y); /* i.e. nan. */
    else if (x == 0.0 && y < 0.0)
    {
        /* This case is not handled correctly in Solaris. It always
           returns -infinity. */
        int iy = (int)floor(y);
        /* If x is -0.0 and y is an odd integer the result is -infinity. */
        if (copysign(1.0, x) < 0.0 && (double)iy == y && (iy & 1))
            return real_result(mdTaskData, negInf); /* -infinity. */
        else return real_result(mdTaskData, posInf); /* +infinity. */
    }
    return real_result(mdTaskData, pow(x, y));
}

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

            if (k == 2) {
#pragma omp single
                {
                    for (j = 0; j < k; ++j) {
                        w[0] = *Q(0,j);
                        w[1] = *Q(1,j);

                        i = indx[0] - 1;
                        *Q(0,j) = w[i];
                        i = indx[1] - 1;
                        *Q(1,j) = w[i];
                    }
                }
            }
            else if (k != 1) {
                // Compute updated W.
                blasf77_scopy( &ik, &w[ib], &ione, &s[ib], &ione);

                // Initialize W(I) = Q(I,I)
                tmp = ldq + 1;
                blasf77_scopy( &ik, Q(ib,ib), &tmp, &w[ib], &ione);

                for (j = 0; j < k; ++j) {
                    magma_int_t i_tmp = min(j, ie);
                    for (i = ib; i < i_tmp; ++i)
                        w[i] = w[i] * ( *Q(i, j) / ( dlamda[i] - dlamda[j] ) );
                    i_tmp = max(j+1, ib);
                    for (i = i_tmp; i < ie; ++i)
                        w[i] = w[i] * ( *Q(i, j) / ( dlamda[i] - dlamda[j] ) );
                }

                for (i = ib; i < ie; ++i)
                    w[i] = copysign( sqrt( -w[i] ), s[i]);

#pragma omp barrier

                //reduce the number of used threads to have enough S workspace
                tot = min(n1, omp_get_num_threads());

                if (id < tot) {
                    ib = (  id   * rk) / tot + iil - 1;
                    ie = ((id+1) * rk) / tot + iil - 1;
                    ik = ie - ib;
                }
                else {
                    ib = -1;
                    ie = -1;
                    ik = -1;
                }

                // Compute eigenvectors of the modified rank-1 modification.
                for (j = ib; j < ie; ++j) {
                    for (i = 0; i < k; ++i)
                        s[id*k + i] = w[i] / *Q(i,j);
                    temp = magma_cblas_snrm2( k, s+id*k, 1 );
                    for (i = 0; i < k; ++i) {
                        magma_int_t iii = indx[i] - 1;
                        *Q(i,j) = s[id*k + iii] / temp;
                    }
                }
            }
        }
    }
    if (*info != 0)
        return *info;

double FPEnvironmentImpl::copySignImpl(double target, double source)
{
	return (float) copysign(target, source);
}

//.........这里部分代码省略.........
				r = nextafter(r, 0);
			feupdateenv(&env);
			return (r);
		case FE_DOWNWARD:
			if (z > 0.0)
				return (x * y);
			feholdexcept(&env);
			r = x * y;
			if (!fetestexcept(FE_INEXACT))
				r = nextafter(r, -INFINITY);
			feupdateenv(&env);
			return (r);
		default:	/* FE_UPWARD */
			if (z < 0.0)
				return (x * y);
			feholdexcept(&env);
			r = x * y;
			if (!fetestexcept(FE_INEXACT))
				r = nextafter(r, INFINITY);
			feupdateenv(&env);
			return (r);
		}
	}
	if (spread < -DBL_MANT_DIG) {
		feraiseexcept(FE_INEXACT);
		if (!isnormal(z))
			feraiseexcept(FE_UNDERFLOW);
		switch (oround) {
		case FE_TONEAREST:
			return (z);
		case FE_TOWARDZERO:
			if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
				return (z);
			else
				return (nextafter(z, 0));
		case FE_DOWNWARD:
			if (x > 0.0 ^ y < 0.0)
				return (z);
			else
				return (nextafter(z, -INFINITY));
		default:	/* FE_UPWARD */
			if (x > 0.0 ^ y < 0.0)
				return (nextafter(z, INFINITY));
			else
				return (z);
		}
	}

	/*
	 * Use Dekker's algorithm to perform the multiplication and
	 * subsequent addition in twice the machine precision.
	 * Arrange so that x * y = c + cc, and x * y + z = r + rr.
	 */
	fesetround(FE_TONEAREST);

	p = xs * split;
	hx = xs - p;
	hx += p;
	tx = xs - hx;

	p = ys * split;
	hy = ys - p;
	hy += p;
	ty = ys - hy;

	p = hx * hy;
	q = hx * ty + tx * hy;
	c = p + q;
	cc = p - c + q + tx * ty;

	zs = ldexp(zs, -spread);
	r = c + zs;
	s = r - c;
	rr = (c - (r - s)) + (zs - s) + cc;

	spread = ex + ey;
	if (spread + ilogb(r) > -1023) {
		fesetround(oround);
		r = r + rr;
	} else {
		/*
		 * The result is subnormal, so we round before scaling to
		 * avoid double rounding.
		 */
		#ifndef __SYMBIAN32__
		p = ldexp(copysign(0x1p-1022, r), -spread);
		#else

		p = ldexp(copysign(0, r), -spread);


		#endif //__SYMBIAN32__
		c = r + p;
		s = c - r;
		cc = (r - (c - s)) + (p - s) + rr;
		fesetround(oround);
		r = (c + cc) - p;
	}
	return (ldexp(r, spread));
}

PyObject*
_PyCode_ConstantKey(PyObject *op)
{
    PyObject *key;

    /* Py_None and Py_Ellipsis are singleton */
    if (op == Py_None || op == Py_Ellipsis
       || PyLong_CheckExact(op)
       || PyBool_Check(op)
       || PyBytes_CheckExact(op)
       || PyUnicode_CheckExact(op)
          /* code_richcompare() uses _PyCode_ConstantKey() internally */
       || PyCode_Check(op)) {
        key = PyTuple_Pack(2, Py_TYPE(op), op);
    }
    else if (PyFloat_CheckExact(op)) {
        double d = PyFloat_AS_DOUBLE(op);
        /* all we need is to make the tuple different in either the 0.0
         * or -0.0 case from all others, just to avoid the "coercion".
         */
        if (d == 0.0 && copysign(1.0, d) < 0.0)
            key = PyTuple_Pack(3, Py_TYPE(op), op, Py_None);
        else
            key = PyTuple_Pack(2, Py_TYPE(op), op);
    }
    else if (PyComplex_CheckExact(op)) {
        Py_complex z;
        int real_negzero, imag_negzero;
        /* For the complex case we must make complex(x, 0.)
           different from complex(x, -0.) and complex(0., y)
           different from complex(-0., y), for any x and y.
           All four complex zeros must be distinguished.*/
        z = PyComplex_AsCComplex(op);
        real_negzero = z.real == 0.0 && copysign(1.0, z.real) < 0.0;
        imag_negzero = z.imag == 0.0 && copysign(1.0, z.imag) < 0.0;
        /* use True, False and None singleton as tags for the real and imag
         * sign, to make tuples different */
        if (real_negzero && imag_negzero) {
            key = PyTuple_Pack(3, Py_TYPE(op), op, Py_True);
        }
        else if (imag_negzero) {
            key = PyTuple_Pack(3, Py_TYPE(op), op, Py_False);
        }
        else if (real_negzero) {
            key = PyTuple_Pack(3, Py_TYPE(op), op, Py_None);
        }
        else {
            key = PyTuple_Pack(2, Py_TYPE(op), op);
        }
    }
    else if (PyTuple_CheckExact(op)) {
        Py_ssize_t i, len;
        PyObject *tuple;

        len = PyTuple_GET_SIZE(op);
        tuple = PyTuple_New(len);
        if (tuple == NULL)
            return NULL;

        for (i=0; i < len; i++) {
            PyObject *item, *item_key;

            item = PyTuple_GET_ITEM(op, i);
            item_key = _PyCode_ConstantKey(item);
            if (item_key == NULL) {
                Py_DECREF(tuple);
                return NULL;
            }

            PyTuple_SET_ITEM(tuple, i, item_key);
        }

        key = PyTuple_Pack(3, Py_TYPE(op), op, tuple);
        Py_DECREF(tuple);
    }
    else if (PyFrozenSet_CheckExact(op)) {
        Py_ssize_t pos = 0;
        PyObject *item;
        Py_hash_t hash;
        Py_ssize_t i, len;
        PyObject *tuple, *set;

        len = PySet_GET_SIZE(op);
        tuple = PyTuple_New(len);
        if (tuple == NULL)
            return NULL;

        i = 0;
        while (_PySet_NextEntry(op, &pos, &item, &hash)) {
            PyObject *item_key;

            item_key = _PyCode_ConstantKey(item);
            if (item_key == NULL) {
                Py_DECREF(tuple);
                return NULL;
            }

            assert(i < len);
            PyTuple_SET_ITEM(tuple, i, item_key);
            i++;
//.........这里部分代码省略.........

double complex
csqrt(double complex z)
{
	double complex result;
	double a, b;
	double t;
	int scale;

	a = creal(z);
	b = cimag(z);

	/* Handle special cases. */
	if (z == 0)
		return (cpack(0, b));
	if (isinf(b))
		return (cpack(INFINITY, b));
	if (isnan(a)) {
		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
		return (cpack(a, t));	/* return NaN + NaN i */
	}
	if (isinf(a)) {
		/*
		 * csqrt(inf + NaN i)  = inf +  NaN i
		 * csqrt(inf + y i)    = inf +  0 i
		 * csqrt(-inf + NaN i) = NaN +- inf i
		 * csqrt(-inf + y i)   = 0   +  inf i
		 */
		if (signbit(a))
			return (cpack(fabs(b - b), copysign(a, b)));
		else
			return (cpack(a, copysign(b - b, b)));
	}
	/*
	 * The remaining special case (b is NaN) is handled just fine by
	 * the normal code path below.
	 */

	/* Scale to avoid overflow. */
	if (fabs(a) >= THRESH || fabs(b) >= THRESH) {
		a *= 0.25;
		b *= 0.25;
		scale = 1;
	} else {
		scale = 0;
	}

	/* Algorithm 312, CACM vol 10, Oct 1967. */
	if (a >= 0) {
		t = sqrt((a + hypot(a, b)) * 0.5);
		result = cpack(t, b / (2 * t));
	} else {
		t = sqrt((-a + hypot(a, b)) * 0.5);
		result = cpack(fabs(b) / (2 * t), copysign(t, b));
	}

	/* Rescale. */
	if (scale)
		return (result * 2);
	else
		return (result);
}

#if __MINGW32__

#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <wchar.h>
#include <float.h>
#include <assert.h>

static double zero = 0;
double Port::nan = copysign(NAN, 1.0);
double Port::infinity = 1 / zero;
double Port::dbl_max = 1.7976931348623157e308;
double Port::dbl_min = 5e-324;
longdouble Port::ldbl_max = LDBL_MAX;

struct PortInitializer
{
    PortInitializer();
};

static PortInitializer portinitializer;

PortInitializer::PortInitializer()
{
    assert(!signbit(Port::nan));

void test2(double x, double y)
{
  if (-tan(x-y) != tan(y-x))
    link_error ();

  if (-sin(x-y) != sin(y-x))
    link_error ();

  if (cos(-x*y) != cos(x*y))
    link_error ();

  if (cos(x*-y) != cos(x*y))
    link_error ();

  if (cos(-x/y) != cos(x/y))
    link_error ();

  if (cos(x/-y) != cos(x/y))
    link_error ();

  if (cos(-fabs(tan(x/-y))) != cos(tan(x/y)))
    link_error ();

  if (cos(y<10 ? -x : y) != cos(y<10 ? x : y))
    link_error ();

  if (cos(y<10 ? x : -y) != cos(y<10 ? x : y))
    link_error ();

  if (cos(y<10 ? -fabs(x) : tan(x<20 ? -x : -fabs(y)))
      != cos(y<10 ? x : tan(x<20 ? x : y)))
    link_error ();

  if (cos((y*=3, -x)) != cos((y*=3,x)))
    link_error ();

  if (cos((y*=2, -fabs(tan(x/-y)))) != cos((y*=2,tan(x/y))))
    link_error ();

  if (cos(copysign(x,y)) != cos(x))
    link_error ();

  if (cos(copysign(-fabs(x),y*=2)) != cos((y*=2,x)))
    link_error ();

  if (hypot (x, 0) != fabs(x))
    link_error ();

  if (hypot (0, x) != fabs(x))
    link_error ();

  if (hypot (x, x) != fabs(x) * __builtin_sqrt(2))
    link_error ();

  if (hypot (-x, y) != hypot (x, y))
    link_error ();

  if (hypot (x, -y) != hypot (x, y))
    link_error ();

  if (hypot (-x, -y) != hypot (x, y))
    link_error ();

  if (hypot (fabs(x), y) != hypot (x, y))
    link_error ();

  if (hypot (x, fabs(y)) != hypot (x, y))
    link_error ();

  if (hypot (fabs(x), fabs(y)) != hypot (x, y))
    link_error ();

  if (hypot (-fabs(-x), -fabs(fabs(fabs(-y)))) != hypot (x, y))
    link_error ();

  if (hypot (-x, 0) != fabs(x))
    link_error ();

  if (hypot (-x, x) != fabs(x) * __builtin_sqrt(2))
    link_error ();

  if (hypot (pure(x), -pure(x)) != fabs(pure(x)) * __builtin_sqrt(2))
    link_error ();

  if (hypot (tan(-x), tan(-fabs(y))) != hypot (tan(x), tan(y)))
    link_error ();

  if (fmin (fmax(x,y),y) != y)
    link_error ();

  if (fmin (fmax(y,x),y) != y)
    link_error ();

  if (fmin (x,fmax(x,y)) != x)
    link_error ();
  
  if (fmin (x,fmax(y,x)) != x)
    link_error ();
  
  if (fmax (fmin(x,y),y) != y)
//.........这里部分代码省略.........

bool _Stl_is_neg_nan(double x)    { return isnan(x) && ( copysign(1., x) < 0 );  }

// bool _Stl_is_neg_inf(double x)    { return _class(x) == FP_MINUS_INF; }
bool _Stl_is_neg_inf(double x)    { return _Stl_is_inf(x) && ( copysign(1., x) < 0 );  }

int
main (void)
{
  int result = 0;

  float i = INFINITY;
  float m = FLT_MAX;
  feclearexcept (FE_ALL_EXCEPT);
  if (nextafterf (m, i) != i)
    {
      puts ("nextafterf+ failed");
      ++result;
    }
  if (fetestexcept (FE_OVERFLOW) == 0)
    {
      puts ("nextafterf+ did not overflow");
      ++result;
    }
  feclearexcept (FE_ALL_EXCEPT);
  if (nextafterf (-m, -i) != -i)
    {
      puts ("nextafterf- failed");
      ++result;
    }
  if (fetestexcept (FE_OVERFLOW) == 0)
    {
      puts ("nextafterf- did not overflow");
      ++result;
    }

  i = 0;
  m = FLT_MIN;
  feclearexcept (FE_ALL_EXCEPT);
  i = nextafterf (m, i);
  if (i < 0 || i >= FLT_MIN)
    {
      puts ("nextafterf+ failed");
      ++result;
    }
  if (fetestexcept (FE_UNDERFLOW) == 0)
    {
      puts ("nextafterf+ did not underflow");
      ++result;
    }
  i = 0;
  feclearexcept (FE_ALL_EXCEPT);
  i = nextafterf (-m, -i);
  if (i > 0 || i <= -FLT_MIN)
    {
      puts ("nextafterf- failed");
      ++result;
    }
  if (fetestexcept (FE_UNDERFLOW) == 0)
    {
      puts ("nextafterf- did not underflow");
      ++result;
    }
  i = -INFINITY;
  feclearexcept (FE_ALL_EXCEPT);
  m = nextafterf (zero, inf);
  if (m < 0.0 || m >= FLT_MIN)
    {
      puts ("nextafterf+ failed");
      ++result;
    }
  if (fetestexcept (FE_UNDERFLOW) == 0)
    {
      puts ("nextafterf+ did not underflow");
      ++result;
    }
  feclearexcept (FE_ALL_EXCEPT);
  if (nextafterf (m, i) != 0.0)
    {
      puts ("nextafterf+ failed");
      ++result;
    }
  if (fetestexcept (FE_UNDERFLOW) == 0)
    {
      puts ("nextafterf+ did not underflow");
      ++result;
    }
  feclearexcept (FE_ALL_EXCEPT);
  m = nextafterf (copysignf (zero, -1.0), -inf);
  if (m > 0.0 || m <= -FLT_MIN)
    {
      puts ("nextafterf- failed");
      ++result;
    }
  if (fetestexcept (FE_UNDERFLOW) == 0)
    {
      puts ("nextafterf- did not underflow");
      ++result;
    }
  feclearexcept (FE_ALL_EXCEPT);
  if (nextafterf (m, -i) != 0.0)
    {
      puts ("nextafterf- failed");
      ++result;
    }
  if (fetestexcept (FE_UNDERFLOW) == 0)
//.........这里部分代码省略.........

int CDataset::extractFeatures(const CConfig& conf){

    int imgRow = this->img.at(0)->rows, imgCol = this->img.at(0)->cols;
    cv::Mat *integralMat;

    if(conf.learningMode != 1){
        if(conf.rgbFeature == 1){ // if got rgb image only, calc hog feature
            feature.clear();
            feature.resize(32);
            for(int i = 0; i < 32; ++i)
                feature.at(i) = new cv::Mat(imgRow, imgCol, CV_8UC1);

            cv::cvtColor(*img.at(0), *(feature.at(0)), CV_RGB2GRAY);

            cv::Mat I_x(imgRow, imgCol, CV_16SC1);
            cv::Mat I_y(imgRow, imgCol, CV_16SC1);


            cv::Sobel(*(feature.at(0)), I_x, CV_16S, 1, 0);
            cv::Sobel(*(feature.at(0)), I_y, CV_16S, 0, 1);

            cv::convertScaleAbs(I_x, *(feature[3]), 0.25);
            cv::convertScaleAbs(I_y, *(feature[4]), 0.25);

            // Orientation of gradients
            for(int  y = 0; y < img.at(0)->rows; y++)
                for(int  x = 0; x < img.at(0)->cols; x++) {
                    // Avoid division by zero
                    float tx = (float)I_x.at<short>(y, x) + (float)copysign(0.000001f, I_x.at<short>(y, x));
                    // Scaling [-pi/2 pi/2] -> [0 80*pi]
                    feature.at(1)->at<uchar>(y, x) = (uchar)(( atan((float)I_y.at<short>(y, x) / tx) + 3.14159265f / 2.0f ) * 80);
                    //std::cout << "scaling" << std::endl;
                    feature.at(2)->at<uchar>(y, x) = (uchar)sqrt((float)I_x.at<short>(y, x)* (float)I_x.at<short>(y, x) + (float)I_y.at<short>(y, x) * (float)I_y.at<short>(y, x));
                }

            // Magunitude of gradients
            for(int y = 0; y < img.at(0)->rows; y++)
                for(int x = 0; x < img.at(0)->cols; x++ ) {
                    feature.at(2)->at<uchar>(y, x) = (uchar)sqrt(I_x.at<short>(y, x)*I_x.at<short>(y, x) + I_y.at<short>(y, x) * I_y.at<short>(y, x));
                }

            hog.extractOBin(feature[1], feature[2], feature, 7);

            // calc I_xx I_yy
            cv::Sobel(*(feature.at(0)), I_x, CV_16S, 2, 0);
            cv::Sobel(*(feature.at(0)), I_y, CV_16S, 0, 2);

            cv::convertScaleAbs(I_x, *(feature[5]), 0.25);
            cv::convertScaleAbs(I_y, *(feature[6]), 0.25);

            cv::Mat img_Lab;
            cv::cvtColor(*img.at(0), img_Lab, CV_RGB2Lab);
            cv::vector<cv::Mat> tempfeature(3);

            cv::split(img_Lab, tempfeature);

            for(int i = 0; i < 3; ++i)
                tempfeature.at(i).copyTo(*(feature.at(i)));

            // min max filter
            for(int c = 0; c < 16; ++c)
                minFilter(feature[c], feature[c + 16], 5);
            for(int c = 0; c < 16; ++c)
                maxFilter(feature[c], feature[c], 5);

        }else{
            feature.clear();

            // calc gray integral image
            cv::Mat grayImg(imgRow + 1, imgCol, CV_8U);
            cv::cvtColor(*img.at(0), grayImg, CV_RGB2GRAY);
            integralMat = new cv::Mat(imgRow + 1, imgCol + 1, CV_64F);
            cv::integral(grayImg, *integralMat, CV_64F);
            feature.push_back(integralMat);

            // calc r g b integral image
            std::vector<cv::Mat> splittedRgb;
            cv::split(*img.at(0), splittedRgb);
            for(int i = 0; i < splittedRgb.size(); ++i){
                integralMat = new cv::Mat(imgRow + 1, imgCol + 1, CV_64F);
                cv::integral(splittedRgb.at(i), *integralMat, CV_64F);
                feature.push_back(integralMat);
            }



            featureFlag = 1;
        }
    }

    if(img.size() > 1){
        cv::Mat tempDepth = cv::Mat(img.at(0)->rows, img.at(0)->cols, CV_8U);// = *img.at(1);

        if(img.at(1)->type() != CV_8U)
            img.at(1)->convertTo(tempDepth, CV_8U, 255.0 / (double)(conf.maxdist - conf.mindist));
        else
            tempDepth = *img.at(1);
        integralMat = new cv::Mat(imgRow + 1, imgCol + 1, CV_64F);
        cv::integral(tempDepth, *integralMat, CV_64F);
        feature.push_back(integralMat);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        vector3d vW(v.x,v.y+height_map_resolution,0);
        vector3d vS(v.x-height_map_resolution,v.y,0);
        vector3d vE(v.x,v.y-height_map_resolution,0);
        
        auto itE=grid_points.find(vE);
        auto itN=grid_points.find(vN);
        auto itW=grid_points.find(vW);
        auto itS=grid_points.find(vS);
        
        if (itE == grid_points.end()||itN == grid_points.end()||itW == grid_points.end()||itS == grid_points.end()){
            exit(-1);
        }
        
        vE=itE->first;
        vN=itN->first;
        vW=itW->first;
        vS=itS->first;
        
        vector3d Normal = ((((vE-v)^(vN-v))+((vW-v)^(vS-v)))/2).norm();       //normal is the crossproduct of two prependicular differences. Avgd.
        
        
        std::vector<double> horizon_elevation_angles;   //vector to hold elevation angles at theta
        horizon_elevation_angles.assign(sun_angles, 0);
        //    std::vector<vector3d> horizon(sun_angles);
        std::vector<double> dists(sun_angles);
        
        size_t edge_points = 0;
        size_t interior_points = 0;
        
        for(int theta = 0; theta<sun_angles;theta++){        //compute for each angle theta
            double th = theta*M_PI/180;
            {
                double sin_theta = sin(th);   //calculate sin of the angle in radians
                double yend = EQ_DBL(copysign(1,sin_theta),1)?east_y_max:0;     //if the sinus is positive endvalue is east_y_max, if negative 0
                for(double y=v.y+copysign(height_map_resolution,sin_theta);(y >= 0 && y <= east_y_max);y+=copysign(height_map_resolution,sin_theta)){//increase/decrease if theta +/-
                    double x = v.x+(y-v.y)/tan(th);                             //find the x that goes with y for this theta
                    double x1 = floor(x-((int)x%(int)height_map_resolution));       //find the nearest lower gridpoint by subtracting remainder according to height_map_resolution
                    double x2 = x1 + height_map_resolution;                  //Add height_map_resolution to get nearest higher gridpoint
                    double phi = horizon_elevation_angles[theta];
                    double dist = (v-vector3d(x,y,v.z)).length();
                    double phiMaxTheta = atan(5000.0/dist);      //maximal phi at this theta (with height 5000 m)
                    
                    
                    if(x1>=north_x_max||x1<0||x2>=north_x_max||x2<0||phi>phiMaxTheta){
                        break;
                    }
                    if((int)x%(int)height_map_resolution){                       //if x is not a grid point
                        
                        auto it_vec1 = grid_points.find(vector3d(x1,y,0));      //get the two gridpoints from the set
                        auto it_vec2 = grid_points.find(vector3d(x2,y,0));
                        if (it_vec1 == grid_points.end()||it_vec2 == grid_points.end()){
                            edge_points++;
                            continue;
                        }
                        vector3d vec1 = it_vec1->first;
                        vector3d vec2 = it_vec2->first;
                        double height = vec1.z*(x2-x)/height_map_resolution + vec2.z*(x-x1)/height_map_resolution-v.z;       //compute height at x via linear interpolation
                        dist=(v-vector3d(x,y,v.z)).length();
                        phi = atan(height/dist);

                        if(phi>horizon_elevation_angles[theta]){//see if larger
                            horizon_elevation_angles[theta]=phi;
                            dists[theta] = dist/1000;
                        }
                        
                    }

double tgamma ( double x )
      {
      register int n, parity, i;
      register double y, y1, result, fact, IsItAnInt, z, numerator, 
	                denominator, ysquared, sum; 
      hexdouble OldEnvironment;
    
      FEGETENVD( OldEnvironment.d );               // save environment, set default
      FESETENVD( 0.0 );
	
/*******************************************************************************
*     The next switch will decipher what sort of argument we have. If argument *
*     is SNaN then a QNaN has to be returned and the invalid flag signaled.    * 
*******************************************************************************/

      switch ( __fpclassifyd ( x ) )
            {
            case FP_NAN:
                  x *= 2.0;                  /* quiets NaN */
                  FESETENVD( OldEnvironment.d ); //   restore caller's environment
                  return x;
                  
            case FP_ZERO:
                  OldEnvironment.i.lo |= FE_DIVBYZERO;
                  FESETENVD( OldEnvironment.d );
                  return copysign( Huge.d, x);

             case FP_INFINITE:
                  if ( x > 0.0 )
                        x = Huge.d;
                  else
                        {
                        x = nan ( GAMMA_NAN );
                        OldEnvironment.i.lo |= SET_INVALID;
                        }

                  FESETENVD( OldEnvironment.d );
                  return x;
                  
            default:                  /*      NORMALNUM and DENORMALNUM      */
                  break;
            }
      
      parity = 0;
      fact = 1.0;
      n = 0;
      y = x;

/*******************************************************************************
*     The argument is negative.                                                *
*******************************************************************************/

      if ( y <= 0.0 )
            {
            y = - x;
			if ( y < MinimumX )
			      {
                  OldEnvironment.i.lo |= FE_OVERFLOW;
                  FESETENVD( OldEnvironment.d );
                  return MinusHuge.d;
				  }
            y1 = trunc ( y );
            IsItAnInt = y - y1;
            if ( IsItAnInt != 0.0 )                   /* is it an integer?   */
                  {                                   /* is it odd or even?  */
                  if ( y1 != trunc ( y1 * 0.5 ) * 2.0 ) parity = 1;
                  fact = - pi / sin ( pi * IsItAnInt );
                  y += 1.0;
                  }
            else
                  {
                  OldEnvironment.i.lo |= SET_INVALID;
                  FESETENVD( OldEnvironment.d );
                  return nan ( GAMMA_NAN );
                  }
            }

/*******************************************************************************
*     The argument is positive.                                                *
*******************************************************************************/

     if ( y < eps )                         /* argument is less than epsilon. */
            {
            if ( y >= MinimumX )          /* x is in [MinimumX,eps].          */
                  result = 1.0 / y;
            else                          /* othewise, x is in [0,MinimumX).  */
                  {
                  OldEnvironment.i.lo |= FE_OVERFLOW;
                  FESETENVD( OldEnvironment.d );
                  return Huge.d;
                  }
            }
      else if ( y < 12.0 )                 /* argument x is eps < x < 12.0.  */
            {
            y1 = y;
            if ( y < 1.0 )                 /* x is in (eps, 1.0).            */
                  {
                  z = y;
                  y += 1.0;
                  }
//.........这里部分代码省略.........

double complex csinh(double complex z)
{
	double x, y, h;
	int32_t hx, hy, ix, iy, lx, ly;

	x = creal(z);
	y = cimag(z);

	EXTRACT_WORDS(hx, lx, x);
	EXTRACT_WORDS(hy, ly, y);

	ix = 0x7fffffff & hx;
	iy = 0x7fffffff & hy;

	/* Handle the nearly-non-exceptional cases where x and y are finite. */
	if (ix < 0x7ff00000 && iy < 0x7ff00000) {
		if ((iy | ly) == 0)
			return CMPLX(sinh(x), y);
		if (ix < 0x40360000)    /* small x: normal case */
			return CMPLX(sinh(x) * cos(y), cosh(x) * sin(y));

		/* |x| >= 22, so cosh(x) ~= exp(|x|) */
		if (ix < 0x40862e42) {
			/* x < 710: exp(|x|) won't overflow */
			h = exp(fabs(x)) * 0.5;
			return CMPLX(copysign(h, x) * cos(y), h * sin(y));
		} else if (ix < 0x4096bbaa) {
			/* x < 1455: scale to avoid overflow */
			z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
			return CMPLX(creal(z) * copysign(1, x), cimag(z));
		} else {
			/* x >= 1455: the result always overflows */
			h = huge * x;
			return CMPLX(h * cos(y), h * h * sin(y));
		}
	}

	/*
	 * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN.
	 * The sign of 0 in the result is unspecified.  Choice = normally
	 * the same as dNaN.  Raise the invalid floating-point exception.
	 *
	 * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN).
	 * The sign of 0 in the result is unspecified.  Choice = normally
	 * the same as d(NaN).
	 */
	if ((ix | lx) == 0 && iy >= 0x7ff00000)
		return CMPLX(copysign(0, x * (y - y)), y - y);

	/*
	 * sinh(+-Inf +- I 0) = +-Inf + I +-0.
	 *
	 * sinh(NaN +- I 0)   = d(NaN) + I +-0.
	 */
	if ((iy | ly) == 0 && ix >= 0x7ff00000) {
		if (((hx & 0xfffff) | lx) == 0)
			return CMPLX(x, y);
		return CMPLX(x, copysign(0, y));
	}

	/*
	 * sinh(x +- I Inf) = dNaN + I dNaN.
	 * Raise the invalid floating-point exception for finite nonzero x.
	 *
	 * sinh(x + I NaN) = d(NaN) + I d(NaN).
	 * Optionally raises the invalid floating-point exception for finite
	 * nonzero x.  Choice = don't raise (except for signaling NaNs).
	 */
	if (ix < 0x7ff00000 && iy >= 0x7ff00000)
		return CMPLX(y - y, x * (y - y));

	/*
	 * sinh(+-Inf + I NaN)  = +-Inf + I d(NaN).
	 * The sign of Inf in the result is unspecified.  Choice = normally
	 * the same as d(NaN).
	 *
	 * sinh(+-Inf +- I Inf) = +Inf + I dNaN.
	 * The sign of Inf in the result is unspecified.  Choice = always +.
	 * Raise the invalid floating-point exception.
	 *
	 * sinh(+-Inf + I y)   = +-Inf cos(y) + I Inf sin(y)
	 */
	if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
		if (iy >= 0x7ff00000)
			return CMPLX(x * x, x * (y - y));
		return CMPLX(x * cos(y), INFINITY * sin(y));
	}

	/*
	 * sinh(NaN + I NaN)  = d(NaN) + I d(NaN).
	 *
	 * sinh(NaN +- I Inf) = d(NaN) + I d(NaN).
	 * Optionally raises the invalid floating-point exception.
	 * Choice = raise.
	 *
	 * sinh(NaN + I y)    = d(NaN) + I d(NaN).
	 * Optionally raises the invalid floating-point exception for finite
	 * nonzero y.  Choice = don't raise (except for signaling NaNs).
	 */
	return CMPLX((x * x) * (y - y), (x + x) * (y - y));
//.........这里部分代码省略.........

void eqs_todo_bulk_ode<real_t>::condevap(
    const mtx::arr<real_t> &rhod,
    mtx::arr<real_t> &rhod_th,
    mtx::arr<real_t> &rhod_rv,
    mtx::arr<real_t> &rhod_rl,
    mtx::arr<real_t> &rhod_rr,
    const quantity<si::time, real_t> dt
)   
{
#  if !defined(USE_BOOST_ODEINT)
    error_macro("eqs_todo_bulk requires icicle to be compiled with Boost.odeint");
#  else
    // odeint::euler< // TODO: opcja?
    odeint::runge_kutta4<
      quantity<multiply_typeof_helper<si::mass_density, si::temperature>::type, real_t>, // state_type
      real_t, // value_type
      quantity<si::temperature, real_t>, // deriv_type
      quantity<si::mass_density, real_t>, // time_type
      odeint::vector_space_algebra, 
      odeint::default_operations, 
      odeint::never_resizer
    > S; // TODO: would be better to instantiate in the ctor (but what about thread safety! :()
    typename detail::rhs F;

    for (int k = rhod.lbound(mtx::k); k <= rhod.ubound(mtx::k); ++k)
      for (int j = rhod.lbound(mtx::j); j <= rhod.ubound(mtx::j); ++j)
        for (int i = rhod.lbound(mtx::i); i <= rhod.ubound(mtx::i); ++i)
        {

          F.init(
            rhod(i,j,k)    * si::kilograms / si::cubic_metres, 
            rhod_th(i,j,k) * si::kilograms / si::cubic_metres * si::kelvins, 
            rhod_rv(i,j,k) * si::kilograms / si::cubic_metres
          );
          real_t // TODO: quantity<si::mass_density
            rho_eps = .00002, // TODO: as an option?
            vapour_excess;
          real_t drho_rr_max = 0; // TODO: quantity<si::mass_density
          if (F.rs > F.r && rhod_rr(i,j,k) > 0 && opt_revp) 
            drho_rr_max = (dt / si::seconds) * (1 - F.r / F.rs) * (1.6 + 124.9 * pow(1e-3 * rhod_rr(i,j,k), .2046)) * pow(1e-3 * rhod_rr(i,j,k), .525) / 
              (5.4e2 + 2.55e5 * (1. / (F.p / si::pascals) / F.rs)); // TODO: move to phc!!!
          bool incloud;

          // TODO: rethink and document 2*rho_eps!!!
          while ( 
            // condensation of cloud water if supersaturated
            (vapour_excess = rhod_rv(i,j,k) - rhod(i,j,k) * F.rs) > rho_eps 
            || (opt_cevp && vapour_excess < -rho_eps && ( // or if subsaturated
              (incloud = (rhod_rl(i,j,k) > 0)) // cloud evaporation if in cloud
              || (opt_revp && rhod_rr(i,j,k) > 0) // or rain evaportation if in a rain shaft (and out-of-cloud)
            )) 
          ) 
          {
            real_t drho_rv = - copysign(.5 * rho_eps, vapour_excess);
            drho_rv = (vapour_excess > 0 || incloud)
              ?  std::min(rhod_rl(i,j,k), drho_rv)
              :  std::min(drho_rr_max, std::min(rhod_rr(i,j,k), drho_rv)); // preventing negative mixing ratios
            assert(drho_rv != 0); // otherwise it should not pass the while condition!

            // theta is modified by do_step, and hence we cannot pass an expression and we need a temp. var.
            quantity<multiply_typeof_helper<si::mass_density, si::temperature>::type, real_t> 
              tmp = rhod_th(i,j,k) * si::kilograms / si::cubic_metres * si::kelvins;

            // integrating the First Law for moist air
            S.do_step(
              boost::ref(F), 
              tmp,
              rhod_rv(i,j,k) * si::kilograms / si::cubic_metres, 
              drho_rv        * si::kilograms / si::cubic_metres
            );

            // latent heat source/sink due to evaporation/condensation
            rhod_th(i,j,k) = tmp / (si::kilograms / si::cubic_metres * si::kelvins); 

            // updating rhod_rv
            rhod_rv(i,j,k) += drho_rv;
            assert(rhod_rv(i,j,k) >= 0);
            assert(isfinite(rhod_rv(i,j,k)));
            
            if (vapour_excess > 0 || incloud) 
            {
              rhod_rl(i,j,k) -= drho_rv; // cloud water 
              assert(rhod_rl(i,j,k) >= 0);
              assert(isfinite(rhod_rl(i,j,k)));
            }
            else // or rain water
            {
              assert(opt_revp); // should be guaranteed by the while() condition above
              rhod_rr(i,j,k) -= drho_rv;
              assert(rhod_rr(i,j,k) >= 0);
              assert(isfinite(rhod_rr(i,j,k)));
              if ((drho_rr_max -= drho_rv) == 0) break; // but not more than Kessler allows
            }
          }
          // hopefully true for RK4
          assert(F.r == real_t(rhod_rv(i,j,k) / rhod(i,j,k)));
          // double-checking....
          assert(rhod_rl(i,j,k) >= 0);
          assert(rhod_rv(i,j,k) >= 0);
          assert(rhod_rr(i,j,k) >= 0);
//.........这里部分代码省略.........