TEST(Join, JoinLargeDim) {
    using af::constant;
    using af::deviceGC;
    using af::span;

    // const int nx = 32;
    const int nx = 1;
    const int ny = 4 * 1024 * 1024;
    const int nw = 4 * 1024 * 1024;

    deviceGC();
    {
        array in         = randu(nx, ny, u8);
        array joined     = join(0, in, in);
        dim4 in_dims     = in.dims();
        dim4 joined_dims = joined.dims();

        ASSERT_EQ(2 * in_dims[0], joined_dims[0]);
        ASSERT_EQ(0.f, sum<float>((joined(0, span) - joined(1, span)).as(f32)));

        array in2 = constant(1, (dim_t)nx, (dim_t)ny, (dim_t)2, (dim_t)nw, u8);
        joined    = join(3, in, in);
        in_dims   = in.dims();
        joined_dims = joined.dims();
        ASSERT_EQ(2 * in_dims[3], joined_dims[3]);
    }
}

TEST(MatrixManipulation, SNIPPET_matrix_manipulation_join) {
    //! [ex_matrix_manipulation_join]
    float hA[] = {1, 2, 3, 4, 5, 6};
    float hB[] = {10, 20, 30, 40, 50, 60, 70, 80, 90};
    array A    = array(3, 2, hA);
    array B    = array(3, 3, hB);

    af_print(join(1, A, B));  // 3x5 matrix
    // array result = join(0, A, B); // fail: dimension mismatch
    //! [ex_matrix_manipulation_join]

    array out = join(1, A, B);
    vector<float> h_out(out.elements());
    out.host(&h_out.front());
    af_print(out);

    ASSERT_EQ(3, out.dims(0));
    ASSERT_EQ(5, out.dims(1));

    unsigned fdim = out.dims(0);
    unsigned sdim = out.dims(1);
    for (unsigned i = 0; i < sdim; i++) {
        for (unsigned j = 0; j < fdim; j++) {
            if (i < 2) {
                ASSERT_FLOAT_EQ(hA[i * fdim + j], h_out[i * fdim + j])
                    << "At [" << i << ", " << j << "]";
            } else {
                ASSERT_FLOAT_EQ(hB[(i - 2) * fdim + j], h_out[i * fdim + j])
                    << "At [" << i << ", " << j << "]";
            }
        }
    }
}

void selectTest(const dim4 &dims)
{
    if (noDoubleTests<T>()) return;
    dtype ty = (dtype)dtype_traits<T>::af_type;

    array a = randu(dims, ty);
    array b = randu(dims, ty);

    if (a.isinteger()) {
        a = (a % (1 << 30)).as(ty);
        b = (b % (1 << 30)).as(ty);
    }

    array cond = randu(dims, ty) > a;

    array c = select(cond, a, b);

    int num = (int)a.elements();

    vector<T> ha(num);
    vector<T> hb(num);
    vector<T> hc(num);
    vector<char> hcond(num);

    a.host(&ha[0]);
    b.host(&hb[0]);
    c.host(&hc[0]);
    cond.host(&hcond[0]);

    for (int i = 0; i < num; i++) {
        ASSERT_EQ(hc[i], hcond[i] ? ha[i] : hb[i]);
    }
}

void selectScalarTest(const dim4 &dims)
{
    if (noDoubleTests<T>()) return;
    dtype ty = (dtype)dtype_traits<T>::af_type;

    array a = randu(dims, ty);
    array cond = randu(dims, ty) > a;
    double b = 3;

    if (a.isinteger()) {
        a = (a % (1 << 30)).as(ty);
    }

    array c = is_right ? select(cond, a, b) : select(cond, b, a);

    int num = (int)a.elements();

    vector<T> ha(num);
    vector<T> hc(num);
    vector<char> hcond(num);

    a.host(&ha[0]);
    c.host(&hc[0]);
    cond.host(&hcond[0]);

    if (is_right) {
        for (int i = 0; i < num; i++) {
            ASSERT_EQ(hc[i], hcond[i] ? ha[i] : T(b));
        }
    } else {
        for (int i = 0; i < num; i++) {
            ASSERT_EQ(hc[i], hcond[i] ? T(b) : ha[i]);
        }
    }
}

void sparseArithTesterMul(const int m, const int n, int factor, const double eps)
{
    deviceGC();

    if (noDoubleTests<T>()) return;

#if 1
    array A = cpu_randu<T>(dim4(m, n));
    array B = cpu_randu<T>(dim4(m, n));
#else
    array A = randu(m, n, (dtype)dtype_traits<T>::af_type);
    array B = randu(m, n, (dtype)dtype_traits<T>::af_type);
#endif

    A = makeSparse<T>(A, factor);

    array RA = sparse(A, AF_STORAGE_CSR);
    array OA = sparse(A, AF_STORAGE_COO);

    // Forward
    {
        // Arith Op
        array resR = arith_op<af_mul_t>()(RA, B);
        array resO = arith_op<af_mul_t>()(OA, B);

        // We will test this by converting the COO to CSR and CSR to COO and
        // comparing them. In essense, we are comparing the resR and resO
        // TODO: Make a better comparison using dense

        // Check resR against conR
        array conR = sparseConvertTo(resR, AF_STORAGE_CSR);
        sparseCompare<T>(resR, conR, eps);

        // Check resO against conO
        array conO = sparseConvertTo(resR, AF_STORAGE_COO);
        sparseCompare<T>(resO, conO, eps);
    }

    // Reverse
    {
        // Arith Op
        array resR = arith_op<af_mul_t>()(B, RA);
        array resO = arith_op<af_mul_t>()(B, OA);

        // We will test this by converting the COO to CSR and CSR to COO and
        // comparing them. In essense, we are comparing the resR and resO
        // TODO: Make a better comparison using dense

        // Check resR against conR
        array conR = sparseConvertTo(resR, AF_STORAGE_CSR);
        sparseCompare<T>(resR, conR, eps);

        // Check resO against conO
        array conO = sparseConvertTo(resR, AF_STORAGE_COO);
        sparseCompare<T>(resO, conO, eps);
    }
}

static void indexArray(af_array &dest, const af_array &src, const unsigned ndims, const af_seq *index)
{
    using af::toOffset;
    using af::toDims;
    using af::toStride;

    const Array<T> &parent = getArray<T>(src);
    vector<af_seq> index_(index, index+ndims);
    Array<T>* dst =  createSubArray(    parent,
                                        toDims(index_, parent.dims()),
                                        toOffset(index_, parent.dims()),
                                        toStride(index_, parent.dims()) );
    dest = getHandle(*dst);
}

void assign(Array<T>& out, const af_index_t idxrs[], const Array<T>& rhs)
{
    out.eval();
    rhs.eval();

    vector<bool> isSeq(4);
    vector<af_seq> seqs(4, af_span);
    // create seq vector to retrieve output dimensions, offsets & offsets
    for (dim_t x=0; x<4; ++x) {
        if (idxrs[x].isSeq) {
            seqs[x] = idxrs[x].idx.seq;
        }
        isSeq[x] = idxrs[x].isSeq;
    }

    vector< Array<uint> > idxArrs(4, createEmptyArray<uint>(dim4()));
    // look through indexs to read af_array indexs
    for (dim_t x=0; x<4; ++x) {
        if (!isSeq[x]) {
            idxArrs[x] = castArray<uint>(idxrs[x].idx.arr);
            idxArrs[x].eval();
        }
    }

    vector<CParam<uint>> idxParams(idxArrs.begin(), idxArrs.end());
    getQueue().enqueue(kernel::assign<T>, out, out.getDataDims(), rhs,
                       move(isSeq), move(seqs), move(idxParams));
}

void stdevDimIndexTest(string pFileName, dim_t dim=-1)
{
    typedef typename sdOutType<T>::type outType;
    if (noDoubleTests<T>()) return;
    if (noDoubleTests<outType>()) return;

    vector<dim4>      numDims;
    vector<vector<int> >       in;
    vector<vector<float> >  tests;

    readTestsFromFile<int,float>(pFileName, numDims, in, tests);

    dim4 dims = numDims[0];
    vector<T> input(in[0].begin(), in[0].end());

    array a(dims, &(input.front()));
    array b = a(seq(2,6), seq(1,7));

    array c = stdev(b, dim);

    vector<outType> currGoldBar(tests[0].begin(), tests[0].end());

    size_t nElems    = currGoldBar.size();
    vector<outType> outData(nElems);

    c.host((void*)outData.data());

    for (size_t elIter=0; elIter<nElems; ++elIter) {
        ASSERT_NEAR(::real(currGoldBar[elIter]), ::real(outData[elIter]), 1.0e-3)<< "at: " << elIter<< endl;
        ASSERT_NEAR(::imag(currGoldBar[elIter]), ::imag(outData[elIter]), 1.0e-3)<< "at: " << elIter<< endl;
    }
}

TEST(JIT, TransposeBuffers)
{
    const int num = 10;
    array a = randu(1, num);
    array b = randu(1, num);
    array c =  a + b;
    array d = a.T() + b.T();

    vector<float> ha(a.elements());
    a.host(ha.data());

    vector<float> hb(b.elements());
    b.host(hb.data());

    vector<float> hc(c.elements());
    c.host(hc.data());

    vector<float> hd(d.elements());
    d.host(hd.data());

    for (int i = 0; i < num; i++) {
        ASSERT_FLOAT_EQ(ha[i] + hb[i], hc[i]);
        ASSERT_FLOAT_EQ(hc[i], hd[i]);
    }
}

TEST(JIT, NonLinearLargeY)
{
    const int d0 = 2;
    // This needs to be > 2 * (1 << 20) to properly check this.
    const int d1 = 3 * (1 << 20);
    array a = randn(d0);
    array b = randn(1, d1);

    // tile is jit-ted for both the operations
    array c = tile(a, 1, d1) + tile(b, d0, 1);
    eval(c);

    vector<float> ha(d0);
    vector<float> hb(d1);
    vector<float> hc(d0 * d1);

    a.host(ha.data());
    b.host(hb.data());
    c.host(hc.data());

    for (int j = 0; j < d1; j++) {
        for (int i = 0; i < d0; i++) {
            ASSERT_EQ(hc[i + j * d0], ha[i] + hb[j]) << " at " << i << " , " << j;
        }
    }
}

void backendTest()
{
    int backends = getAvailableBackends();

    ASSERT_NE(backends, 0);

    bool cpu    = backends & AF_BACKEND_CPU;
    bool cuda   = backends & AF_BACKEND_CUDA;
    bool opencl = backends & AF_BACKEND_OPENCL;

    printf("\nRunning Default Backend...\n");
    testFunction<float>();

    if(cpu) {
        printf("\nRunning CPU Backend...\n");
        setBackend(AF_BACKEND_CPU);
        testFunction<float>();
    }

    if(cuda) {
        printf("\nRunning CUDA Backend...\n");
        setBackend(AF_BACKEND_CUDA);
        testFunction<float>();
    }

    if(opencl) {
        printf("\nRunning OpenCL Backend...\n");
        setBackend(AF_BACKEND_OPENCL);
        testFunction<float>();
    }
}

TEST(JIT, CPP_Multi_strided)
{
    const int num = 1024;
    gforSet(true);
    array a = randu(num, 1, s32);
    array b = randu(1, num, s32);
    array x = a + b;
    array y = a - b;
    eval(x, y);
    gforSet(false);

    vector<int> ha(num);
    vector<int> hb(num);
    vector<int> hx(num * num);
    vector<int> hy(num * num);

    a.host(&ha[0]);
    b.host(&hb[0]);
    x.host(&hx[0]);
    y.host(&hy[0]);

    for (int j = 0; j < num; j++) {
        for (int i = 0; i < num; i++) {
            ASSERT_EQ((ha[i] + hb[j]), hx[j*num + i]);
            ASSERT_EQ((ha[i] - hb[j]), hy[j*num + i]);
        }
    }
}

void sparseCompare(array A, array B, const double eps)
{
// This macro is used to check if either value is finite and then call assert
// If neither value is finite, then they can be assumed to be equal to either inf or nan
#define ASSERT_FINITE_EQ(V1, V2)                    \
    if(std::isfinite(V1) || std::isfinite(V2)) {    \
        ASSERT_NEAR(V1, V2, eps) << "at : " << i;   \
    }                                               \

    array AValues = sparseGetValues(A);
    array ARowIdx = sparseGetRowIdx(A);
    array AColIdx = sparseGetColIdx(A);

    array BValues = sparseGetValues(B);
    array BRowIdx = sparseGetRowIdx(B);
    array BColIdx = sparseGetColIdx(B);

    // Verify row and col indices
    ASSERT_EQ(0, max<int>(ARowIdx - BRowIdx));
    ASSERT_EQ(0, max<int>(AColIdx - BColIdx));

    T *ptrA = AValues.host<T>();
    T *ptrB = BValues.host<T>();
    for(int i = 0; i < AValues.elements(); i++) {
        ASSERT_FINITE_EQ(real(ptrA[i]), real(ptrB[i]));

        if(A.iscomplex()) {
            ASSERT_FINITE_EQ(imag(ptrA[i]), imag(ptrB[i]));
        }
    }
    freeHost(ptrA);
    freeHost(ptrB);

#undef ASSERT_FINITE_EQ
}

TEST(MatrixMultiply, RhsBroadcastBatched)
{
    const int M = 512;
    const int K = 512;
    const int N = 10;
    const int D2 = 2;
    const int D3 = 3;

    for (int d3 = 1; d3 <= D3; d3 *= D3) {
        for (int d2 = 1; d2 <= D2; d2 *= D2) {
            array a = randu(M, K, d2, d3);
            array b = randu(K, N);
            array c = matmul(a, b);

            for (int j = 0; j < d3; j++) {
                for (int i = 0; i < d2; i++) {
                    array a_ij = a(span, span, i, j);
                    array c_ij = c(span, span, i, j);
                    array res = matmul(a_ij, b);
                    EXPECT_LT(max<float>(abs(c_ij - res)), 1E-3)
                        << " for d2 = " << d2 << " for d3 = " << d3;
                }
            }
        }
    }
}

void cppMatMulCheck(string TestFile)
{
    if (noDoubleTests<T>()) return;

    vector<dim4> numDims;

    vector<vector<T> > hData;
    vector<vector<T> > tests;
    readTests<T,T,int>(TestFile, numDims, hData, tests);

    array a(numDims[0], &hData[0].front());
    array b(numDims[1], &hData[1].front());

    dim4 atdims = numDims[0];
    {
        dim_t f  =    atdims[0];
        atdims[0]   =    atdims[1];
        atdims[1]   =    f;
    }
    dim4 btdims = numDims[1];
    {
        dim_t f = btdims[0];
        btdims[0] = btdims[1];
        btdims[1] = f;
    }

    array aT = moddims(a, atdims.ndims(), atdims.get());
    array bT = moddims(b, btdims.ndims(), btdims.get());

    vector<array> out(tests.size());
    if(isBVector) {
        out[0] = matmul(aT, b,    AF_MAT_NONE,    AF_MAT_NONE);
        out[1] = matmul(bT, a,   AF_MAT_NONE,    AF_MAT_NONE);
        out[2] = matmul(b, a,    AF_MAT_TRANS,       AF_MAT_NONE);
        out[3] = matmul(bT, aT,   AF_MAT_NONE,    AF_MAT_TRANS);
        out[4] = matmul(b, aT,    AF_MAT_TRANS,       AF_MAT_TRANS);
    }
    else {
        out[0] = matmul(a, b, AF_MAT_NONE,   AF_MAT_NONE);
        out[1] = matmul(a, bT, AF_MAT_NONE,   AF_MAT_TRANS);
        out[2] = matmul(a, bT, AF_MAT_TRANS,      AF_MAT_NONE);
        out[3] = matmul(aT, bT, AF_MAT_TRANS,      AF_MAT_TRANS);
    }

    for(size_t i = 0; i < tests.size(); i++) {
        dim_t elems = out[i].elements();
        vector<T> h_out(elems);
        out[i].host((void*)&h_out.front());

        if (false == equal(h_out.begin(), h_out.end(), tests[i].begin())) {

            cout << "Failed test " << i << "\nCalculated: " << endl;
            copy(h_out.begin(), h_out.end(), ostream_iterator<T>(cout, ", "));
            cout << "Expected: " << endl;
            copy(tests[i].begin(), tests[i].end(), ostream_iterator<T>(cout, ", "));
            FAIL();
        }
    }
}

///////////////////////////////////// CPP ////////////////////////////////
//
TEST(Transform, CPP)
{
    if (noImageIOTests()) return;

    vector<dim4>   inDims;
    vector<string> inFiles;
    vector<dim_t>  goldDim;
    vector<string> goldFiles;

    vector<dim4> HDims;
    vector<vector<float> >   HIn;
    vector<vector<float> >   HTests;
    readTests<float, float, float>(TEST_DIR"/transform/tux_tmat.test",HDims,HIn,HTests);

    readImageTests(string(TEST_DIR"/transform/tux_nearest.test"), inDims, inFiles, goldDim, goldFiles);

    inFiles[0].insert(0,string(TEST_DIR"/transform/"));
    inFiles[1].insert(0,string(TEST_DIR"/transform/"));

    goldFiles[0].insert(0,string(TEST_DIR"/transform/"));

    array H = array(HDims[0][0], HDims[0][1], &(HIn[0].front()));
    array IH = array(HDims[0][0], HDims[0][1], &(HIn[0].front()));

    array scene_img = loadImage(inFiles[1].c_str(), false);

    array gold_img = loadImage(goldFiles[0].c_str(), false);

    array out_img = transform(scene_img, IH, inDims[0][0], inDims[0][1], AF_INTERP_NEAREST, false);

    dim4 outDims = out_img.dims();
    dim4 goldDims = gold_img.dims();

    vector<float> h_out_img(outDims[0] * outDims[1]);
    out_img.host(&h_out_img.front());
    vector<float> h_gold_img(goldDims[0] * goldDims[1]);
    gold_img.host(&h_gold_img.front());

    const dim_t n = gold_img.elements();
    const float thr = 1.0f;

    // Maximum number of wrong pixels must be <= 0.01% of number of elements,
    // this metric is necessary due to rounding errors between different
    // backends for AF_INTERP_NEAREST and AF_INTERP_LOWER
    const size_t maxErr = n * 0.0001f;
    size_t err = 0;

    for (dim_t elIter = 0; elIter < n; elIter++) {
        err += fabs((int)h_out_img[elIter] - h_gold_img[elIter]) > thr;
        if (err > maxErr) {
            ASSERT_LE(err, maxErr) << "at: " << elIter << endl;
        }
    }
}

TEST(Where, MaxDim) {
    const size_t largeDim = 65535 * 32 + 2;

    array input  = range(dim4(1, largeDim), 1);
    array output = where(input % 2 == 0);
    array gold   = 2 * range(largeDim / 2);
    ASSERT_ARRAYS_EQ(gold.as(u32), output);

    input  = range(dim4(1, 1, 1, largeDim), 3);
    output = where(input % 2 == 0);
    ASSERT_ARRAYS_EQ(gold.as(u32), output);
}

TEST(JIT, CPP_common_node)
{
    array r = seq(-3, 3, 0.5);

    int n = r.dims(0);

    array x = tile(r, 1, r.dims(0));
    array y = tile(r.T(), r.dims(0), 1);
    x.eval();
    y.eval();


    vector<float> hx(x.elements());
    vector<float> hy(y.elements());
    vector<float> hr(r.elements());

    x.host(&hx[0]);
    y.host(&hy[0]);
    r.host(&hr[0]);

    for (int j = 0; j < n; j++) {
        for (int i = 0; i < n; i++) {
            ASSERT_EQ(hx[j * n + i], hr[i]);
            ASSERT_EQ(hy[j * n + i], hr[j]);
        }
    }
}

TEST(Susan, InvalidEdge) {
    try {
        array a      = randu(128, 128);
        features out = susan(a, 3, 32, 10, 1.3f, 129);
        EXPECT_TRUE(false);
    } catch (exception &e) { EXPECT_TRUE(true); }
}