 }
 
 
-/*
-__global__
-void device_layer_create_sparse_exprs(expr_t** expr_array, const double* weights, const double* bias, const size_t* dims, const size_t num_out_neurons, const size_t num_in_neurons)
-{
-    size_t i = blockIdx.x;
-
-    if(i < num_out_neurons)
-    {
-        const double* weight_i = weights + i*num_in_neurons;
-        const size_t* dims_i = dims + i*num_in_neurons;
-        const double bias_i = bias[i];
-        expr_array[i] = create_sparse_expr(weight_i, bias_i, dims_i, num_in_neurons);
-        sort_sparse_expr(expr_array[i]);
-    }
-}
-
-
-void layer_create_sparse_exprs(expr_t** expr_array, const double* weights, const double* bias, const size_t* dims, const size_t num_out_neurons, const size_t num_in_neurons)
-{
-    double* tmp_weights;
-    cudaMalloc((void**) &tmp_weights, num_out_neurons*num_in_neurons*sizeof(double));
-    size_t* tmp_dims;
-    cudaMalloc((void**) &tmp_dims   , num_out_neurons*num_in_neurons*sizeof(size_t));
-
-    cudaMemcpy(tmp_weights, weights, num_out_neurons*num_in_neurons*sizeof(double), cudaMemcpyHostToDevice);
-    cudaMemcpy(tmp_dims   , dims   , num_out_neurons*num_in_neurons*sizeof(size_t), cudaMemcpyHostToDevice);
-
-    double* tmp_bias;
-    cudaMalloc((void**) &tmp_bias, num_out_neurons*sizeof(double));
-    cudaMemcpy(tmp_bias, bias, num_out_neurons*sizeof(double), cudaMemcpyHostToDevice);
-
-    device_layer_create_sparse_exprs<<<num_out_neurons, 1>>>(expr_array, tmp_weights, tmp_bias, tmp_dims, num_out_neurons, num_in_neurons);
-
-    cudaFree(tmp_weights);
-    cudaFree(tmp_dims);
-    cudaFree(tmp_bias);
-}
-
-
-elina_abstract0_t* fppoly_from_network_input_poly(elina_manager_t* man, const size_t intdim, const size_t realdim, const double* inf_array,
-                                                  const double* sup_array, const double* lexpr_weights, const double* lexpr_cst, const size_t* lexpr_dim,
-                                                  const double* uexpr_weights, const double* uexpr_cst, const size_t* uexpr_dim, const size_t expr_size)
-{
-    fppoly_t* res = (fppoly_t*) malloc(sizeof(fppoly_t));
-    fppoly_from_network_input_box(res, intdim, realdim, inf_array, sup_array);
-
-    size_t num_pixels = intdim + realdim;
-
-    cudaMalloc((void**) &res->input_lexpr, num_pixels*sizeof(expr_t*));
-    cudaMalloc((void**) &res->input_uexpr, num_pixels*sizeof(expr_t*));
-
-    layer_create_sparse_exprs(res->input_lexpr, lexpr_weights, lexpr_cst, lexpr_dim, num_pixels, expr_size);
-    layer_create_sparse_exprs(res->input_uexpr, uexpr_weights, uexpr_cst, uexpr_dim, num_pixels, expr_size);
-
-    return abstract0_of_fppoly(man, res);
-}
-*/
-
-
 void fppoly_add_new_layer(fppoly_t* const fp, const size_t num_out_neurons, const size_t num_in_neurons, const layertype_t type, const activation_type_t activation)
 {
     const size_t numlayers = fp->numlayers;

    }
}

/*
void compute_chord_slope(double *slope_inf, double *slope_sup, double f_sup_l, double f_sup_u,
             double f_inf_l, double f_inf_u, double inf_l, double inf_u, double sup_l, double sup_u){
    double num_l =  f_sup_l + f_inf_u;
    double num_u =  f_sup_u + f_inf_l;

    double den_l = sup_l + inf_u;
    double den_u = sup_u + inf_l;

        elina_double_interval_div(slope_inf, slope_sup, num_l, num_u, den_l, den_u);
}
*/

/*
void compute_derivative(double *slope_inf, double *slope_sup, double s_curve_l, double s_curve_u, double sq_l, double sq_u, bool is_sigmoid){
    double sq_den_sup_l, sq_den_sup_u;
    elina_double_interval_mul(&sq_den_sup_l, &sq_den_sup_u, sq_l, sq_u, sq_l, sq_u);
    if(is_sigmoid){
            elina_double_interval_div(slope_inf, slope_sup, s_curve_l, s_curve_u, sq_den_sup_l, sq_den_sup_u);
    }
    else{
            *slope_inf = -1 + sq_den_sup_u;
            *slope_sup = 1 + sq_den_sup_l;
    }

}
*/

/*
expr_t * lexpr_replace_s_curve_bounds(fppoly_internal_t * pr, expr_t * expr, neuron_t ** neurons, bool is_sigmoid){
    size_t num_neurons = expr->size;
    size_t i,k;
    expr_t * res = alloc_expr();
    res->inf_coeff = (double *)malloc(num_neurons*sizeof(double));
    res->sup_coeff = (double *)malloc(num_neurons*sizeof(double));
    res->inf_cst = expr->inf_cst;
    res->sup_cst = expr->sup_cst;
    res->type = expr->type;
    res->size = num_neurons;

    for(i = 0; i < num_neurons; i++){
        if(expr->type==DENSE){
            k = i;
        }
        else{
            k = expr->dim[i];
        }
        double lb = neurons[k]->lb;
        double ub = neurons[k]->ub;
        //fesetround(FE_DOWNWARD);
        double e_sup_l = is_sigmoid ? -exp(ub) : -tanh(ub);
        double e_inf_l = is_sigmoid ? -exp(-lb) : -tanh(-lb);

        //fesetround(FE_UPWARD);
        double e_sup_u = is_sigmoid ? exp(ub) : tanh(ub);
        double e_inf_u = is_sigmoid ? exp(-lb) : tanh(-lb);

        double f_sup_l, f_sup_u;
        double f_inf_l, f_inf_u;
        double den_sup_l, den_sup_u;
        double den_inf_l, den_inf_u;


        if(is_sigmoid){
            den_sup_l = -1 + e_sup_l;
            den_sup_u = 1 + e_sup_u;
            den_inf_l = -1 + e_inf_l;
            den_inf_u = 1 + e_inf_u;
            elina_double_interval_div(&f_sup_l, &f_sup_u, e_sup_l, e_sup_u, den_sup_l, den_sup_u);
            elina_double_interval_div(&f_inf_l, &f_inf_u, e_inf_l, e_inf_u, den_inf_l, den_inf_u);
        }
        else{
            f_inf_l = e_inf_l;
            f_inf_u = e_inf_u;
            f_sup_l = e_sup_l;
            f_sup_u = e_sup_u;
            den_inf_l = e_inf_l;
            den_inf_u = e_inf_u;
            den_sup_l = e_sup_l;
            den_sup_u = e_sup_u;
        }

        if((-lb==ub)|| (-f_inf_l==f_sup_u)){

            res->inf_coeff[i] = 0.0;
            res->sup_coeff[i] = 0.0;
            double tmp1, tmp2;
            elina_double_interval_mul(&tmp1,&tmp2,expr->inf_coeff[i],expr->sup_coeff[i],f_inf_l,f_sup_u);
            res->inf_cst = res->inf_cst + tmp1;
            res->sup_cst = res->sup_cst - tmp1;

        }

        else if(expr->sup_coeff[i]<0 || expr->inf_coeff[i] < 0){
            double slope_inf, slope_sup;
            double intercept_inf, intercept_sup;
            double add_inf, add_sup;
                    double mul_inf, mul_sup;
                    double x_l, x_u;
                    double f_x_l, f_x_u;
                    bool boxify = false;
                    if(expr->sup_coeff[i] < 0){
                        if(ub<0){

                           compute_chord_slope(&slope_inf, &slope_sup, f_sup_l, f_sup_u, f_inf_l, f_inf_u, lb, -lb, -ub, ub);
                       //if(slope_inf>0){
                    //boxify=true;
                       //}
                       //else{
                    //compute_derivative( &slope_inf, &slope_sup, e_inf_l, e_inf_u, den_inf_l, den_inf_u, is_sigmoid);
                    //printf("slope2: %.30f %.30f\n",slope_inf*(ub+lb),slope_sup*(ub+lb));
                       // }    //fflush(stdout);

                    x_l = ub;
                    x_u = -ub;
                    f_x_l = f_sup_l;
                    f_x_u = f_sup_u;
                    //elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,-lb, lb, slope_inf,slope_sup);
                    //elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_inf_l, f_inf_u, intercept_inf, intercept_sup);
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                    double tmp1, tmp2, tmp3, tmp4;
                           elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,lb, -lb, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);
                    if(tmp4<f_inf_u){

                                boxify = true;
                                }

                   //}
                //printf("slope: %.30f %.30f %.30f %.30f\n",slope_inf*(ub+lb),slope_sup*(ub+lb),- tmp3,tmp4);
                //    fflush(stdout);
                        //elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);

                        }
                         else if(lb<=0){

                            compute_derivative( &slope_inf, &slope_sup, e_sup_l, e_sup_u, den_sup_l, den_sup_u, is_sigmoid);
                           x_l = ub;
                            x_u = -ub;
                            f_x_l = f_sup_l;
                            f_x_u = f_sup_u;
                            elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                            elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                    double tmp1, tmp2, tmp3, tmp4;
                            elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,lb, -lb, slope_inf,slope_sup);
                            elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);
                    if(-tmp3<f_inf_u){

                                boxify = true;
                            }
                      }
                    else{

                            if(lb<=ub){
                                //double slope_inf1, slope_sup1;
                                //double slope_inf2, slope_sup2;
                                compute_derivative( &slope_inf, &slope_sup, e_sup_l, e_sup_u, den_sup_l, den_sup_u, is_sigmoid);

                            }
                            else{
                                compute_derivative( &slope_inf, &slope_sup, e_inf_l, e_inf_u, den_inf_l, den_inf_u, is_sigmoid);
                            }
                   // if(slope_inf1>=slope_sup2){
                     //   slope_inf = slope_inf2;
                       // slope_sup = slope_sup2;
                    //}
                    //else if(slope_inf2>=slope_sup1){
                      //  slope_inf = slope_inf1;
                        //slope_sup = slope_sup1;
                    //}
                    //else{
                      //  boxify = true;
                    //}
                    x_l = ub;
                    x_u = -ub;
                    f_x_l = f_sup_l;
                    f_x_u = f_sup_u;
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                    double tmp1, tmp2, tmp3, tmp4;
                    elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,lb, -lb, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);
                    if(-tmp3<f_inf_u){

                        boxify = true;
                    }
                    }




                }
                else{
                    if(ub < 0){

                    compute_derivative( &slope_inf, &slope_sup, e_inf_l, e_inf_u, den_inf_l, den_inf_u, is_sigmoid);
                    x_l = -lb;
                    x_u = lb;
                    f_x_l = f_inf_l;
                    f_x_u = f_inf_u;
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                double tmp1, tmp2, tmp3, tmp4;
                    elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,-ub, ub, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);
                    if(tmp4>-f_sup_l){

                        boxify = true;
                    }
                        }
                else if(lb<=0){
                    compute_chord_slope(&slope_inf, &slope_sup, f_sup_l, f_sup_u, f_inf_l, f_inf_u, lb, -lb, -ub, ub);
                //if(slope_inf>0){
                    //boxify=true;
                //}
               // compute_derivative( &slope_inf, &slope_sup, e_sup_l, e_sup_u, den_sup_l, den_sup_u, is_sigmoid);
                //}
                //else{
                    x_l = -lb;
                    x_u = lb;
                    f_x_l = f_inf_l;
                    f_x_u = f_inf_u;
                    //elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_mul(&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    //add_inf = f_x_l + intercept_inf;
                    //add_sup = f_x_u + intercept_sup;

                           elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                       double tmp1, tmp2, tmp3, tmp4;
                           elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,-ub, ub, slope_inf,slope_sup);
                           elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);

                   if(-tmp3>f_sup_u){

                            boxify = true;
                           }
                   //}
                }
                else{

                    //double slope_inf1, slope_sup1;
                    //double slope_inf2, slope_sup2;
                    if(lb<=ub){
                    compute_derivative( &slope_inf, &slope_sup, e_sup_l, e_sup_u, den_sup_l, den_sup_u, is_sigmoid);

                    }
                    else{
                    compute_derivative( &slope_inf, &slope_sup, e_inf_l, e_inf_u, den_inf_l, den_inf_u, is_sigmoid);
                    }
                    //if(slope_inf1>=slope_sup2){
                      //  slope_inf = slope_inf2;
                        //slope_sup = slope_sup2;
                   // }
                    //else if(slope_inf2>=slope_sup1){
                      //  slope_inf = slope_inf1;
                       // slope_sup = slope_sup1;
                    //}
                    //else{
                      //  boxify = true;
                    //}
                    x_l = -lb;
                    x_u = lb;
                    f_x_l = f_inf_l;
                    f_x_u = f_inf_u;
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                    double tmp1, tmp2, tmp3, tmp4;
                    elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,-ub, ub, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);
                    if(tmp4>-f_sup_l){

                        boxify = true;
                    }
                }


            }
            if(boxify){
                //printf("boxify lexpr\n");
            //fflush(stdout);
                res->inf_coeff[i] = 0.0;
                res->sup_coeff[i] = 0.0;
                double tmp1, tmp2;
                elina_double_interval_mul(&tmp1,&tmp2,expr->inf_coeff[i],expr->sup_coeff[i],f_inf_l,f_sup_u);
                res->inf_cst = res->inf_cst + tmp1;
                res->sup_cst = res->sup_cst - tmp1;

            }
            else{
                elina_double_interval_mul_expr_coeff(pr,&res->inf_coeff[i],&res->sup_coeff[i],slope_inf,slope_sup,expr->inf_coeff[i],expr->sup_coeff[i]);
                    elina_double_interval_mul_cst_coeff(pr, &mul_inf, &mul_sup, add_inf, add_sup, expr->inf_coeff[i], expr->sup_coeff[i] );
                    elina_double_interval_add_cst_coeff(pr,&res->inf_cst,&res->sup_cst,mul_inf, mul_sup, res->inf_cst, res->sup_cst);
            }

        }


        else{

            res->inf_coeff[i] = 0.0;
            res->sup_coeff[i] = 0.0;
            double tmp1, tmp2;
            elina_double_interval_mul(&tmp1,&tmp2,expr->inf_coeff[i],expr->sup_coeff[i],f_inf_l,f_sup_u);
            res->inf_cst = res->inf_cst + tmp1;
            res->sup_cst = res->sup_cst - tmp1;
        }
    }
    if(expr->type==SPARSE){
        res->dim = (size_t*)malloc(num_neurons*sizeof(size_t));
        for(i=0; i < num_neurons; i++){
            res->dim[i] = expr->dim[i];
        }
    }
    return res;
}
*/

/*
expr_t * uexpr_replace_s_curve_bounds(fppoly_internal_t *pr, expr_t * expr, neuron_t ** neurons, bool is_sigmoid){
    size_t num_neurons = expr->size;
    size_t i,k;
    expr_t * res = alloc_expr();
    res->inf_coeff = (double *)malloc(num_neurons*sizeof(double));
    res->sup_coeff = (double *)malloc(num_neurons*sizeof(double));
    res->inf_cst = expr->inf_cst;
    res->sup_cst = expr->sup_cst;
    res->type = expr->type;
    res->size = num_neurons;

    for(i = 0; i < num_neurons; i++){
        if(expr->type==DENSE){
            k = i;
        }
        else{
            k = expr->dim[i];
        }
        double lb = neurons[k]->lb;
        double ub = neurons[k]->ub;
        //fesetround(FE_DOWNWARD);
        double e_sup_l = is_sigmoid ? -exp(ub) : -tanh(ub);
        double e_inf_l = is_sigmoid ? -exp(-lb) : -tanh(-lb);

        //fesetround(FE_UPWARD);
        double e_sup_u = is_sigmoid ? exp(ub) : tanh(ub);
        double e_inf_u = is_sigmoid ? exp(-lb) : tanh(-lb);

        double f_sup_l, f_sup_u;
        double f_inf_l, f_inf_u;
        double den_sup_l, den_sup_u;
        double den_inf_l, den_inf_u;

        if(is_sigmoid){
            den_sup_l = -1 + e_sup_l;
            den_sup_u = 1 + e_sup_u;
            den_inf_l = -1 + e_inf_l;
            den_inf_u = 1 + e_inf_u;
            elina_double_interval_div(&f_sup_l, &f_sup_u, e_sup_l, e_sup_u, den_sup_l, den_sup_u);
            elina_double_interval_div(&f_inf_l, &f_inf_u, e_inf_l, e_inf_u, den_inf_l, den_inf_u);
        }
        else{
            f_inf_l = e_inf_l;
            f_inf_u = e_inf_u;
            f_sup_l = e_sup_l;
            f_sup_u = e_sup_u;
            den_inf_l = e_inf_l;
            den_inf_u = e_inf_u;
            den_sup_l = e_sup_l;
            den_sup_u = e_sup_u;
        }

        if((-lb==ub) || (-f_inf_l==f_sup_u)){
            //printf("boxify %.30f %.30f %.30f %.30f\n",-f_inf_l,f_inf_u,-f_sup_l,f_sup_u);
            res->inf_coeff[i] = 0.0;
            res->sup_coeff[i] = 0.0;
            double tmp1, tmp2;
            elina_double_interval_mul(&tmp1,&tmp2,expr->inf_coeff[i],expr->sup_coeff[i],f_inf_l,f_sup_u);
            res->inf_cst = res->inf_cst - tmp2;
            res->sup_cst = res->sup_cst + tmp2;
        }
       // else if(lb+ub<0.000000000001){
    //    printf("boxify %g %g %g %g\n",lb,ub,f_inf_u,f_sup_u);
    //    fflush(stdout);
          //  res->inf_coeff[i] = 0.0;
            //res->sup_coeff[i] = 0.0;
            //double tmp1, tmp2;
            //elina_double_interval_mul(&tmp1,&tmp2,expr->inf_coeff[i],expr->sup_coeff[i],f_inf_l,f_sup_u);
            //res->inf_cst = res->inf_cst + tmp1;
            //res->sup_cst = res->sup_cst - tmp1;
       // }
        else if(expr->sup_coeff[i]<0 || expr->inf_coeff[i] < 0){
            double slope_inf, slope_sup;
            double intercept_inf, intercept_sup;
            double add_inf, add_sup;
            double mul_inf, mul_sup;
            double x_l, x_u;
            double f_x_l, f_x_u;
            bool boxify = false;
            if(expr->sup_coeff[i] < 0){
                if(ub<0){

                    compute_derivative( &slope_inf, &slope_sup, e_inf_l, e_inf_u, den_inf_l, den_inf_u, is_sigmoid);
                    x_l = -lb;
                    x_u = lb;
                    f_x_l = f_inf_l;
                    f_x_u = f_inf_u;
                    //elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,-lb, lb, slope_inf,slope_sup);
                    //elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_inf_l, f_inf_u, intercept_inf, intercept_sup);
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
            double tmp1, tmp2, tmp3, tmp4;
                    elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,-ub, ub, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);
                    if(tmp4>-f_sup_l){
                        boxify = true;
                    }
                }
                else if(lb<=0){

                    compute_chord_slope(&slope_inf, &slope_sup, f_sup_l, f_sup_u, f_inf_l, f_inf_u, lb, -lb, -ub, ub);
                //if(slope_inf>0){
                //boxify = true;
                //}
                //else{
                        //compute_derivative( &slope_inf, &slope_sup, e_sup_l, e_sup_u, den_sup_l, den_sup_u, is_sigmoid);
               //}
                    x_l = -lb;
                    x_u = lb;
                    f_x_l = f_inf_l;
                    f_x_u = f_inf_u;
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                    double tmp1, tmp2, tmp3, tmp4;
                                elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,-ub, ub, slope_inf,slope_sup);
                                elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);
                                if(-tmp3>f_sup_u){
                                boxify = true;
                                }
                    //}
                }
                else{
                    //double slope_inf1, slope_sup1;
                    //double slope_inf2, slope_sup2;
                    if(lb<=ub){
                    compute_derivative( &slope_inf, &slope_sup, e_sup_l, e_sup_u, den_sup_l, den_sup_u, is_sigmoid);

                    }
                    else{
                    compute_derivative( &slope_inf, &slope_sup, e_inf_l, e_inf_u, den_inf_l, den_inf_u, is_sigmoid);
                    }
                    //if(slope_inf1>=slope_sup2){
                      //  slope_inf = slope_inf2;
                        //slope_sup = slope_sup2;
                    //}
                    //else if(slope_inf2>=slope_sup1){
                      //  slope_inf = slope_inf1;
                        //slope_sup = slope_sup1;
                    //}
                    //else{
                      //  boxify = true;
                    //}
                    x_l = -lb;
                    x_u = lb;
                    f_x_l = f_inf_l;
                    f_x_u = f_inf_u;
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                    double tmp1, tmp2, tmp3, tmp4;
                    elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,-ub, ub, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);
                    if(tmp4>-f_sup_l){
                        boxify = true;
                    }
                }

            }
            else{
                if(ub < 0){

                    compute_chord_slope(&slope_inf, &slope_sup, f_sup_l, f_sup_u, f_inf_l, f_inf_u, lb, -lb, -ub, ub);
                    //if(slope_inf>0){
                //boxify = true;
                //}
                //else{
                    //compute_derivative( &slope_inf, &slope_sup, e_inf_l, e_inf_u, den_inf_l, den_inf_u, is_sigmoid);
                //}
                    x_l = ub;
                    x_u = -ub;
                    f_x_l = f_sup_l;
                    f_x_u = f_sup_u;
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                    double tmp1, tmp2, tmp3, tmp4;
                                elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,lb, -lb, slope_inf,slope_sup);
                                elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);

                                if(tmp4<f_inf_u){
                                boxify = true;
                                }
                   // }


                }
                else if(lb<=0){

                    compute_derivative( &slope_inf, &slope_sup, e_sup_l, e_sup_u, den_sup_l, den_sup_u, is_sigmoid);

                    x_l = ub;
                    x_u = -ub;
                    f_x_l = f_sup_l;
                    f_x_u = f_sup_u;
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                    double tmp1, tmp2, tmp3, tmp4;
                    elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,lb, -lb, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);

                    if(-tmp3<f_inf_u){
                        boxify = true;
                    }


                }
                else{
                    //double slope_inf1, slope_sup1;
                    //double slope_inf2, slope_sup2;
                    if(lb<=ub){
                    compute_derivative( &slope_inf, &slope_sup, e_sup_l, e_sup_u, den_sup_l, den_sup_u, is_sigmoid);

                    }
                    else{
                    compute_derivative( &slope_inf, &slope_sup, e_inf_l, e_inf_u, den_inf_l, den_inf_u, is_sigmoid);
                    }
                    //if(slope_inf1>=slope_sup2){
                      //  slope_inf = slope_inf2;
                        //slope_sup = slope_sup2;
                    //}
                    //else if(slope_inf2>=slope_sup1){
                      //  slope_inf = slope_inf1;
                        //slope_sup = slope_sup1;
                    //}
                    //else{
                      //  boxify = true;
                    //}

                    x_l = ub;
                    x_u = -ub;
                    f_x_l = f_sup_l;
                    f_x_u = f_sup_u;
                    elina_double_interval_mul_cst_coeff(pr,&intercept_inf,&intercept_sup,x_l, x_u, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&add_inf,&add_sup,f_x_l, f_x_u, intercept_inf, intercept_sup);
                    double tmp1, tmp2, tmp3, tmp4;
                    elina_double_interval_mul_cst_coeff(pr,&tmp1,&tmp2,lb, -lb, slope_inf,slope_sup);
                    elina_double_interval_add_cst_coeff(pr,&tmp3,&tmp4,add_inf, add_sup, tmp1, tmp2);
                    if(-tmp3<f_inf_u){
                        boxify = true;
                    }

                }


            }
       // printf("slope inf: %g slope sup: %g\n",slope_inf,slope_sup);
       // fflush(stdout);
            if(boxify){
                //printf("boxify uexpr\n");
                //fflush(stdout);
                res->inf_coeff[i] = 0.0;
                res->sup_coeff[i] = 0.0;
                double tmp1, tmp2;
                elina_double_interval_mul(&tmp1,&tmp2,expr->inf_coeff[i],expr->sup_coeff[i],f_inf_l,f_sup_u);
                res->inf_cst = res->inf_cst - tmp2;
                res->sup_cst = res->sup_cst + tmp2;

            }
            else{
            elina_double_interval_mul_expr_coeff(pr,&res->inf_coeff[i],&res->sup_coeff[i],slope_inf,slope_sup,expr->inf_coeff[i],expr->sup_coeff[i]);

            elina_double_interval_mul_cst_coeff(pr, &mul_inf, &mul_sup, add_inf, add_sup, expr->inf_coeff[i], expr->sup_coeff[i] );
            elina_double_interval_add_cst_coeff(pr,&res->inf_cst,&res->sup_cst,mul_inf, mul_sup, res->inf_cst, res->sup_cst);
            }
           // printf("cst: %g %g\n",res->inf_cst,res->sup_cst);
        //fflush(stdout);
        }


        else{

            res->inf_coeff[i] = 0.0;
            res->sup_coeff[i] = 0.0;
            double tmp1, tmp2;
            elina_double_interval_mul(&tmp1,&tmp2,expr->inf_coeff[i],expr->sup_coeff[i],f_inf_l,f_sup_u);
            res->inf_cst = res->inf_cst - tmp2;
            res->sup_cst = res->sup_cst + tmp2;
        }
    }
    if(expr->type==SPARSE){
        res->dim = (size_t*)malloc(num_neurons*sizeof(size_t));
        for(i=0; i < num_neurons; i++){
            res->dim[i] = expr->dim[i];
        }
    }
    return res;
}
*/


/*
expr_t * uexpr_replace_sigmoid_bounds(fppoly_internal_t *pr, expr_t * expr, neuron_t ** neurons){
    return uexpr_replace_s_curve_bounds(pr, expr, neurons, true);
}

expr_t * uexpr_replace_tanh_bounds(fppoly_internal_t *pr, expr_t * expr, neuron_t ** neurons){
    return uexpr_replace_s_curve_bounds(pr, expr, neurons, false);
}

expr_t * lexpr_replace_sigmoid_bounds(fppoly_internal_t *pr, expr_t * expr, neuron_t ** neurons){
    return lexpr_replace_s_curve_bounds(pr, expr, neurons, true);
}

expr_t * lexpr_replace_tanh_bounds(fppoly_internal_t *pr, expr_t * expr, neuron_t ** neurons){
    return lexpr_replace_s_curve_bounds(pr, expr, neurons, false);
}
*/


__global__
void lexpr_replace_maxpool_bounds(expr_t** expr_array, expr_t** aux_maxpool_lexpr_array, expr_t** aux_maxpool_uexpr_array, double* aux_lb_array, double* aux_ub_array, const size_t num_out_neurons)
{
    size_t n = blockIdx.x;

    if(n < num_out_neurons)
    {
        expr_t* expr = expr_array[n];

        const size_t num_neurons = expr->size;

        expr_t* res;

        size_t k;

        if(expr->type == DENSE)
        {
            k = 0;
        }
        else
        {
            k = expr->dim[0];
        }

        const double lb = aux_lb_array[k];
        const double ub = aux_ub_array[k];

        expr_t* maxpool_lexpr = aux_maxpool_lexpr_array[k];
        expr_t* maxpool_uexpr = aux_maxpool_lexpr_array[k];

        if(expr->sup_coeff[0] < 0)
        {
            //expr_print(maxpool_uexpr);
            if(maxpool_uexpr == nullptr)
            {
                res = (expr_t*) malloc(sizeof(expr_t));
                res->inf_coeff = res->sup_coeff =  nullptr;
                res->dim = nullptr;
                res->size = 0;
                res->type = SPARSE;
                elina_double_interval_mul_cst_coeff(&res->inf_cst, &res->sup_cst, lb, ub, expr->inf_coeff[0], expr->sup_coeff[0]);
            }
            else
            {
                res = multiply_expr(maxpool_uexpr, expr->inf_coeff[0], expr->sup_coeff[0]);
            }
            //printf("multiply end %zu \n",k);
            //expr_print(res);
            //fflush(stdout);
        }
        else if(expr->inf_coeff[0] < 0)
        {
            //expr_print(maxpool_lexpr);
            if(maxpool_lexpr == nullptr)
            {
                res = (expr_t*) malloc(sizeof(expr_t));
                res->inf_coeff = nullptr;
                res->sup_coeff = nullptr;
                res->dim = nullptr;
                res->size = 0;
                res->type = SPARSE;
                elina_double_interval_mul_cst_coeff(&res->inf_cst, &res->sup_cst, lb, ub, expr->inf_coeff[0], expr->sup_coeff[0]);
            }
            else
            {
                res = multiply_expr(maxpool_lexpr, expr->inf_coeff[0], expr->sup_coeff[0]);
            }
            //printf("multiply end %zu \n",k);
            //expr_print(res);
            //fflush(stdout);
        }
        else
        {
            //printf("WTF1\n");
            //fflush(stdout);
            double tmp1, tmp2;
            elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, lb, ub, expr->inf_coeff[0], expr->sup_coeff[0]);
            double coeff[1];
            size_t dim[1];
            coeff[0] = 0;
            dim[0] = 0;
            res = create_sparse_expr(coeff, -tmp1, dim, 1);
        }
        //printf("middle\n");
        //fflush(stdout);
        expr_t* mul_expr;

        for(size_t i = 1; i < num_neurons; i++)
        {
            if(expr->type == DENSE)
            {
                k = i;
            }
            else
            {
                k = expr->dim[i];
            }

            const double lb = aux_lb_array[k];
            const double ub = aux_ub_array[k];

            expr_t* maxpool_lexpr = aux_maxpool_lexpr_array[k];
            expr_t* maxpool_uexpr = aux_maxpool_lexpr_array[k];

            if(expr->sup_coeff[i] < 0)
            {
                //expr_print(maxpool_uexpr);
                //printf("add start %zu %zu\n",k,i);

                    //expr_print(res);

                if(maxpool_uexpr == nullptr)
                {
                    mul_expr = (expr_t*) malloc(sizeof(expr_t));
                    mul_expr->inf_coeff = nullptr;
                    mul_expr->sup_coeff = nullptr;
                    mul_expr->dim = nullptr;
                    mul_expr->size = 0;
                    mul_expr->type = SPARSE;
                    //printf("lb: %g %g\n");
                    elina_double_interval_mul_cst_coeff(&mul_expr->inf_cst, &mul_expr->sup_cst, lb, ub, expr->inf_coeff[i], expr->sup_coeff[i]);
                    res->inf_cst += mul_expr->inf_cst;
                    res->sup_cst += mul_expr->sup_cst;
                }
                else
                {
                    mul_expr = multiply_expr(maxpool_uexpr, expr->inf_coeff[i], expr->sup_coeff[i]);
                    add_expr(res, mul_expr);
                }
                //expr_print(mul_expr);
                    //fflush(stdout);
                //printf("add finish\n");
                    //expr_print(res);
                    //fflush(stdout);
                free_expr(mul_expr);
            }
            else if (expr->inf_coeff[i] < 0)
            {
                //expr_print(maxpool_lexpr);
                //printf("add start %zu %zu\n",k,i);

                    //expr_print(res);

                if(maxpool_lexpr == nullptr)
                {
                    mul_expr = (expr_t*) malloc(sizeof(expr_t));
                    mul_expr->inf_coeff = nullptr;
                    mul_expr->sup_coeff = nullptr;
                    mul_expr->dim = nullptr;
                    mul_expr->size = 0;
                    mul_expr->type = SPARSE;
                    elina_double_interval_mul_cst_coeff(&mul_expr->inf_cst, &mul_expr->sup_cst, lb, ub, expr->inf_coeff[i], expr->sup_coeff[i]);
                    res->inf_cst += mul_expr->inf_cst;
                    res->sup_cst += mul_expr->sup_cst;
                }
                else
                {
                    mul_expr = multiply_expr(maxpool_lexpr, expr->inf_coeff[i], expr->sup_coeff[i]);
                    //printf("add start1 %zu %zu\n",k,i);
                    //expr_print(res);
                    //expr_print(mul_expr);
                    //fflush(stdout);
                    add_expr(res, mul_expr);
                }
                //expr_print(mul_expr);
                //    fflush(stdout);
                //printf("add finish1\n");
                //expr_print(res);
                //fflush(stdout);
                free_expr(mul_expr);
            }
            else
            {
                //printf("WTF2\n");
                //fflush(stdout);
                double tmp1, tmp2;
                elina_double_interval_mul_expr_coeff(&tmp1, &tmp2, lb, ub, expr->inf_coeff[i], expr->sup_coeff[i]);
                res->inf_cst = res->inf_cst + tmp1;
                res->sup_cst = res->sup_cst - tmp1;
            }
        }
        //printf("finish\n");
        //fflush(stdout);
        res->inf_cst = res->inf_cst + expr->inf_cst;
        res->sup_cst = res->sup_cst + expr->sup_cst;

        free_expr(expr);
        expr_array[n] = res;
    }
}


__global__
void uexpr_replace_maxpool_bounds(expr_t** expr_array, expr_t** aux_maxpool_lexpr_array, expr_t** aux_maxpool_uexpr_array, double* aux_lb_array, double* aux_ub_array, const size_t num_out_neurons)
{
    size_t n = blockIdx.x;

    if(n < num_out_neurons)
    {
        expr_t* expr = expr_array[n];

        const size_t num_neurons = expr->size;

        expr_t* res;

        size_t k;

        if(expr->type == DENSE)
        {
            k = 0;
        }
        else
        {
            k = expr->dim[0];
        }

        const double lb = aux_lb_array[k];
        const double ub = aux_ub_array[k];

        expr_t* maxpool_lexpr = aux_maxpool_lexpr_array[k];
        expr_t* maxpool_uexpr = aux_maxpool_lexpr_array[k];

        if(expr->sup_coeff[0] < 0)
        {
            res = multiply_expr(maxpool_lexpr, expr->inf_coeff[0], expr->sup_coeff[0]);
        }
        else if(expr->inf_coeff[0] < 0)
        {
            res = multiply_expr(maxpool_uexpr, expr->inf_coeff[0], expr->sup_coeff[0]);
        }
        else
        {
            double tmp1, tmp2;
            elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, lb, ub, expr->inf_coeff[0], expr->sup_coeff[0]);
            double coeff[1];
            size_t dim[1];
            coeff[0] = 0;
            dim[0] = 0;
            res = create_sparse_expr(coeff, tmp2, dim, 1);
        }

        for(size_t i = 1; i < num_neurons; i++)
        {
            if(expr->type == DENSE)
            {
                k = i;
            }
            else
            {
                k = expr->dim[i];
            }

            const double lb = aux_lb_array[k];
            const double ub = aux_ub_array[k];

            expr_t* maxpool_lexpr = aux_maxpool_lexpr_array[k];
            expr_t* maxpool_uexpr = aux_maxpool_lexpr_array[k];

            if(expr->sup_coeff[i] < 0)
            {
                expr_t* mul_expr = multiply_expr(maxpool_lexpr, expr->inf_coeff[i], expr->sup_coeff[i]);
                add_expr(res, mul_expr);
                free_expr(mul_expr);
            }
            else if (expr->inf_coeff[i] < 0)
            {
                expr_t* mul_expr = multiply_expr(maxpool_uexpr, expr->inf_coeff[i], expr->sup_coeff[i]);
                add_expr(res, mul_expr);
                free_expr(mul_expr);
            }
            else
            {
                double tmp1, tmp2;
                elina_double_interval_mul_cst_coeff(&tmp1, &tmp2, lb, ub, expr->inf_coeff[i], expr->sup_coeff[i]);
                res->inf_cst = res->inf_cst - tmp2;
                res->sup_cst = res->sup_cst + tmp2;
            }
        }

        res->inf_cst = res->inf_cst + expr->inf_cst;
        res->sup_cst = res->sup_cst + expr->sup_cst;

        free_expr(expr);
        expr_array[n] = res;
    }
}


__global__
void expr_from_previous_layer(expr_t** expr_array, expr_t** aux_expr_array, const size_t num_out_neurons)
{
    size_t n = blockIdx.x;

    if(n < num_out_neurons)
    {
        expr_t* expr = expr_array[n];

        if(expr->size == 0)
        {
            return;
        }

        if((expr->inf_coeff == nullptr) || (expr->sup_coeff == nullptr))
        {
            return;
        }

        const size_t in_num_neurons = expr->size;
        size_t k;
        expr_t* res;

        if(expr->type == DENSE)
        {
            k = 0;
        }
        else
        {
            k = expr->dim[0];
        }

        if(aux_expr_array[k]->size == 0)
        {
            res = multiply_cst_expr(aux_expr_array[k], expr->inf_coeff[0], expr->sup_coeff[0]);
        }
        else
        {
            res = multiply_expr(aux_expr_array[k], expr->inf_coeff[0], expr->sup_coeff[0]);
        }

        for(size_t i = 1; i < in_num_neurons; i++)
        {
            if(expr->type == DENSE)
            {
                k = i;
            }
            else
            {
                k = expr->dim[i];
            }

            expr_t* mul_expr;

            if(aux_expr_array[k]->size == 0)
            {
                mul_expr = multiply_cst_expr(aux_expr_array[k], expr->inf_coeff[i], expr->sup_coeff[i]);
                add_cst_expr(res, mul_expr);
                free_expr(mul_expr);
            }
            else if((expr->inf_coeff[i] != 0) || (expr->sup_coeff[i] != 0))
            {
                mul_expr = multiply_expr(aux_expr_array[k], expr->inf_coeff[i], expr->sup_coeff[i]);
                add_expr(res, mul_expr);
                free_expr(mul_expr);
            }
        }

        res->inf_cst = res->inf_cst + expr->inf_cst;
        res->sup_cst = res->sup_cst + expr->sup_cst;

        free_expr(expr);
        expr_array[n] = res;
    }
}


__global__
void expr_from_previous_layer_specialized(expr_t** expr_array, expr_t** res_array, expr_t** aux_expr_array, const size_t num_out_neurons_last_layer, const size_t num_out_neurons_current_layer, const size_t num_in_neurons_current_layer)
{
    size_t n = blockIdx.x;
    size_t j = blockIdx.y*blockDim.x + threadIdx.x;

        std::cout << "num_threads" << num_threads << " num_blocks_relu " << num_blocks_relu.y << " num_blocks_linear " << num_blocks_linear.y << std::endl;

        expr_t** aux_expr_array = fp->layers[k]->expr_array;

        double* aux_lb_array = fp->layers[k]->lb_array;
        double* aux_ub_array = fp->layers[k]->ub_array;

        if(fp->layers[k]->activation == RELU)
        {
            lexpr_replace_relu_bounds<<<num_blocks_relu, num_threads>>>(lexpr_array, aux_lb_array, aux_ub_array, num_out_neurons_last_layer, num_out_neurons_current_layer);
            uexpr_replace_relu_bounds<<<num_blocks_relu, num_threads>>>(uexpr_array, aux_lb_array, aux_ub_array, num_out_neurons_last_layer, num_out_neurons_current_layer);
        }

        layer_allocate_exprs<<<num_out_neurons_last_layer, 1>>>(lexpr_array_tmp, num_out_neurons_last_layer, num_in_neurons_current_layer);
        layer_allocate_exprs<<<num_out_neurons_last_layer, 1>>>(uexpr_array_tmp, num_out_neurons_last_layer, num_in_neurons_current_layer);

        expr_from_previous_layer<<<num_blocks_linear, num_threads>>>(lexpr_array, lexpr_array_tmp, aux_expr_array, num_out_neurons_last_layer, num_out_neurons_current_layer, num_in_neurons_current_layer);
        expr_from_previous_layer<<<num_blocks_linear, num_threads>>>(uexpr_array, uexpr_array_tmp, aux_expr_array, num_out_neurons_last_layer, num_out_neurons_current_layer, num_in_neurons_current_layer);

        std::swap(lexpr_array, lexpr_array_tmp);
        std::swap(uexpr_array, uexpr_array_tmp);

        free_expr_array<<<num_out_neurons_last_layer, 1>>>(lexpr_array_tmp, num_out_neurons_current_layer);
        free_expr_array<<<num_out_neurons_last_layer, 1>>>(uexpr_array_tmp, num_out_neurons_current_layer);
    }

    compute_lb_from_expr<<<num_out_neurons_last_layer, 1>>>(lb_array, lexpr_array, fp->input_inf, fp->input_sup, num_out_neurons_last_layer);
