     if(fp->layers[k]->activation == RELU)
     {
         expr_replace_relu_bounds_conv_sparse<<<dim3(fp->layers[layerno]->output_size[0], fp->layers[layerno]->output_size[1], fp->layers[layerno]->output_size[2]/num_chunks), fp->layers[k]->output_size[2]>>>(*linf_coeff, *lsup_coeff, *uinf_coeff, *usup_coeff, *linf_cst, *lsup_cst, *uinf_cst, *usup_cst, aux_lb_array, aux_ub_array, fp->layers[k]->output_size[0], fp->layers[k]->output_size[1], fp->layers[k]->output_size[2], offset_x, offset_y, length_x, length_y, shift_x, shift_y, use_area_heuristic);
-        //lexpr_replace_relu_bounds_conv_sparse<<<dim3(fp->layers[layerno]->output_size[0], fp->layers[layerno]->output_size[1], fp->layers[layerno]->output_size[2]/num_chunks), fp->layers[k]->output_size[2]>>>(*linf_coeff, *lsup_coeff, *linf_cst, *lsup_cst, aux_lb_array, aux_ub_array, fp->layers[k]->output_size[0], fp->layers[k]->output_size[1], fp->layers[k]->output_size[2], offset_x, offset_y, length_x, length_y, shift_x, shift_y, use_area_heuristic);
-        //uexpr_replace_relu_bounds_conv_sparse<<<dim3(fp->layers[layerno]->output_size[0], fp->layers[layerno]->output_size[1], fp->layers[layerno]->output_size[2]/num_chunks), fp->layers[k]->output_size[2]>>>(*uinf_coeff, *usup_coeff, *uinf_cst, *usup_cst, aux_lb_array, aux_ub_array, fp->layers[k]->output_size[0], fp->layers[k]->output_size[1], fp->layers[k]->output_size[2], offset_x, offset_y, length_x, length_y, shift_x, shift_y, use_area_heuristic);
     }
 
     size_t missing_length = ((length_x - 1)*fp->layers[k]->strides[0] + fp->layers[k]->filter_size[0])*((length_y - 1)*fp->layers[k]->strides[1] + fp->layers[k]->filter_size[1])*fp->layers[k]->input_size[2];



__global__
void lexpr_replace_relu_bounds(float_type* __restrict__ inf_coeff, float_type* __restrict__ sup_coeff, float_type* __restrict__ inf_cst, float_type* __restrict__ sup_cst, const float_type* __restrict__ lb_array, const float_type* __restrict__ ub_array, const size_t num_out_neurons_current_layer, const bool use_area_heuristic)
{
    const size_t n = blockIdx.x;

    size_t i = threadIdx.x;

    float_type res_inf_cst = 0;
    float_type res_sup_cst = 0;

    while(i < num_out_neurons_current_layer)
    {
        const size_t a = n*num_out_neurons_current_layer + i;

        const float_type lb = lb_array[i];
        const float_type ub = ub_array[i];
        const float_type width = ub - lb;
        const float_type lambda_inf = ub/width;
        const float_type lambda_sup = ub/width;

        const float_type old_inf_coeff = inf_coeff[a];
        const float_type old_sup_coeff = sup_coeff[a];

        if((old_sup_coeff == 0) && (old_inf_coeff == 0))
        {
            inf_coeff[a] = 0.0;
            sup_coeff[a] = 0.0;
        }
        else if(ub <= 0)
        {
            inf_coeff[a] = 0.0;
            sup_coeff[a] = 0.0;
        }
        else if(lb > 0)
        {
            inf_coeff[a] = old_inf_coeff;
            sup_coeff[a] = old_sup_coeff;
        }
        else if(old_sup_coeff < 0)
        {
            const float_type mu_inf = -lambda_inf*lb;
            const float_type mu_sup = -lambda_sup*lb;
            elina_double_interval_mul_expr_coeff(&inf_coeff[a], &sup_coeff[a], lambda_inf, lambda_sup, old_inf_coeff, old_sup_coeff);
            float_type tmp1, tmp2;
            elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, mu_inf, mu_sup, old_inf_coeff, old_sup_coeff);

            res_inf_cst += tmp1 - min_denormal;
            res_sup_cst += tmp2 + min_denormal;
        }
        else if (old_inf_coeff > 0)
        {
            if(use_area_heuristic)
            {
                const float_type area1 = -lb*ub;
                const float_type area2 = 0.5*ub*width;
                const float_type area3 = -0.5*lb*width;

                if((area1 < area2) && (area1 < area3))
                {
                    elina_double_interval_mul_expr_coeff(&inf_coeff[a], &sup_coeff[a], lambda_inf, lambda_sup, old_inf_coeff, old_sup_coeff);
                }
                else if((area2 < area1) && (area2 < area3))
                {
                    inf_coeff[a] = 0.0;
                    sup_coeff[a] = 0.0;
                }
                else
                {
                    inf_coeff[a] = old_inf_coeff;
                    sup_coeff[a] = old_sup_coeff;
                }
            }
            else
            {
                inf_coeff[a] = 0.0;
                sup_coeff[a] = 0.0;
            }
        }
        else
        {
            inf_coeff[a] = 0.0;
            sup_coeff[a] = 0.0;
            float_type tmp1, tmp2;
            elina_double_interval_mul2(&tmp1, &tmp2, old_inf_coeff, old_sup_coeff, 0, ub);

            res_inf_cst += tmp1;
            res_sup_cst += tmp1;
        }

        i += blockDim.x;
    }

    block_reduce_sum(res_inf_cst, blockDim.x);
    block_reduce_sum(res_sup_cst, blockDim.x);

    if(threadIdx.x == 0)
    {
        inf_cst[n] = res_inf_cst + inf_cst[n];
        sup_cst[n] = res_sup_cst + sup_cst[n];
    }
}


__global__
void uexpr_replace_relu_bounds(float_type* __restrict__ inf_coeff, float_type* __restrict__ sup_coeff, float_type* __restrict__ inf_cst, float_type* __restrict__ sup_cst, const float_type* __restrict__ lb_array, const float_type* __restrict__ ub_array, const size_t num_out_neurons_current_layer, const bool use_area_heuristic)
{
    const size_t n = blockIdx.x;

    size_t i = threadIdx.x;

    float_type res_inf_cst = 0;
    float_type res_sup_cst = 0;

    while(i < num_out_neurons_current_layer)
    {
        const size_t a = n*num_out_neurons_current_layer + i;

        const float_type lb = lb_array[i];
        const float_type ub = ub_array[i];
        const float_type width = ub - lb;
        const float_type lambda_inf = ub/width;
        const float_type lambda_sup = ub/width;

        const float_type old_inf_coeff = inf_coeff[a];
        const float_type old_sup_coeff = sup_coeff[a];

        if((old_sup_coeff == 0) && (old_inf_coeff == 0))
        {
            inf_coeff[a] = 0.0;
            sup_coeff[a] = 0.0;
        }
        else if(ub <= 0)
        {
            inf_coeff[a] = 0.0;
            sup_coeff[a] = 0.0;
        }
        else if(lb > 0)
        {
            inf_coeff[a] = old_inf_coeff;
            sup_coeff[a] = old_sup_coeff;
        }
        else if(old_inf_coeff > 0)
        {
            const float_type mu_inf = -lambda_inf*lb;
            const float_type mu_sup = -lambda_sup*lb;
            elina_double_interval_mul_expr_coeff(&inf_coeff[a], &sup_coeff[a], lambda_inf, lambda_sup, old_inf_coeff, old_sup_coeff);
            float_type tmp1, tmp2;
            elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, mu_inf, mu_sup, old_inf_coeff, old_sup_coeff);

            res_inf_cst += tmp1 - min_denormal;
            res_sup_cst += tmp2 + min_denormal;
        }
        else if(old_sup_coeff < 0)
        {
            if(use_area_heuristic)
            {
                const float_type area1 = -lb*ub;
                const float_type area2 = 0.5*ub*width;
                const float_type area3 = -0.5*lb*width;

                if((area1 < area2) && (area1 < area3))
                {
                    elina_double_interval_mul_expr_coeff(&inf_coeff[a], &sup_coeff[a], lambda_inf, lambda_sup, old_inf_coeff, old_sup_coeff);
                }
                else if((area2 < area1) && (area2 < area3))
                {
                    inf_coeff[a] = 0.0;
                    sup_coeff[a] = 0.0;
                }
                else
                {
                    inf_coeff[a] = old_inf_coeff;
                    sup_coeff[a] = old_sup_coeff;
                }
            }
            else
            {
                inf_coeff[a] = 0.0;
                sup_coeff[a] = 0.0;
            }
        }
        else
        {
            inf_coeff[a] = 0.0;
            sup_coeff[a] = 0.0;
            float_type tmp1, tmp2;
            elina_double_interval_mul2(&tmp1, &tmp2, old_inf_coeff, old_sup_coeff, 0, ub);

            res_inf_cst += tmp2;
            res_sup_cst += tmp2;
        }

        i += blockDim.x;
    }

    block_reduce_sum(res_inf_cst, blockDim.x);
    block_reduce_sum(res_sup_cst, blockDim.x);

    if(threadIdx.x == 0)
    {
        inf_cst[n] = res_inf_cst + inf_cst[n];
        sup_cst[n] = res_sup_cst + sup_cst[n];
    }
}


__global__
void expr_replace_relu_bounds(float_type* __restrict__ linf_coeff, float_type* __restrict__ lsup_coeff, float_type* __restrict__ uinf_coeff, float_type* __restrict__ usup_coeff, float_type* __restrict__ linf_cst, float_type* __restrict__ lsup_cst, float_type* __restrict__ uinf_cst, float_type* __restrict__ usup_cst, const float_type* __restrict__ lb_array, const float_type* __restrict__ ub_array, const size_t num_out_neurons_current_layer, const bool use_area_heuristic)
{
    const size_t n = blockIdx.x;

        res_lsup_cst[n] = lsup_cst + expr_lsup_cst[n];
        res_uinf_cst[n] = uinf_cst + expr_uinf_cst[n];
        res_usup_cst[n] = usup_cst + expr_usup_cst[n];
    }
}

