Trig Heuristic Method - Alpha V1.1 - C++ version - 11/18/15

Under The Artistic License 2.0, https://opensource.org/licenses/Artistic-2.0

This alpha version enables an increased number of contraints. If more than 100 constraints are required, please increase the array bound definition.

//Project starts here: special case when objective variable distance to the origin intersects the k-th constraints is not yet address here
//Two variable optimization project, which includes more than three constrains
#include <iostream> // std::cout, std::cin
#include <cmath> // std:: sqrt, std::pow, std::abs
using namespace std;
int i=0, j=0; // incremental operators
int a=0; // incremental operator in an array that is used to record the smallest angle for the shortest distance gamma, among equally small angles
int b=0; // incremental operator used to organize constraints bound by the k-th constraint; these constraints are below a threshold that requires them to considered in the set of primary constraints
int f=0; // incremental operators
int r=0; // operator set to sum the number of primary reference constraints, which excludes erased secondary set of reference constraints, and system constraints, which are established based on the k-th bounds
int q=0; // operator used to identify second reference constraint after the k-th one; this constraint is included in a primary set of constraint that are not bounded by the k-th constraint
int k=0; // operator used to identify first reference constraint, which yields the shortest distance gamma and smallest angle.
int kk=0; //resequenced k-th constraint
long double p=0; // operator used to subtract one array from another, and operator used towards computing the power function
int num=2; // number of variables in the objective function
int t=6; // number of constraints
int u=0; // increment operator used to assign value to array beyond the 6th referenced constraint
int g=0; // operator used to replace the primary set of reference constraints and erase secondary set of constraints
long double newconstr[100][2]; // matrix containing six of the original constraints and two project generated constraints
long double formatconstr[100][2]; // matrix containing set of primary reference constraints, which excludes erased secondary set of reference constraints
 
long double epi[2][2]; // matrix inverse of the matrix containing two constraints meant to be used for obtaining the quantity for each variable
long double x=0; // operator used to sum the squares, to measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
int y=0; // identifier when the angle for the k-th constraint is null
int s=0; // identifier for goto loop, when angle for th k-th constraint is null
long double z=0; // operator used to sum the squares, to measure distance between variable coefficients and origin
int LC[100] = {}; 
long double LK[100] = {}; // measure distance between the k-th constraint and second likely constraint
long double L7[100] = {};  // array used to store angles
long double LA[100] = {};  // measure distance between constraint coordinates and origin; there are six constraints specified by the user and two system generated
// array used to store the distance between system constraint coordinates and origin; there are six constraints specified by the user and two system generated
long double L6[100] = {}; // measure distance between constraint coordinates and variable coefficients; there are six constraints specified by the user and two system generated
// array used to store the distance between system constraint coordinates and variable coefficients; there are six constraints specified by the user and two system generated
long double L8[100] = {}; // compute cosine, and thereafter record the angle, value of the triangle formed between the variable coefficients, a given constraint, and the origin
long double L9[100] = {};  // multiply distance between the variable coefficients and a given constraints and the cosine value; this provides distance between the two points if they were aligned between the variable coefficients and the origin
long double LB[100] = {};  // array that includes all constraints that will be excluded from the set of primary constraints, since the array includes all constraints bound by the k-th constraint below a target threshold
long double LE[100] = {};  // array used record the angle of the constraint yielding the smallest distance
int LF[100] = {};  // array used store which constraint yields the smallest distance
// array used to store the distance between system constraint coordinates and k-th constraint
long double w=0, v=0; // angles recorded from triangles formed, which are used to determine which plane formed by the binding a system constraint and the k-th constraint intersects the line formed by binding the coefficient variables and the origin
long double sol=0; // operator used to sum parameters for each constraint; this is equivalent to a sumproduct of the constraint with a unity vector
long double SOLL[2]; // array used to store the solution set for the maximization project
long double N[2] = { 18.3333333333333333, 20 }; // array of coefficients for variables in the objective function, which will be used in the maximization project within the prescribed constraint
// parameters for each of the six constraints are assigned to the variable to be optimized
long double constraints[][2]= {
{1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}, {15, 12}, {12, 15},
{13, 16}, {14, 17}, {15, 18}, {16, 19}, {17, 20}, {18, 21}, {19, 22},
{20, 23}
};
// limiting parameters for each of the six constraints
long double limit[]= {
10, 20, 30, 40, 50, 60, 100, 100, 100, 100, 100, 100, 120, 120, 120, 120
};
//list of functions
void cosine_and_angle(long double array1[], long double array2[], int upper_bound); //compute cosine value and angle of the triangle formed between the variable coefficients, a given constraint, and the origin
void min(long double array3[], int upper_bound); //function used to determine the minimum within an array and below an upper bound
void matrix_inverse(long double array4[2][2], int amount6, int amoutn7); //function used to compute the inverse of a 2 x 2 matrix
long double cosine(long double amount1, long double amount2, long double amount3); //function used to compute the cosine 
void dist_constr_var (long double array5[][2], int amount5, int amount4, long double return_array[], int upper_bound, int increment); //measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
//int main begins
int main () 
{
    
    // defined coefficients are provided for each variable of the objective function for the maximization problem
    for(i=0;i<t;i++){
        for(j=0;j<num;j++){
            newconstr[i][j] = constraints[i][j]/ limit[i];
        }
    }
    
    //measure distance between constraint coordinates and origin, by taking the square root for the sum of squares
    dist_constr_var (newconstr, 0, 0, LA, t, 0) ;
    
    //measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
    dist_constr_var (newconstr, 0, 1, L6, t, 0) ;
    //measure distance between variable coefficients and origin, by taking the square root for the sum of squares
    for(j=0;j<num;j++){
        z += pow(N[j],2);
    }
    z = sqrt(z);   //measure distance between variable coefficients and origin
        
    //compute cosine value and angle of the triangle formed between the variable coefficients, a given constraint, and the origin    
    cosine_and_angle(LA,L6,t) ; 
    
    //find the smallest distance gamma recorded in array L9
    min(L9, t);
    //record all constraints for which the smallest distance gamma is equal
    for(i=0;i<t;i++){
        if (L9[i] == x) {

LF[a]=i; //array used record which constraint yields the smallest distance gamma

            LE[a]=L8[i]; // array used record the angle of the constraint yielding the smallest distance gamma
            a++;
        }
    }
    //record the smallest angle for the shortest distance gamma, which is recorded in the set of ‘’a’’ for equally small angles
    min(LE, a);
    
    //record the k-th constraint that yields the shortest distance gamma and angle, which is the first reference constraint
    for(i=0;i<a;i++){
        if (LE[i] == x) {
            k = LF[i];
            x=0; 
        }
    }
    
    //determnine if the anlge for the k-th constraint is null
    for(j=0;j<num;j++){
        LE[j] = N[j]/ newconstr[k][j];
        x +=  LE[j];
    }
    
    for(j=0;j<num;j++){
        if ( x/num - LE[j] < 1e-10) { // y==1, this case applies when the angle for the k-th constraint is null, within rounding erro 
            p++;
        }
    }
    
    if (p == num){
        y = 1;
    }
    
    p = 0;
    
    //mapping of back-up constraints by using the k-th constraint as binding
    newconstr[6][0] = newconstr[k][0];
    newconstr[6][1] = 0;
    newconstr[7][0] = 0;
    newconstr[7][1] = newconstr[k][1];
    if (y != 1) {
        //record the distance between the newly created system constraint, by taking the square root of the sum of squares
        dist_constr_var (newconstr, 0, 0, LA, num, t) ;
        dist_constr_var (newconstr, 0, 1, L6, num, t) ;
        //measure distance between k-th constraint and system generated constraint
        dist_constr_var (newconstr, 1, 0, LK, num, t);
        //back-up constraints are now available
        //measure the angle for each triangle and verify whether the angle formed by the k-th constraint is the difference between the former two
        for(i=0;i<num;i++){
            u = i+t;

v = cosine(LK[u], L6[u], L6[k]);

            v = acos(v);
    
            w = cosine(LA[u], z, L6[u]);
            w = acos(w);
            //determine which of the two system constraints should be used as the second reference L-th constraint with the k-th one
    
            //computational error prevents the certitude of additional for angles; an error of 1e-17 is chosen as an acceptable threshold
            if (abs( L8[k] + w - v )< 1e-10 ){

q=i+t;

}

        }
        for(j=0;j<num;j++){
            if (newconstr[q][j] == newconstr[k][j]){
                a = j;
                v = newconstr[k][j];
            }
        }
    }
    kk=k;
    repeat:
    k=kk;
    //verify whether all constraints, which exclude the k-th constraint, are bound by the k-th primary constraint and its shared axis with the q-th constraint
    //repeat:
    b=0;
    g=0;
   
    u = 2+t;
    for(i=0;i<u;i++) {
    f=0;
    x=0;    

// if and only if the angle for the k-th constraint is not null, then all constraint that are below a threshold, which is defined by the segment binding the k-th constraint and the q-th constraint, will chosen for removal

// review all constraints not excluded from the set of primary constraints, since the array includes all constraints bound by the k-th constraint below a target threshold

if (y != 1) {

if (i != k) {

if (newconstr[i][a] < v) {

                    f++;
                }
                if (newconstr[i][a] >= v) {

x++;

p++;

}

}

if (i == k) {

x++;

}

}

if (y == 1) {

if (i != k) {

if (newconstr[i][s] < newconstr[k][s]) {//repeat loop for each "s" coordinate of the k-th constraint

                    f++;
                }

if (newconstr[i][s] >= newconstr[k][s]) {//repeat loop for each "s" coordinate of the k-th constraint

                    x++;
                    p++;

}

}

if (i == k) {

x++;

}

}

if ( f != 0 ) {

LB[b] = i; //i-th constraint will be recorded in the LB array, which includes all constraints that will be excluded from the set of primary constraints

           b++;// ”b” operator indicates the total number of constraints bound by the k-th constraint
        }
        if ( x != 0 ) {// array used to sort the reference set of primary constraints

LF[g] = i;

            if (i == k) {
                k=g;
            }
            g++;

}

    }
    f=0;
    // if the number of constraints excluded from the set of primary constraints amounts to the total number constraints minus one, which is the k-th constraint, then the k-th constraint is the only binding one
    a=0;
    g=0;
    
    if (p == 0) { a=1; } // a is an operator used to determine whether the first reference constraint is the only binding constraint
    p=0;
    
    if (a == 1) {
        
        //first constraint reference provides bounds for solution; the L-th and k-th constraints are used to determine the solution
        //the following provide the inverse of a 2x2 matrix and output of solution
        if (y != 1) {
            matrix_inverse(newconstr, q, k);
        }
        
        if (y == 1) {
            for(i=0;i<num;i++){
                q=i+t;
                matrix_inverse(newconstr, q, k);
                cout<< "\n";
            }
        }
        
        s++;
        if (s < num) {
            goto repeat;
        }
        
        return 0; //the project ends here, if a=1
    }
    
    //matrix containing set of primary reference constraints, which excludes erased secondary set of reference constraints
    r = (u-b);
    for (i=0;i<r;i++){
        a=LF[i];

LE[i]=LF[i];

for (j=0;j<num;j++){

            formatconstr[i][j] = newconstr[a][j];
        }
    }
    
    //measure distance between constraint coordinates and origin, by taking the square root for the sum of squares
    dist_constr_var (formatconstr, 0, 0, LA, r, 0) ;
    //measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
    dist_constr_var (formatconstr, 0, 1, L6, r, 0) ;
    //recalculating cosine values
    cosine_and_angle(LA,L6,r) ; 
    //begin to locate second reference constraint
    dist_constr_var (formatconstr, 1, 0, LK, r, 0);
    
    f=0;
    for(i=0;i<r;i++){

if (i != k) {

            // angle recorded from the reference k-th constraint and likely second constraint
            v = cosine(LK[i], L6[k], L6[i]);
            v = acos(v); // determine angle v            
        
            //sum of angles should amount to the addition of the k-th constraint angle and likely second reference constraint;
            //this provides an array of distances between the k-th constraints and other constraints; these distances must intersect with the distance from the origin to the objective variable coordinates


            //computational error prevents the certitude of additional for angles; an error of 1e-12 is chosen as an acceptable threshold rather than 1e-17
      
            // L6 array used to store the second reference constraint after the k-th one; this constraint is included in a primary set of constraints that are not bounded by the k-th constraint
        
            if (abs( L8[k] + L8[i] - v ) < 1e-10) {

L6[f] = i;

f++;

            }
        }
    }
    
    //determine angle formed between a slope orthogonal to gamma and the distance between the k-th constraint and the likely second reference constraint
    for(i=0;i<f;i++){
        a = L6[i];
        L7[i] = acos( cosine( L9[a], L9[k], LK[a] ) );
    }
    
    //determine the minimum
    min(L7, f);


    for(i=0;i<f;i++){

if (L7[i]==x){

            q=L6[i];
            a=q;
        }
    }
    
    if (y != 1) {
        matrix_inverse(formatconstr, q, k);    //the following provide the inverse of a 2x2 matrix and output solution     
        x=0;
        p=0;
        for(i=0;i<r;i++) {
            if (i != q) {
                if (i != k) {
                    for (j=0;j<num;j++) {
                        x += formatconstr[i][j] * SOLL[j];
                        if (1 - x < 1e-10) {
                            LF[i]=i;
                            p++;
                        }
                    }
                    x=0;
                }
            }
        }
        
        for(g=0;g<p;g++) {
            q=LF[g];
            matrix_inverse(formatconstr, q, k);
        }
    }
    if (y == 1) {//case when the angle for the k-th constraint is null
        matrix_inverse(formatconstr, q, k);    //the following provide the inverse of a 2x2 matrix and output solution     
        x=0;
        p=0;
        for(i=0;i<r;i++) {
            if (i != q) {
                if (i != k) {
                    for (j=0;j<num;j++) {
                        x += formatconstr[i][j] * SOLL[j];
                            if (1 - x < 1e-10) {
                                LF[i]=i;
                                p++;
                            }
                    }
                    x=0;
                }
            }
        }
           
        for(g=0;g<p;g++) {
            q=LF[g];
           
            matrix_inverse(formatconstr, q, k);
        }
           
        s++;
        if (s < num) {
            goto repeat;
        }
    }  
    
    return 0;
    //This is the end of the project
}
void cosine_and_angle(long double array1[], long double array2[], int upper_bound)
{
    for(i=0; i < upper_bound ;i++){
        //compute cosine value of the triangle formed between the variable coefficients, a given constraint, and the origin
        L8[i] = cosine(array1[i], array2[i], z); 
        
        //multiply distance, which is measured between the variable coefficients and a given constraints, and the cosine value; this provides distance between the two points if they were aligned between the variable coefficients and the origin; this distance will be referenced as gamma
        L9[i]=L8[i]*array2[i];
        
        //the angle for the cosine
        L8[i] = acos(L8[i]);
    }
}
void min(long double array3[], int upper_bound)
{
    //x will be assigned the minimum value
    x = array3[0];
    for(i=0 ; i < upper_bound ; i++) {
        if (array3[i] < x) {
            x = array3[i];
        }
    } 
}
void matrix_inverse(long double array4[2][2], int amount6, int amount7) // inverse of 2 x 2 matrix and output of solution
{
    x=0;
    b=amount6;
    a=amount7;
    
    // matrix containing two constraints meant to be used for obtaining the quantity for each variable
   
    //matrix inverse
    //determinant
    v = array4[a][0] * array4[b][1] - array4[a][1] * array4[b][0];
    
    epi[0][0] = array4[b][1] / v; 
    epi[1][1] = array4[a][0] / v; 
    epi[0][1] = - array4[a][1] / v; 
    epi[1][0] = - array4[b][0] / v; 
    
    //output of solution
    for(i=0;i<num;i++){
        for(j=0;j<num;j++){
            x += epi[i][j];
        }
        SOLL[i]=x;
        x=0;
    }
    sol=0;
    for(i=0;i<num;i++){

sol += N[i] * SOLL[i];

    }
    
    //optimal quantities for each variable are
    cout<< "\noptimal quantities for each decision variable are respectively "<<SOLL[0]<<"  "<< SOLL[1]<< "\n";
    //optimal solution is 
    cout<< "optimal solution is " <<sol << "\n";
    
    //set of optimal constraints which can be combined to determine the optimal quantity for decision variables
    cout<< "\nthis set of constraints yields the optimal quantity for the decision variables \n";
        
    b = LE[b]; //warning: convert double to int
    a = LE[a]; //warning: convert double to int
       
    cout<< "this k-th constraint is pre-requisit for any combination with the set \n[";    
    for(j=0;j<num;j++){
        cout<<constraints[a][j]<<" ";
    }
    cout<<"   <=  "<<limit[a]<<"\n\n";
    
    if ( b > t-1 ) {    
        cout<< "all other constraints within the set are listed \n";
        for(j=0;j<num;j++){
            cout<<newconstr[b][j]*limit[a]<<" ";
        }       
        cout<<"]    <=  "<<limit[a]<<"\n\n";
    }
    
    if ( b <= t-1 ) {    
        cout<< "all other constraints within the set are listed \n";
        for(j=0;j<num;j++){
            cout<<constraints[b][j]<<" ";
        }       
        cout<<"]    <=  "<<limit[b]<<"\n\n";
    } 
    
}
long double cosine(long double amount1, long double amount2, long double amount3) //function used to compute the cosine 
{
    return (pow(amount1,2) - pow(amount2,2) - pow(amount3,2)) / (-2 * amount2 * amount3);    
}
//measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares
void dist_constr_var (long double array5[][2], int amount5, int amount4, long double return_array[], int upper_bound, int increment)
{
    x=0;
    
    for(i=0; i < upper_bound ;i++){
        
        u = i + increment; 
        
        for(j=0;j<num;j++){
            //measure distance between constraint coordinates and variable coefficients, by taking the square root for the sum of squares

x += pow((array5[u][j] - amount5*array5[k][j] - amount4*N[j]), 2);

        }    
        // measure distance between constraint coordinates and variable coefficients

return_array[u] = sqrt(x);

x = 0;

    }
}