{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":1370,"sourceType":"datasetVersion","datasetId":738},{"sourceId":9243,"sourceType":"datasetVersion","datasetId":2243},{"sourceId":1553068,"sourceType":"datasetVersion","datasetId":916586},{"sourceId":2085889,"sourceType":"datasetVersion","datasetId":1250649},{"sourceId":4494999,"sourceType":"datasetVersion","datasetId":2628764},{"sourceId":6144787,"sourceType":"datasetVersion","datasetId":3523739}],"dockerImageVersionId":30615,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\ndef warn(*args, **kwargs):\n    pass\nimport warnings\nwarnings.warn = warn\n%pip install gurobipy \nimport gurobipy as gp\nfrom gurobipy import GRB \nfrom gurobipy import  *\n\nfrom math import sqrt\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport time\n#import Pytorch\nfrom sklearn.ensemble import *\nfrom sklearn import tree\nfrom sklearn.metrics import mean_squared_error, log_loss, accuracy_score\nfrom sklearn.neural_network import MLPRegressor, MLPClassifier\nfrom sklearn.preprocessing import LabelEncoder, OneHotEncoder\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\nGLOBAL_RANDOM=42\nN=10\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-12-19T20:32:35.699138Z","iopub.execute_input":"2023-12-19T20:32:35.699601Z","iopub.status.idle":"2023-12-19T20:32:50.777830Z","shell.execute_reply.started":"2023-12-19T20:32:35.699568Z","shell.execute_reply":"2023-12-19T20:32:50.776173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Load Data\nall_train = np.genfromtxt('/kaggle/input/fashionmnist/fashion-mnist_train.csv', delimiter=',')\nall_test = np.genfromtxt('/kaggle/input/fashionmnist/fashion-mnist_test.csv', delimiter=',')\ntrain_data=all_train[1:,1:]\ntrain_label=all_train[1:,0]\ntest_data=all_test[1:,1:]\ntest_label=all_test[1:,0]\nEvery_Dataset=[]\nEvery_Dataset.append([[train_data,train_label],[test_data,test_label],'fashion-mnist'])\n","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:03:03.237404Z","iopub.execute_input":"2023-12-19T20:03:03.237980Z","iopub.status.idle":"2023-12-19T20:04:20.150848Z","shell.execute_reply.started":"2023-12-19T20:03:03.237942Z","shell.execute_reply":"2023-12-19T20:04:20.149829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all =pd.read_csv('/kaggle/input/loan-approval-prediction-dataset/loan_approval_dataset.csv')#, delimiter=',')\n#print(all.columns)\nall = all.drop('loan_id', axis=1)\n\n#print(all.iloc[1,:])\none_hot = pd.get_dummies(all[' education'],dtype=float)\nall = all.drop(' education', axis=1)\nall = all.join(one_hot)\none_hot = pd.get_dummies(all[' self_employed'],dtype=float)\nall = all.drop(' self_employed', axis=1)\nall = all.join(one_hot)\none_hot = pd.get_dummies(all[' loan_status'],dtype=float)\nall = all.drop(' loan_status', axis=1)\nall = all.join(one_hot)\nprint(all.iloc[1,:])\nall=np.array(all)\n\ntrain_data,train_label=all[1:int(.85*all.shape[0]),0:-1],all[1:int(.85*all.shape[0]),-1:].flatten()\ntest_data,test_label=all[1+int(.85*all.shape[0]):,0:-1],all[1+int(.85*all.shape[0]):,-1:].flatten()\nEvery_Dataset.append([[train_data,train_label],[test_data,test_label],'loan'])","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:04:20.152235Z","iopub.execute_input":"2023-12-19T20:04:20.153128Z","iopub.status.idle":"2023-12-19T20:04:20.237424Z","shell.execute_reply.started":"2023-12-19T20:04:20.153092Z","shell.execute_reply":"2023-12-19T20:04:20.235901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"files=[['/kaggle/input/fetal-health-classification/fetal_health.csv','fetal_health'],\n      #['/kaggle/input/glass/glass.csv','glass'],\n      #['/kaggle/input/oil-spill-detection/oil_spill.csv','oil_spill'],\n      ['/kaggle/input/pulsar-classification-for-class-prediction/Pulsar.csv','Pulsar']  \n      ]\nfor file in files:\n    all =np.genfromtxt(file[0], delimiter=',')\n    np.random.shuffle(all)\n    train_data,train_label=all[1:int(.85*all.shape[0]),0:-1],all[1:int(.85*all.shape[0]),-1:].flatten()\n    test_data,test_label=all[1+int(.85*all.shape[0]):,0:-1],all[1+int(.85*all.shape[0]):,-1:].flatten()\n    d=np.isnan(train_label)\n    t=np.isnan(test_label)\n    \n    train_label=train_label[np.logical_not(d)]\n    train_data=train_data[np.logical_not(d)]\n    test_label=test_label[np.logical_not(t)]\n    test_data=test_data[np.logical_not(t)]\n    \n    \n    \n    \n    Every_Dataset.append([[train_data,train_label],[test_data,test_label],file[1]])","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:04:20.241208Z","iopub.execute_input":"2023-12-19T20:04:20.241709Z","iopub.status.idle":"2023-12-19T20:04:20.531381Z","shell.execute_reply.started":"2023-12-19T20:04:20.241667Z","shell.execute_reply":"2023-12-19T20:04:20.529992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nfor item in Every_Dataset:\n    #print(len(item))\n    #print(item)\n    print(item[0][0].shape,item[0][1].shape,item[1][0].shape,item[1][1].shape,item[2])","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:04:20.533304Z","iopub.execute_input":"2023-12-19T20:04:20.533674Z","iopub.status.idle":"2023-12-19T20:04:20.541077Z","shell.execute_reply.started":"2023-12-19T20:04:20.533641Z","shell.execute_reply":"2023-12-19T20:04:20.539728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"print(train_data.shape)\nprint(train_label.shape)\nprint(test_data.shape)\nprint(test_label.shape)","metadata":{}},{"cell_type":"code","source":"#FNN\n#Sklearn\nresults = np.array([[\"Iter\",\"Time\",\"Accuracy\",'Loss',\"Dataset\"]])\nn=2*N#5*N\nfor dataset in Every_Dataset:\n    print(dataset[2])\n    train_data=dataset[0][0]\n    train_label=dataset[0][1]\n    test_data=dataset[1][0]\n    test_label=dataset[1][1]\n    labels=np.unique(test_label)\n    \n    print(train_data.shape)\n    print(test_data.shape)  \n    labels=np.append(labels,np.unique(train_label), axis=0)\n    labels=np.unique(labels)\n      \n    print(labels)\n    for i in range(1, n):\n        #FNN=MLPRegressor(random_state=GLOBAL_RANDOM, max_iter=i)\n        FNN=MLPClassifier(random_state=GLOBAL_RANDOM, max_iter=i,shuffle=False,hidden_layer_sizes=(15,))\n        start = time.time()\n        FNN.fit(train_data, train_label)\n        end = time.time()\n        y_pred = FNN.predict(test_data)\n        #print(np.unique(y_pred),np.unique(test_label))\n        y_pred_log=FNN.predict_proba(test_data)\n        #print(y_pred_log,np.unique(y_pred),np.unique(test_label))\n        #print(np.unique(y_pred),)\n        loss = log_loss(test_label, y_pred_log,)#labels=labels)\n        acc  = accuracy_score(test_label, y_pred)\n        # Calculate the time elapsed\n        elapsed = end - start\n        # Append the accuracy and time to the results array\n       \n        \n        results = np.append(results, [[i,elapsed,acc,loss,dataset[2]]], axis=0)\n    # save the dataframe as a csv file\n#print(results.shape)\n#print(results[1:,:])\ndf = pd.DataFrame(results[1:,:], columns=[str(x) for x in results[0,:]])\n#print(df.head())\n#DF = pd.DataFrame(results)\n \n# save the dataframe as a csv file\ndf.to_csv(\"FNN_\"+str(n)+\".csv\")\n","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:32:50.781199Z","iopub.execute_input":"2023-12-19T20:32:50.781742Z","iopub.status.idle":"2023-12-19T20:36:18.326966Z","shell.execute_reply.started":"2023-12-19T20:32:50.781686Z","shell.execute_reply":"2023-12-19T20:36:18.324611Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#RF\n#USE TREE COUNT instead of Epochs\n#results = np.empty((0, 2))\nresults = np.array([[\"Iter\",\"Time\",\"Accuracy\",'Loss',\"Dataset\"]])\nn=2*N\nfor dataset in Every_Dataset:\n    print(dataset[2])\n    train_data=dataset[0][0]\n    train_label=dataset[0][1]\n    test_data=dataset[1][0]\n    test_label=dataset[1][1]\n    labels=np.unique(test_label)\n    \n    print(train_data.shape)\n    print(test_data.shape)  \n    labels=np.append(labels,np.unique(train_label), axis=0)\n    labels=np.unique(labels)\n      \n    print(labels)\n    for i in range(1, n):\n        # Create a random forest classifier with i trees\n        clf = RandomForestClassifier(n_estimators=i, random_state=GLOBAL_RANDOM)\n        # Record the start time\n        start = time.time()\n        # Fit the classifier on the train data\n\n        clf.fit(train_data, np.ravel(train_label))\n        # Predict on the test data    \n        # Record the end time\n        end = time.time()\n        y_pred = clf.predict(test_data)\n        y_pred_log=clf.predict_proba(test_data)\n        \n        loss = log_loss(test_label, y_pred_log,)#labels=labels)\n        acc  = accuracy_score(test_label, y_pred)\n        # Calculate the time elapsed\n        elapsed = end - start\n        # Append the accuracy and time to the results array\n                # Append the accuracy and time to the results array\n\n        #print(dataset[2])\n        results = np.append(results, [[i,elapsed,acc,loss,dataset[2]]], axis=0)\n        # save the dataframe as a csv file\n\n#print(results.shape)\n#print(results[1:,:])\ndf = pd.DataFrame(results[1:,:], columns=[str(x) for x in results[0,:]])\n#print(df.head())\n#DF = pd.DataFrame(results)\n \n# save the dataframe as a csv file\ndf.to_csv(\"RF_\"+str(n)+\".csv\")\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-19T20:36:18.330899Z","iopub.execute_input":"2023-12-19T20:36:18.331876Z","iopub.status.idle":"2023-12-19T20:40:24.981075Z","shell.execute_reply.started":"2023-12-19T20:36:18.331790Z","shell.execute_reply":"2023-12-19T20:40:24.979543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#GBT\n#USE TREE COUNT instead of Epochs\n#results = np.empty((0, 2))\nresults = np.array([[\"Iter\",\"Time\",\"Accuracy\",'Loss',\"Dataset\"]])\nn=N\nfor dataset in Every_Dataset:\n    print(dataset[2])\n    train_data=dataset[0][0]\n    train_label=dataset[0][1]\n    test_data=dataset[1][0]\n    test_label=dataset[1][1]\n    labels=np.unique(test_label)\n    \n    print(train_data.shape)\n    print(test_data.shape)  \n    labels=np.append(labels,np.unique(train_label), axis=0)\n    labels=np.unique(labels)\n      \n    print(labels)\n    for i in range(1,n):\n        # Create a random forest classifier with i trees\n        GBT = HistGradientBoostingClassifier(max_iter=i, random_state=GLOBAL_RANDOM)\n        # Record the start time\n        start = time.time()\n        # Fit the classifier on the train data\n        GBT.fit(train_data, train_label)\n        # Predict on the test data\n\n        # Record the end time\n        end = time.time()\n        y_pred = GBT.predict(test_data)\n        y_pred_log=GBT.predict_proba(test_data)\n        \n        loss = log_loss(test_label, y_pred_log,)#labels=labels)\n        acc  = accuracy_score(test_label, y_pred)\n        # Calculate the time elapsed\n        elapsed = end - start\n        # Append the accuracy and time to the results array\n                # Append the accuracy and time to the results array\n\n        #print(dataset[2])\n        results = np.append(results, [[i,elapsed,acc,loss,dataset[2]]], axis=0)\n\n# Print the results array\n#print(results.shape)\n#print(results[1:,:])\ndf = pd.DataFrame(results[1:,:], columns=[str(x) for x in results[0,:]])\n#print(df.head())\n#DF = pd.DataFrame(results)\n \n# save the dataframe as a csv file\ndf.to_csv(\"GBT_\"+str(n)+\".csv\")\n","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:40:24.982683Z","iopub.execute_input":"2023-12-19T20:40:24.983068Z","iopub.status.idle":"2023-12-19T20:43:16.890583Z","shell.execute_reply.started":"2023-12-19T20:40:24.983028Z","shell.execute_reply":"2023-12-19T20:43:16.889059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#BAY\n","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:06:58.239262Z","iopub.execute_input":"2023-12-19T20:06:58.244527Z","iopub.status.idle":"2023-12-19T20:06:58.250246Z","shell.execute_reply.started":"2023-12-19T20:06:58.244477Z","shell.execute_reply":"2023-12-19T20:06:58.248774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#SVM","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:06:58.254988Z","iopub.execute_input":"2023-12-19T20:06:58.255392Z","iopub.status.idle":"2023-12-19T20:06:58.264185Z","shell.execute_reply.started":"2023-12-19T20:06:58.255360Z","shell.execute_reply":"2023-12-19T20:06:58.262853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#REG","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:06:58.265598Z","iopub.execute_input":"2023-12-19T20:06:58.265983Z","iopub.status.idle":"2023-12-19T20:06:58.278009Z","shell.execute_reply.started":"2023-12-19T20:06:58.265952Z","shell.execute_reply":"2023-12-19T20:06:58.276793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#CNN","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:06:58.279452Z","iopub.execute_input":"2023-12-19T20:06:58.280188Z","iopub.status.idle":"2023-12-19T20:06:58.291991Z","shell.execute_reply.started":"2023-12-19T20:06:58.280152Z","shell.execute_reply":"2023-12-19T20:06:58.290647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#RNN","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:06:58.293551Z","iopub.execute_input":"2023-12-19T20:06:58.293959Z","iopub.status.idle":"2023-12-19T20:06:58.303535Z","shell.execute_reply.started":"2023-12-19T20:06:58.293926Z","shell.execute_reply":"2023-12-19T20:06:58.302118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#N=5\nfiles=['/kaggle/working/FNN_'+str(2*N)+'.csv',\n       '/kaggle/working/RF_'+str(2*N)+'.csv',\n       '/kaggle/working/GBT_'+str(N)+'.csv'\n      ]\nall_epochs=np.array([[0,0,0,0]])\nall_epochs_time=np.array([[0,0,0,0]])\nnames=[]\nc=0\nfor file in files:\n    data=pd.read_csv(file)\n    #print(data)\n    #data=data.iloc[:,1:]\n    nam=file\n    for i in range(1,3*N):\n        c+=1\n        #print(df)\n       #np.array(data['Accuracy'])\n        #print(all_epochs.shape)\n        try:\n            df=data.loc[data['Iter'] == i]\n            #names.append(str(file[-10:-7]+'_'+str(i)))\n            #print(df)\n            #print(i,'time',all_epochs_time.shape)\n            #X = np.insert(X\n            #a=np.append(np.array(df['Accuracy']),values=file[-10:-7]+'_'+str(i), axis=0)\n            all_epochs=np.vstack((all_epochs, np.array(df['Accuracy'])))\n            all_epochs_time=np.vstack((all_epochs_time, np.array(df['Time'])))\n            names.append(str(file[-10:-7]+'_'+str(i)))\n            #print('time',all_epochs_time.shape)\n            #print(i,'time',all_epochs_time.shape)\n            #print(df['Accuracy'])\n            #print(np.array(df['Time']))\n            #print(np.array(df['Accuracy']))\n            #print(np.array(df['Dataset']))\n            #all_epochs=np.append(all_epochs,values=np.array(data['Accuracy']), axis=1)\n        except Exception as e:\n            print(c,e)\n            print('time',all_epochs_time.shape)\n\n            break\n\n#all_epochs=all_epochs.transpose()\n#names=np.array(names)\nall_epochs=np.delete(all_epochs, 0, axis=0)\nall_epochs_time=np.delete(all_epochs_time, 0, axis=0)\n\n#print('acc',all_epochs)\n\nprint('time',all_epochs_time.shape)\ncovar_matrix=np.cov(all_epochs)\nmean_acc=np.array([x.mean() for x in all_epochs]) #for a given model and epoch, over all datasets\nmean_var=np.array([x.var() for x in all_epochs]) #for a given model and epoch, over all datasets\nmean_time=np.array([x.mean() for x in all_epochs_time])\n#print(covar_matrix.shape)\n#print(covar_matrix)\n#print(mean_acc.shape)\n#print('mean time',mean_time)\n#print('part',np.sort(mean_time)[-2])\nMAX_TIME=np.sort(mean_time)[-2]#(mean_time.max()+np.sort(mean_time)[-2])/2\nprint('max',MAX_TIME)\nnorm_mean_time=np.array([x/MAX_TIME for x in mean_time])\n#print(all_epochs[0])\nprint(norm_mean_time)\nprint(norm_mean_time.shape)\nprint(names)\nprint(len(names))\nmean_shift=mean_acc\nmean_shift=np.array(list(mean_shift[-1:]) + list(mean_shift[:-1]))#np.array(list(mean_acc).rotate(1))    \nprint(mean_shift)","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:47:04.078363Z","iopub.execute_input":"2023-12-19T20:47:04.078804Z","iopub.status.idle":"2023-12-19T20:47:04.138564Z","shell.execute_reply.started":"2023-12-19T20:47:04.078769Z","shell.execute_reply":"2023-12-19T20:47:04.137480Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Create all constraints for only one epoch per model\nimport itertools\n\ndef gen_comb():\n    all_const=[]#list(itertools.combinations([x for x in range(0,N-1)], r=2))]\n    #for i in range(0,3):  #number of models\n    for i in range(0,3):  #number of models\n        if i==0:\n            #perms=list(itertools.permutations([x for x in range(0,5*N)], r=2))\n            perms=list(itertools.combinations([x for x in range(1,2*N)], r=2))\n\n            all_const.append(perms)\n        if i==1:\n            perms=list(itertools.combinations([x for x in range(2*N,4*N-1)], r=2))\n            all_const.append(perms)\n        if i==2:\n            perms=list(itertools.combinations([x for x in range(4*N-1,5*N-1)], r=2))\n            all_const.append(perms)\n    #print(np.array(all_const).shape)\n    #print(all_const)\n    return(all_const)\na=gen_comb()\nprint(len(a))\nfor elm in a:\n    print(len(elm))\nprint(a[2])","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:49:59.457055Z","iopub.execute_input":"2023-12-19T20:49:59.457542Z","iopub.status.idle":"2023-12-19T20:49:59.469551Z","shell.execute_reply.started":"2023-12-19T20:49:59.457505Z","shell.execute_reply":"2023-12-19T20:49:59.468158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Put time porportion constraints on each epoch","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:06:58.431172Z","iopub.execute_input":"2023-12-19T20:06:58.431628Z","iopub.status.idle":"2023-12-19T20:06:58.444471Z","shell.execute_reply.started":"2023-12-19T20:06:58.431593Z","shell.execute_reply":"2023-12-19T20:06:58.443168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def epoch_exlusiviy(m):\n    for model in gen_comb():\n        for epoch in model:\n            #print((epoch[0]-1,epoch[1]-1))\n            m.addConstr(x[epoch[0]-1]*x[epoch[1]-1] ==0,\"exlusive\"+str(epoch[0])+str(epoch[1]))\n\ndef epoch_bounds(m):\n    for i in range(0,2*N-1): #NN\n        if i==0:\n            #m.addConstr(0<=(x[0]*(x[0]-,'0l')\n            m.addConstr(x[0]*(x[0]-norm_mean_time[1])<=0,'0u')\n            print('lower:',0,' x ','upper:',norm_mean_time[1])\n        else:\n            m.addConstr(x[i]*(x[i]-(norm_mean_time[i-1]+epsilon))>=0,str(i)+'l')\n            m.addConstr(x[i]*(x[i]-norm_mean_time[i])<=0,str(i)+'u')\n            print('lower:',norm_mean_time[i-1], ' x ','upper:',norm_mean_time[i])\n\n    for i in range(2*N-1,4*N-2): #NN\n        if i==2*N-1:\n            #m.addConstr(0<=(x[0]*(x[0]-,'0l')\n            m.addConstr(x[2*N-1]*(x[2*N-1]-norm_mean_time[2*N])<=0,'0u')\n            print('lower:',0,' x ','upper:',norm_mean_time[2*N])\n        else:\n            m.addConstr(x[i]*(x[i]-(norm_mean_time[i-1]+epsilon))>=0,str(i)+'l')\n            m.addConstr(x[i]*(x[i]-norm_mean_time[i])<=0,str(i)+'u')\n            print('lower:',norm_mean_time[i-1], ' x ','upper:',norm_mean_time[i])\n    \n    for i in range(4*N-2,5N-3): #NN\n        if i==4*N-2:\n            #m.addConstr(0<=(x[0]*(x[0]-,'0l')\n            m.addConstr(x[4*N-2]*(x[4*N-2]-norm_mean_time[4*N-1])<=0,'0u')\n            print('lower:',0,' x ','upper:',norm_mean_time[4*N-1])\n        else:\n            m.addConstr(x[i]*(x[i]-(norm_mean_time[i-1]+epsilon))>=0,str(i)+'l')\n            m.addConstr(x[i]*(x[i]-norm_mean_time[i])<=0,str(i)+'u')\n            print('lower:',norm_mean_time[i-1], ' x ','upper:',norm_mean_time[i])\n    ","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:06:58.446310Z","iopub.execute_input":"2023-12-19T20:06:58.446860Z","iopub.status.idle":"2023-12-19T20:06:58.462430Z","shell.execute_reply.started":"2023-12-19T20:06:58.446797Z","shell.execute_reply":"2023-12-19T20:06:58.461391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def XOR_test(m):\n    m.addConstr(e1.sum() == 1, 'exclusive_sum_1')\n    m.addConstr(e2.sum() == 1, 'exclusive_sum_2')\n    m.addConstr(e3.sum() == 1, 'exclusive_sum_3')\n    for i in range (0,2*N-1):\n        m.addConstr(x[i]<=e1[i],'exclusive_var_'+str((1,i))) #x<=0 or x<=1\n    \n        #m.addConstr(x[i]*e1[i]<=norm_mean_time[i]*e1[i])\n    m.addConstr(x[0]<=norm_mean_time[0]*e1[0])#*e1[0]\n    for i in range (1,2*N-1):\n        m.addConstr(x[i]>=norm_mean_time[i-1]*e1[i])#*e1[i]\n        m.addConstr(x[i]<=norm_mean_time[i]*e1[i])#*e1[i]\n        print('Epoch ',i+1, 'lower:',norm_mean_time[i-1], ' x ','upper:',norm_mean_time[i])\n    \n   ##################################\n    for i in range (2*N-1,4*N-2):\n        m.addConstr(x[i]<=e2[i-(2*N-1)],'exclusive_var_'+str((2,i)))\n        \n    m.addConstr(x[2*N-1]<=norm_mean_time[2*N-1]*e2[0]) #*e2[0]\n    #print('Epoch ',i+1, ' x ','upper:',norm_mean_time[i])\n    for i in range (2*N,4*N-2):\n        print('Epoch ',i+1, 'lower:',norm_mean_time[i-1], ' x ','upper:',norm_mean_time[i])\n\n        m.addConstr(x[i]>=norm_mean_time[i-1]*e2[i-2*N]) #*e2[i-9]\n        m.addConstr(x[i]<=norm_mean_time[i]*e2[i-2*N])#*e2[i-9]\n    ##################################\n    for i in range (4*N-2,5*N-3):\n        m.addConstr(x[i]<=e3[i-(4*N-2)],'exclusive_var_'+str((3,i)))\n        \n    m.addConstr(x[4*N-2]<=norm_mean_time[4*N-2]*e3[0]) #*e2[0]\n    #print('Epoch ',i+1, ' x ','upper:',norm_mean_time[i])\n    for i in range (4*N-1,5*N-3):\n        print('Epoch ',i+1, 'lower:',norm_mean_time[i-1], ' x ','upper:',norm_mean_time[i])\n\n        m.addConstr(x[i]>=norm_mean_time[i-1]*e3[i-(4*N-1)]) #*e2[i-9]\n        m.addConstr(x[i]<=norm_mean_time[i]*e3[i-(4*N-1)])#*e2[i-9]\n\n        ","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:57:59.050736Z","iopub.execute_input":"2023-12-19T20:57:59.051250Z","iopub.status.idle":"2023-12-19T20:57:59.072976Z","shell.execute_reply.started":"2023-12-19T20:57:59.051215Z","shell.execute_reply":"2023-12-19T20:57:59.070795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#https://colab.research.google.com/github/Gurobi/modeling-examples/blob/master/portfolio_selection_optimization/portfolio_selection_optimization.ipynb#scrollTo=qPMdJ5Vdzlxp\n\nwith gp.Env(empty=True) as env:\n    env.setParam('OutputFlag', 0)\n    env.start()\n    with gp.Model(env=env) as m:\n            # Create an empty model\n        m = gp.Model('portfolio')\n        print(len(covar_matrix[0]))\n\n        # Add matrix variable for the epochs\n        x = m.addMVar(len(covar_matrix[0]))\n        z = m.addMVar(len(covar_matrix[0]))\n\n        #Exclusive\n        e1 = m.addMVar(19, vtype=GRB.BINARY) #XOR NN\n        e2 = m.addMVar(19, vtype=GRB.BINARY) #XOR RF\n        e3 = m.addMVar(9, vtype=GRB.BINARY)\n        portfolio_risk = x @ covar_matrix @ x\n        #portfolio_risk = gp.quicksum([z[i] * covar_matrix[i] * x[i] for i in range(len(covar_matrix[0]))])\n        m.setObjective(portfolio_risk, GRB.MINIMIZE)\n\n\n        # Fix budget with a constraint\n        #for i in range(len(covar_matrix[0])):\n            #m.addConstrs(z[i] == discount(x[i]),\"z_def_\"+str(i))\n            #m.addConstr((z[i]+1)*(1-x[i]) == 1,\"z_def_\"+str(i))\n\n        #m.addConstr(x.sum() >=.5, 'Lbudget')\n        #m.addConstr(x.sum() <=.6, 'Ubudget')\n        m.addConstr(x.sum() <=1, 'Ubudget')\n        epsilon=0.01\n        #or model in all_const: #try to add exclusivisty constrains\n         #   for i in range(len(x)):\n          #      m.addGenConstrIndicator(x[i], 1, t[i] >=  + epsilon, name=f\"lower_bound_{i}\")\n        #for i in range(0,18):\n            #m.addConstr(x[i]<=.6)\n        XOR_test(m)\n        #m.setParam(\"Presolve\", 0)\n        #make epochs exclusive\n        #epoch_exlusiviy(m)\n        #Add investment bounds on each epoch\n        #epoch_bounds(m)\n\n        # Verify model formulation\n        m.write('portfolio_selection_optimization.lp')\n\n        # Optimize model to find the minimum risk portfolio\n        m.optimize()","metadata":{"execution":{"iopub.status.busy":"2023-12-19T21:01:19.361088Z","iopub.execute_input":"2023-12-19T21:01:19.361547Z","iopub.status.idle":"2023-12-19T21:01:19.635031Z","shell.execute_reply.started":"2023-12-19T21:01:19.361514Z","shell.execute_reply":"2023-12-19T21:01:19.633888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def calc_return(e1,e2,x):\n    returns=0\n    for i in range(0,N*2-1):\n        if e1[i]==1 and x[i]!=0:\n            returns+=mean_acc[i-1]+x[i]*(mean_acc[i]-mean_acc[i-1])\n            #print(returns)\n    for i in range(N*2-1,4*N-2):\n        if e2[i-N*2-1]==1 and x[i]!=0:\n            returns+=mean_acc[i-1]+x[i]*(mean_acc[i]-mean_acc[i-1])\n            #print(mean_acc[i-1],mean_acc[i]-mean_acc[i-1])\n    for i in range(4*N-2,5*N-3):\n        if e2[i-4*N-2]==1 and x[i]!=0:\n            returns+=mean_acc[i-1]+x[i]*(mean_acc[i]-mean_acc[i-1])\n\n    return returns\n            ","metadata":{"execution":{"iopub.status.busy":"2023-12-19T21:01:25.902720Z","iopub.execute_input":"2023-12-19T21:01:25.903317Z","iopub.status.idle":"2023-12-19T21:01:25.916671Z","shell.execute_reply.started":"2023-12-19T21:01:25.903272Z","shell.execute_reply":"2023-12-19T21:01:25.914830Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"minrisk_volatility = sqrt(m.ObjVal)\n#minrisk_return = mean_acc @ x.X\n#minrisk_return = np.array(mean_acc[0:9]) @ e1 + np.array(mean_acc[9:]) @ e2\n#minrisk_return = mean_acc[0:9] @ e1.X + mean_acc[9:] @ e2.X\nminrisk_return = calc_return(e1.X,e2.X,x.X)\n\n#minrisk_return = aa(e1,e2,x,x1,x2)\n#print(e1.X)\n#print(e2.X)\n#print(x.X)\n#print(mean_acc)\nprint(minrisk_return)\n#print(x.X)\nprint(names)\n#print(norm_mean_time)\n#minrisk_return= mean_acc @ np.array([1 if x > 0 else 0 for x in x.X])\npd.DataFrame(data=np.append(x.X, [minrisk_volatility, minrisk_return]),\n             index=names + ['Volatility', 'Expected Return'],\n             columns=['Minimum Risk Portfolio'])","metadata":{"execution":{"iopub.status.busy":"2023-12-19T21:01:39.005265Z","iopub.execute_input":"2023-12-19T21:01:39.005731Z","iopub.status.idle":"2023-12-19T21:01:39.035146Z","shell.execute_reply.started":"2023-12-19T21:01:39.005691Z","shell.execute_reply":"2023-12-19T21:01:39.033937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def aa(e1,e2,x1,x2):\n    u1=mean_acc[0:9] @ e1\n    u2=mean_acc[9:] @ e2\n    l1=mean_shift[0:9] @ e1\n    l2=mean_shift[9:] @ e2\n    \n    #diff1 = np.full(9,u1-l1) # create the array\n    #diff2 = np.full(u2-l2,9)\n    a1=u1-l1\n    a2=u2-l2\n    s1 = gp.quicksum([a1*x1[i] for i in range(0,9)])\n    s2 = gp.quicksum([a2*x2[i] for i in range(0,9)])\n    \n    #print(s1)\n    returns=l1+l2+s1+s2\n    print(returns)\n    return returns","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:18:40.587718Z","iopub.execute_input":"2023-12-19T20:18:40.588190Z","iopub.status.idle":"2023-12-19T20:18:40.597072Z","shell.execute_reply.started":"2023-12-19T20:18:40.588155Z","shell.execute_reply":"2023-12-19T20:18:40.595880Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#portfolio_return = mean_acc @ x\nprint(mean_acc[0:9])\nportfolio_return = mean_acc[0:2*N-1] @ e1 + mean_acc[2*N-1:4*N-2] @ e2 + mean_acc[4*N-2:] @ e3\n#portfolio_return = aa(e1,e2,x1,x2)\n#portfolio_return = mean_acc @ var.setAttr([1 if x > 0 else 0 for x in x.X])\ntarget = m.addConstr(portfolio_return <= minrisk_return+.05, 'target')\ntarget = m.addConstr(portfolio_return >= minrisk_return-.05, 'target')\n#target = m.addConstr(portfolio_return == minrisk_return, 'target')\n\n# Solve for efficient frontier by varying target return\nfeasible=[]\ninfeasible=[]\nfrontier = np.empty((2,0))\nprint (mean_acc.max())\nfor r in np.linspace(mean_acc.min(),mean_acc.max(), 50):#mean_acc.max(), 25): #1.1, 25):#mean_acc.max(), 25):\n    #print(r,'##############')\n    target.rhs = r\n    #m.addConstr(portfolio_return == r, 'target')\n    m.optimize()\n    if m.status == GRB.OPTIMAL:\n        \n        frontier = np.append(frontier, [[sqrt(m.ObjVal)],[r]], axis=1)\n        feasible.append(r)\n    else:\n        infeasible.append(r)\nprint('feas',feasible)\nprint('infeas',infeasible)\n","metadata":{"execution":{"iopub.status.busy":"2023-12-19T20:59:30.539144Z","iopub.execute_input":"2023-12-19T20:59:30.539642Z","iopub.status.idle":"2023-12-19T20:59:32.321857Z","shell.execute_reply.started":"2023-12-19T20:59:30.539605Z","shell.execute_reply":"2023-12-19T20:59:32.319203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n#plt.figure(figsize=(10,10))\nprint(frontier.shape)\nprint(frontier)\nfig, ax = plt.subplots(figsize=(10,8))\n\n# Plot volatility versus expected return for individual stocks\nax.scatter(x=mean_var, y=mean_acc,\n           color='Blue', label='Individual Epochs')\nfor i, name in enumerate(names):\n    ax.annotate(name, (mean_var[i], mean_acc[i]))\n\n# Plot volatility versus expected return for minimum risk portfolio\nax.scatter(x=minrisk_volatility, y=minrisk_return, color='DarkGreen')\nax.annotate('Minimum\\nRisk\\nPortfolio', (minrisk_volatility, minrisk_return),\n            horizontalalignment='right')\n\n# Plot efficient frontier\nax.plot(frontier[0], frontier[1], label='Efficient Frontier', color='DarkGreen')\n\n# Format and display the final plot\nax.axis([frontier[0].min()*0, frontier[0].max()*5, mean_acc.min()*1, mean_acc.max()*5])\nax.set_xlabel('Volatility (standard deviation)')\nax.set_ylabel('Expected Return')\nax.legend()\nax.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-19T21:00:10.826846Z","iopub.execute_input":"2023-12-19T21:00:10.827250Z","iopub.status.idle":"2023-12-19T21:00:11.498607Z","shell.execute_reply.started":"2023-12-19T21:00:10.827220Z","shell.execute_reply":"2023-12-19T21:00:11.496409Z"},"trusted":true},"execution_count":null,"outputs":[]}]}