 
         std::cout << "shift_x " << shift_x << " shift_y " << shift_y << std::endl;
 
-        const long int mat_min_min_0 =  x_min_0     *fp->layers[layerno]->input_size[1]*fp->layers[layerno]->input_size[2] + y_min_0*fp->layers[layerno]->input_size[2];
-        const long int mat_min_max_0 =  x_min_0     *fp->layers[layerno]->input_size[1]*fp->layers[layerno]->input_size[2] + y_max_0*fp->layers[layerno]->input_size[2];
-        const long int mat_max_min_0 = (x_max_0 - 1)*fp->layers[layerno]->input_size[1]*fp->layers[layerno]->input_size[2] + y_min_0*fp->layers[layerno]->input_size[2];
-        const long int mat_max_max_0 = (x_max_0 - 1)*fp->layers[layerno]->input_size[1]*fp->layers[layerno]->input_size[2] + y_max_0*fp->layers[layerno]->input_size[2];
+        const long int mat_min_min_0 =  x_min_0     *fp->layers[k]->input_size[1]*fp->layers[k]->input_size[2] + y_min_0*fp->layers[k]->input_size[2];
+        const long int mat_min_max_0 =  x_min_0     *fp->layers[k]->input_size[1]*fp->layers[k]->input_size[2] + y_max_0*fp->layers[k]->input_size[2];
+        const long int mat_max_min_0 = (x_max_0 - 1)*fp->layers[k]->input_size[1]*fp->layers[k]->input_size[2] + y_min_0*fp->layers[k]->input_size[2];
+        const long int mat_max_max_0 = (x_max_0 - 1)*fp->layers[k]->input_size[1]*fp->layers[k]->input_size[2] + y_max_0*fp->layers[k]->input_size[2];
 
         std::cout << "(x_min_0, y_min_0) " << mat_min_min_0 << std::endl;
         std::cout << "(x_min_0, y_max_0) " << mat_min_max_0 << std::endl;
         std::cout << "(x_max_0, y_min_0) " << mat_max_min_0 << std::endl;
         std::cout << "(x_max_0, y_max_0) " << mat_max_max_0 << std::endl;
 
-        const long int mat_min_min_1 =  x_min_1     *fp->layers[layerno]->input_size[1]*fp->layers[layerno]->input_size[2] + y_min_1*fp->layers[layerno]->input_size[2];
-        const long int mat_min_max_1 =  x_min_1     *fp->layers[layerno]->input_size[1]*fp->layers[layerno]->input_size[2] + y_max_1*fp->layers[layerno]->input_size[2];
-        const long int mat_max_min_1 = (x_max_1 - 1)*fp->layers[layerno]->input_size[1]*fp->layers[layerno]->input_size[2] + y_min_1*fp->layers[layerno]->input_size[2];
-        const long int mat_max_max_1 = (x_max_1 - 1)*fp->layers[layerno]->input_size[1]*fp->layers[layerno]->input_size[2] + y_max_1*fp->layers[layerno]->input_size[2];
+        const long int mat_min_min_1 =  x_min_1     *fp->layers[k]->input_size[1]*fp->layers[k]->input_size[2] + y_min_1*fp->layers[k]->input_size[2];
+        const long int mat_min_max_1 =  x_min_1     *fp->layers[k]->input_size[1]*fp->layers[k]->input_size[2] + y_max_1*fp->layers[k]->input_size[2];
+        const long int mat_max_min_1 = (x_max_1 - 1)*fp->layers[k]->input_size[1]*fp->layers[k]->input_size[2] + y_min_1*fp->layers[k]->input_size[2];
+        const long int mat_max_max_1 = (x_max_1 - 1)*fp->layers[k]->input_size[1]*fp->layers[k]->input_size[2] + y_max_1*fp->layers[k]->input_size[2];
 
         std::cout << "(x_min_1, y_min_1) " << mat_min_min_1 << std::endl;
         std::cout << "(x_min_1, y_max_1) " << mat_min_max_1 << std::endl;


    std::cout << "num_out_neurons_last " << num_out_neurons_last_layer << std::endl;

    float_type* coeffs;
    float_type* csts;

    if(fp->layers[layerno]->type == CONV)
    {
        std::cout << "INITIAL" << std::endl;

        std::cout << "Parameters: " << std::endl
            << "x_pad = " << fp->layers[layerno]->pad[0] << " x_filter_size = " << fp->layers[layerno]->filter_size[0] << " x_stride = " << fp->layers[layerno]->strides[0] << std::endl;

        const size_t out_x_0 = 0;
        const size_t out_y_0 = 0;

        const long int x_min_0 = fp->layers[layerno]->strides[0]*out_x_0 - fp->layers[layerno]->pad[0];
        const long int x_max_0 = fp->layers[layerno]->strides[0]*out_x_0 + fp->layers[layerno]->filter_size[0] - fp->layers[layerno]->pad[0];

        std::cout << "xmin_0 " << x_min_0 << " xmax_0 " << x_max_0 << std::endl;

        const long int y_min_0 = fp->layers[layerno]->strides[1]*out_y_0 - fp->layers[layerno]->pad[1];
        const long int y_max_0 = fp->layers[layerno]->strides[1]*out_y_0 + fp->layers[layerno]->filter_size[1] - fp->layers[layerno]->pad[1];

        std::cout << "ymin_0 " << y_min_0 << " ymax_0 " << y_max_0 << std::endl;

        const long int out_x_1 = 1;
        const long int out_y_1 = 1;

        const long int x_min_1 = fp->layers[layerno]->strides[0]*out_x_1 - fp->layers[layerno]->pad[0];
        const long int x_max_1 = fp->layers[layerno]->strides[0]*out_x_1 + fp->layers[layerno]->filter_size[0] - fp->layers[layerno]->pad[0];

        std::cout << "xmin_1 " << x_min_1 << " xmax_1 " << x_max_1 << std::endl;

        const long int y_min_1 = fp->layers[layerno]->strides[1]*out_y_1 - fp->layers[layerno]->pad[1];
        const long int y_max_1 = fp->layers[layerno]->strides[1]*out_y_1 + fp->layers[layerno]->filter_size[1] - fp->layers[layerno]->pad[1];

        std::cout << "ymin_1 " << y_min_1 << " ymax_1 " << y_max_1 << std::endl;



        std::cout << std::endl;

        std::swap(linf_coeff, linf_coeff_tmp);
        std::swap(lsup_coeff, lsup_coeff_tmp);
        std::swap(linf_cst, linf_cst_tmp);
        std::swap(lsup_cst, lsup_cst_tmp);

        std::swap(uinf_coeff, uinf_coeff_tmp);
        std::swap(usup_coeff, usup_coeff_tmp);
        std::swap(uinf_cst, uinf_cst_tmp);
        std::swap(usup_cst, usup_cst_tmp);

        cudaFree(linf_coeff_tmp);
        cudaFree(lsup_coeff_tmp);
        cudaFree(uinf_coeff_tmp);
        cudaFree(usup_coeff_tmp);
    }

    compute_lb_from_expr_conv_sparse<<<dim3(fp->layers[layerno]->output_size[0], fp->layers[layerno]->output_size[1], fp->layers[layerno]->output_size[2]), 1>>>(lb_array, linf_coeff, lsup_coeff, linf_cst, fp->input_inf, fp->input_sup, num_out_neurons_last_layer, fp->layers[0]->input_size[0], fp->layers[0]->input_size[1], fp->layers[0]->input_size[2], offset_x, offset_y, length_x, length_y, shift_x, shift_y);
    compute_ub_from_expr_conv_sparse<<<dim3(fp->layers[layerno]->output_size[0], fp->layers[layerno]->output_size[1], fp->layers[layerno]->output_size[2]), 1>>>(ub_array, uinf_coeff, usup_coeff, usup_cst, fp->input_inf, fp->input_sup, num_out_neurons_last_layer, fp->layers[0]->input_size[0], fp->layers[0]->input_size[1], fp->layers[0]->input_size[2], offset_x, offset_y, length_x, length_y, shift_x, shift_y);

    cudaFree(linf_coeff);
    cudaFree(lsup_coeff);
    cudaFree(linf_cst);
    cudaFree(lsup_cst);

    cudaFree(uinf_coeff);
    cudaFree(usup_coeff);
    cudaFree(uinf_cst);
    cudaFree(usup_cst);

    cudaFree(linf_cst_tmp);
    cudaFree(lsup_cst_tmp);

    cudaFree(uinf_cst_tmp);
    cudaFree(usup_cst_tmp);

    cudaDeviceSynchronize();

    auto end = std::chrono::system_clock::now();

    std::chrono::duration<double> elapsed_seconds = end-start;
    std::cout << "elapsed time: " << elapsed_seconds.count() << "s" << std::endl << std::endl;
}


void ffn_handle_intermediate_layer(elina_manager_t* man, elina_abstract0_t* element, const double** weights, const double* bias, const size_t num_out_neurons, const size_t num_in_neurons, const activation_type_t activation)
{
    fppoly_t* fp = fppoly_of_abstract0(element);
    ffn_add_layer(fp, num_out_neurons, num_in_neurons, FFN, activation);

    float_type* coeffs = fp->layers[fp->numlayers - 1]->coeffs;
    float_type* csts = fp->layers[fp->numlayers - 1]->csts;

    layer_create_dense_exprs(coeffs, csts, weights, bias, num_out_neurons, num_in_neurons);

    update_state_using_previous_layers(man, fp, fp->numlayers - 1);
}


void ffn_handle_intermediate_relu_layer(elina_manager_t* man, elina_abstract0_t* element, const double** weights, const double* bias, const size_t num_out_neurons, const size_t num_in_neurons)
{
    ffn_handle_intermediate_layer(man, element, weights, bias, num_out_neurons, num_in_neurons, RELU);
}

void ffn_handle_intermediate_sigmoid_layer(elina_manager_t* man, elina_abstract0_t* element, const double** weights, const double* bias, const size_t num_out_neurons, const size_t num_in_neurons)
{
    //ffn_handle_intermediate_layer(man, element, weights, bias, num_out_neurons, num_in_neurons, SIGMOID);
}

void ffn_handle_intermediate_tanh_layer(elina_manager_t* man, elina_abstract0_t* element, const double** weights, const double* bias, const size_t num_out_neurons, const size_t num_in_neurons)
{
    //ffn_handle_intermediate_layer(man, element, weights, bias, num_out_neurons, num_in_neurons, TANH);
}


__global__
void print_bounds(const float_type* __restrict__ bounds_array, const size_t num_out_neurons)
{
    for(size_t i = 0; i < num_out_neurons; i++)
    {
        printf("out inf number %i is: %g\n", i, bounds_array[i]);
    }
}


void ffn_handle_last_layer(elina_manager_t* man, elina_abstract0_t* element, const double** weights, const double* bias, const size_t num_out_neurons, const size_t num_in_neurons, const bool has_activation, const activation_type_t activation)
{
    fppoly_t* fp = fppoly_of_abstract0(element);
    fppoly_internal_t* pr = fppoly_init_from_manager(man, ELINA_FUNID_ASSIGN_LINEXPR_ARRAY);

    if(has_activation)
    {
        ffn_add_layer(fp, num_out_neurons, num_in_neurons, FFN, activation);
    }
    else
    {
        ffn_add_layer(fp, num_out_neurons, num_in_neurons, FFN, NONE);
    }

    float_type* coeffs = fp->layers[fp->numlayers - 1]->coeffs;
    float_type* csts = fp->layers[fp->numlayers - 1]->csts;

    layer_create_dense_exprs(coeffs, csts, weights, bias, num_out_neurons, num_in_neurons);

    update_state_using_previous_layers(man, fp, fp->numlayers - 1);

    float_type* lb_array = fp->layers[fp->numlayers - 1]->lb_array;
    float_type* ub_array = fp->layers[fp->numlayers - 1]->ub_array;

    print_bounds<<<1, 1>>>(lb_array, num_out_neurons);
    print_bounds<<<1, 1>>>(ub_array, num_out_neurons);

    std::cout << std::endl;
}

void ffn_handle_last_relu_layer(elina_manager_t* man, elina_abstract0_t* element, const double** weights, const double* bias, const size_t num_out_neurons, const size_t num_in_neurons, const bool has_relu)
{
    ffn_handle_last_layer(man, element, weights, bias, num_out_neurons, num_in_neurons, has_relu, RELU);
}

void ffn_handle_last_sigmoid_layer(elina_manager_t* man, elina_abstract0_t* element, const double** weights, const double* bias, const size_t num_out_neurons, const size_t num_in_neurons, const bool has_sigmoid)
{
    //ffn_handle_last_layer(man, element, weights, bias, num_out_neurons, num_in_neurons, has_sigmoid, SIGMOID);
}

void ffn_handle_last_tanh_layer(elina_manager_t* man, elina_abstract0_t* element, const double** weights, const double* bias, const size_t num_out_neurons, const size_t num_in_neurons, const bool has_tanh)
{
    //ffn_handle_last_layer(man, element, weights, bias, num_out_neurons, num_in_neurons, has_tanh, TANH);
}


__global__
void create_sub_expr(float_type* __restrict__ inf_coeff, float_type* __restrict__ sup_coeff, float_type* __restrict__ inf_cst, float_type* __restrict__ sup_cst, const size_t index, const elina_dim_t y, const elina_dim_t x)
{
    inf_cst[index] = 0;
    sup_cst[index] = 0;

    for(size_t i = 0; i < 10; i++)
    {
        inf_coeff[index*10 + i] = 0.;
        sup_coeff[index*10 + i] = 0.;
    }

    inf_coeff[index*10 + y] = -1.;
    sup_coeff[index*10 + y] = 1.;

    inf_coeff[index*10 + x] = 1.;
    sup_coeff[index*10 + x] = -1.;
}


void get_lb_using_previous_layers(elina_manager_t* man, const fppoly_t* const fp)
{
    const size_t numlayers = fp->numlayers;
    fppoly_internal_t* pr = fppoly_init_from_manager(man, ELINA_FUNID_ASSIGN_LINEXPR_ARRAY);

    const size_t num_out_neurons_last_layer = 90;

    const size_t num_in_neurons_first_layer = fp->layers[0]->num_in_neurons;

    float_type* lb_dev;
    cudaMalloc((void**) &lb_dev, num_out_neurons_last_layer*sizeof(float_type));

    float_type* linf_coeff;
    float_type* lsup_coeff;
    float_type* linf_cst;
    float_type* lsup_cst;

    cudaMalloc((void**) &linf_coeff, num_out_neurons_last_layer*10*sizeof(float_type*));
    cudaMalloc((void**) &lsup_coeff, num_out_neurons_last_layer*10*sizeof(float_type*));
    cudaMalloc((void**) &linf_cst, num_out_neurons_last_layer*10*sizeof(float_type));
    cudaMalloc((void**) &lsup_cst, num_out_neurons_last_layer*10*sizeof(float_type));

    size_t index = 0;

    for(elina_dim_t y = 0; y < 10; y++)
    {
        for(elina_dim_t x = 0; x < 10; x++)
        {
            if(y != x)
            {
                create_sub_expr<<<1, 1>>>(linf_coeff, lsup_coeff, linf_cst, lsup_cst, index, y, x);
                index++;
            }
        }
    }

    float_type* linf_coeff_tmp;
