diff --git a/Feasibility_Problem/src/OrderConstraints_Real.hpp b/Feasibility_Problem/src/OrderConstraints_Real.hpp index 99777ed..7793d69 100644 --- a/Feasibility_Problem/src/OrderConstraints_Real.hpp +++ b/Feasibility_Problem/src/OrderConstraints_Real.hpp @@ -54,6 +54,7 @@ If r_i is real, this simplifies to + 1/(r_i r_j* r_k) + 1/(r_i r_j* r_k*) = [4 Re(r_j) Re(r_k) + Re(r_k^2) + Im(r_k)^2] / [r_i Radius(r_j)^2 Radius(r_k)^2] + */ /// SECOND ORDER /// diff --git a/Feasibility_Problem/src/OrderConstraints_RealImag.hpp b/Feasibility_Problem/src/OrderConstraints_RealImag.hpp index 29c85b6..255d0d9 100644 --- a/Feasibility_Problem/src/OrderConstraints_RealImag.hpp +++ b/Feasibility_Problem/src/OrderConstraints_RealImag.hpp @@ -20,12 +20,41 @@ Since we have that for complex-conjugated roots 1/r + 1/r* = 2 * Re(r)/Radius(r) 0.5 = - sum_i^(S/2) 2 * Re(r_i)/Radius(r_i)^2 => 0.25 = - sum_i^(S/2) Re(r_i)/Radius(r_i)^2 + For third order, we demand 1/6 = sum_{i,j i != j}^S 1/(r_i r_j) For complex conjugated roots we have 1/(r_i r_i*) + 1/(r_j r_j*) + 1/(r_i r_j) + 1/(r_i r_j*) + 1/(r_i* r_j) + 1/(r_i* r_j*) = 1/Radius(r_i) + 1/Radius(r_j) + 4 * Re(r_i) * Re(r_j) / (Radius(r_i)^2 + Radius(r_j)^2) => 1.0/24.0 = sum_{i}^(S/2) 0.25/Radius(r_i) + sum_{j; i != j}^(S/2) Re(r_i) * Re(r_j) / (Radius(r_i)^2 + Radius(r_j)^2) +If r_i is real, this simplifies to +1/(r_i r_j) + 1/(r_i r_j*) = 2 * Re(r_j) / ( r_i Radius(r_j)^2 ) + + +For fourth order, we demand 1/24 = sum_{i,j,k i != j != k}^S 1/(r_i r_j r_k) +For complex conjugated roots we have + 1/(r_i r_i* r_j) + 1/(r_i r_i* r_j*) + 1/(r_i r_i* r_k) + 1/(r_i r_i* r_k*) ++ 1/(r_i r_j r_j*) + 1/(r_i r_j r_k) + 1/(r_i r_j r_k*) ++ 1/(r_i r_j* r_k) + 1/(r_i r_j* r_k*) + 1/(r_i r_k r_k*) ++ 1/(r_i* r_j r_j*) + 1/(r_i* r_j r_k) + 1/(r_i* r_j r_k* ++ 1/(r_i* r_j* r_k) + 1/(r_i* r_j* r_k*) + 1/(r_i* r_k r_k*) + ++ 1/(r_j r_j* r_k) + 1/(r_j r_j* r_k*) + 1/(r_j r_k r_k*) + 1/(r_j* r_k r_k*) + +The first 5 lines can be rewritten as +2 * [Re(r_i) * ( Re(r_j)^2 + 4Re(r_j)Re(r_k) + Im(r_j)^2 + Re(r_k)^2 + Im(r_k)^2) + + Re(r_k) * ( Re(r_j) * ( Re(r_j) + Re(r_k) ) + Im(r_j)^2 ) + + Re(r_j) * Im(r_k)^2] +/ ( Radius(r_i)^2 * Radius(r_j)^2 * Radius(r_k)^2 ) + +and this is what will be implemented, similar to the third order constraint. + +If r_i is real, this simplifies to + 1/(r_i r_j r_j*) + 1/(r_i r_j r_k) + 1/(r_i r_j r_k*) ++ 1/(r_i r_j* r_k) + 1/(r_i r_j* r_k*) += [4 Re(r_j) Re(r_k) + Re(r_k^2) + Im(r_k)^2] + / [r_i Radius(r_j)^2 Radius(r_k)^2] + */ /// SECOND ORDER /// diff --git a/Optimization_Problem/Roots_Real.opt b/Optimization_Problem/Roots_Real.opt index a6b74e0..d341e28 100644 --- a/Optimization_Problem/Roots_Real.opt +++ b/Optimization_Problem/Roots_Real.opt @@ -1,4 +1,4 @@ -max_iter 500 +max_iter 2000 # Objective is linear grad_f_constant yes diff --git a/Optimization_Problem/Roots_RealImag.opt b/Optimization_Problem/Roots_RealImag.opt index 80afd26..f20098b 100644 --- a/Optimization_Problem/Roots_RealImag.opt +++ b/Optimization_Problem/Roots_RealImag.opt @@ -1,4 +1,4 @@ -max_iter 392 +max_iter 2000 # Objective is linear grad_f_constant yes diff --git a/Optimization_Problem/src/Main_Roots_Real.cpp b/Optimization_Problem/src/Main_Roots_Real.cpp index 96d1bf0..12b90b5 100644 --- a/Optimization_Problem/src/Main_Roots_Real.cpp +++ b/Optimization_Problem/src/Main_Roots_Real.cpp @@ -23,8 +23,8 @@ int main(int argc, char** argv) { const int NumStagesRef = std::stoi(argv[3]); const Number dtRef = std::stod(argv[4]); - // Currently only order 1 to 3 are implemented - assert(ConsOrder == 1 || ConsOrder == 2 || ConsOrder == 3); + // Currently only order 1 to 4 are implemented + assert(ConsOrder == 1 || ConsOrder == 2 || ConsOrder == 3 || ConsOrder == 4); // Create a new instance of your nlp (use Ipopt::SmartPtr) SmartPtr mynlp; diff --git a/Optimization_Problem/src/Main_Roots_RealImag.cpp b/Optimization_Problem/src/Main_Roots_RealImag.cpp index 4634614..2e58ce4 100644 --- a/Optimization_Problem/src/Main_Roots_RealImag.cpp +++ b/Optimization_Problem/src/Main_Roots_RealImag.cpp @@ -23,8 +23,8 @@ int main(int argc, char** argv) { const int NumStagesRef = std::stoi(argv[3]); const Number dtRef = std::stod(argv[4]); - // Currently only order 1 to 3 are implemented - assert(ConsOrder == 1 || ConsOrder == 2 || ConsOrder == 3); + // Currently only order 1 to 4 are implemented + assert(ConsOrder == 1 || ConsOrder == 2 || ConsOrder == 3 || ConsOrder == 4); // Create a new instance of your nlp (use Ipopt::SmartPtr) SmartPtr mynlp; diff --git a/Optimization_Problem/src/OrderConstraints_Real.hpp b/Optimization_Problem/src/OrderConstraints_Real.hpp index e132a60..c72ae85 100644 --- a/Optimization_Problem/src/OrderConstraints_Real.hpp +++ b/Optimization_Problem/src/OrderConstraints_Real.hpp @@ -20,12 +20,41 @@ Since we have that for complex-conjugated roots 1/r + 1/r* = 2 * Re(r)/Radius(r) 0.5 = - sum_i^(S/2) 2 * Re(r_i)/Radius(r_i)^2 => 0.25 = - sum_i^(S/2) Re(r_i)/Radius(r_i)^2 + For third order, we demand 1/6 = sum_{i,j i != j}^S 1/(r_i r_j) For complex conjugated roots we have 1/(r_i r_i*) + 1/(r_j r_j*) + 1/(r_i r_j) + 1/(r_i r_j*) + 1/(r_i* r_j) + 1/(r_i* r_j*) = 1/Radius(r_i) + 1/Radius(r_j) + 4 * Re(r_i) * Re(r_j) / (Radius(r_i)^2 + Radius(r_j)^2) => 1.0/24.0 = sum_{i}^(S/2) 0.25/Radius(r_i) + sum_{j; i != j}^(S/2) Re(r_i) * Re(r_j) / (Radius(r_i)^2 + Radius(r_j)^2) +If r_i is real, this simplifies to +1/(r_i r_j) + 1/(r_i r_j*) = 2 * Re(r_j) / ( r_i Radius(r_j)^2 ) + + +For fourth order, we demand 1/24 = sum_{i,j,k i != j != k}^S 1/(r_i r_j r_k) +For complex conjugated roots we have + 1/(r_i r_i* r_j) + 1/(r_i r_i* r_j*) + 1/(r_i r_i* r_k) + 1/(r_i r_i* r_k*) ++ 1/(r_i r_j r_j*) + 1/(r_i r_j r_k) + 1/(r_i r_j r_k*) ++ 1/(r_i r_j* r_k) + 1/(r_i r_j* r_k*) + 1/(r_i r_k r_k*) ++ 1/(r_i* r_j r_j*) + 1/(r_i* r_j r_k) + 1/(r_i* r_j r_k* ++ 1/(r_i* r_j* r_k) + 1/(r_i* r_j* r_k*) + 1/(r_i* r_k r_k*) + ++ 1/(r_j r_j* r_k) + 1/(r_j r_j* r_k*) + 1/(r_j r_k r_k*) + 1/(r_j* r_k r_k*) + +The first 5 lines can be rewritten as +2 * [Re(r_i) * ( Re(r_j)^2 + 4Re(r_j)Re(r_k) + Im(r_j)^2 + Re(r_k)^2 + Im(r_k)^2) + + Re(r_k) * ( Re(r_j) * ( Re(r_j) + Re(r_k) ) + Im(r_j)^2 ) + + Re(r_j) * Im(r_k)^2] +/ ( Radius(r_i)^2 * Radius(r_j)^2 * Radius(r_k)^2 ) + +and this is what will be implemented, similar to the third order constraint. + +If r_i is real, this simplifies to + 1/(r_i r_j r_j*) + 1/(r_i r_j r_k) + 1/(r_i r_j r_k*) ++ 1/(r_i r_j* r_k) + 1/(r_i r_j* r_k*) += [4 Re(r_j) Re(r_k) + Re(r_k^2) + Im(r_k)^2] + / [r_i Radius(r_j)^2 Radius(r_k)^2] + */ /// SECOND ORDER /// @@ -442,4 +471,414 @@ T ThirdOrder(const std::vector& x, const int NumRoots, return g; } + +/// FOURTH ORDER /// + + +// For Odd Base Polynom => Even Lower Degree Polynomial +template +void FourthOrder(const T* x, T* g, const int NumRoots, const int NumEigVals, + const std::vector& RealRange, const std::vector& ImagRange) +{ + g[NumEigVals+2] = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + for(size_t j = 0; j < NumRoots; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } +} + +// For Even Base Polynomial => Odd Lower degree polynomial +template +void FourthOrder(const T* x, T* g, const int NumRoots, const int NumEigVals, + const std::vector& RealRange, const std::vector& ImagRange, + const std::vector& ImagDiff_over_RealDiff, const size_t i_min) +{ + g[NumEigVals+2] = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + // Handle combination: Real with complex conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*x[j]*x[k] + x[k]*x[k] + b2*b2)/(x[i_min] * Radius1 * Radius2); + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*x[j]*x[k] + x[k]*x[k] + b2*b2)/(x[i_min] * Radius1 * Radius2); + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*x[j]*x[k] + x[k]*x[k] + b2*b2)/(x[i_min] * Radius1 * Radius2); + } + } + + // Now: Complex Conjugated <-> Complex Conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = l+1; l < i_min; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + + for(size_t l = i_min+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } +} + +// For Odd Base Polynom => Even Lower Degree Polynomial +template // T: for dco types PT: Passive Type: For usual real types +void FourthOrder(const std::vector& x, std::vector& g, const int NumRoots, const int NumEigVals, + const std::vector& RealRange, const std::vector& ImagRange) +{ + g[NumEigVals+2] = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + for(size_t j = 0; j < NumRoots; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } +} + +// For Even Base Polynom => Odd Lower Degree Polynomial +template // T: for dco types PT: Passive Type: For usual real types +void FourthOrder(const std::vector& x, std::vector& g, const int NumRoots, const int NumEigVals, + const std::vector& RealRange, const std::vector& ImagRange, + const std::vector& ImagDiff_over_RealDiff, const size_t i_min) +{ + g[NumEigVals+2] = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + // Handle combination: Real with complex conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*x[j]*x[k] + x[k]*x[k] + b2*b2)/(x[i_min] * Radius1 * Radius2); + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*x[j]*x[k] + x[k]*x[k] + b2*b2)/(x[i_min] * Radius1 * Radius2); + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*x[j]*x[k] + x[k]*x[k] + b2*b2)/(x[i_min] * Radius1 * Radius2); + } + } + + // Now: Complex Conjugated <-> Complex Conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = l+1; l < i_min; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + + for(size_t l = i_min+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g[NumEigVals+2] += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } +} + +// For Odd Base Polynom => Even Lower Degree Polynomial +template // T: for dco types PT: Passive Type: For usual real types +T FourthOrder(const std::vector& x, const int NumRoots, + const std::vector& RealRange, const std::vector& ImagRange) +{ + T g = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + for(size_t j = 0; j < NumRoots; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange); + Radius3 = x[l]*x[l] + b3*b3; + + g += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + return g; +} + +// For Even Base Polynomial => Odd Lower degree polynomial +template // T: for dco types PT: Passive Type: For usual real types +T FourthOrder(const std::vector& x, const int NumRoots, + const std::vector& RealRange, const std::vector& ImagRange, + const std::vector& ImagDiff_over_RealDiff, const size_t i_min) +{ + T g = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + // Handle combination: Real with complex conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + g += 0.5 * (4.0*x[j]*x[k] + x[k]*x[k] + b2*b2)/(x[i_min] * Radius1 * Radius2); + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + g += 0.5 * (4.0*x[j]*x[k] + x[k]*x[k] + b2*b2)/(x[i_min] * Radius1 * Radius2); + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + g += 0.5 * (4.0*x[j]*x[k] + x[k]*x[k] + b2*b2)/(x[i_min] * Radius1 * Radius2); + } + } + + // Now: Complex Conjugated <-> Complex Conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = l+1; l < i_min; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + + for(size_t l = i_min+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius1 = x[j]*x[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius2 = x[k]*x[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(x[l], RealRange, ImagRange, ImagDiff_over_RealDiff); + Radius3 = x[l]*x[l] + b3*b3; + + g += (x[j] * (b2*b2 + 4 * x[j] * x[k] + x[l]*x[l] + b3*b3) + + x[l] * (x[k] * (x[k] + x[l]) + b2*b2) + + x[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + return g; +} + #endif \ No newline at end of file diff --git a/Optimization_Problem/src/OrderConstraints_RealImag.hpp b/Optimization_Problem/src/OrderConstraints_RealImag.hpp index 362a3bd..c6617ba 100644 --- a/Optimization_Problem/src/OrderConstraints_RealImag.hpp +++ b/Optimization_Problem/src/OrderConstraints_RealImag.hpp @@ -20,12 +20,41 @@ Since we have that for complex-conjugated roots 1/r + 1/r* = 2 * Re(r)/Radius(r) 0.5 = - sum_i^(S/2) 2 * Re(r_i)/Radius(r_i)^2 => 0.25 = - sum_i^(S/2) Re(r_i)/Radius(r_i)^2 + For third order, we demand 1/6 = sum_{i,j i != j}^S 1/(r_i r_j) For complex conjugated roots we have 1/(r_i r_i*) + 1/(r_j r_j*) + 1/(r_i r_j) + 1/(r_i r_j*) + 1/(r_i* r_j) + 1/(r_i* r_j*) = 1/Radius(r_i) + 1/Radius(r_j) + 4 * Re(r_i) * Re(r_j) / (Radius(r_i)^2 + Radius(r_j)^2) => 1.0/24.0 = sum_{i}^(S/2) 0.25/Radius(r_i) + sum_{j; i != j}^(S/2) Re(r_i) * Re(r_j) / (Radius(r_i)^2 + Radius(r_j)^2) +If r_i is real, this simplifies to +1/(r_i r_j) + 1/(r_i r_j*) = 2 * Re(r_j) / ( r_i Radius(r_j)^2 ) + + +For fourth order, we demand 1/24 = sum_{i,j,k i != j != k}^S 1/(r_i r_j r_k) +For complex conjugated roots we have + 1/(r_i r_i* r_j) + 1/(r_i r_i* r_j*) + 1/(r_i r_i* r_k) + 1/(r_i r_i* r_k*) ++ 1/(r_i r_j r_j*) + 1/(r_i r_j r_k) + 1/(r_i r_j r_k*) ++ 1/(r_i r_j* r_k) + 1/(r_i r_j* r_k*) + 1/(r_i r_k r_k*) ++ 1/(r_i* r_j r_j*) + 1/(r_i* r_j r_k) + 1/(r_i* r_j r_k* ++ 1/(r_i* r_j* r_k) + 1/(r_i* r_j* r_k*) + 1/(r_i* r_k r_k*) + ++ 1/(r_j r_j* r_k) + 1/(r_j r_j* r_k*) + 1/(r_j r_k r_k*) + 1/(r_j* r_k r_k*) + +The first 5 lines can be rewritten as +2 * [Re(r_i) * ( Re(r_j)^2 + 4Re(r_j)Re(r_k) + Im(r_j)^2 + Re(r_k)^2 + Im(r_k)^2) + + Re(r_k) * ( Re(r_j) * ( Re(r_j) + Re(r_k) ) + Im(r_j)^2 ) + + Re(r_j) * Im(r_k)^2] +/ ( Radius(r_i)^2 * Radius(r_j)^2 * Radius(r_k)^2 ) + +and this is what will be implemented, similar to the third order constraint. + +If r_i is real, this simplifies to + 1/(r_i r_j r_j*) + 1/(r_i r_j r_k) + 1/(r_i r_j r_k*) ++ 1/(r_i r_j* r_k) + 1/(r_i r_j* r_k*) += [4 Re(r_j) Re(r_k) + Re(r_k^2) + Im(r_k)^2] + / [r_i Radius(r_j)^2 Radius(r_k)^2] + */ /// SECOND ORDER /// @@ -443,4 +472,414 @@ T ThirdOrder(const std::vector& xy, const int NumRoots, return g; } + +/// FOURTH ORDER /// + + +// For Odd Base Polynom => Even Lower Degree Polynomial +template +void FourthOrder(const T* xy, T* g, const int NumRoots, const int NumEigVals, + const std::vector& RealRange, const std::vector& ImagRange) +{ + g[NumEigVals+2] = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + for(size_t j = 0; j < NumRoots; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } +} + +// For Even Base Polynomial => Odd Lower degree polynomial +template +void FourthOrder(const T* xy, T* g, const int NumRoots, const int NumEigVals, + const std::vector& RealRange, const std::vector& ImagRange, + const std::vector& ImagDiff_over_RealDiff, const size_t i_min) +{ + g[NumEigVals+2] = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + // Handle combination: Real with complex conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*xy[j]*xy[k] + xy[k]*xy[k] + b2*b2)/(xy[i_min] * Radius1 * Radius2); + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*xy[j]*xy[k] + xy[k]*xy[k] + b2*b2)/(xy[i_min] * Radius1 * Radius2); + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*xy[j]*xy[k] + xy[k]*xy[k] + b2*b2)/(xy[i_min] * Radius1 * Radius2); + } + } + + // Now: Complex Conjugated <-> Complex Conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = l+1; l < i_min; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + + for(size_t l = i_min+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } +} + +// For Odd Base Polynom => Even Lower Degree Polynomial +template // T: for dco types PT: Passive Type: For usual real types +void FourthOrder(const std::vector& xy, std::vector& g, const int NumRoots, const int NumEigVals, + const std::vector& RealRange, const std::vector& ImagRange) +{ + g[NumEigVals+2] = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + for(size_t j = 0; j < NumRoots; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } +} + +// For Even Base Polynom => Odd Lower Degree Polynomial +template // T: for dco types PT: Passive Type: For usual real types +void FourthOrder(const std::vector& xy, std::vector& g, const int NumRoots, const int NumEigVals, + const std::vector& RealRange, const std::vector& ImagRange, + const std::vector& ImagDiff_over_RealDiff, const size_t i_min) +{ + g[NumEigVals+2] = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + // Handle combination: Real with complex conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*xy[j]*xy[k] + xy[k]*xy[k] + b2*b2)/(xy[i_min] * Radius1 * Radius2); + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*xy[j]*xy[k] + xy[k]*xy[k] + b2*b2)/(xy[i_min] * Radius1 * Radius2); + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + g[NumEigVals+2] += 0.5 * (4.0*xy[j]*xy[k] + xy[k]*xy[k] + b2*b2)/(xy[i_min] * Radius1 * Radius2); + } + } + + // Now: Complex Conjugated <-> Complex Conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = l+1; l < i_min; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + + for(size_t l = i_min+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g[NumEigVals+2] += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } +} + +// For Odd Base Polynom => Even Lower Degree Polynomial +template // T: for dco types PT: Passive Type: For usual real types +T FourthOrder(const std::vector& xy, const int NumRoots, + const std::vector& RealRange, const std::vector& ImagRange) +{ + T g = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + for(size_t j = 0; j < NumRoots; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + return g; +} + +// For Even Base Polynomial => Odd Lower degree polynomial +template // T: for dco types PT: Passive Type: For usual real types +T FourthOrder(const std::vector& xy, const int NumRoots, + const std::vector& RealRange, const std::vector& ImagRange, + const std::vector& ImagDiff_over_RealDiff, const size_t i_min) +{ + T g = 1.0/48.0; + + T b1, b2, b3, Radius1, Radius2, Radius3; + + // Handle combination: Real with complex conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + g += 0.5 * (4.0*xy[j]*xy[k] + xy[k]*xy[k] + b2*b2)/(xy[i_min] * Radius1 * Radius2); + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + g += 0.5 * (4.0*xy[j]*xy[k] + xy[k]*xy[k] + b2*b2)/(xy[i_min] * Radius1 * Radius2); + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + g += 0.5 * (4.0*xy[j]*xy[k] + xy[k]*xy[k] + b2*b2)/(xy[i_min] * Radius1 * Radius2); + } + } + + // Now: Complex Conjugated <-> Complex Conjugated + for(size_t j = 0; j < i_min; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < i_min; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = l+1; l < i_min; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + + for(size_t l = i_min+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + + for(size_t k = i_min+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + for(size_t j = i_min + 1; j < NumRoots; j++) { + b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots]; + Radius1 = xy[j]*xy[j] + b1*b1; + + for(size_t k = j+1; k < NumRoots; k++) { + b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots]; + Radius2 = xy[k]*xy[k] + b2*b2; + + for(size_t l = k+1; l < NumRoots; l++) { + b3 = Lin_IntPol(xy[l], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[l + NumRoots]; + Radius3 = xy[l]*xy[l] + b3*b3; + + g += (xy[j] * (b2*b2 + 4 * xy[j] * xy[k] + xy[l]*xy[l] + b3*b3) + + xy[l] * (xy[k] * (xy[k] + xy[l]) + b2*b2) + + xy[k] * b3*b3) / (Radius1 * Radius2 * Radius3); + } + } + } + + return g; +} + #endif \ No newline at end of file diff --git a/Optimization_Problem/src/Roots_Real.cpp b/Optimization_Problem/src/Roots_Real.cpp index 9ffcea4..324c3af 100644 --- a/Optimization_Problem/src/Roots_Real.cpp +++ b/Optimization_Problem/src/Roots_Real.cpp @@ -392,6 +392,13 @@ bool Roots_Real::eval_g( ThirdOrder(x, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled); else ThirdOrder(x, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled); + + if(ConsOrder == 4){ + if(UseHull) + FourthOrder(x, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled); + else + FourthOrder(x, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled); + } } } } @@ -420,6 +427,13 @@ bool Roots_Real::eval_g( ThirdOrder(x, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled, ImagDiff_over_RealDiff, i_min); else ThirdOrder(x, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled, ImagDiff_over_RealDiff, i_min); + + if(ConsOrder == 4) { + if(UseHull) + FourthOrder(x, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled, ImagDiff_over_RealDiff, i_min); + else + FourthOrder(x, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled, ImagDiff_over_RealDiff, i_min); + } } } } @@ -573,6 +587,36 @@ bool Roots_Real::eval_jac_g( values[ind] = dco::derivative(x_dco[j]); ind++; } + + if(ConsOrder == 4) { + DCO_M::global_tape->zero_adjoints(); + + if(OddDegree) { + if(UseHull) + FourthOrder(x_dco, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled); + else + FourthOrder(x_dco, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled); + } + else { + if(UseHull) + FourthOrder(x_dco, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled, + ImagDiff_over_RealDiff, i_min); + else + FourthOrder(x_dco, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled, + ImagDiff_over_RealDiff, i_min); + + } + + dco::derivative(g)[NumEigVals + 2] = 1.; // Seed component + + DCO_M::global_tape->interpret_adjoint(); // Interpret (stored) tape + + // Harvest + for (size_t j = 0; j < NumUnknowns; j++) { + values[ind] = dco::derivative(x_dco[j]); + ind++; + } + } } } @@ -826,6 +870,46 @@ bool Roots_Real::eval_h( DCO_BM::global_tape->zero_adjoints(); // Unseed } + + if(ConsOrder == 4) { + for (size_t j = 0; j < NumUnknowns; j++) { + DCO_BM::global_tape->register_variable(dco::value(x_dco[j]) ); // record active input + DCO_BM::global_tape->register_variable(dco::derivative(x_dco[j]) ); // record active input + + dco::value(dco::value(x_dco[j]) ) = x[j]; + + DCO_M::global_tape->register_variable(x_dco[j]); + } + + if(OddDegree) { + if(UseHull) + g = FourthOrder(x_dco, NumRoots, HullRealScaled, HullImagScaled); + else + g = FourthOrder(x_dco, NumRoots, RealEigValsScaled, ImagEigValsScaled); + } + else { + if(UseHull) + g = FourthOrder(x_dco, NumRoots, HullRealScaled, HullImagScaled, ImagDiff_over_RealDiff, i_min); + else + g = FourthOrder(x_dco, NumRoots, RealEigValsScaled, ImagEigValsScaled, ImagDiff_over_RealDiff, i_min); + } + + dco::value(dco::derivative(g) ) = 1.; // Seed + DCO_M::global_tape->interpret_adjoint(); // Back-propagate from output/adjoint + + ind = 0; + for(size_t k = 0; k < NumUnknowns; k++) { + dco::derivative(dco::derivative(x_dco[k]) ) = 1.; // Seed + + DCO_BM::global_tape->interpret_adjoint(); // Back-propagate from output/adjoint + for(size_t j = 0; j <= k; j++) { + values[ind] += lambda[NumEigVals + 2] * dco::derivative(dco::value(x_dco[j]) ); + ind++; + } + + DCO_BM::global_tape->zero_adjoints(); // Unseed + } + } } } @@ -1000,6 +1084,24 @@ void Roots_Real::finalize_solution( ImagDiff_over_RealDiff, i_min); } std::cout << "Final value of 3rd order constraint: " << Constr[NumEigVals+1] << std::endl; + + if(ConsOrder == 4) { + if(OddDegree) { + if(UseHull) + FourthOrder(xMaxdt, Constr, NumRoots, NumEigVals, HullRealScaled, HullImagScaled); + else + FourthOrder(xMaxdt, Constr, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled); + } + else { + if(UseHull) + FourthOrder(xMaxdt, Constr, NumRoots, NumEigVals, HullRealScaled, HullImagScaled, + ImagDiff_over_RealDiff, i_min); + else + FourthOrder(xMaxdt, Constr, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled, + ImagDiff_over_RealDiff, i_min); + } + std::cout << "Final value of 4th order constraint: " << Constr[NumEigVals+2] << std::endl; + } } } diff --git a/Optimization_Problem/src/Roots_RealImag.cpp b/Optimization_Problem/src/Roots_RealImag.cpp index 11d275b..bc64793 100644 --- a/Optimization_Problem/src/Roots_RealImag.cpp +++ b/Optimization_Problem/src/Roots_RealImag.cpp @@ -350,6 +350,13 @@ bool Roots_RealImag::eval_g( ThirdOrder(xy, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled); else ThirdOrder(xy, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled); + + if(ConsOrder == 4) { + if(UseHull) + FourthOrder(xy, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled); + else + FourthOrder(xy, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled); + } } } } @@ -382,6 +389,15 @@ bool Roots_RealImag::eval_g( else ThirdOrder(xy, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled, ImagDiff_over_RealDiff, i_min); + + if(ConsOrder == 4) { + if(UseHull) + FourthOrder(xy, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled, + ImagDiff_over_RealDiff, i_min); + else + FourthOrder(xy, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled, + ImagDiff_over_RealDiff, i_min); + } } } } @@ -533,6 +549,36 @@ bool Roots_RealImag::eval_jac_g( values[ind] = dco::derivative(xy_dco[j]); ind++; } + + if(ConsOrder == 4){ + // Reset adjoints + DCO_M::global_tape->zero_adjoints(); + + if(OddDegree) { + if(UseHull) + FourthOrder(xy_dco, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled); + else + FourthOrder(xy_dco, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled); + } + else { + if(UseHull) + FourthOrder(xy_dco, g, NumRoots, NumEigVals, HullRealScaled, HullImagScaled, + ImagDiff_over_RealDiff, i_min); + else + FourthOrder(xy_dco, g, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled, + ImagDiff_over_RealDiff, i_min); + } + + dco::derivative(g)[NumEigVals+2] = 1.; // Seed component + + DCO_M::global_tape->interpret_adjoint(); // Interpret (stored) tape + + // Harvest + for (size_t j = 0; j < NumUnknowns; j++) { + values[ind] = dco::derivative(xy_dco[j]); + ind++; + } + } } } @@ -792,6 +838,48 @@ bool Roots_RealImag::eval_h( DCO_BM::global_tape->zero_adjoints(); // Unseed } + + if(ConsOrder == 4){ + for (size_t j = 0; j < NumUnknowns; j++) { + DCO_BM::global_tape->register_variable(dco::value(xy_dco[j]) ); // record active input + DCO_BM::global_tape->register_variable(dco::derivative(xy_dco[j]) ); // record active input + + dco::value(dco::value(xy_dco[j]) ) = xy[j]; + + DCO_M::global_tape->register_variable(xy_dco[j]); + } + + if(OddDegree) { + if(UseHull) + g = FourthOrder(xy_dco, NumRoots, HullRealScaled, HullImagScaled); + else + g = FourthOrder(xy_dco, NumRoots, RealEigValsScaled, ImagEigValsScaled); + } + else { + if(UseHull) + g = FourthOrder(xy_dco, NumRoots, HullRealScaled, HullImagScaled, + ImagDiff_over_RealDiff, i_min); + else + g = FourthOrder(xy_dco, NumRoots, RealEigValsScaled, ImagEigValsScaled, + ImagDiff_over_RealDiff, i_min); + } + + dco::value(dco::derivative(g) ) = 1.; // Seed + DCO_M::global_tape->interpret_adjoint(); // Back-propagate from output/adjoint + + ind = 0; + for(size_t k = 0; k < NumUnknowns; k++) { + dco::derivative(dco::derivative(xy_dco[k]) ) = 1.; // Seed + + DCO_BM::global_tape->interpret_adjoint(); // Back-propagate from output/adjoint + for(size_t j = 0; j <= k; j++) { + values[ind] += lambda[NumEigVals+2] * dco::derivative(dco::value(xy_dco[j]) ); + ind++; + } + + DCO_BM::global_tape->zero_adjoints(); // Unseed + } + } } } @@ -954,7 +1042,7 @@ void Roots_RealImag::finalize_solution( } std::cout << std::endl << "Final value of 2nd order constraint: " << Constr[NumEigVals] << std::endl; - if(ConsOrder >= 3) { + if(ConsOrder >= 3) { if(OddDegree) { if(UseHull) ThirdOrder(xy, Constr, NumRoots, NumEigVals, HullRealScaled, HullImagScaled); @@ -970,6 +1058,24 @@ void Roots_RealImag::finalize_solution( ImagDiff_over_RealDiff, i_min); } std::cout << "Final value of 3rd order constraint: " << Constr[NumEigVals+1] << std::endl; + + if(ConsOrder == 4){ + if(OddDegree) { + if(UseHull) + FourthOrder(xy, Constr, NumRoots, NumEigVals, HullRealScaled, HullImagScaled); + else + FourthOrder(xy, Constr, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled); + } + else { + if(UseHull) + FourthOrder(xy, Constr, NumRoots, NumEigVals, HullRealScaled, HullImagScaled, + ImagDiff_over_RealDiff, i_min); + else + FourthOrder(xy, Constr, NumRoots, NumEigVals, RealEigValsScaled, ImagEigValsScaled, + ImagDiff_over_RealDiff, i_min); + } + std::cout << "Final value of 4th order constraint: " << Constr[NumEigVals+2] << std::endl; + } } }