#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN

#include <stdio.h>

#include <iostream>

#include <string>

#include <functional>

#include <stdio.h>

#include <cstring>

#include <fstream>

#include <stdio.h>

#include <stdlib.h>

#include <gravity/solver.h>

#include <gravity/doctest.h>

#include <PowerNet.h>

using namespace std;

using namespace gravity;

TEST_CASE("Add in row") {

Model<> M("Test");

/* Symbolic variable x in R^8 */

param<> x_lb("x_lb"), x_ub("x_ub");

x_lb = {100,1000,1000,10,10,10,10,10};

x_ub = {10000,10000,10000,1000,1000,1000,1000,1000};

var<> x("x",x_lb,x_ub);

/* Symbolic variable x in R^8 */

param<> y_lb("y_lb"), y_ub("y_ub");

y_lb = {-100,-1000,-1000,-10,-10,-10,-10,-10};

y_ub = {2,2,0,10,12,-1,-2,-3};

var<> y("y",y_lb,y_ub);

auto x_ids = indices("x_ids"), y_ids = indices("y_ids");

x_ids = range(1,8);

y_ids = range(1,8);

M.add(x.in(x_ids), y.in(y_ids));

indices y_mat("y_mat");/* Declaring a matrix index set based on the indices of y */

y_mat = y_ids; /* y_mat is based on y_ids */

for(int i = 0; i<4; i++){

y_mat.add_empty_row();

for(int j = 0; j<7; j++){

if(i%2==0 && j%(i+1)==0)

y_mat.add_in_row(i, to_string(j+1));

}

if(i%2!=0)

y_mat.add_in_row(i, "1");

}

Constraint<> LinCons("LinCons");

LinCons = x.in(range(1,4)) - x.in(range(2,5)) - y.in(y_mat);

M.add(LinCons.in(range(1,4)) == 0);

M.print();

CHECK(M.get_nb_vars()==16);

CHECK(M.get_nb_cons()==4);

}

TEST_CASE("Add in row, weighted sum") {

gravity::Model<> M("Test");

/* Symbolic variable x in R^8 */

gravity::param<> x_lb("x_lb"), x_ub("x_ub");

x_lb = {100,1000,1000,10,10,10,10,10};

x_ub = {10000,10000,10000,1000,1000,1000,1000,1000};

gravity::var<> x("x",x_lb,x_ub);

/* Symbolic variable x in R^8 */

gravity::param<> y_lb("y_lb"), y_ub("y_ub");

y_lb = {-100,-1000,-1000,-10,-10,-10,-10,-10};

y_ub = {2,2,0,10,12,-1,-2,-3};

gravity::var<> y("y",y_lb,y_ub);

auto x_ids = gravity::indices("x_ids"), y_ids = gravity::indices("y_ids");

x_ids = gravity::range(1,8);

y_ids = gravity::range(1,8);

M.add(x.in(x_ids), y.in(y_ids));

gravity::indices y_mat("y_mat");/* Declaring a matrix index set based on the indices of y */

y_mat = y_ids; /* y_mat is based on y_ids */

gravity::param<double> w("w");

w = {0.1,0.2,0.3,0.4,0.5,0.6,0.7};

w.in(range(1,7));

gravity::indices w_mat("weight mat");

w_mat = gravity::range(1,7);

y_mat.add_empty_row();

y_mat.add_in_row(0,"1");

y_mat.add_in_row(0,"2");

y_mat.add_in_row(0,"3");

w_mat.add_empty_row();

w_mat.add_in_row(0,"1");

w_mat.add_in_row(0,"2");

w_mat.add_in_row(0,"3");

y_mat.add_empty_row();

y_mat.add_in_row(1,"1");

y_mat.add_in_row(1,"2");

w_mat.add_empty_row();

w_mat.add_in_row(1,"4");

w_mat.add_in_row(1,"5");

y_mat.add_empty_row();

y_mat.add_in_row(2,"1");

w_mat.add_empty_row();

w_mat.add_in_row(2,"6");

y_mat.add_empty_row();

y_mat.add_in_row(3,"1");

w_mat.add_empty_row();

w_mat.add_in_row(3,"7");

gravity::Constraint<> LinCons("LinCons");

LinCons = x.in(range(1,4)) - x.in(range(2,5)) - w.in(w_mat) * y.in(y_mat);

M.add(LinCons.in(range(1,4)) == 0);

M.print();

}

#ifdef USE_MP

TEST_CASE("testing readNL() function on ex4.nl") {

Model<> M;

string NL_file = string(prj_dir)+"/data_sets/NL/ex4.nl";

int status = M.readNL(NL_file);

M.print();

CHECK(status==0);

CHECK(M.get_nb_vars()==36);

CHECK(M.get_nb_vars()==36);

CHECK(M.get_nb_cons()==30);

CHECK(M.is_convex());

M.restructure();

DebugOn("Done Restructuring!");

M.write();

}

TEST_CASE("testing readNL() function on waterlund32.nl") {

Model<> M;

string NL_file = string(prj_dir)+"/data_sets/NL/waterund32.nl";

int status = M.readNL(NL_file);

M.print();

CHECK(status==0);

CHECK(M.get_nb_vars()==660);

CHECK(M.get_nb_cons()==200);

CHECK(!M.is_convex());

M.restructure();

DebugOn("Done Restructuring!");

}

#endif

TEST_CASE("testing param, var anf func copy operators") {

indices ids("ids");

ids.add("key1", "key2");

param<int> ip("ip");

ip.in(ids);

ip.print();

ip.set_val("key1",2);

ip.set_val("key2",3);

ip.print();

CHECK(ip._range->first==2);

CHECK(ip._range->second==3);

param<int> ip2 = ip("key1");/* ip2 is a masked copy of ip, they will be sharing the same _val vector */

ip2.print();

CHECK(ip2.is_param());

CHECK(ip2.get_intype()==integer_);

CHECK(ip2.is_positive());

param<int> ip3 = ip.deep_copy();/* ip3 will not be sharing the same _val vector with ip or ip2 */

ip.set_val("key1",-1);

CHECK(ip._range->first==-1);/* The range of ip should be updated */

CHECK(ip2._range->first==-1);/* So is the range of ip2 */

CHECK(ip3.eval("key1")==2);/* The _val of ip3 was not modified */

var<> v1("v1", -2,3), v2("v2", -10,40);

v1.in(ids);v2.in(ids);

func<> f = v1 + v1.get_lb()*v2;

f.print();

auto fcpy = f;

func<> deep_cpy;

deep_cpy.deep_copy(f);

fcpy.print();

v1.set_lb("key1", -3);

f.uneval();

f.eval_all();

fcpy.uneval();

fcpy.eval_all();

f.print();

fcpy.print();

deep_cpy.print();

CHECK(fcpy.to_str(0,5)=="v1[key1] - 3v2[key1]");

CHECK(deep_cpy.to_str(0,5)=="v1[key1] - 2v2[key1]");

}

TEST_CASE("Variable Scaling") {

Model<> M("Test");

param<> lb("x_lb");

lb = {100,1000,1000,10,10,10,10,10};

param<> ub("x_ub");

ub = {10000,10000,10000,1000,1000,1000,1000,1000};

var<> x("x",lb,ub);

M.add(x.in(range(1,8)));

Constraint<> C1("C1");

C1 = 0.0025*(x[4] + x[6]);

M.add(C1 <= 1);

Constraint<> C2("C2");

C2 = 0.0025*(x[5] - x[4] + x[7]);

M.add(C2 <= 1);

Constraint<> C3("C3");

C3 = 0.01*(x[8]-x[5]);

M.add(C3 <= 1);

Constraint<> C4("C4");

C4 = 100*x[1] - x[1]*x[6] + 833.33252*x[4];

M.add(C4 <= 83333.333);

Constraint<> C5("C5");

C5 = x[2]*x[4] - x[2]*x[7] - 1250*x[4] + 1250*x[5];

M.add(C5 <= 0);

Constraint<> C6("C6");

C6 = x[3]*x[5] - x[3]*x[8] - 2500*x[5] + 1250000;

M.add(C6 <= 0);

M.min(x[1]+x[2]+x[3]);

M.print();

solver<> s(M,ipopt);

s.run(5,1e-6);

auto obj_val = M.get_obj_val();

auto M_scale = M;

M_scale.reset();

M_scale.scale_vars(10);

M_scale.print();

double coef_scale = 100;

M_scale.scale_coefs(coef_scale);

M_scale.print();

M_scale.initialize_midpoint();

solver<> s2(M_scale,ipopt);

s2.run(5,1e-6);

CHECK(std::abs(M_scale.get_obj_val()-obj_val) < 1e-3);

M_scale.write_solution();

}

TEST_CASE("Model.relax()") {

Model<> M("Test");

param<> lb("x_lb");

lb = {100,1000,1000,10,10,10,10,10};

param<> ub("x_lb");

ub = {10000,10000,10000,1000,1000,1000,1000,1000};

var<> x("x",lb,ub);

auto x_ids = indices("x_ids");

x_ids = range(1,8);

M.add(x.in(x_ids));

indices x_mat("x_mat"), x_mat_RLT1("x_mat_RLT1"), x_mat_RLT2("x_mat_RLT2"), x_mat_RLT3("x_mat_RLT3");

param<> A("A");

A.in(x_ids);

A._indices->add_in_row(0, "1");

A._indices->add_in_row(0, "2");

x_mat.add_in_row(0, "4");

x_mat.add_in_row(0, "6");

x_mat_RLT1.add_in_row(0, "1");

x_mat_RLT1.add_in_row(0, "1");

x_mat_RLT2.add_in_row(0, "2");

x_mat_RLT2.add_in_row(0, "2");

x_mat_RLT3.add_in_row(0, "3");

x_mat_RLT3.add_in_row(0, "3");

A.set_val("1", 0.0025);A.set_val("2", 0.0025);

A._indices->add_in_row(1, "3");

A._indices->add_in_row(1, "4");

A._indices->add_in_row(1, "5");

x_mat.add_in_row(1, "5");

x_mat.add_in_row(1, "4");

x_mat.add_in_row(1, "7");

x_mat_RLT1.add_in_row(1, "1");

x_mat_RLT1.add_in_row(1, "1");

x_mat_RLT1.add_in_row(1, "1");

x_mat_RLT2.add_in_row(1, "2");

x_mat_RLT2.add_in_row(1, "2");

x_mat_RLT2.add_in_row(1, "2");

x_mat_RLT3.add_in_row(1, "3");

x_mat_RLT3.add_in_row(1, "3");

x_mat_RLT3.add_in_row(1, "3");

A.set_val("3", 0.0025);A.set_val("4", -0.0025);A.set_val("5", 0.0025);

A._indices->add_in_row(2, "6");

A._indices->add_in_row(2, "7");

x_mat.add_in_row(2, "8");

x_mat.add_in_row(2, "5");

x_mat_RLT1.add_in_row(2, "1");

x_mat_RLT1.add_in_row(2, "1");

x_mat_RLT2.add_in_row(2, "2");

x_mat_RLT2.add_in_row(2, "2");

x_mat_RLT3.add_in_row(2, "3");

x_mat_RLT3.add_in_row(2, "3");

A.set_val("6", 0.01);A.set_val("7", -0.01);

Constraint<> LinCons("LinCons");

LinCons = A*x.in(x_mat);

M.add(LinCons.in(range(1,3)) <= 1);

Constraint<> RLT1("RLT1");

RLT1 = A*x.in(x_mat_RLT1)*x.in(x_mat) - x.in(x.repeat_id(3,0));

M.add(RLT1.in(range(1,3)) <= 0);

Constraint<> RLT2("RLT2");

RLT2 = A*x.in(x_mat_RLT2)*x.in(x_mat) - x.in(x.repeat_id(3,1));

M.add(RLT2.in(range(1,3)) <= 0);

Constraint<> RLT3("RLT3");

RLT3 = A*x.in(x_mat_RLT3)*x.in(x_mat) - x.in(x.repeat_id(3,2));

M.add(RLT3.in(range(1,3)) <= 0);

/*

Constraint<> C4("C4");

C4 = 100*x[1] - x[1]*x[6] + 833.33252*x[4];

M.add(C4 <= 83333.333);

Constraint<> C5("C5");

C5 = x[2]*x[4] - x[2]*x[7] - 1250*x[4] + 1250*x[5];

M.add(C5 <= 0);

Constraint<> C6("C6");

C6 = x[3]*x[5] - x[3]*x[8] - 2500*x[5] + 1250000;

M.add(C6 <= 0);

M.min(x[1]+x[2]+x[3]);

M.print();

auto determinant_level = 1;

bool add_Kim_Kojima = false, add_SDP_3d = false;

auto LB = M.relax(determinant_level,add_Kim_Kojima, add_SDP_3d);

LB->print();

double max_time = 54000,ub_solver_tol=1e-6, lb_solver_tol=1e-6, range_tol=1e-4, opt_rel_tol=1e-2, opt_abs_tol=1e6;

unsigned max_iter=30, nb_threads = thread::hardware_concurrency();

SolverType ub_solver_type = ipopt, lb_solver_type = ipopt;

M.run_obbt(LB, max_time, max_iter, opt_rel_tol, opt_abs_tol, nb_threads=1, ub_solver_type, lb_solver_type, ub_solver_tol, lb_solver_tol, range_tol);

LB->print_constraints_stats(1e-6);

LB->print();

auto final_gap = 100*(M.get_obj_val() - LB->get_obj_val())/std::abs(M.get_obj_val());

CHECK(final_gap<1);

}

TEST_CASE("hard nlp") {

Model<> M("Test");

param<> lb("x_lb");

lb = {100,1000,1000,10,10,10,10,10};

param<> ub("x_lb");

ub = {10000,10000,10000,1000,1000,1000,1000,1000};

var<> x("x",lb,ub);

M.add(x.in(range(1,8)));

Constraint<> C1("C1");

C1 = 0.0025*(x[4] + x[6]);

M.add(C1 <= 1);

Constraint<> RLT1_1("RLT1_1");

RLT1_1 = x[1]*(0.0025*(x[4] + x[6]) - 1);

M.add(RLT1_1 <= 0);

Constraint<> RLT1_2("RLT1_2");

RLT1_2 = x[2]*(0.0025*(x[4] + x[6]) - 1);

M.add(RLT1_2 <= 0);

Constraint<> RLT1_3("RLT1_3");

RLT1_3 = x[3]*(0.0025*(x[4] + x[6]) - 1);

M.add(RLT1_3 <= 0);

Constraint<> C2("C2");

C2 = 0.0025*(x[5] - x[4] + x[7]);

M.add(C2 <= 1);

Constraint<> RLT2_1("RLT2_1");

RLT2_1 = x[1]*(0.0025*(x[5] - x[4] + x[7])-1);

M.add(RLT2_1 <= 0);

Constraint<> RLT2_2("RLT2_2");

RLT2_2 = x[2]*(0.0025*(x[5] - x[4] + x[7])-1);

M.add(RLT2_2 <= 0);

Constraint<> RLT2_3("RLT2_3");

RLT2_3 = x[3]*(0.0025*(x[5] - x[4] + x[7])-1);

M.add(RLT2_3 <= 0);

Constraint<> C3("C3");

C3 = 0.01*(x[8]-x[5]);

M.add(C3 <= 1);

Constraint<> RLT3_1("RLT3_1");

RLT3_1 = x[1]*(0.01*(x[8]-x[5])-1);

M.add(RLT3_1 <= 0);

Constraint<> RLT3_2("RLT3_2");

RLT3_2 = x[2]*(0.01*(x[8]-x[5])-1);

M.add(RLT3_2 <= 0);

Constraint<> RLT3_3("RLT3_3");

RLT3_3 = x[3]*(0.01*(x[8]-x[5])-1);

M.add(RLT3_3 <= 0);

Constraint<> C4("C4");

C4 = 100*x[1] - x[1]*x[6] + 833.33252*x[4];

M.add(C4 <= 83333.333);

Constraint<> C5("C5");

C5 = x[2]*x[4] - x[2]*x[7] - 1250*x[4] + 1250*x[5];

M.add(C5 <= 0);

Constraint<> C6("C6");

C6 = x[3]*x[5] - x[3]*x[8] - 2500*x[5] + 1250000;

M.add(C6 <= 0);

M.min(x[1]+x[2]+x[3]);

M.print();

auto determinant_level = 1;

bool add_Kim_Kojima = false, add_SDP_3d = false;

auto LB = M.relax(determinant_level,add_Kim_Kojima, add_SDP_3d);

LB->print();

double max_time = 54000,ub_solver_tol=1e-6, lb_solver_tol=1e-6, range_tol=1e-4, opt_rel_tol=1e-2, opt_abs_tol=1e6;

unsigned max_iter=1, nb_threads = thread::hardware_concurrency();

SolverType ub_solver_type = ipopt, lb_solver_type = ipopt;

M.run_obbt(LB, max_time, max_iter, opt_rel_tol, opt_abs_tol, nb_threads=1, ub_solver_type, lb_solver_type, ub_solver_tol, lb_solver_tol, range_tol);

LB->print_constraints_stats(1e-6);

LB->print();

}

TEST_CASE("testing projection1") {

indices buses("buses");

buses.insert("1", "2", "3", "4");

indices node_pairs("bpairs");

node_pairs.insert("1,2", "1,3", "3,4", "4,1");

Model<> Mtest("Mtest");

var<> R_Wij("R_Wij", -1, 1);

/* Imaginary part of Wij = ViVj */

var<> Im_Wij("Im_Wij", -1, 1);

var<> Wii("Wii", 0.8, 1.21);

Mtest.add(R_Wij.in(node_pairs), Im_Wij.in(node_pairs), Wii.in(buses));

Constraint<> SOC("SOC");

SOC = 2*R_Wij + pow(Im_Wij, 2) - 4*Wii.from(node_pairs);

Mtest.add(SOC.in(node_pairs) == 0);

auto subset = node_pairs.exclude("4,1");

Constraint<> PAD("PAD");

PAD = 2*R_Wij.in(subset) - Im_Wij.in(subset);

Mtest.add(PAD.in(subset)<=2);

Constraint<> PAD2("PAD2");

PAD2 = R_Wij("4,1") - 2*Im_Wij("4,1");

Mtest.add(PAD2>=1);

Mtest.min(sum(R_Wij));

Mtest.print();

CHECK(Mtest.get_nb_cons() == 8);

Mtest.project();

Mtest.print();

CHECK(Mtest.get_nb_eq() == 0);

CHECK(Mtest.get_nb_ineq() == 12);

}

TEST_CASE("testing projection2") {

indices buses("buses");

buses.insert("1", "2", "3", "4");

indices node_pairs("bpairs");

node_pairs.insert("1,2", "1,3", "3,4", "4,1");

auto subset = node_pairs.exclude("4,1");

Model<> Mtest("Mtest");

var<> R_Wij("R_Wij");

/* Imaginary part of Wij = ViVj */

var<> Im_Wij("Im_Wij", -1, 1);

var<> Wii("Wii", 0.8, 1.21);

Mtest.add(R_Wij.in(node_pairs), Im_Wij.in(node_pairs), Wii.in(buses));

Constraint<> SOC("SOC");

SOC = 2*R_Wij.in(subset) + pow(Im_Wij.in(subset), 2) - 4*Wii.from(subset);

Mtest.add(SOC.in(subset) == 0);

Constraint<> PAD("PAD");

PAD = 2*R_Wij.in(node_pairs) - Im_Wij.in(node_pairs);

Mtest.add(PAD.in(node_pairs)<=2);

Mtest.min(sum(R_Wij));

Mtest.print();

CHECK(Mtest.get_nb_cons() == 7);

Mtest.project();

Mtest.print();

CHECK(Mtest.get_nb_eq() == 0);

CHECK(Mtest.get_nb_ineq() == 4);

}

TEST_CASE("testing constants") {

constant<> c0;

CHECK(c0.is_double());

c0 = 3.5;

CHECK(c0==3.5);

CHECK(c0.is_positive());

CHECK(!c0.is_negative());

constant<Cpx> cx0;

cx0 = Cpx(-1,1);

CHECK(cx0.is_complex());

CHECK(cx0==Cpx(-1,1));

constant<Cpx> cx1;

cx1 = Cpx(-1,-2);

auto cx2 = cx0 + cx1;

cx2.print();

CHECK(cx2.is_complex());

CHECK(cx2.is_negative());

auto mag0 = sqrmag(cx2);

CHECK(std::abs(mag0.eval()-5)<1e-8);

auto ang0 = angle(cx2);

ang0.println();

CHECK(std::abs(ang0.eval()-(-2.677945045))<1e-8);

auto cx3 = conj(cx1);

CHECK(cx3.eval().real()==cx1.eval().real());

CHECK(cx3.eval().imag()==-1*cx1.eval().imag());

CHECK(real(cx3) == cx3.eval().real());

CHECK(imag(cx3) == cx3.eval().imag());

CHECK(cx3.get_dim() == 1);

try{

cx3.get_dim(2);

}

catch(invalid_argument& arg){

cout << "Error successfully caught: "<< endl;

cout << arg.what() << endl;

}

CHECK(cx3.is_number());

}

TEST_CASE("testing parameters") {

param<> c0("c0");

CHECK(c0.is_double());

c0 = 3.5;

CHECK(c0.is_positive());

CHECK(!c0.is_negative());

param<Cpx> cx0("cx0");

cx0 = Cpx(-1,1);

CHECK(cx0.is_complex());

param<Cpx> cx1("cx1");

cx1 = Cpx(-1,-2);

auto mag0 = sqrmag(cx1);

auto ang0 = ang(cx1);

mag0.print();

ang0.print();

auto cx2 = conj(cx1);

CHECK(cx2.eval().real()==cx1.eval().real());

CHECK(cx2.get_dim() == 1);

param<Cpx> cx3("cx3");

cx3.in(C(6));

CHECK(cx3.get_dim() == 6);

}

TEST_CASE("testing param copy operator") {

param<int> ip("ip");

ip.print();

ip.add_val(2);

ip.add_val(3);

ip.print();

param<int> ip2(ip);

ip2.print();

CHECK(ip==ip2);

CHECK(ip2.is_param());

CHECK(ip2.get_intype()==integer_);

CHECK(ip2.is_positive());

CHECK(ip2._range->first==2);

CHECK(ip2._range->second==3);

}

TEST_CASE("testing param indexing, add_val() and set_val() functions") {

param<> ip("ip");

ip.print();

ip.add_val(2);

ip.add_val(-1.3);

CHECK(ip._range->first==-1.3);

CHECK(ip._range->second==2);

ip.print();

ip.set_val(1, 1.5);

ip.print();

CHECK(ip.eval(1)==1.5);

CHECK(ip._range->first==1.5);

CHECK(ip._range->second==2);

indices ids("index_set");

ids.add({"id1", "id2", "key3"});

param<> dp("dp");

dp.in(ids);

dp.print();

dp("id1") = 1.5;

dp("id2") = -231.5;

dp.print();

CHECK(dp.eval("id1")==1.5);

CHECK(dp.eval("id2")==-231.5);

CHECK(dp.eval("key3")==0);

CHECK(dp._range->first==-231.5);

CHECK(dp._range->second==1.5);

REQUIRE_THROWS_AS(dp("unexisting_key").eval(), invalid_argument);

auto ndp = dp.in(ids.exclude("id2"));

ndp.print();

CHECK(ndp.get_dim()==2);

}

TEST_CASE("testing matrix params") {

param<> mat("M");

mat.set_size(3,4);

for (auto i = 0; i<3;i++) {

for (auto j = 0; j<4;j++) {

mat.set_val(i, j, 10*i+j);

}

}

mat.print();

CHECK(mat.eval(1,2)==12);

CHECK(mat(1,2).eval()==12);

auto tr_mat = mat.tr();

tr_mat.print();

CHECK(tr_mat.eval(2,1)==12);

CHECK(tr_mat(2,1).eval()==12);

var<> v("v",0,1);

v.in(R(4));

Constraint<> Mv("Mv");

Mv = product(mat,v);

Mv.print();

/* Complex matrices */

param<Cpx> Cmat("Cmat");

Cmat.set_size(3,3);

Cpx cval = Cpx(-1,2);

Cmat.set_val(0, 1, cval);

Cmat.print();

CHECK(Cmat.eval(0,1)==cval);

CHECK(Cmat(0,1).eval()==cval);

}

TEST_CASE("testing param sign functions") {

param<int> ip("ip");

ip.print();

ip.add_val(2);

ip.add_val(3);

ip.print();

CHECK(ip.is_non_negative());

CHECK(ip.is_positive());

CHECK(!ip.is_negative());

CHECK(!ip.is_non_positive());

param<> dp("dp");

dp.add_val(0);

dp.add_val(-3);

dp.print();

CHECK(!dp.is_positive());

CHECK(!dp.is_non_negative());

CHECK(!dp.is_negative());

CHECK(dp.is_non_positive());

}

TEST_CASE("testing variables indexing") {

indices ids("index_set");

ids.add({"id1", "id2", "key3"});

var<> iv("iv",-2, 5);

iv.in(ids);

iv.print();

param<Cpx> lb("lb"), ub("ub");

lb.in(ids);

lb = Cpx(0,-1);

ub.in(ids);

ub = Cpx(1,1);

var<Cpx> cv("cv", lb, ub);

cv.in(ids.exclude("id2"));

cv.print();

CHECK(cv.get_dim()==2);

var<Cpx> cv1("cv", Cpx(0,-1),Cpx(1,1));

cv1.in(ids.exclude("id2"));

cv1.print();

CHECK(cv1.get_dim()==2);

CHECK(cv==cv1);

CHECK(cv.get_indices()!=ids);

}

TEST_CASE("testing vector dot product"){

param<> a("a");

a.set_size(3);

a.set_val(0, 1);

a.set_val(1, -1);

a.set_val(2, 2);

param<int> b("b");

b=5;

CHECK(b.get_dim()==1);

CHECK(b.is_integer());

CHECK(b.eval(0)==5);

b=0;

b=2;

b=5;

b=4;

CHECK(b.get_dim()==5);

CHECK(b.eval(1)==0);

b.set_val(1, 3);

CHECK(b.eval(1)==3);

b.print();

var<> z("z",-1,1);

z.in(R(4));

CHECK(z.get_dim()==4);

z.print();

z.in(R(3));

CHECK(z.get_dim()==3);

z.print();

param<Cpx> cp("cp");

cp.in(C(3));

cp.set_val(Cpx(1,1));

auto lin = cp+z;

CHECK(!lin.is_evaluated());

lin.print_symbolic();

lin.print();

CHECK(lin._range->first==Cpx(0,1));

CHECK(lin._range->second==Cpx(2,1));

CHECK(lin.is_complex());

CHECK(lin.is_linear());

auto lin2 = cp*z;

auto lin3 = a.tr()*z - b;

CHECK(lin3.is_double());

CHECK(lin3.get_dim()==5);

lin3.print();

var<Cpx> y("y");

y.in(C(5));

auto lin4 = b*y*y;

lin4.print_symbolic();

lin4.print();

CHECK(lin4.is_complex());

CHECK(lin4.is_convex());

CHECK(lin4.get_dim()==5);

auto lin5 = lin4 - lin3;

lin5.print_symbolic();

lin5.print();

CHECK(lin5.is_complex());

CHECK(lin5.get_dim()==5);

CHECK(lin5.get_nb_vars()==4);

CHECK(lin5.is_convex());

CHECK(lin5.to_str()=="(b)y² - ([a]ᵀ)[z] + b");

auto lin6 = (2*a-exp(a)+1).tr()*z;

lin6.print_symbolic();

CHECK(lin6.to_str()=="([2a + 1 + -exp(a)]ᵀ)[z]");

CHECK(lin6.is_linear());

lin6.print();

}

TEST_CASE("testing complex numbers") {

indices ids("index_set");

ids.add({"id1", "id2", "key3", "key4"});

var<> iv("x",-2, 5);

iv.in(ids);

var<Cpx> cv("y", Cpx(0,-1),Cpx(1,1));

constant<Cpx> cx(Cpx(-1,-2));

auto cx_conj = conj(cx);

CHECK(real(cx_conj)==real(cx));

CHECK(imag(cx_conj).eval()==-1*imag(cx).eval());

param<Cpx> px("px");

px = Cpx(-1,-2);

px.print();

auto px_conj = conj(px);

CHECK(real(px_conj)._is_real);

CHECK(imag(px_conj)._is_imag);

}

TEST_CASE("testing complex functions") {

var<> z("z",-1,1);

z.in(R(3));

param<Cpx> cpx("cpx");

cpx = Cpx(1,1);

cpx = Cpx(2,2);

cpx = Cpx(3,1);

auto cpx_f = 2*exp(cpx)*z;

cpx_f.print_symbolic();

CHECK(cpx_f.to_str()=="((2,0)exp(cpx))z");

cpx_f.print();

}

TEST_CASE("testing ReLU") {

var<> A("A",-1,1), B("B",-1,1);

A.in(R(4));

B.in(R(1));

param<> X("X");

X.in(R(4));

X(0) = 0.2;

X(1) = 0.3;

X(2) = 0.4;

X(3) = -0.2;

auto f = ReLU(A.tr()*X + B);

f.print_symbolic();

f.print();

auto dfdA = f.get_derivative(B);

dfdA.print_symbolic();

CHECK(dfdA.to_str()=="UnitStep(B + (Xᵀ)[A])");

dfdA.print();

}

TEST_CASE("2d Polynomial") {

var<> x1("x1",-1,1), x2("x2",-1,1);

Model<> M("test");

M.add(x1.in(R(1)),x2.in(R(1)));

M.initialize_uniform();

M.min(pow(x1,2) - 2*pow(x1,4) + x1*x2 - 4*pow(x2,2) + 4*pow(x2,8));

solver<> s(M,ipopt);

s.set_option("max_iter", 2);

s.run(5,1e-6);

M.print();

M.print_solution();

}

TEST_CASE("testing range propagation") {

indices ids("index_set");

ids.add({"id1", "id2", "key3", "key4"});

var<> x("x",-2, 5);

x.in(ids);

auto f = 2*x + 2;

f.print_symbolic();

CHECK(f.is_linear());

CHECK(f.is_convex());

CHECK(f._range->first==-2);

CHECK(f._range->second==12);

f+= pow(x,2);

CHECK(f._range->first==-2);

CHECK(f._range->second==37);

f+=2;

CHECK(f._range->first==0);

CHECK(f._range->second==39);

CHECK(f.to_str()=="x² + 2x + 4");

f.print_symbolic();

f.print();

CHECK(f.is_quadratic());

CHECK(f.is_convex());

auto dfx = f.get_derivative(x);

dfx.print();

CHECK(dfx.is_linear());

CHECK(dfx.get_dim()==4);

param<> a("a");

a.in(R(4));

a = 2;

a.print();

var<> x1("x1",-2, 2);

x1.in(ids);

var<> x2("x2",0, 3);

x2.in(ids);

auto f2 = a*pow(x1,4) * pow(x2,2);

f2.print_symbolic();

f2.print();

CHECK(f2.is_polynomial());

CHECK(f2._range->first==0);

CHECK(f2._range->second==2*pow(2,4)*pow(3,2));

var<> x3("x3",-2, 2);

x3.in(ids);

f2 *= pow(x3,3);

f2.print_symbolic();

f2.print();

CHECK(f2._range->first==-2*pow(2,4)*pow(3,2)*pow(2,3));

CHECK(f2._range->second==2*pow(2,4)*pow(3,2)*pow(2,3));

auto dfdx1 = f2.get_derivative(x1);

dfdx1.print_symbolic();

dfdx1.print();

CHECK(dfdx1._range->first==-8*pow(2,3)*pow(3,2)*pow(2,3));

CHECK(dfdx1._range->second==8*pow(2,3)*pow(3,2)*pow(2,3));

}

TEST_CASE("testing polynomial functions") {

param<int> a("a");

a.set_size(3);

a.set_val(0,1);

a.set_val(1,-1);

a.set_val(2,2);

param<int> b("b");

b=5;

b=0;

b=2;

var<> x("x",-1,1);

x.in(R(3));