Merge pull request #28 from khalil-research/portfolio

Add Portfolio dataset

pkg/meta.yaml

@@ -3,7 +3,7 @@
 
 package:
   name: pyepo
-  version: 0.3.7
+  version: 0.3.8
 
 source:
   path: ./

pkg/pyepo/data/__init__.py

@@ -4,4 +4,4 @@
 Synthetic data generation
 """
 
-from pyepo.data import dataset, shortestpath, knapsack, tsp
+from pyepo.data import dataset, shortestpath, knapsack, tsp, portfolio

pkg/pyepo/data/portfolio.py

@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+# coding: utf-8
+"""
+Synthetic data for portfolio
+"""
+
+import numpy as np
+
+
+def genData(num_data, num_features, num_assets, deg=1, noise_level=1, seed=135):
+    """
+    A function to generate synthetic data and features for travelling salesman
+
+    Args:
+        num_data (int): number of data points
+        num_features (int): dimension of features
+        num_assets (int): number of assets
+        deg (int): data polynomial degree
+        noise_level (float): level of data random noise
+        seed (int): random seed
+
+    Returns:
+        tuple: data features (np.ndarray), costs (np.ndarray)
+    """
+    # positive integer parameter
+    if type(deg) is not int:
+        raise ValueError("deg = {} should be int.".format(deg))
+    if deg <= 0:
+        raise ValueError("deg = {} should be positive.".format(deg))
+    # set seed
+    rnd = np.random.RandomState(seed)
+    # number of data points
+    n = num_data
+    # dimension of features
+    p = num_features
+    # number of assets
+    m = num_assets
+    # random matrix parameter B
+    B = rnd.binomial(1, 0.5, (m, p))
+    # random matrix parameter L
+    L = rnd.uniform(-2.5e-3*noise_level, 2.5e-3*noise_level, (num_assets, num_features))
+    # feature vectors
+    x = rnd.normal(0, 1, (n, p))
+    # value of items
+    r = np.zeros((n, m))
+    for i in range(n):
+        # mean return of assets
+        r[i] = (0.05 * np.dot(B, x[i].reshape(p, 1)).T / np.sqrt(p) + \
+                0.1 ** (1 / deg)) ** deg
+        # random noise
+        f = rnd.randn(num_features)
+        eps = rnd.randn(num_assets)
+        r[i] += L @ f + 0.01 * noise_level * eps
+    # covariance matrix of the returns
+    cov = L @ L.T + (1e-2 * noise_level) * np.eye(num_assets)
+    return cov, x, r

pkg/pyepo/func/utlis.py

@@ -42,6 +42,9 @@ def _solve_in_pass(cp, optmodel, processes, pool):
             solp, objp = optmodel.solve()
             sol.append(solp)
             obj.append(objp)
+        # to numpy
+        sol = np.array(sol)
+        obj = np.array(obj)
     # multi-core
     else:
         # get class

pkg/pyepo/model/grb/__init__.py

@@ -8,3 +8,4 @@
 from pyepo.model.grb.shortestpath import shortestPathModel
 from pyepo.model.grb.knapsack import knapsackModel
 from pyepo.model.grb.tsp import tspGGModel, tspDFJModel, tspMTZModel
+from pyepo.model.grb.portfolio import portfolioModel

pkg/pyepo/model/grb/grbmodel.py

@@ -34,6 +34,13 @@ def __init__(self):
     def __repr__(self):
         return "optGRBModel " + self.__class__.__name__
 
+    @property
+    def num_cost(self):
+        """
+        number of cost to be predicted
+        """
+        return self.x.size if isinstance(self.x, gp.MVar) else len(self.x)
+
     def setObj(self, c):
         """
         A method to set objective function
@@ -43,7 +50,12 @@ def setObj(self, c):
         """
         if len(c) != self.num_cost:
             raise ValueError("Size of cost vector cannot match vars.")
-        obj = gp.quicksum(c[i] * self.x[k] for i, k in enumerate(self.x))
+        # mvar
+        if isinstance(self.x, gp.MVar):
+            obj = c @ self.x
+        # vars
+        else:
+            obj = gp.quicksum(c[i] * self.x[k] for i, k in enumerate(self.x))
         self._model.setObjective(obj)
 
     def solve(self):
@@ -55,7 +67,14 @@ def solve(self):
         """
         self._model.update()
         self._model.optimize()
-        return [self.x[k].x for k in self.x], self._model.objVal
+        # solution
+        if isinstance(self.x, gp.MVar):
+            sol = self.x.x
+        else:
+            sol = [self.x[k].x for k in self.x]
+        # objective value
+        obj = self._model.objVal
+        return sol, obj
 
     def copy(self):
         """
@@ -71,7 +90,7 @@ def copy(self):
         new_model._model = self._model.copy()
         # variables for new model
         x = new_model._model.getVars()
-        new_model.x = {key: x[i] for i, key in enumerate(self.x)}
+        new_model.x = {key: x[i] for i, key in enumerate(x)}
         return new_model
 
     def addConstr(self, coefs, rhs):

pkg/pyepo/model/grb/knapsack.py

@@ -30,7 +30,7 @@ def __init__(self, weights, capacity):
         """
         self.weights = np.array(weights)
         self.capacity = np.array(capacity)
-        self.items = list(range(self.weights.shape[1]))
+        self.items = self.weights.shape[1]
         super().__init__()
 
     def _getModel(self):
@@ -43,13 +43,11 @@ def _getModel(self):
         # ceate a model
         m = gp.Model("knapsack")
         # varibles
-        x = m.addVars(self.items, name="x", vtype=GRB.BINARY)
+        x = m.addMVar(self.items, name="x", vtype=GRB.BINARY)
         # sense
         m.modelSense = GRB.MAXIMIZE
         # constraints
-        for i in range(len(self.capacity)):
-            m.addConstr(gp.quicksum(self.weights[i,j] * x[j]
-                        for j in self.items) <= self.capacity[i])
+        m.addConstr(self.weights @ x <= self.capacity)
         return m, x
 
     def relax(self):
@@ -75,32 +73,29 @@ def _getModel(self):
         # turn off output
         m.Params.outputFlag = 0
         # varibles
-        x = m.addVars(self.items, name="x", ub=1)
+        x = m.addMVar(self.items, name="x", ub=1)
         # sense
         m.modelSense = GRB.MAXIMIZE
         # constraints
-        for i in range(len(self.capacity)):
-            m.addConstr(gp.quicksum(self.weights[i,j] * x[j]
-                        for j in self.items) <= self.capacity[i])
+        m.addConstr(self.weights @ x <= self.capacity)
         return m, x
 
     def relax(self):
         """
         A forbidden method to relax MIP model
         """
         raise RuntimeError("Model has already been relaxed.")
-        
-        
+
+
 if __name__ == "__main__":
-
-    import random
+
     # random seed
-    random.seed(42)
+    np.random.seed(42)
     # set random cost for test
-    cost = [random.random() for _ in range(16)]
+    cost = np.random.random(16)
     weights = np.random.choice(range(300, 800), size=(2,16)) / 100
     capacity = [20, 20]
-    
+
     # solve model
     optmodel = knapsackModel(weights=weights, capacity=capacity) # init model
     optmodel = optmodel.copy()
@@ -111,7 +106,7 @@ def relax(self):
     for i in range(16):
         if sol[i] > 1e-3:
             print(i)
-            
+
     # relax
     optmodel = optmodel.relax()
     optmodel.setObj(cost) # set objective function
@@ -121,7 +116,7 @@ def relax(self):
     for i in range(16):
         if sol[i] > 1e-3:
             print(i)
-            
+
     # add constraint
     optmodel = optmodel.addConstr([weights[0,i] for i in range(16)], 10)
     optmodel.setObj(cost) # set objective function

pkg/pyepo/model/grb/portfolio.py

@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+# coding: utf-8
+"""
+Portfolio problem
+"""
+
+import numpy as np
+import gurobipy as gp
+from gurobipy import GRB
+
+from pyepo.model.grb.grbmodel import optGrbModel
+
+
+class portfolioModel(optGrbModel):
+    """
+    This class is optimization model for portfolio problem
+
+    Attributes:
+        _model (GurobiPy model): Gurobi model
+        num_assets (int): number of assets
+        cov (numpy.ndarray): covariance matrix of the returns
+        risk_level (float): risk level
+    """
+
+    def __init__(self, num_assets, covariance, gamma=2.25):
+        """
+        Args:
+            num_assets (int): number of assets
+            covariance (numpy.ndarray): covariance matrix of the returns
+            gamma (float): risk level parameter
+        """
+        self.num_assets = num_assets
+        self.cov = covariance
+        self.risk_level = self._getRiskLevel(gamma)
+        super().__init__()
+
+    def _getRiskLevel(self, gamma):
+        """
+        A method to calculate risk level
+
+        Returns:
+            float: risk level
+        """
+        risk_level = gamma * np.mean(self.cov)
+        return risk_level
+
+    def _getModel(self):
+        """
+        A method to build Gurobi model
+
+        Returns:
+            tuple: optimization model and variables
+        """
+        # ceate a model
+        m = gp.Model("portfolio")
+        # varibles
+        x = m.addMVar(self.num_assets, ub=1, name="x")
+        # sense
+        m.modelSense = GRB.MAXIMIZE
+        # constraints
+        m.addConstr(x.sum() == 1, "budget")
+        m.addConstr(x.T @ self.cov @ x <= self.risk_level, "risk_limit")
+        return m, x
+
+
+if __name__ == "__main__":
+
+    import random
+    from pyepo.data.portfolio import genData
+    # random seed
+    random.seed(42)
+    # set random cost for test
+    cov, _, revenue = genData(num_data=100, num_features=4, num_assets=50, deg=2)
+
+    # solve model
+    optmodel = portfolioModel(num_assets=50, covariance=cov) # init model
+    optmodel = optmodel.copy()
+    optmodel.setObj(revenue[0]) # set objective function
+    sol, obj = optmodel.solve() # solve
+    # print res
+    print('Obj: {}'.format(obj))
+    for i in range(50):
+        if sol[i] > 1e-3:
+            print("Asset {}: {:.2f}%".format(i, 100*sol[i]))
+
+    # add constraint
+    optmodel = optmodel.addConstr([1]*50, 30)
+    optmodel.setObj(revenue[0]) # set objective function
+    sol, obj = optmodel.solve() # solve
+    # print res
+    print('Obj: {}'.format(obj))
+    for i in range(50):
+        if sol[i] > 1e-3:
+            print("Asset {}: {:.2f}%".format(i, 100*sol[i]))

pkg/setup.py

@@ -17,7 +17,7 @@
     long_description=long_description,
     long_description_content_type="text/markdown",
     # version
-    version = "0.3.7",
+    version = "0.3.8",
     # Github repo
     url = "https://github.com/khalil-research/PyEPO",
     # author name