{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"02eaf587-125d-a4fe-62c7-3dac63b6c75f"},"source":"Following a similar recipe to lewis' R script (https://www.kaggle.com/cartographic/bosch-production-line-performance/bish-bash-xgboost), sampling the data to select features before running on the full set in order to stay within kaggle's memory limits. Here I add in the train_date data too.\n\nPlease feel free to fork and improve."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"b8ff0786-470f-f3de-e826-2334abc61cda"},"outputs":[],"source":"import numpy as np\nimport pandas as pd\nfrom xgboost import XGBClassifier\nfrom sklearn.metrics import matthews_corrcoef, roc_auc_score\nfrom sklearn.cross_validation import cross_val_score, StratifiedKFold\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n%matplotlib inline"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"eb108407-22f9-9716-1905-d6ebfd6a5b59"},"outputs":[],"source":"# I'm limited by RAM here and taking the first N rows is likely to be\n# a bad idea for the date data since it is ordered.\n# Sample the data in a roundabout way:\ndate_chunks = pd.read_csv(\"../input/train_date.csv\", index_col=0, chunksize=100000, dtype=np.float32)\nnum_chunks = pd.read_csv(\"../input/train_numeric.csv\", index_col=0,\n                         usecols=list(range(969)), chunksize=100000, dtype=np.float32)\nX = pd.concat([pd.concat([dchunk, nchunk], axis=1).sample(frac=0.05)\n               for dchunk, nchunk in zip(date_chunks, num_chunks)])\ny = pd.read_csv(\"../input/train_numeric.csv\", index_col=0, usecols=[0,969], dtype=np.float32).loc[X.index].values.ravel()\nX = X.values"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"7fa0923e-ddfa-d046-f7a8-842667883110"},"outputs":[],"source":"clf = XGBClassifier(base_score=0.005)\nclf.fit(X, y)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c6159cc6-47b5-e632-b0b6-2f8e0578876d"},"outputs":[],"source":"# threshold for a manageable number of features\nplt.hist(clf.feature_importances_[clf.feature_importances_>0])\nimportant_indices = np.where(clf.feature_importances_>0.005)[0]\nprint(important_indices)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"d2a8fa06-d180-d316-2f92-aef473339726"},"outputs":[],"source":"# load entire dataset for these features. \n# note where the feature indices are split so we can load the correct ones straight from read_csv\nn_date_features = 1156\nX = np.concatenate([\n    pd.read_csv(\"../input/train_date.csv\", index_col=0, dtype=np.float32,\n                usecols=np.concatenate([[0], important_indices[important_indices < n_date_features] + 1])).values,\n    pd.read_csv(\"../input/train_numeric.csv\", index_col=0, dtype=np.float32,\n                usecols=np.concatenate([[0], important_indices[important_indices >= n_date_features] + 1 - 1156])).values\n], axis=1)\ny = pd.read_csv(\"../input/train_numeric.csv\", index_col=0, dtype=np.float32, usecols=[0,969]).values.ravel()"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"041f6652-0318-ea87-79c6-8f1ecae047f4"},"outputs":[],"source":"from numpy import random as rd\ndef GA(pop, cf, ng=100, ns=100):\n    for gen in range(ng):\n        fits = np.array([[],[]]) # [ [params], [fitnesses]]\n        for specimen in ns:\n            params = [p * rd.normal( 1, .5*p) for p in pop]\n            for p in params:\n                fits[0].append(params[0])\n                fits[1].append(cf(params))\n            sorted_index = sort(range(len(fits[1])), key=lambda x:fits[1][x])\n            fits[0] = fits[0][sorted_index]\n            fits[1] = fits[1][sorted_index]\n        "},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"b693e51d-a933-c475-425e-1988236205e8"},"outputs":[],"source":"import operator\nimport random\n\nimport numpy\n\nfrom deap import base\nfrom deap import benchmarks\nfrom deap import creator\nfrom deap import tools\n\ncreator.create(\"FitnessMax\", base.Fitness, weights=(-1.0,))\ncreator.create(\"Particle\", list, fitness=creator.FitnessMax, speed=list, \n    smin=None, smax=None, best=None)\n\ndef generate(size, pmin, pmax, smin, smax):\n    part = creator.Particle(random.uniform(pmin, pmax) for _ in range(size)) \n    part.speed = [random.uniform(smin, smax) for _ in range(size)]\n    part.smin = smin\n    part.smax = smax\n    return part\n\ndef updateParticle(part, best, phi1, phi2):\n    u1 = (random.uniform(0, phi1) for _ in range(len(part)))\n    u2 = (random.uniform(0, phi2) for _ in range(len(part)))\n    v_u1 = map(operator.mul, u1, map(operator.sub, part.best, part))\n    v_u2 = map(operator.mul, u2, map(operator.sub, best, part))\n    part.speed = list(map(operator.add, part.speed, map(operator.add, v_u1, v_u2)))\n    for i, speed in enumerate(part.speed):\n        if speed < part.smin:\n            part.speed[i] = part.smin\n        elif speed > part.smax:\n            part.speed[i] = part.smax\n    part[:] = list(map(operator.add, part, part.speed))\n\ntoolbox = base.Toolbox()\ntoolbox.register(\"particle\", generate, size=2, pmin=-6, pmax=6, smin=-3, smax=3)\ntoolbox.register(\"population\", tools.initRepeat, list, toolbox.particle)\ntoolbox.register(\"update\", updateParticle, phi1=2.0, phi2=2.0)\ntoolbox.register(\"evaluate\", cost)\n\ndef main():\n    pop = toolbox.population(n=5)\n    stats = tools.Statistics(lambda ind: ind.fitness.values)\n    stats.register(\"avg\", numpy.mean)\n    stats.register(\"std\", numpy.std)\n    stats.register(\"min\", numpy.min)\n    stats.register(\"max\", numpy.max)\n\n    logbook = tools.Logbook()\n    logbook.header = [\"gen\", \"evals\"] + stats.fields\n\n    GEN = 1000\n    best = None\n\n    for g in range(GEN):\n        for part in pop:\n            part.fitness.values = toolbox.evaluate(part)\n            if not part.best or part.best.fitness < part.fitness:\n                part.best = creator.Particle(part)\n                part.best.fitness.values = part.fitness.values\n            if not best or best.fitness < part.fitness:\n                best = creator.Particle(part)\n                best.fitness.values = part.fitness.values\n        for part in pop:\n            toolbox.update(part, best)\n\n        # Gather all the fitnesses in one list and print the stats\n        logbook.record(gen=g, evals=len(pop), **stats.compile(pop))\n        print(logbook.stream)\n        print \"  Best so far: %s - %s\" % (best, best.fitness)\n\n    return pop, logbook, best\n\npop,logbook,best= main()"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0dd2d440-bb20-82c5-7519-d1eb24a7e510"},"outputs":[],"source":"clf = XGBClassifier(max_depth=5, base_score=0.005)\ncv = StratifiedKFold(y, n_folds=3)\npreds = np.ones(y.shape[0])\nfor i, (train, test) in enumerate(cv):\n    preds[test] = clf.fit(X[train], y[train]).predict_proba(X[test])[:,1]\n    print(\"fold {}, ROC AUC: {:.3f}\".format(i, roc_auc_score(y[test], preds[test])))\nprint(roc_auc_score(y, preds))\n\ndef cost(input_pop):\n\n   max_depth = input_pop[0] \n   base_score = input_pop[1]\n   learning_rate = input_pop[2] \n   n_estimators = input_pop[3] \n   subsample = input_pop[4]\n    \n   clf = XGBClassifier(max_depth=max_depth, base_score=base_score, learning_rate=learning_rate, n_estimators=n_estimators, subsample=subsample)\n   cv = StratifiedKFold(y, n_folds=3)\n   preds = np.ones(y.shape[0])\n   for i, (train, test) in enumerate(cv):\n       preds[test] = clf.fit(X[train], y[train]).predict_proba(X[test])[:,1]\n       #print(\"fold {}, ROC AUC: {:.3f}\".format(i, roc_auc_score(y[test], preds[test])))\n   return(roc_auc_score(y, preds))"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0f1e63c3-19fe-0a22-87c7-70eb54108694"},"outputs":[],"source":"\ndef cost(input_pop):\n\n   max_depth = input_pop[0] \n   base_score = input_pop[1]\n   learning_rate = input_pop[2] \n   n_estimators = input_pop[3] \n   subsample = input_pop[4]\n    \n   clf = XGBClassifier(max_depth=max_depth, base_score=base_score, learning_rate=learning_rate, n_estimators=n_estimators, subsample=subsample)\n   cv = StratifiedKFold(y, n_folds=3)\n   preds = np.ones(y.shape[0])\n   for i, (train, test) in enumerate(cv):\n       preds[test] = clf.fit(X[train], y[train]).predict_proba(X[test])[:,1]\n       #print(\"fold {}, ROC AUC: {:.3f}\".format(i, roc_auc_score(y[test], preds[test])))\n   return(roc_auc_score(y, preds))"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"be34f9d0-86b1-b0a8-257d-83561254b49c"},"outputs":[],"source":"# Thank donbeao http://stackoverflow.com/questions/22403982/how-to-minimize-a-function-using-deap\n\nimport operator\nimport random\n\nimport numpy\n\nfrom deap import base\nfrom deap import benchmarks\nfrom deap import creator\nfrom deap import tools\n\ncreator.create(\"FitnessMax\", base.Fitness, weights=(-1.0,))\ncreator.create(\"Particle\", list, fitness=creator.FitnessMax, speed=list, \n    smin=None, smax=None, best=None)\n\ndef generate(size, pmin, pmax, smin, smax):\n    part = creator.Particle(random.uniform(pmin, pmax) for _ in range(size)) \n    part.speed = [random.uniform(smin, smax) for _ in range(size)]\n    part.smin = smin\n    part.smax = smax\n    return part\n\ndef updateParticle(part, best, phi1, phi2):\n    u1 = (random.uniform(0, phi1) for _ in range(len(part)))\n    u2 = (random.uniform(0, phi2) for _ in range(len(part)))\n    v_u1 = map(operator.mul, u1, map(operator.sub, part.best, part))\n    v_u2 = map(operator.mul, u2, map(operator.sub, best, part))\n    part.speed = list(map(operator.add, part.speed, map(operator.add, v_u1, v_u2)))\n    for i, speed in enumerate(part.speed):\n        if speed < part.smin:\n            part.speed[i] = part.smin\n        elif speed > part.smax:\n            part.speed[i] = part.smax\n    part[:] = list(map(operator.add, part, part.speed))\n\ntoolbox = base.Toolbox()\ntoolbox.register(\"particle\", generate, size=2, pmin=-6, pmax=6, smin=-3, smax=3)\ntoolbox.register(\"population\", tools.initRepeat, list, toolbox.particle)\ntoolbox.register(\"update\", updateParticle, phi1=2.0, phi2=2.0)\ntoolbox.register(\"evaluate\", cost)\n\ndef main():\n    pop = toolbox.population(n=5)\n    stats = tools.Statistics(lambda ind: ind.fitness.values)\n    stats.register(\"avg\", numpy.mean)\n    stats.register(\"std\", numpy.std)\n    stats.register(\"min\", numpy.min)\n    stats.register(\"max\", numpy.max)\n\n    logbook = tools.Logbook()\n    logbook.header = [\"gen\", \"evals\"] + stats.fields\n\n    GEN = 3# 1000\n    best = None\n\n    for g in range(GEN):\n        for part in pop:\n            part.fitness.values = toolbox.evaluate(part)\n            if not part.best or part.best.fitness < part.fitness:\n                part.best = creator.Particle(part)\n                part.best.fitness.values = part.fitness.values\n            if not best or best.fitness < part.fitness:\n                best = creator.Particle(part)\n                best.fitness.values = part.fitness.values\n        for part in pop:\n            toolbox.update(part, best)\n\n        # Gather all the fitnesses in one list and print the stats\n        logbook.record(gen=g, evals=len(pop), **stats.compile(pop))\n        print(logbook.stream)\n        print (best, best.fitness)\n\n    return pop, logbook, best\n\npop,logbook,best= main()"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"eca4bd3d-586a-233d-6dcd-4380939e2b1e"},"outputs":[],"source":"# pick the best threshold out-of-fold\nthresholds = np.linspace(0.01, 0.99, 50)\nmcc = np.array([matthews_corrcoef(y, preds>thr) for thr in thresholds])\nplt.plot(thresholds, mcc)\nbest_threshold = thresholds[mcc.argmax()]\nprint(mcc.max())"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0e9aceeb-3b06-74af-3238-c500a33d5a71"},"outputs":[],"source":"# load test data\nX = np.concatenate([\n    pd.read_csv(\"../input/test_date.csv\", index_col=0, dtype=np.float32,\n                usecols=np.concatenate([[0], important_indices[important_indices<1156]+1])).values,\n    pd.read_csv(\"../input/test_numeric.csv\", index_col=0, dtype=np.float32,\n                usecols=np.concatenate([[0], important_indices[important_indices>=1156] +1 - 1156])).values\n], axis=1)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"a796c10a-d149-0608-350b-6adb5c733e8b"},"outputs":[],"source":"# generate predictions at the chosen threshold\npreds = (clf.predict_proba(X)[:,1] > best_threshold).astype(np.int8)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"e0574c4a-ece0-6545-fbc9-cafe21d0844c"},"outputs":[],"source":"# and submit\nsub = pd.read_csv(\"../input/sample_submission.csv\", index_col=0)\nsub[\"Response\"] = preds\nsub.to_csv(\"submission.csv.gz\", compression=\"gzip\")"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"6a735779-2db7-947a-944e-f67a54a63826"},"outputs":[],"source":""}],"metadata":{"_change_revision":0,"_is_fork":false,"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.5.2"}},"nbformat":4,"nbformat_minor":0}