```python
from gio import GIO
import numpy as np

##Data from Chan et al. (2018)'s Example 1##
#Notice that we have numpy arrays here.
A = np.array([[2,5],[2,-3],[2,1],[-2,-1]])
b = np.array([[10],[-6],[4],[-10]])
x0 = np.array([[2.5],[3]])

##Creating the GIO object and Generating all of the x^0-ep* for all GIO Models## 
#.GIO_all_measures() calls all of the GIO methods
gio_testing = GIO(A,b,x0)
gio_testing.GIO_all_measures()
print("This is x^0-ep* for p=1: ",gio_testing.x0_epsilon_p[0])
print("This is x^0-ep* for p=2: ",gio_testing.x0_epsilon_p[1])
print("This is x^0-ep* for p='inf': ",gio_testing.x0_epsilon_p[2])
print("This is x^0-ep* for absolute duality: ",gio_testing.x0_epsilon_a)
print("This is x^0-ep* for relative duality: ",gio_testing.x0_epsilon_r)
```

    This is x^0-ep* for p=1:  [[ 2.5       ]
     [ 3.66666667]]
    This is x^0-ep* for p=2:  [[ 2.19230769]
     [ 3.46153846]]
    This is x^0-ep* for p='inf':  [[ 2.1]
     [ 3.4]]
    This is x^0-ep* for absolute duality:  [array([[ 2.1],
           [ 3.4]])]
    This is x^0-ep* for relative duality:  [array([[ 3.16666667],
           [ 3.66666667]])]



```python
##Because we are storing the calculated attributes in lists, we see that, 
##for the last two print statements in the previous cell, a list was returned 
##with an array inside of it.  We can solve this problem by add an index element [0] 
print("This is x^0-ep* for absolute duality: ",gio_testing.x0_epsilon_a[0])
print("This is x^0-ep* for relative duality: ",gio_testing.x0_epsilon_r[0])
```

    This is x^0-ep* for absolute duality:  [[ 2.1]
     [ 3.4]]
    This is x^0-ep* for relative duality:  [[ 3.16666667]
     [ 3.66666667]]



```python
##The GIO class also handles cases when you might minimally project onto multiple
##hyperplanes. The class will put the multiple istar indices into the .istar_multi 
##attribute but will choose the first istar index as the one to continue calculations.

##The class will also output a message to the user to let the user know that 
##this has occurred. 

A_1 = np.array([[-1,-1],[1,-1],[1,1],[-1,1]]) 
b_1 = np.array([[-1],[-1],[-1],[-1]])
x0_1 = np.array([[0],[0.7]])

```

The corresponding feasible region set up for the above $A$ and $b$ is as follows:

\begin{equation}
-x_1 - x_2 \geq -1
\end{equation}
\begin{equation}
x_1 - x_2 \geq -1
\end{equation}
\begin{equation}
x_1 + x_2 \geq -1
\end{equation}
\begin{equation}
-x_1 + x_2 \geq -1
\end{equation}

We can then demonstrate the feasible region using the image below (with approximate placement of $x^0$ below). The numbering corresponds to the order of the constraints below (and then the 0 indexing of Python).




```python
##We will create another GIO object to find the optimal epsilon* under the p=1 norm.
##Due to the position of x0_1 and the properties of the p=1 norm, we know that this
##will result into a projection onto multiple hyperplanes

#We will need to bring in everything again because we had a break in the code 
#cells (with the markup cell above)

from gio import GIO
import numpy as np

A_1 = np.array([[-1,-1],[1,-1],[1,1],[-1,1]]) 
b_1 = np.array([[-1],[-1],[-1],[-1]])
x0_1 = np.array([[0],[0.7]])

gio_multi_project_testing = GIO(A_1,b_1,x0_1)
gio_multi_project_testing.GIO_p(1,'F') 
```

    Under the inf dual norm, x^0 has been projected onto multiple hyperplanes.  For now, we will choose the first i index and will put the rest of the indices in the istar_multi attribute.



```python
##As we can see, a message is printed to the user when x^0 is projected onto multiple
##hyperplanes.
##We can see the multiple istar if we print the .istar_multi attribute
print(gio_multi_project_testing.istar_multi)
```

    [array([0, 1], dtype=int64)]



```python
##To obtain the raw indices, index into the list
print("Raw Indices: ", gio_multi_project_testing.istar_multi[0])
##To obtain only one of the indices, do a double index
print("Single index: ",gio_multi_project_testing.istar_multi[0][0])
```

    Raw Indices:  [0 1]
    Single index:  0



```python
##We see though that only one epsilon* is calculated and thus only one x^0-epsilon* 
##is also calculated
print("This is the epsilon* for p=1: ",gio_multi_project_testing.epsilon_p)
print("This is the (x^0 - epsilon*): ",gio_multi_project_testing.x0_epsilon_p)
```

    This is the epsilon* for p=1:  [array([[-0.3],
           [-0. ]])]
    This is the (x^0 - epsilon*):  [array([[ 0.3],
           [ 0.7]])]



```python
###Calculating rho values###
from gio import GIO
import numpy as np

##Chan et al. (2018)'s Example 1##
##We have changed x0 to demonstrate the rho_p_approx <= rho_p relationship
A = np.array([[2,5],[2,-3],[2,1],[-2,-1]])
b = np.array([[10],[-6],[4],[-10]])
x0 = np.array([[4],[1]])

GIO_rho_test = GIO(A,b,x0)
GIO_rho_test.calculate_rho_p(2) #notice that we did not
                #have to specify if_append='F' because
                #we set a default value in the code
GIO_rho_test.calculate_rho_p_approx(2)
print("This is rho_p exact:",GIO_rho_test.rho_p)
print("This is rho_p - rho_p_approx:",(GIO_rho_test.rho_p[0] - GIO_rho_test.rho_p_approx[0]))

GIO_rho_test.calculate_rho_a()
print("This is rho_a:",GIO_rho_test.rho_a)
GIO_rho_test.calculate_rho_r()
print("This is rho_r:",GIO_rho_test.rho_r)

```

    This is rho_p exact: [0.73886235873548278]
    This is rho_p - rho_p_approx: 0.0232038316283
    This is rho_a: [0.7119341563786008]
    This is rho_r: [0.8851674641148326]



```python
###Testing some values for Validating rho_r#####
from gio import GIO
import numpy as np

##Chan et al. (2018)'s Example 1##
A = np.array([[2,5],[2,-3],[2,1],[-2,-1]])
b = np.array([[10],[-6],[4],[-10]])
x0_1 = np.array([[3],[2]])
x0_2 = np.array([[4],[1]])
x0_3 = np.array([[2],[2]])

validate_rho_r = GIO(A,b,x0_1)
validate_rho_r.calculate_rho_r()
print("This is rho_r for x0_1:",validate_rho_r.rho_r)

validate_rho_r = GIO(A,b,x0_2)
validate_rho_r.calculate_rho_r()
print("This is rho_r for x0_2:",validate_rho_r.rho_r)

validate_rho_r = GIO(A,b,x0_3)
validate_rho_r.calculate_rho_r()
print("This is rho_r for x0_3:",validate_rho_r.rho_r)


```

    This is rho_r for x0_1: [0.71428571428571441]
    This is rho_r for x0_2: [0.8851674641148326]
    Under the b dual norm, x^0 has been projected onto multiple hyperplanes.  For now, we will choose the first i index and will put the rest of the indices in the istar_multi attribute.
    This is rho_r for x0_3: [0.18644067796610142]



```python
###Structural Considerations for GIO Models###
from gio import GIO
import numpy as np
import pyomo.environ as pyo

A = np.array([[2,5],[2,-3],[2,1],[-2,-1]])
b = np.array([[10],[-6],[4],[-10]])
x0 = np.array([[2.5],[3]])

structural_ep_gio_test = GIO(A,b,x0) #instigating the GIO object

###Step 1: Create the Base Model###
#The attribute where the model is stored is called GIO_struc_ep, which is
#the Pyomo model.
structural_ep_gio_test.GIO_structural_epsilon_setup() 

### Step 2: If desired, Add constraints on Epsilon ###
### to GIO_struc_ep ###
### The variable for epsilon is called ep ###
### The index set for variables is varindex ###
def ep_constraint(model):  #should provide the details of the index sets and the numvar parameters
    return model.ep[1] <= model.ep[2] #specifically did not ID the epsilon as nonnegative in gio.py
structural_ep_gio_test.GIO_struc_ep.constraint_ep = pyo.Constraint(rule=ep_constraint)

def neg_ep(model,i):
    return model.ep[i] <= 0
structural_ep_gio_test.GIO_struc_ep.neg_ep = pyo.Constraint(\
            structural_ep_gio_test.GIO_struc_ep.varindex,rule=neg_ep)

### Step 3: Solve the Model ###
#You can specify p as p=1,2,infty
#The solver will generate some output letting you know about infeasibility
#but we also produce output too, letting you know which constraint number 
#in terms of Python indexing that, when forced to be an equality, 
#results in infeasibility
structural_ep_gio_test.GIO_structural_epsilon_solve(2)

print("*************************************************")
print("The c vector found by this workflow is:",structural_ep_gio_test.c_p[0])

###########################################################################################
#We can also generate the same results as the close form solution
structural_ep_no_added_constraints_test = GIO(A,b,x0)
structural_ep_no_added_constraints_test.GIO_structural_epsilon_setup()
structural_ep_no_added_constraints_test.GIO_structural_epsilon_solve(2)

print("The c vector found by this workflow (not adding epsilon constraints) is:",\
     structural_ep_no_added_constraints_test.c_p[0])


```

    WARNING: Loading a SolverResults object with a warning status into
        model=unknown;
            message from solver=Ipopt 3.11.1\x3a Converged to a locally infeasible
            point. Problem may be infeasible.
    We have infeasibility for constraint= 0 .  Putting FLAG in the  container_for_obj_vals vector
    WARNING: Loading a SolverResults object with a warning status into
        model=unknown;
            message from solver=Ipopt 3.11.1\x3a Converged to a locally infeasible
            point. Problem may be infeasible.
    We have infeasibility for constraint= 1 .  Putting FLAG in the  container_for_obj_vals vector
    WARNING: Loading a SolverResults object with a warning status into
        model=unknown;
            message from solver=Ipopt 3.11.1\x3a Converged to a locally infeasible
            point. Problem may be infeasible.
    We have infeasibility for constraint= 2 .  Putting FLAG in the  container_for_obj_vals vector
    *************************************************
    The c vector found by this workflow is: [[-0.66666667]
     [-0.33333333]]
    The c vector found by this workflow (not adding epsilon constraints) is: [[ 0.4]
     [-0.6]]



```python

```

91851df778b0d78b2b02fdf6ed2a625b81adaeb5
ipynb
Jupyter Notebook
GIO_files/jupyter_notebook_demo/demo_gio_class.ipynb
sallen7/inverse_optimization
ef7187a3f9831483abcd7b3d584a9ff0f5f443c2
BSD-3-Clause
2019-07-13T03:44:28.000Z
2019-07-13T03:44:28.000Z
GIO_files/jupyter_notebook_demo/demo_gio_class.ipynb
sallen7/inverse_optimization
ef7187a3f9831483abcd7b3d584a9ff0f5f443c2
BSD-3-Clause
GIO_files/jupyter_notebook_demo/demo_gio_class.ipynb
sallen7/inverse_optimization
ef7187a3f9831483abcd7b3d584a9ff0f5f443c2
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES