diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..c35ceb65ccfbc6a9ffcc42a3674ae0db278ab8ee --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*.swp +*.swo +.DS_Store +build \ No newline at end of file diff --git a/cpp_frep/CMakeLists.txt b/cpp_frep/CMakeLists.txt index 16627a6f1642a54631cf932fb80ae60903225d97..223013a50946bc276c84a91e7a5873fdd6152744 100755 --- a/cpp_frep/CMakeLists.txt +++ b/cpp_frep/CMakeLists.txt @@ -1,5 +1,7 @@ add_library(cpp_frep_lib frep.cpp + frep_tree.cpp + nelder_mead.h ) target_include_directories(cpp_frep_lib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) target_link_libraries(cpp_frep_lib shared_settings Eigen3::Eigen) diff --git a/cpp_frep/frep.cpp b/cpp_frep/frep.cpp index da3ac736590e6224fff19beb293e1a3fcb37608e..72476d520d02f5117795fbbb3bd2d71e0e399797 100755 --- a/cpp_frep/frep.cpp +++ b/cpp_frep/frep.cpp @@ -1,5 +1,6 @@ #include <iostream> #include <math.h> +#include <vector> namespace cpp_frep { using scalar=float; @@ -9,6 +10,7 @@ namespace cpp_frep { scalar X; scalar Y; std::string name; + std::vector<scalar> params; Frep():X(0.0f),Y(0.0f){}; Frep(std::string name_,scalar X_, scalar Y_) { @@ -28,6 +30,9 @@ namespace cpp_frep { virtual scalar eval(){ return 0.0f; } + virtual void setParameters(std::vector<scalar> params_){ + params=params_; + } // private: @@ -67,6 +72,8 @@ namespace cpp_frep { }; scalar eval(){ + // std::cout<<"circleee\n"; + return ((radius)*(radius)-((X-(center_x))*(X-(center_x))+(Y-(center_y))*(Y-(center_y)))); } @@ -112,15 +119,22 @@ namespace cpp_frep { //.......................OPERATORS..................................// class FrepOperators: public Frep { public: - + Frep* shape1; + Frep* shape2; FrepOperators(){}; + + void setShape1(Frep* shape1_){ + shape1=shape1_; + } + void setShape2(Frep* shape2_){ + shape2=shape2_; + } + }; class Add: public FrepOperators { public: - Frep* shape1; - Frep* shape2; - + Add(){}; Add(Frep* shape1_,Frep* shape2_) { shape1=shape1_; @@ -131,9 +145,10 @@ namespace cpp_frep { //change x and y shape1->setXY(X,Y); shape2->setXY(X,Y); + // std::cout<<"addddd\n"; return std::max(shape1->eval(),shape2->eval()); } - + Add get(){ return *this; } @@ -143,9 +158,6 @@ namespace cpp_frep { class Subtract: public FrepOperators { public: - Frep* shape1; - Frep* shape2; - Subtract(){}; Subtract(Frep* shape1_,Frep* shape2_) { shape1=shape1_; @@ -158,6 +170,7 @@ namespace cpp_frep { shape2->setXY(X,Y); return std::min(shape1->eval(),-shape2->eval()); } + Subtract get(){ return *this; @@ -168,9 +181,7 @@ namespace cpp_frep { class Intersect: public FrepOperators { public: - Frep* shape1; - Frep* shape2; - + Intersect(){}; Intersect(Frep* shape1_,Frep* shape2_) { shape1=shape1_; @@ -184,6 +195,7 @@ namespace cpp_frep { return std::min(shape1->eval(),shape2->eval()); } + Intersect get(){ return *this; } @@ -192,11 +204,10 @@ namespace cpp_frep { class Clearance: public FrepOperators { public: - Frep* shape1; - Frep* shape2; + scalar offset; - Clearance(){}; + Clearance():offset(0.5f){}; Clearance(Frep* shape1_,Frep* shape2_,scalar offset_) { shape1=shape1_; shape2=shape2_; @@ -209,6 +220,9 @@ namespace cpp_frep { shape2->setXY(X,Y); return std::min(shape1->eval(),-(offset+shape2->eval())); } + void setParameters(std::vector<scalar> offset_){ + offset=offset_[0]; + } Clearance get(){ return *this; @@ -217,11 +231,10 @@ namespace cpp_frep { class Blend: public FrepOperators { public: - Frep* shape1; - Frep* shape2; + scalar offset; - Blend(){}; + Blend():offset(0.5f){}; Blend(Frep* shape1_,Frep* shape2_,scalar offset_) { shape1=shape1_; shape2=shape2_; @@ -232,6 +245,11 @@ namespace cpp_frep { return *this; } + void setParameters(std::vector<scalar> offset_){ + offset=offset_[0]; + } + + float eval(){ //change x and y shape1->setXY(X,Y); @@ -244,25 +262,30 @@ namespace cpp_frep { class Morph: public FrepOperators { public: - Frep* shape1; - Frep* shape2; + scalar weight; - Morph(){}; + Morph():weight(0.5f){}; Morph(Frep* shape1_,Frep* shape2_,scalar weight_) { shape1=shape1_; shape2=shape2_; weight=weight_; }; + void setParameters(std::vector<scalar> weight_){ + weight=weight_[0]; + } + Morph get(){ return *this; } float eval(){ //change x and y + // std::cout<<"xv cv\n"; shape1->setXY(X,Y); shape2->setXY(X,Y); + // std::cout<<"moorrrphhhh\n"; return (weight*shape1->eval()+(1-weight)*shape2->eval()); } @@ -272,18 +295,22 @@ namespace cpp_frep { //.......................TRANSFORMS.................................// class FrepTransforms: public Frep { public: - + Frep* shape; //Frep(std::string const& name): name_(std::move(name)) {} FrepTransforms(){}; + + void setShape(Frep* shape_){ + shape=shape_; + } + }; class Translate: public FrepTransforms { public: - Frep* shape; scalar dx; scalar dy; - Translate():dx(0.0f),dy(0.0f){}; + Translate():dx(1.0f),dy(1.0f){}; Translate(Frep* shape_,scalar dx_,scalar dy_) { shape= shape_; dx=dx_; @@ -296,6 +323,11 @@ namespace cpp_frep { return shape->eval(); } + void setParameters(std::vector<scalar> params_){ + dx=params_[0]; + dy=params_[1]; + } + Translate get(){ return *this; } @@ -303,13 +335,12 @@ namespace cpp_frep { class Scale: public FrepTransforms { public: - Frep* shape; scalar sx; scalar sy; scalar ox; scalar oy; - Scale(){}; + Scale():sx(0.5f),sy(0.5f),ox(0.0f),oy(0.0f){}; Scale(Frep* shape_,scalar ox_, scalar oy_,scalar sx_,scalar sy_) { shape=shape_; sx=sx_; @@ -324,6 +355,13 @@ namespace cpp_frep { return shape->eval(); } + void setParameters(std::vector<scalar> params_){ + sx=params_[0]; + sy=params_[1]; + ox=params_[2]; + oy=params_[3]; + } + Scale get(){ return *this; } @@ -331,13 +369,12 @@ namespace cpp_frep { class Rotate: public FrepTransforms { public: - Frep* shape; scalar angle; scalar ox; scalar oy; const float PI_F=3.14159265358979f; - Rotate(){}; + Rotate():angle(45.0f),ox(0.0f),oy(0.0f){}; Rotate(Frep* shape_,scalar ox_, scalar oy_, scalar angle_) { shape=shape_; angle=angle_; @@ -349,6 +386,12 @@ namespace cpp_frep { return *this; } + void setParameters(std::vector<scalar> params_){ + angle=params_[0]; + ox=params_[1]; + oy=params_[2]; + } + float eval(){ //change x and y @@ -356,6 +399,7 @@ namespace cpp_frep { scalar newX=((ox)+(X-(ox))*cos(r)+(Y-(oy))*sin(r)); scalar newY=((oy)-(X-(ox))*sin(r)+(Y-(oy))*cos(r)); shape->setXY(newX,newY); + // std::cout<<"rooottaattte\n"; return shape->eval(); } @@ -364,14 +408,13 @@ namespace cpp_frep { class LinearArrayX: public FrepTransforms { public: - Frep* shape; scalar spacing; scalar number; scalar size; const float PI_F=3.14159265358979f; - LinearArrayX(){}; + LinearArrayX():spacing(0.5f),number(2.0f),size(1.0f){}; LinearArrayX(Frep* shape_,scalar spacing_,scalar number_,scalar size_) { shape=shape_; spacing=spacing_; @@ -383,6 +426,13 @@ namespace cpp_frep { return *this; } + void setParameters(std::vector<scalar> params_){ + spacing=params_[0]; + number=params_[1]; + size=params_[2]; + } + + float eval(){ //change x and y scalar xmin=-size; @@ -403,14 +453,13 @@ namespace cpp_frep { class LinearArrayY: public FrepTransforms { public: - Frep* shape; scalar spacing; scalar number; scalar size; const float PI_F=3.14159265358979f; - LinearArrayY(){}; + LinearArrayY():spacing(0.5f),number(2.0f),size(1.0f){}; LinearArrayY(Frep* shape_,scalar spacing_,scalar number_,scalar size_) { shape=shape_; spacing=spacing_; @@ -422,6 +471,12 @@ namespace cpp_frep { return *this; } + void setParameters(std::vector<scalar> params_){ + spacing=params_[0]; + number=params_[1]; + size=params_[2]; + } + float eval(){ //change x and y scalar ymin=-size; @@ -442,26 +497,33 @@ namespace cpp_frep { class PolarArray: public FrepTransforms { public: - Frep* shape; scalar angle; scalar radius; scalar ox; scalar oy; const float PI_F=3.14159265358979f; - PolarArray(){}; + PolarArray():angle(45.0f),radius(1.0f),ox(0.0f),oy(0.0f){}; PolarArray(Frep* shape_,scalar ox_, scalar oy_, scalar angle_,scalar radius_) { shape=shape_; angle=angle_; + radius=radius_; ox=ox_; oy=oy_; - radius=radius_; + }; PolarArray get(){ return *this; } + void setParameters(std::vector<scalar> params_){ + angle=params_[0]; + radius=params_[1]; + ox=params_[2]; + oy=params_[3]; + } + float eval(){ //change x and y diff --git a/cpp_frep/frep_tree.cpp b/cpp_frep/frep_tree.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4889fa1c7e26b2a14d6d555ef01982b70ff6b364 --- /dev/null +++ b/cpp_frep/frep_tree.cpp @@ -0,0 +1,491 @@ +#include <iostream> +#include <math.h> +#include "frep.cpp" +#include <vector> + +namespace cpp_frep { + using scalar=float; + class Frep_tree { + public: + // std::vector<Frep*> inputs; + int maxLength; + std::vector<std::vector<scalar>> genome; + scalar X; + scalar Y; + Frep* parsedTree; + scalar nodeLength=12.0f; + + Frep_tree():X(0.0f),Y(0.0f){}; + Frep_tree(std::vector<std::vector<scalar>> genome_) { + // inputs=inputs_; + genome=genome_; + X=0.0f; + Y=0.0f; + std::vector<scalar> g,g1; + g.push_back(10.0f/nodeLength); + g1.push_back(11.0f/nodeLength); + genome.push_back(g); + genome.push_back(g1); + maxLength = genome.size(); + parseTree(); + }; + void updateGenome(std::vector<std::vector<scalar>> genome_){ + genome=genome_; + std::vector<scalar> g,g1; + g.push_back(10.0f/nodeLength); + g1.push_back(11.0f/nodeLength); + genome.push_back(g); + genome.push_back(g1); + maxLength=genome.size(); + parseTree(); + } + + Frep_tree get(){ + return *this; + } + + void parseTree(){ + parsedTree=getNode(0); + // parsedTree=morphNodes(0); + std::cout<<"\n"; + } + + Frep* getNode(int index){ + + int nodeIndex=(int)(genome[index][0]*nodeLength); + int shape1Index,shape2Index,shapeIndex; + + // std::cout<<nodeIndex <<"\n"; + + //0 add, 1 subtract, 2 intersect, 3 morph, 4 translate, 5 scale, 6 rotate, 7 linearArrayX, 8 linearArrayX, 9 polar array + //TODO: Add parameters + switch (nodeIndex) + { + case 0: { + static Add f; + std::cout<<"Add("; + shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape1(getNode(shape1Index)); + std::cout<<","; + shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + f.setShape2(getNode(shape2Index)); + std::cout<<")"; + return &f; + break; + } + case 1: { + static Subtract f; + std::cout<<"Subtract("; + shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape1(getNode(shape1Index)); + std::cout<<","; + shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + f.setShape2(getNode(shape2Index)); + std::cout<<")"; + return &f; + break; + } + case 2: { + static Intersect f; + std::cout<<"Intersect("; + shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape1(getNode(shape1Index)); + std::cout<<","; + shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + f.setShape2(getNode(shape2Index)); + std::cout<<")"; + return &f; + break; + } + case 3: { + static Morph f; + std::cout<<"Morph("; + shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape1(getNode(shape1Index)); + std::cout<<","; + shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + f.setShape2(getNode(shape2Index)); + std::cout<<")"; + return &f; + break; + } + case 4: { + static Translate f; + std::cout<<"Translate("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(getNode(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 5: { + static Scale f; + std::cout<<"Scale("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(getNode(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 6: { + static Rotate f; + std::cout<<"Rotate("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(getNode(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 7: { + static LinearArrayX f; + std::cout<<"LinearArrayX("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(getNode(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 8: { + static LinearArrayY f; + std::cout<<"LinearArrayY("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(getNode(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 9: { + static PolarArray f; + std::cout<<"PolarArray("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(getNode(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 10: { + static Circle f(1.0f,0.0f,0.0f); + std::cout<<"Circle"; + return &f; + break; + } + case 11: { + static Rectangle f(-1.0f,1.0f,-1.0f,1.0f); + std::cout<<"Rectangle"; + return &f; + break; + } + default: { + static Rectangle f(-1.0f,1.0f,-1.0f,1.0f); + std::cout<<"rectangle1"; + return &f; + break; + } + } + // static Rectangle f1(-1.0f,1.0f,-1.0f,1.0f); + // std::cout<<"rectangle2"; + // return &f1; + } + + Frep* morphNodes(int index){ + + int nodeIndex1,nodeIndex2; + nodeIndex1=(int)(genome[index][0]*nodeLength); + nodeIndex2=nodeIndex1==nodeLength-1?10.0:nodeIndex1+1; + float weight=genome[index][0]*nodeLength-(float)nodeIndex1; + + static Morph f1; + std::cout<<"Morph("; + f1.setShape1(getNodeSharp( index ,nodeIndex1)); + std::cout<<","; + f1.setShape2(getNodeSharp(index ,nodeIndex2)); + std::cout<<")"; + f1.weight=1.0f-weight; + return &f1; + + } + + + Frep* getNodeSharp(int index, int nodeIndex){ + + int shape1Index,shape2Index,shapeIndex; + + //0 add, 1 subtract, 2 intersect, 3 morph, 4 translate, 5 scale, 6 rotate, 7 linearArrayX, 8 linearArrayX, 9 polar array + //TODO: Add parameters + switch (nodeIndex) + { + case 0: { + static Add f; + std::cout<<"Add("; + shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape1(morphNodes(genome[index][1])); + std::cout<<","; + shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + f.setShape2(morphNodes(shape2Index)); + std::cout<<")"; + return &f; + break; + } + case 1: { + static Subtract f; + std::cout<<"Subtract("; + shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape1(morphNodes(shape1Index)); + std::cout<<","; + shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + f.setShape2(morphNodes(shape2Index)); + std::cout<<")"; + return &f; + break; + } + case 2: { + static Intersect f; + std::cout<<"Intersect("; + shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape1(morphNodes(shape1Index)); + std::cout<<","; + shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + f.setShape2(morphNodes(shape2Index)); + std::cout<<")"; + return &f; + break; + } + case 3: { + static Morph f; + std::cout<<"Morph("; + shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape1(morphNodes(shape1Index)); + std::cout<<","; + shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + f.setShape2(morphNodes(shape2Index)); + std::cout<<")"; + return &f; + break; + } + case 4: { + static Translate f; + std::cout<<"Translate("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(morphNodes(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 5: { + static Scale f; + std::cout<<"Scale("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(morphNodes(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 6: { + static Rotate f; + std::cout<<"Rotate("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(morphNodes(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 7: { + static LinearArrayX f; + std::cout<<"LinearArrayX("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(morphNodes(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 8: { + static LinearArrayY f; + std::cout<<"LinearArrayY("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(morphNodes(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 9: { + static PolarArray f; + std::cout<<"PolarArray("; + shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + f.setShape(morphNodes(shapeIndex)); + std::cout<<")"; + return &f; + break; + } + case 10: { + static Circle f(1.0f,0.0f,0.0f); + std::cout<<"Circle"; + return &f; + break; + } + case 11: { + static Rectangle f(-1.0f,1.0f,-1.0f,1.0f); + std::cout<<"Rectangle"; + return &f; + break; + } + default: { + static Rectangle f(-1.0f,1.0f,-1.0f,1.0f); + std::cout<<"rectangle1"; + std::cout<<nodeIndex; + return &f; + break; + } + } + // static Rectangle f1(-1.0f,1.0f,-1.0f,1.0f); + // std::cout<<"rectangle2"; + // return &f1; + } + + // Frep* getNodeSharp(int index, int nodeIndex){ + + // int shape1Index,shape2Index,shapeIndex; + + // //0 add, 1 subtract, 2 intersect, 3 morph, 4 translate, 5 scale, 6 rotate, 7 linearArrayX, 8 linearArrayX, 9 polar array + // //TODO: Add parameters + // switch (nodeIndex) + // { + // case 0: { + // static Add f; + // std::cout<<"Add("; + // shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape1(getNode(shape1Index)); + // std::cout<<","; + // shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + // f.setShape2(getNode(shape2Index)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 1: { + // static Subtract f; + // std::cout<<"Subtract("; + // shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape1(getNode(shape1Index)); + // std::cout<<","; + // shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + // f.setShape2(getNode(shape2Index)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 2: { + // static Intersect f; + // std::cout<<"Intersect("; + // shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape1(getNode(shape1Index)); + // std::cout<<","; + // shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + // f.setShape2(getNode(shape2Index)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 3: { + // static Morph f; + // std::cout<<"Morph("; + // shape1Index=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape1(getNode(shape1Index)); + // std::cout<<","; + // shape2Index=index+(int)ceilf(genome[index][2]* (maxLength-1-index)); + // f.setShape2(getNode(shape2Index)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 4: { + // static Translate f; + // std::cout<<"Translate("; + // shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape(getNode(shapeIndex)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 5: { + // static Scale f; + // std::cout<<"Scale("; + // shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape(getNode(shapeIndex)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 6: { + // static Rotate f; + // std::cout<<"Rotate("; + // shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape(getNode(shapeIndex)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 7: { + // static LinearArrayX f; + // std::cout<<"LinearArrayX("; + // shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape(getNode(shapeIndex)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 8: { + // static LinearArrayY f; + // std::cout<<"LinearArrayY("; + // shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape(getNode(shapeIndex)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 9: { + // static PolarArray f; + // std::cout<<"PolarArray("; + // shapeIndex=index+(int)ceilf(genome[index][1]* (maxLength-1-index)); + // f.setShape(getNode(shapeIndex)); + // std::cout<<")"; + // return &f; + // break; + // } + // case 10: { + // static Circle f(1.0f,0.0f,0.0f); + // std::cout<<"Circle"; + // return &f; + // break; + // } + // case 11: { + // static Rectangle f(-1.0f,1.0f,-1.0f,1.0f); + // std::cout<<"Rectangle"; + // return &f; + // break; + // } + // default: { + // static Rectangle f(-1.0f,1.0f,-1.0f,1.0f); + // std::cout<<"rectangle1"; + // std::cout<<nodeIndex; + // return &f; + // break; + // } + // } + // // static Rectangle f1(-1.0f,1.0f,-1.0f,1.0f); + // // std::cout<<"rectangle2"; + // // return &f1; + // } + + void setXY(scalar X_,scalar Y_){ + X=X_; + Y=Y_; + } + scalar eval(){ + parsedTree->setXY(X,Y); + return parsedTree->eval(); + } + + // private: + + }; +} \ No newline at end of file diff --git a/cpp_frep/main.cpp b/cpp_frep/main.cpp index 534f16b6926b95c23d5add79fd45445817dc677a..f8428bb822e0d5310666e19ced3c6a735c300328 100755 --- a/cpp_frep/main.cpp +++ b/cpp_frep/main.cpp @@ -1,5 +1,10 @@ -#include "frep.cpp" +// #include "frep.cpp" +#include "frep_tree.cpp" +#include "nelder_mead.h" #include <iostream> +#include <vector> +#include <stdio.h> +#include <stdlib.h> using namespace cpp_frep; @@ -23,47 +28,299 @@ void drawFrep(Frep* frep,scalar minX,scalar maxX,scalar minY,scalar maxY,scalar } +void drawFrep(Frep_tree frep,scalar minX,scalar maxX,scalar minY,scalar maxY,scalar dx,scalar dy){ + for(scalar i=minX;i<maxX;i+=dx){ + for(scalar j=minY;j<maxY;j+=dy){ + frep.setXY(i,j); + scalar val=frep.eval(); + if(val>0.0f){ + std::cout << "+" ; + + }else{ + std::cout << "o" ; + } + } + std::cout << '\n'; + } + std::cout << '\n'; + +} + +void drawFrep(Frep_tree frep,scalar minX,scalar maxX,scalar minY,scalar maxY,scalar dx,scalar dy,std::string name){ + + FILE *imageFile; + int height=(int)((maxX-minX)/dx); + int width=(int)((maxY-minY)/dy); + int size=height*width*3; + // int x,y,pixel; + + std::string fileName=name+".ppm"; + const char* c = fileName.c_str(); + imageFile=fopen(c,"wb"); + if(imageFile==NULL){ + perror("ERROR: Cannot open output file"); + exit(EXIT_FAILURE); + } + + fprintf(imageFile,"P6\n"); // P6 filetype + fprintf(imageFile,"%d %d\n",width,height); // dimensions + fprintf(imageFile,"255\n"); // Max pixel + + + int stride=0; + unsigned char pix[size]; + + + for(scalar i=minX;i<maxX-dx;i+=dx){ + for(scalar j=minY;j<maxY-dy;j+=dy){ + frep.setXY(i,j); + scalar val=frep.eval(); + if(val>0.0f){ + std::cout << "+" ; + pix[stride+0]=0; + pix[stride+1]=0; + pix[stride+2]=0; + + }else{ + std::cout << "o" ; + pix[stride+0]=255; + pix[stride+1]=255; + pix[stride+2]=255; + } + stride+=3; + } + std::cout << '\n'; + } + std::cout << '\n'; + + + // unsigned char pix[]={200,200,200, 100,100,100, 0,0,0, 255,0,0, 0,255,0, 0,0,255}; + fwrite(pix,1,size,imageFile); + fclose(imageFile); + +} + +float f(Vector v) { + float x = v[0]; + float y = v[1]; + return ((-x*x*x*x+4.5*x*x+2)/pow(2.71828,2*y*y)); +} + +float f1(Vector v) { + float x = v[0]; + return -(2*x*x+6*x-3); +} + +float constraint(float num){ + num=num>1.0f?1.0f:num; + num=num<0.0001f?0.000001f:num; + return num; +} + +float compareL1(Vector v, Frep* frep,Frep_tree* frep_tree,scalar minX,scalar maxX,scalar minY,scalar maxY,scalar dx,scalar dy,std::string name){ + + FILE *imageFile; + int height=(int)((maxX-minX)/dx); + int width=(int)((maxY-minY)/dy); + int size=height*width*3; + // int x,y,pixel; + + std::string fileName=name+".ppm"; + const char* c = fileName.c_str(); + imageFile=fopen(c,"wb"); + if(imageFile==NULL){ + perror("ERROR: Cannot open output file"); + exit(EXIT_FAILURE); + } + + fprintf(imageFile,"P6\n"); // P6 filetype + fprintf(imageFile,"%d %d\n",width,height); // dimensions + fprintf(imageFile,"255\n"); // Max pixel + + + int stride=0; + unsigned char pix[size]; + int same=0; + + + std::cout<< constraint(v.at(0))<<","<<constraint(v.at(1))<<","<<constraint(v.at(2))<<"\n"; + + frep_tree->genome[0][0]=constraint(v.at(0)); + frep_tree->genome[0][1]=constraint(v.at(1)); + frep_tree->genome[0][2]=constraint(v.at(2)); + frep_tree->updateGenome(frep_tree->genome); + + + for(scalar i=minX;i<maxX-dx;i+=dx){ + for(scalar j=minY;j<maxY-dy;j+=dy){ + frep_tree->setXY(i,j); + scalar val_tree=frep_tree->eval(); + + + frep->setXY(i,j); + scalar val=frep->eval(); + + if((val>0.0&&val_tree>0.0) || (val<0.0&&val_tree<0.0)){ + same++; + } + + if(val_tree>0.0f){ + // std::cout << "+" ; + pix[stride+0]=0; + pix[stride+1]=0; + pix[stride+2]=0; + + }else{ + // std::cout << "o" ; + pix[stride+0]=255; + pix[stride+1]=255; + pix[stride+2]=255; + } + stride+=3; + } + // std::cout << '\n'; + } + // std::cout << '\n'; + float result= (float)same/((float)size/3.0f); + std::cout<<result <<"\n"; + // unsigned char pix[]={200,200,200, 100,100,100, 0,0,0, 255,0,0, 0,255,0, 0,0,255}; + fwrite(pix,1,size,imageFile); + fclose(imageFile); + return result; + +} + int main(int const, char const**) { scalar minX=-3.0f; scalar maxX=3.0f; - scalar minY=-5.0f; - scalar maxY=5.0f; + scalar minY=-3.0f; + scalar maxY=3.0f; scalar dx=0.1f; scalar dy=0.1f; + + std::vector<std::vector<scalar>> genome; - Circle c(1.0f,0.0f,0.0f); - drawFrep(&c,minX, maxX, minY, maxY, dx, dy); + std::vector<scalar> node1; + node1.push_back(0.55); + node1.push_back(0.32); + node1.push_back(0.67); + genome.push_back(node1); + + // std::vector<scalar> node2; + // node2.push_back(0.55); + // node2.push_back(0.32); + // node2.push_back(0.67); + // genome.push_back(node2); + + Frep_tree tree(genome); + // genome[0][0]=0.1; + // tree.updateGenome(genome); + // drawFrep(tree,minX, maxX, minY, maxY, dx, dy,"image"); - Translate t(&c,0.0f,1.5f); - drawFrep(&t,minX, maxX, minY, maxY, dx, dy); + + float precision = 0.01; + int dimension = 3; + NelderMeadOptimizer o(dimension, precision); - Scale sc(&t,0.0f,0.0f,2.0f,2.0f); - drawFrep(&sc,minX, maxX, minY, maxY, dx, dy); + // request a simplex to start with + Vector v(0.6,0.1,0.7667); + o.insert(v); + o.insert(Vector(0.3,0.1,0.7)); + o.insert(Vector(0.8,0.4,0.9)); + Circle c(1.0f,0.0f,0.0f); Rectangle r(-1.0f,1.0f,-1.0f,1.0f); - drawFrep(&r,minX, maxX, minY, maxY, dx, dy); + Subtract s(&r,&c); + drawFrep(&s,minX, maxX, minY, maxY, dx, dy); - Rotate rt(&r,0.0f,0.0f,15.0f); - drawFrep(&rt,minX, maxX, minY, maxY, dx, dy); + int count=0; + + float score; + while (!o.done()&& !o.done(score)) { + std::string name1 = "i"; + std::string name=name1+std::to_string(count); + score=compareL1( v, &c, &tree, minX, maxX, minY, maxY, dx, dy, name); + v = o.step(v,score ); + count++; + } - LinearArrayY la(&r,0.2f,4.0f,1.0); - drawFrep(&la,minX, maxX, minY, maxY, dx, dy); - Add a(&t,&r); - drawFrep(&a,minX, maxX, minY, maxY, dx, dy); - Subtract s(&r,&t); - drawFrep(&s,minX, maxX, minY, maxY, dx, dy); - Intersect i(&t,&r); - drawFrep(&i,minX, maxX, minY, maxY, dx, dy); + // Circle c(1.0f,0.0f,0.0f); + // drawFrep(&c,minX, maxX, minY, maxY, dx, dy); + + // Translate t(&c,0.0f,1.5f); + // drawFrep(&t,minX, maxX, minY, maxY, dx, dy); + + // Scale sc(&t,0.0f,0.0f,2.0f,2.0f); + // drawFrep(&sc,minX, maxX, minY, maxY, dx, dy); + + // Rectangle r(-1.0f,1.0f,-1.0f,1.0f); + // drawFrep(&r,minX, maxX, minY, maxY, dx, dy); + + // Rotate rt(&r,0.0f,0.0f,15.0f); + // drawFrep(&rt,minX, maxX, minY, maxY, dx, dy); + + // LinearArrayY la(&r,0.2f,4.0f,1.0); + // drawFrep(&la,minX, maxX, minY, maxY, dx, dy); + + // Add a(&t,&r); + // drawFrep(&a,minX, maxX, minY, maxY, dx, dy); + + // Subtract s(&r,&t); + // drawFrep(&s,minX, maxX, minY, maxY, dx, dy); + + // Intersect i(&t,&r); + // drawFrep(&i,minX, maxX, minY, maxY, dx, dy); + + // Clearance cl(&t,&r,1.0); + // drawFrep(&cl,minX, maxX, minY, maxY, dx, dy); + + // PolarArray pa(&r,0.0f, 0.0f, 45.0f,-1.0f); + // drawFrep(&pa,minX, maxX, minY, maxY, dx, dy); + + //0 add, 1 subtract, 2 intersect, 3 morph, 4 translate, 5 scale, 6 rotate, 7 linearArray, 8 polar array + + + // float precision = 0.01; + // int dimension = 2; + // NelderMeadOptimizer o(dimension, precision); + + // // request a simplex to start with + // Vector v(0.5, 0.5); + // o.insert(v); + // o.insert(Vector(0.1, 0.1)); + // o.insert(Vector(0.2, 0.7)); + + + // std::cout<< v.at(0)<<","<<v.at(1)<<"\n"; + + // while (!o.done()) { + // v = o.step(v, f(v)); + // std::cout<< v.at(0)<<","<<v.at(1)<<"\n"; + // std::cout<< f(v)<<"\n"; + // } + + // float precision = 0.001; + // int dimension = 1; + // NelderMeadOptimizer o(dimension, precision); + + // // request a simplex to start with + // Vector v(0.5); + // o.insert(v); + // o.insert(Vector(0.1)); + // o.insert(Vector(0.2)); + - Clearance cl(&t,&r,1.0); - drawFrep(&cl,minX, maxX, minY, maxY, dx, dy); + // std::cout<< v.at(0)<<"\n"; - PolarArray pa(&r,0.0f, 0.0f, 45.0f,-1.0f); - drawFrep(&pa,minX, maxX, minY, maxY, dx, dy); + // while (!o.done()) { + // v = o.step(v, f1(v)); + // std::cout<< f1(v)<<"\n"; + // } + // std::cout<< f1(v)<<"\n"; return 0; diff --git a/cpp_frep/nelder_mead.h b/cpp_frep/nelder_mead.h new file mode 100644 index 0000000000000000000000000000000000000000..67f79b8667204fe21a5f271a2364d6834b9d5294 --- /dev/null +++ b/cpp_frep/nelder_mead.h @@ -0,0 +1,258 @@ +//from https://github.com/blinry/nelder-mead-optimizer +#include <ctime> +#include <vector> +#include <map> +#include <cmath> +#include <algorithm> +using namespace std; + +// Float vector with standard operations +class Vector { + public: + Vector() { + } + Vector(float c0) { + coords.push_back(c0); + } + Vector(float c0, float c1) { + coords.push_back(c0); + coords.push_back(c1); + } + Vector(float c0, float c1, float c2) { + coords.push_back(c0); + coords.push_back(c1); + coords.push_back(c2); + } + + // add more constructors when N gets > 3 + + float& operator[](int i) { + return coords[i]; + } + float at(int i) const { + return coords[i]; + } + int dimension() const { + return coords.size(); + } + void prepare(int size) { + for (int i=0; i<size; i++) { + coords.push_back(0); + } + } + Vector operator+(Vector other) { + Vector result; + result.prepare(dimension()); + for (int i=0; i<dimension(); i++) { + result[i] = coords[i] + other[i]; + } + return result; + } + void operator+=(Vector other) { + for (int i=0; i<dimension(); i++) { + coords[i] += other[i]; + } + } + Vector operator-(Vector other) { + Vector result; + result.prepare(dimension()); + for (int i=0; i<dimension(); i++) { + result[i] = coords[i] - other[i]; + } + return result; + } + bool operator==(Vector other) { + if (dimension() != other.dimension()) { + return false; + } + for (int i=0; i<dimension(); i++) { + if (other[i] != coords[i]) { + return false; + } + } + return true; + } + Vector operator*(float factor) { + Vector result; + result.prepare(dimension()); + for (int i=0; i<dimension(); i++) { + result[i] = coords[i]*factor; + } + return result; + } + Vector operator/(float factor) { + Vector result; + result.prepare(dimension()); + for (int i=0; i<dimension(); i++) { + result[i] = coords[i]/factor; + } + return result; + } + void operator/=(float factor) { + for (int i=0; i<dimension(); i++) { + coords[i] /= factor; + } + } + bool operator<(const Vector other) const { + for (int i=0; i<dimension(); i++) { + if (at(i) < other.at(i)) + return false; + else if (at(i) > other.at(i)) + return true; + } + return false; + } + float length() { + float sum = 0; + for (int i=0; i<dimension(); i++) { + sum += coords[i]*coords[i]; + } + return pow(sum, 0.5f); + } + private: + vector<float> coords; +}; + +// This class stores known values for vectors. It throws unknown vectors. +class ValueDB { + public: + ValueDB() { + } + float lookup(Vector vec) { + if (!contains(vec)) { + throw vec; + } else { + return values[vec]; + } + } + void insert(Vector vec, float value) { + values[vec] = value; + } + private: + bool contains(Vector vec) { + map<Vector, float>::iterator it = values.find(vec); // TODO add tolerance + return it != values.end(); + } + map<Vector, float> values; +}; + +class NelderMeadOptimizer { + public: + NelderMeadOptimizer(int dimension, float termination_distance=0.001) { + this->dimension = dimension; + srand(time(NULL)); + alpha = 1; + gamma = 2; + rho = -0.5; + sigma = 0.5; + this->termination_distance = termination_distance; + } + // used in `step` to sort the vectors + bool operator()(const Vector& a, const Vector& b) { + return db.lookup(a) < db.lookup(b); + } + // termination criteria: each pair of vectors in the simplex has to + // have a distance of at most `termination_distance` + bool done() { + if (vectors.size() < dimension) { + return false; + } + for (int i=0; i<dimension+1; i++) { + for (int j=0; j<dimension+1; j++) { + if (i==j) continue; + if ((vectors[i]-vectors[j]).length() > termination_distance) { + return false; + } + } + } + return true; + } + + bool done(float score) { + if (score==1.0f) { + return true; + } + return false; + + } + void insert(Vector vec) { + if (vectors.size() < dimension+1) { + vectors.push_back(vec); + } + } + Vector step(Vector vec, float score) { + db.insert(vec, score); + try { + if (vectors.size() < dimension+1) { + vectors.push_back(vec); + } + + // otherwise: optimize! + if (vectors.size() == dimension+1) { + while(!done(score)) { + sort(vectors.begin(), vectors.end(), *this); + Vector cog; // center of gravity + cog.prepare(dimension); + for (int i = 1; i<=dimension; i++) { + cog += vectors[i]; + } + cog /= dimension; + Vector best = vectors[dimension]; + Vector worst = vectors[0]; + Vector second_worst = vectors[1]; + // reflect + Vector reflected = cog + (cog - worst)*alpha; + if (f(reflected) > f(second_worst) && f(reflected) < f(best)) { + vectors[0] = reflected; + } else if (f(reflected) > f(best)) { + // expand + Vector expanded = cog + (cog - worst)*gamma; + if (f(expanded) > f(reflected)) { + vectors[0] = expanded; + } else { + vectors[0] = reflected; + } + } else { + // contract + Vector contracted = cog + (cog - worst)*rho; + if (f(contracted) > f(worst)) { + vectors[0] = contracted; + } else { + for (int i=0; i<dimension; i++) { + vectors[i] = best + (vectors[i] - best)*sigma; + } + } + } + } + + // algorithm is terminating, output: simplex' center of gravity + Vector cog; + for (int i = 0; i<=dimension; i++) { + cog += vectors[i]; + } + return cog/(dimension+1); + } else { + // as long as we don't have enough vectors, request random ones, + // with coordinates between 0 and 1. If you want other start vectors, + // simply ignore these and use `step` on the vectors you want. + Vector result; + result.prepare(dimension); + for (int i = 0; i<dimension; ++i) { + result[i] = 0.001*(rand()%1000); + } + return result; + } + } catch (Vector v) { + return v; + } + } + private: + float f(Vector vec) { + return db.lookup(vec); + } + int dimension; + float alpha, gamma, rho, sigma; + float termination_distance; + vector<Vector> vectors; + ValueDB db; +}; \ No newline at end of file