{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":false},"cell_type":"markdown","source":"## About Simple\n[Simple](https://github.com/chrisstroemel/Simple) is very fast optimizer which likely to converge to global optimum.\n\n    Simple is a radically more scalable alternative to Bayesian Optimization. Like Bayesian Optimization, it is highly sample-efficient, converging to the global optimum in as few samples as possible. Unlike Bayesian Optimization, it has a runtime performance of O(log(n)) instead of O(n^3) (or O(n^2) with approximations), as well as a constant factor that is roughly three orders of magnitude smaller.\n\n## Optimizing ensemble weights\n\nThis kernel shows an example of ensemble weights optimization by using [Simple](https://github.com/chrisstroemel/Simple).\nYou might think that we already have `scipy.optimize.minimize`, and I was also using it in my last kernel ([Example of weighted ensemble](https://www.kaggle.com/daisukelab/example-of-weighted-ensemble)).\n\nBut I know that we also need to find best parameter for `scipy.optimize.minimize` to find the _better_ optimized results, or the results may not be the best.\n\nUnlike conventional optimization algorithm, I found Simlpe works fine without tricky tuning."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"trusted":true},"cell_type":"code","source":"# This is copy of Simple, just wanted to import..\n# Visit https://github.com/chrisstroemel/Simple/blob/master/Simple.py\n#  for find the latest code.\n\nfrom heapq import heappush, heappop, heappushpop\nimport numpy\nimport math\nimport time\nimport matplotlib.pyplot as plotter\n\nCAPACITY_INCREMENT = 1000\n\n\nclass _Simplex:\n\tdef __init__(self, pointIndices, testCoords, contentFractions, objectiveScore, opportunityCost, contentFraction, difference):\n\t\tself.pointIndices = pointIndices\n\t\tself.testCoords = testCoords\n\t\tself.contentFractions = contentFractions\n\t\tself.contentFraction = contentFraction\n\t\tself.__objectiveScore = objectiveScore\n\t\tself.__opportunityCost = opportunityCost\n\t\tself.update(difference)\n\n\tdef update(self, difference):\n\t\tself.acquisitionValue = -(self.__objectiveScore + (self.__opportunityCost * difference))\n\t\tself.difference = difference\n\n\tdef __eq__(self, other):\n\t\treturn self.acquisitionValue == other.acquisitionValue\n\n\tdef __lt__(self, other):\n\t\treturn self.acquisitionValue < other.acquisitionValue\n\nclass SimpleTuner:\n\tdef __init__(self, cornerPoints, objectiveFunction, exploration_preference=0.15):\n\t\tself.__cornerPoints = cornerPoints\n\t\tself.__numberOfVertices = len(cornerPoints)\n\t\tself.queue = []\n\t\tself.capacity = self.__numberOfVertices + CAPACITY_INCREMENT\n\t\tself.testPoints = numpy.empty((self.capacity, self.__numberOfVertices))\n\t\tself.objective = objectiveFunction\n\t\tself.iterations = 0\n\t\tself.maxValue = None\n\t\tself.minValue = None\n\t\tself.bestCoords = []\n\t\tself.opportunityCostFactor = exploration_preference #/ self.__numberOfVertices\n\t\t\t\n\n\tdef optimize(self, maxSteps=10):\n\t\tfor step in range(maxSteps):\n\t\t\t#print(self.maxValue, self.iterations, self.bestCoords)\n\t\t\tif len(self.queue) > 0:\n\t\t\t\ttargetSimplex = self.__getNextSimplex()\n\t\t\t\tnewPointIndex = self.__testCoords(targetSimplex.testCoords)\n\t\t\t\tfor i in range(0, self.__numberOfVertices):\n\t\t\t\t\ttempIndex = targetSimplex.pointIndices[i]\n\t\t\t\t\ttargetSimplex.pointIndices[i] = newPointIndex\n\t\t\t\t\tnewContentFraction = targetSimplex.contentFraction * targetSimplex.contentFractions[i]\n\t\t\t\t\tnewSimplex = self.__makeSimplex(targetSimplex.pointIndices, newContentFraction)\n\t\t\t\t\theappush(self.queue, newSimplex)\n\t\t\t\t\ttargetSimplex.pointIndices[i] = tempIndex\n\t\t\telse:\n\t\t\t\ttestPoint = self.__cornerPoints[self.iterations]\n\t\t\t\ttestPoint.append(0)\n\t\t\t\ttestPoint = numpy.array(testPoint, dtype=numpy.float64)\n\t\t\t\tself.__testCoords(testPoint)\n\t\t\t\tif self.iterations == (self.__numberOfVertices - 1):\n\t\t\t\t\tinitialSimplex = self.__makeSimplex(numpy.arange(self.__numberOfVertices, dtype=numpy.intp), 1)\n\t\t\t\t\theappush(self.queue, initialSimplex)\n\t\t\tself.iterations += 1\n\n\tdef get_best(self):\n\t\treturn (self.maxValue, self.bestCoords[0:-1])\n\n\tdef __getNextSimplex(self):\n\t\ttargetSimplex = heappop(self.queue)\n\t\tcurrentDifference = self.maxValue - self.minValue\n\t\twhile currentDifference > targetSimplex.difference:\n\t\t\ttargetSimplex.update(currentDifference)\n\t\t\t# if greater than because heapq is in ascending order\n\t\t\tif targetSimplex.acquisitionValue > self.queue[0].acquisitionValue:\n\t\t\t\ttargetSimplex = heappushpop(self.queue, targetSimplex)\n\t\treturn targetSimplex\n\t\t\n\tdef __testCoords(self, testCoords):\n\t\tobjectiveValue = self.objective(testCoords[0:-1])\n\t\tif self.maxValue == None or objectiveValue > self.maxValue: \n\t\t\tself.maxValue = objectiveValue\n\t\t\tself.bestCoords = testCoords\n\t\t\tif self.minValue == None: self.minValue = objectiveValue\n\t\telif objectiveValue < self.minValue:\n\t\t\tself.minValue = objectiveValue\n\t\ttestCoords[-1] = objectiveValue\n\t\tif self.capacity == self.iterations:\n\t\t\tself.capacity += CAPACITY_INCREMENT\n\t\t\tself.testPoints.resize((self.capacity, self.__numberOfVertices))\n\t\tnewPointIndex = self.iterations\n\t\tself.testPoints[newPointIndex] = testCoords\n\t\treturn newPointIndex\n\n\n\tdef __makeSimplex(self, pointIndices, contentFraction):\n\t\tvertexMatrix = self.testPoints[pointIndices]\n\t\tcoordMatrix = vertexMatrix[:, 0:-1]\n\t\tbarycenterLocation = numpy.sum(vertexMatrix, axis=0) / self.__numberOfVertices\n\n\t\tdifferences = coordMatrix - barycenterLocation[0:-1]\n\t\tdistances = numpy.sqrt(numpy.sum(differences * differences, axis=1))\n\t\ttotalDistance = numpy.sum(distances)\n\t\tbarycentricTestCoords = distances / totalDistance\n\n\t\teuclideanTestCoords = vertexMatrix.T.dot(barycentricTestCoords)\n\t\t\n\t\tvertexValues = vertexMatrix[:,-1]\n\n\t\ttestpointDifferences = coordMatrix - euclideanTestCoords[0:-1]\n\t\ttestPointDistances = numpy.sqrt(numpy.sum(testpointDifferences * testpointDifferences, axis=1))\n\n\n\n\t\tinverseDistances = 1 / testPointDistances\n\t\tinverseSum = numpy.sum(inverseDistances)\n\t\tinterpolatedValue = inverseDistances.dot(vertexValues) / inverseSum\n\n\n\t\tcurrentDifference = self.maxValue - self.minValue\n\t\topportunityCost = self.opportunityCostFactor * math.log(contentFraction, self.__numberOfVertices)\n\n\t\treturn _Simplex(pointIndices.copy(), euclideanTestCoords, barycentricTestCoords, interpolatedValue, opportunityCost, contentFraction, currentDifference)\n\n\tdef plot(self):\n\t\tif self.__numberOfVertices != 3: raise RuntimeError('Plotting only supported in 2D')\n\t\tmatrix = self.testPoints[0:self.iterations, :]\n\n\t\tx = matrix[:,0].flat\n\t\ty = matrix[:,1].flat\n\t\tz = matrix[:,2].flat\n\n\t\tcoords = []\n\t\tacquisitions = []\n\n\t\tfor triangle in self.queue:\n\t\t\tcoords.append(triangle.pointIndices)\n\t\t\tacquisitions.append(-1 * triangle.acquisitionValue)\n\n\n\t\tplotter.figure()\n\t\tplotter.tricontourf(x, y, coords, z)\n\t\tplotter.triplot(x, y, coords, color='white', lw=0.5)\n\t\tplotter.colorbar()\n\n\n\t\tplotter.figure()\n\t\tplotter.tripcolor(x, y, coords, acquisitions)\n\t\tplotter.triplot(x, y, coords, color='white', lw=0.5)\n\t\tplotter.colorbar()\n\n\t\tplotter.show()","execution_count":2,"outputs":[]},{"metadata":{"_uuid":"ad57f05ea53073c53e8ca8a12950a4e6df9b597e","_cell_guid":"9cae9b1f-c356-4cb8-9fd4-cd2386afa364","trusted":false},"cell_type":"markdown","source":"Here we have following results from 4 models:\n\n- predXX4valid.npy - validation results for estimating ensembled score and calculating weight balance.\n- predXX.npy - test results for building final submission results.\n\nAll prediction results are predicted probability for all classes:\n\n- Shape of one model prediction results for validation set (1232, 15)\n- Shape of one model prediction results for test set (1500, 15)"},{"metadata":{"_uuid":"3d873c58e5d24561dd7eb35d709ef2f8a1ed4757","_cell_guid":"a4582b44-545c-44f2-bfc8-ab93da228d2d","trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport keras\n\n# We are ensembling these four models.\nmodel_id_list_to_ensemble = ['D6', 'E', 'E4', 'H2']\n\n# And each results are stored here.\nexampledir = \"../input/tutasc-ensemble-example/\"\n\n# Load the training data, then make list of labels and dictionaries for converting labels and its numbers.\nraw_y_train = pd.read_csv('../input/acoustic-scene-2018/y_train.csv', sep=',')['scene_label'].tolist()\nlabels = sorted(list(set(raw_y_train)))\nlabel2int = {l:i for i, l in enumerate(labels)}\nint2label = {i:l for i, l in enumerate(labels)}\n\n# Make reference answers for valid as y_valid_ref, split from the raw training answers.\ny_train_int = [label2int[l] for l in raw_y_train]\nsplitlist = pd.read_csv('../input/acoustic-scene-2018/crossvalidation_train.csv', sep=',')['set'].tolist()\ny_valid_ref = keras.utils.to_categorical(np.array([y for i, y in enumerate(y_train_int) if splitlist[i] == 'test']))","execution_count":6,"outputs":[]},{"metadata":{"_uuid":"3021eb3f12c4f683ebc67db3e2a8cd502628610a","_cell_guid":"0fc0b2b3-b8dd-4a76-a76b-b8ad74f49982","trusted":true},"cell_type":"code","source":"# Load all predicted answers, and convert them to each result as class number.\nraw_valid_preds = [np.load(exampledir + 'preds%s4valid.npy' % e) for e in model_id_list_to_ensemble]\nref_valid_cls = [np.argmax(y) for y in y_valid_ref]\n\n# f() calculates combined prediction results by weighted ensembling.\ndef f(weights):\n    norm_weights = weights / np.sum(weights)\n    valid_preds = np.average(raw_valid_preds, axis=0, weights=norm_weights)\n    y_valid_pred_cls = [np.argmax(pred) for pred in valid_preds]\n    return y_valid_pred_cls\n\n# Calculates accuracy for particular weights\ndef acc_function(weights):\n    y_valid_pred_cls = f(weights)\n    n_eq = [result == ref for result, ref in zip(y_valid_pred_cls, ref_valid_cls)]\n    return np.sum(n_eq) / len(y_valid_pred_cls)\n\n# This defines the search area, and other optimization parameters.\noptimization_domain_vertices = [[0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]]\nnumber_of_iterations = 3000\nexploration = 0.01 # optional, default 0.15\n\n# Optimize weights\ntuner = SimpleTuner(optimization_domain_vertices, acc_function, exploration_preference=exploration)\ntuner.optimize(number_of_iterations)\nbest_objective_value, best_weights = tuner.get_best()\n\nprint('Best objective value =', best_objective_value)\nprint('Optimum weights =', best_weights)\nprint('Ensembled Accuracy (same as best objective value) =', acc_function(best_weights))","execution_count":4,"outputs":[]},{"metadata":{"_uuid":"22835f0e934ff86b4d597f64586eca670edb8eaf","_cell_guid":"82de8c2c-3494-4b3e-abff-d7546ba457f7","collapsed":true,"trusted":false},"cell_type":"markdown","source":"### Here's what if we use `scipy.optimize.minimize`\n\nWeights looks not to be optimized, they are all almost close to the initial 0.25."},{"metadata":{"trusted":true,"_uuid":"8d7c60634ef989eb7ad49b1e98371c30eef81d43"},"cell_type":"code","source":"from scipy import optimize\n\n# We need loss function to minimize\ndef loss_function(weights):\n    y_valid_pred_cls = f(weights)\n    n_lost = [result != ref for result, ref in zip(y_valid_pred_cls, ref_valid_cls)]\n    return np.sum(n_lost) / len(y_valid_pred_cls)\n\nopt_weights = optimize.minimize(loss_function,\n                                [1/len(model_id_list_to_ensemble)] * len(model_id_list_to_ensemble),\n                                constraints=({'type': 'eq','fun': lambda w: 1-sum(w)}),\n                                method= 'Nelder-Mead', #'SLSQP',\n                                bounds=[(0.0, 1.0)] * len(model_id_list_to_ensemble),\n                                options = {'ftol':1e-10},\n                            )['x']\n\nprint('Optimum weights =', opt_weights, 'with loss', loss_function(opt_weights))\nprint('Ensembled Accuracy =', acc_function(opt_weights))","execution_count":5,"outputs":[]},{"metadata":{"_uuid":"2da56bb934d4f7ed8cd8b5b77a682e35c32976f0"},"cell_type":"markdown","source":"## Summary\n\nSimple is better, at least in this case.\n\n| Optimizer | Ensembled Accuracy |\n|---------|----------------------------|\n| Simple | 0.8376623376623377 |\n| `scipy.optimize.minimize` | 0.823051948051948 |\n\n"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"cc47e4f91cf8efa43b9b317cb16e6ba6e203b9d0"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}},"nbformat":4,"nbformat_minor":1}