{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Quasar Classifier ","metadata":{}},{"cell_type":"code","source":"__author__ = \"Xinyue Sheng\"\n__copyright__ = \"Copyright (C) 2020 Xinyue Sheng\"\n__license__ = \"Public Domain\"\n__version__ = \"1.0\"","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### input data building & feature construction","metadata":{}},{"cell_type":"code","source":"from numpy.random import seed\n# seed(1)\nimport tensorflow\n# tensorflow.random.set_seed(1)\nimport numpy as np\nimport pandas as pd\nimport sys\nimport tensorflow.keras as ks \nfrom sklearn.gaussian_process import GaussianProcessRegressor\nfrom sklearn.gaussian_process.kernels import RBF, ConstantKernel as C\n\n\n\ndef difference_data(id_list = list, data = list, features = list, group_bound = list):\n\t'''\n\tthe difference between the neighboring magnitudes\n\tInputs:\n\t- id_list: the id list of all objects\n\t- data: the data from the preprocessed file\n\t- features: the targeted feature list. Other features will be emitted\n\t- group_bound: group mjd intervals after setting each group to the same size\n\tReturns:\n\t- new_data: modified data\n\n\t'''\n\n\tnew_data = data.copy()\n\tfor i in id_list:\n\t\tind_list = data[data.id == i].index.tolist()\n\t\tmjd_list = data['mjd'][ind_list].tolist()\n\t\tfor bound in group_bound:\n\t\t\tbound_right = bound[1]\n\t\t\tbound_left = bound[0]\n\t\t\tgroup_mjd = [x for x in mjd_list if bound_left<=x<=bound_right]\n\t\t\tif len(group_mjd)>1:\n\t\t\t\tn = 0\n\t\t\t\twhile n<len(group_mjd):\n\t\t\t\t\tif n == 0:\n\t\t\t\t\t\tidx = ind_list[0]+mjd_list.index(group_mjd[n])\n\t\t\t\t\t\tfor feature in features:\n\t\t\t\t\t\t\tnew_data.loc[idx,feature] = 0\n\t\t\t\t\telse:\n\t\t\t\t\t\tidx = ind_list[0]+mjd_list.index(group_mjd[n])\n\t\t\t\t\t\tfor feature in features:\n\t\t\t\t\t\t\tnew_data.loc[idx, feature] = data.loc[idx, feature] - data.loc[idx-1, feature]\t\n\n\t\t\t\t\tn +=1\n\n\treturn new_data \n\n\t\n\n\ndef normalization(id_list = list, data = list, features = list):\n\t'''\n\tRescaling (min-max normalization)\n\tInputs:\n\t- id_list: the id list of all objects\n\t- data: the data from the preprocessed file\n\t- features: the targeted feature list. Other features will be emitted\n\tReturns:\n\t- new_data: modified data\n\t'''\n\n\tfor i in id_list:\n\t\tind_list = data[data.id == i].index.tolist()\n\t\tnormalize_dict = {}\n\t\tfor feature in features:\n\t\t\tvalue_list = data[feature][ind_list].tolist()\n\t\t\t_range = np.max(value_list) - np.min(value_list)\n\t\t\t_min = np.min(value_list)\n\t\t\tdata.loc[ind_list, feature] = data.loc[ind_list, feature].map(lambda x: (x-_min)/_range)\n\t\t\t\n\n\treturn data\n\ndef standardization(id_list = list, data = list, features = list):\n\t'''\n\tStandardization(Z-score normalization)\n\tInputs:\n\t- id_list: the id list of all objects\n\t- data: the data from the preprocessed file\n\t- features: the targeted feature list. Other features will be emitted\n\tReturns:\n\t- new_data: modified data\n\t'''\n\n\tfor i in id_list:\n\t\tind_list = data[data.id == i].index.tolist()\n\t\tnormalize_dict = {}\n\t\tfor feature in features:\n\t\t\tvalue_list = data[feature][ind_list].tolist()\n\t\t\t_mean = np.nanmean(value_list)\n\t\t\t_std = np.nanstd(value_list)\n\t\t\tdata.loc[ind_list, feature] = data.loc[ind_list, feature].map(lambda x: (x-_mean)/_std)\n\n\treturn data\n\n\ndef GPR(id_list = list, data = list, group_bound = list, features = list, withErr = False):\n\t'''\n\tGuassian Process regression.\n\tInputs:\n\t- id_list: the id list of all objects\n\t- data: the data from the preprocessed file\n\t- group_bound: group mjd intervals after setting each group to the same size\n\t- features: the targeted feature list. Other features will be emitted\n\t- withErr: True means add the magnitude errors into GP regression calculation\n\tReturns:\n\t- new_data: modified data\n\t'''\n\tprint(\"start Gaussian process regression...\")\n\n\tnew_data = pd.DataFrame(columns = ['id','mjd']+features+['type'])\n\n\tfor i in id_list:\n\t\tind_list = data[data.id == i].index.tolist()\n\t\tmjd_list = data['mjd'][ind_list].tolist()\n\t\tfill_mjd_list = []\n\t\tobj_type = data['type'][ind_list[0]]\n\t\tinsert_dataframe = pd.DataFrame(columns = ['id','mjd']+features+['type'])\n\t\tpred_dict = {}\n\t\tfor f in features:\n\t\t\tpred_dict[f] = []\n\t\tfor bound in group_bound:\n\t\t\tbound_left = bound[0]\n\t\t\tbound_right = bound[1]\n\t\t\tgroup_mjd = [x for x in mjd_list if bound_left<=x<=bound_right]\n\t\t\tgroup_mjd_idx = [mjd_list.index(x)+ind_list[0] for x in mjd_list if bound_left<=x<=bound_right]\n\t\t\tif len(group_mjd)>0:\n\t\t\t\tX = np.atleast_2d(group_mjd).T\n\t\t\t\tfill_mjd = np.atleast_2d(np.array([t for t in range(group_mjd[0],group_mjd[-1]+1,1)])).T\n\t\t\t\tgroup_fill_mjd = [t for t in range(group_mjd[0],group_mjd[-1]+1,1)]\n\t\t\t\tfill_mjd_list += group_fill_mjd\n\n\t\t\t\tfor f in features:\n\t\t\t\t\ty = np.array(data[f][group_mjd_idx].tolist()).ravel()\n\t\t\t\t\ty_given = data[f][group_mjd_idx].tolist()\n\t\t\t\t\tnoise = np.array(data[f+'_error'][group_mjd_idx].tolist())\n\t\t\t\t\tif withErr == True:\n\t\t\t\t\t\ty += noise\n\t\t\t\t\tkernel = C(constant_value=1.0, constant_value_bounds=(1e-3, 1e3)) * RBF(length_scale=10, length_scale_bounds=(1e-1, 1e1))\n\t\t\t\t\tgp = GaussianProcessRegressor(kernel=kernel, alpha = noise)\n\t\t\t\t\tgp.fit(X,y)\n\t\t\t\t\ty_pred = list(gp.predict(fill_mjd))\n\t\t\t\t\tnew_filled_values = []\n\t\t\t\t\t\n\t\t\t\t\tn = 0\n\t\t\t\t\tm = 0\n\t\t\t\t\twhile(n<len(group_fill_mjd)):\n\t\t\t\t\t\tif group_fill_mjd[n] not in group_mjd:\n\t\t\t\t\t\t\tnew_filled_values.append(y_pred[n])\n\t\t\t\t\t\telse:\n\t\t\t\t\t\t\tnew_filled_values.append(y_given[m])\n\t\t\t\t\t\t\tm +=1\n\t\t\t\t\t\tn +=1\n\t\t\t\t\tpred_dict[f] += new_filled_values\n\n\t\tinsert_dataframe['id'] = [i for t in fill_mjd_list]\n\t\tinsert_dataframe['mjd'] = fill_mjd_list\n\t\tfor f in features:\n\t\t\tinsert_dataframe[f] = pred_dict[f]\n\t\tinsert_dataframe['type'] = [obj_type for t in fill_mjd_list]\n\t\tnew_data = new_data.append(insert_dataframe, ignore_index=True)\n\n\tprint(\"GPR is Done.\")\n\n\treturn new_data\n\n\n\ndef remove_alone_mjd(id_list = list, data = list, check_delta = 300, min_size = 3):\n\t'''\n\tThis function is used to remove the data points whose mjd is far away from other data points.\n\tThese alone data pionts is not useful for group format and season format input data.\n\tInput:\n\t- id_list: the id list of all objects\n\t- data: the data from the preprocessed file\n\t- check_delta: if the difference between two neighboring data points' mjd is larger than this value, the earlier one will be removed.\n\t- min_size: the minimal size of the group before padding\n\tReturns:\n\t- data: the modified data\n\t'''\n\n\trecord_row = []\n\n\tfor i in id_list:\n\t\tind_list = data[data.id == i].index.tolist()\n\t\tmjd_exit_list = data['mjd'][ind_list].tolist()\n\t\twarn_index = []\n\t\tn = 0\n\t\twhile n+1<len(mjd_exit_list):\n\t\t\tdelta = mjd_exit_list[n+1]-mjd_exit_list[n]\n\n\t\t\tif delta>check_delta:\n\t\t\t\twarn_index.append(ind_list[n])\n\t\t\tif n == 0 and delta>260:\n\t\t\t\trecord_row.append(ind_list[0])\n\t\t\tn +=1\n\t\tif len(warn_index)!=0:\n\t\t\tt = 0\n\t\t\twhile t+1<len(warn_index):\n\t\t\t\tif t == 0:\n\t\t\t\t\tif warn_index[0]-ind_list[0]<=min_size: record_row+=list(range(ind_list[0],warn_index[0]+1))\n\n\t\t\t\tif warn_index[t+1]-warn_index[t]<=min_size: \n\t\t\t\t\trecord_row += list(range(warn_index[t]+1, warn_index[t+1]+1))\n\t\t\t\tt +=1\n\t\t\tif len(warn_index)==1:\n\t\t\t\tif warn_index[0]-ind_list[0]<=min_size: record_row+=list(range(ind_list[0],warn_index[0]+1))\n\n\n\tnew_data = data.drop(index=list(set(record_row)), axis=0).reset_index(drop=True)\n\n\treturn new_data\n\n\n\n\ndef combine_narrow_mjd(id_list = list, data=list, check_delta = 0.7):\n\t'''\n\tThis function is used for combine the data points whose mjds are closed. \n\tFor example, 2 data points in one mjd.\n\tInput:\n\t- id_list: the id list of all objects\n\t- data: the data from the preprocessed file\n\t- check_delta: if the difference between two neighboring data points' mjd is smaller than this value, the later one will be removed.\n\tReturns:\n\t- data: the modified data\n\n\t'''\n\trecord_row = []\n\tfor i in id_list:\n\t\tind_list = data[data.id == i].index.tolist()\n\t\tmjd_exit_list = data['mjd'][ind_list].tolist()\n\t\tn = 0\n\t\twhile n+1<len(mjd_exit_list):\n\t\t\tdelta = mjd_exit_list[n+1]-mjd_exit_list[n]\n\t\t\tif delta<check_delta:\n\t\t\t\trecord_row.append(ind_list[n])\n\t\t\tn +=1\n\tnew_data = data.drop(index=list(set(record_row)), axis=0).reset_index(drop=True)\n\t\t\n\treturn new_data\n\n\ndef make_int_mjd(data=list):\n\t'''\n\tchange the mjd's format to integer\n\t'''\n\t\n\tdata.loc[:,'mjd'] = data.loc[:,'mjd'].map(lambda x: int(x))\n\treturn data\n\n\ndef group_observations(data = list, remove_check_delta = 233, combine_check_delta = 0.7, sequence_length = None, group_size = None, min_size = 3):\n\t'''\n\tThis function is used to extract the main features of the data in order to group the observations\n\tIt is called in format_group function.\n\tInputs:\n\t- data: the data from the preprocessed file\n\t- remove_check_delta: used in remove_alone_mjd function. if the difference between two neighboring data points' mjd is larger than this value, the earlier one will be removed.\n\t- combine_check_delta: used in combine_narrow_mjd function. if the difference between two neighboring data points' mjd is smaller than this value, the later one will be removed.\n\t- sequence_length: the number of groups in one sequence\n\t- group_size: the fixed size of one group after padding\n\t- min_size: the minimal size of the group before padding\n\tReturns:\n\t- data: modified data\n\t- max_delta: the maximal gap between neighboring mjd data points\n\t- num_group: the fixed/designed number of group among all objects\n\t- max_group: the largest number of groups among all objects\n\t- min_group: the smallest number of groups among all objects\n\t- group_size: the fixed/designed size of group\n\t- modifid_group_bound: the modified group mjd intervals after setting each group to the same size\n\t'''\n\n\tcombine = data.groupby(data['id']) \n\tid_list = []\n\tfor _id, group in combine:\n\t    id_list.append(_id)\n\n\tdata =remove_alone_mjd(id_list, data, remove_check_delta, min_size)\n\n\tdata = combine_narrow_mjd(id_list, data, combine_check_delta)\n\n\tdelta_time = []\n\tnum_group = []\n\tlast_mjd = []\n\tsize_group = []\n\tmjd_group = []\n\n\tfor i in id_list:\n\t\tind_list = data[data.id == i].index.tolist()\n\t\tmjd_exit_list = data['mjd'][ind_list].tolist()\n\t\tm = 0\n\t\tobj_delta_time = []\n\t\twhile m<len(mjd_exit_list)-1:\n\t\t\tdelta = mjd_exit_list[m+1]-mjd_exit_list[m]\n\t\t\tobj_delta_time.append(delta)\n\t\t\tm +=1\n\t\tdelta_time.append(obj_delta_time)\n\n\t\t# the real time length of observations for each object\n\t\tlast_mjd.append(int(mjd_exit_list[-1]))\n\n\t\tobs_time = mjd_exit_list[-1]-mjd_exit_list[0]\n\n\t\t# note: the num of groups is not equal to the observation year\n\t\tobj_num_group = len([x for x in obj_delta_time if x>=remove_check_delta])+1\n\t\tnum_group.append(obj_num_group) \n\t\t\n\t\t# calculate the mjd in each group for each object\n\t\tidx_group_gap = [obj_delta_time.index(x) for x in obj_delta_time if x>=remove_check_delta]\n\n\t\t# print(idx_group_gap)\n\n\t\tn = 0\n\t\tobj_group_size = []\n\t\tobj_group_mjd = []\n\t\tif len(idx_group_gap)>1:\n\t\t\twhile n<len(idx_group_gap):\n\t\t\t\tif n == 0:\n\t\t\t\t\tobj_group_size.append(sum(obj_delta_time[:idx_group_gap[n]]))\n\t\t\t\t\tobj_group_mjd.append(mjd_exit_list[:idx_group_gap[n]+1])\n\t\t\t\telif n == len(idx_group_gap)-1:\n\t\t\t\t\tobj_group_size.append(sum(obj_delta_time[idx_group_gap[n]+1:]))\n\t\t\t\t\tobj_group_mjd.append(mjd_exit_list[idx_group_gap[n]+1:])\n\t\t\t\telse:\n\t\t\t\t\tobj_group_size.append(sum(obj_delta_time[idx_group_gap[n]+1:idx_group_gap[n+1]]))\n\t\t\t\t\tobj_group_mjd.append(mjd_exit_list[idx_group_gap[n]+1:idx_group_gap[n+1]+1])\n\n\t\t\t\tn +=1\n\t\telif len(idx_group_gap)==1:\n\t\t\tobj_group_size += [sum(obj_delta_time[:idx_group_gap[0]]), sum(obj_delta_time[idx_group_gap[0]+1:])]\n\t\t\tobj_group_mjd +=[mjd_exit_list[:idx_group_gap[0]+1], mjd_exit_list[idx_group_gap[0]+1:]]\n\n\t\tif len(idx_group_gap)==0:\n\t\t\tobj_group_size.append(sum(obj_delta_time))\n\t\t\tobj_group_mjd.append(mjd_exit_list)\n\n\n\t\tsize_group.append(obj_group_size)\n\t\tmjd_group.append(obj_group_mjd)\n\t\n\tmax_delta = int(np.max([np.max(x) for x in delta_time]))\n\tmax_group = np.max(num_group)\n\tmin_group = np.min(num_group)\n\tmax_group_size = int(np.max([np.max(x) for x in size_group]))+1\n\tmean_group_size = int(np.mean([np.mean(x) for x in size_group]))+1\n\n\tprint(\"maximal group number: \",max_group)\n\tprint(\"minimal group number: \",min_group)\n\tprint(\"maximal group size: \", max_group_size+1)\n\tprint(\"mean group size: \", mean_group_size+1)\n\n\t# this mjd is regarded as the last mjd for every object\n\tbase_mjd = np.max(last_mjd)\n\tprint(\"the latest mjd is \", int(base_mjd)+1)\n\n\n\t#calculate the gap time between neigboring groups\n\tn = 0\n\twhile n<len(id_list):\n\t\tgroup = num_group[n]\n\t\tobj_delta_time = delta_time[n]\n\n\t\tif group > 1 : \n\t\t\twhile len([x for x in obj_delta_time if x>=max_delta]) < group-1 and max_delta > 0:\n\t\t\t\tmax_delta -= 1\n\t\tn +=1\n\n\tsuit_delta = int(max_delta-1)\n\tprint(\"the suitable gap time between each group: \"+str(suit_delta))\n\n\t\n\n\t# calculate the boundary mjd of each group\n\tn = 1\n\tgroup_bound = []\n\tif sequence_length == None:\n\t\tsequence_length = max_group\n\tbound_left  = base_mjd - mean_group_size - suit_delta\n\tbound_right = base_mjd\n\n\n\twhile n<=sequence_length:\n\t\tobj_bound_pre = []\n\t\tobj_bound_post = []\n\t\t\n\t\tfor x in mjd_group:\n\t\t\tmatch = [t for t in x if bound_left<=t[0]<=bound_right] \n\t\t\tif len(match)==1:\n\t\t\t\tobj_bound_pre.append(match[0][0])\n\t\t\t\tobj_bound_post.append(match[0][-1])\n\n\t\tif len(obj_bound_pre)>0 and len(obj_bound_post)>0:\n\t\t\tgroup_bound.append([np.min(obj_bound_pre),np.max(obj_bound_post)])\n\t\t\tgap = np.max(obj_bound_post)-np.min(obj_bound_pre)\n\t\t\t\n\t\t\tbound_left = bound_left - gap - suit_delta\n\t\t\tbound_right = bound_right - gap - suit_delta\n\t\telse:\n\t\t\tbound_left = bound_left - mean_group_size - suit_delta\n\t\t\tbound_right = bound_right - mean_group_size - suit_delta\n\n\t\tn +=1\n\t\n\t# print(group_bound)\n\n\t# modify the group size, and make it equal\n\tmodified_group_bound = []\n\tif group_size == None:\n\t\tgroup_size = max_group_size\n\tfor x in group_bound:\n\t\tnew_bound = int(x[1]-group_size)\n\t\tmodified_group_bound.append([new_bound,int(x[1])])\n\n\n\t\n\tmodified_group_bound.reverse() \n\tprint(modified_group_bound)\n\n\n\treturn data, max_delta, num_group, max_group, min_group, group_size, modified_group_bound\n\n\ndef format_season(data = list, train_id = list, test_id = list, features = list, sequence_length = None, remove_check_delta = 233, combine_check_delta = 0.7, group_size = None, min_size = 3, preprocess = 's', set_GPR = False):\n\t'''\n\tThis function regard one season of observations for one object as one sequence.\n\tTherefore, one objects can offer several sequences.\n\tIn this way, we could figure out the season trend of the variable object.\n\tThe vector in each timestep is the combination of targeted bands.\n\tInputs:\n\t- data: the data from the preprocessed file\n\t- train_id: the id of train objects\n\t- test_id: the id of test objects\n\t- features: the targeted feature list. Other features will be emitted.\n\t- sequence_length: the number of data points in one sequence\n\t- remove_check_delta: used in remove_alone_mjd function. if the difference between two neighboring data points' mjd is larger than this value, the earlier one will be removed.\n\t- combine_check_delta: used in combine_narrow_mjd function. if the difference between two neighboring data points' mjd is smaller than this value, the later one will be removed.\n\t- group_size: the fixed size of one group after padding\n\t- min_size: the minimal size of the group before padding\n\t- preprocess: choose the method of preprocessing: \n\t\t\t\t  's': standardization\n\t\t\t\t  'n': normalization\n\t\t\t\t  'd': difference_data - set the difference between neighboring values as the preprocessed value\n\t- set_GPR: True means using Gaussian process regression to the magnitude data.\n\tReturns:\n\t- X_train/X_test: the train/test data with 3 dimensions: \n\t\tobjects, the number of observations in an objects, the vector in 1 observation\n\t- Y_train/Y_test: the label of train/test data. \n\t'''\n\n\n\tdata, max_delta, num_group, max_group, min_group, group_size, group_bound = group_observations(data = data, remove_check_delta = remove_check_delta, \n\t\tcombine_check_delta = combine_check_delta, sequence_length = sequence_length, group_size = group_size,  min_size = min_size)\n\t\n\tdata = make_int_mjd(data)\n\n\tcombine = data.groupby(data['id']) \n\tid_list = []\n\tfor _id, group in combine:\n\t    id_list.append(_id)\n\n\tif set_GPR == True:\n\t\tdata = GPR(id_list, data, group_bound, features, withErr = False)\n\n\tif preprocess == 's':\n\t\tdata = standardization(id_list, data, features)\n\telif preprocess == 'n':\n\t\tdata = normalization(id_list, data, features)\n\telif preprocess == 'd':\n\t\tdata = difference_data(id_list, data, features, group_bound)\n\telse:\n\t\tprint('ERROR: wrong preprocess method input!')\n\n\tobj_data = []\n\tX_train = []\n\tY_train = []\n\tX_test = []\n\tY_test = []\n\n\n\tfor i in id_list:\n\t\tind_list = data[data.id == i].index.tolist()\n\t\tmjd_list = data['mjd'][ind_list].tolist()\n\t\tn = 0\n\t\t\n\t\tfor bound in group_bound:\n\t\t\tbound_left = bound[0]\n\t\t\tbound_right = bound[1]\n\t\t\tgroup_mjd = [x for x in mjd_list if bound_left<=x<=bound_right]\n\t\t\tif len(group_mjd)>=10:\n\t\t\t\tband_data = []\n\t\t\t\twhile bound_left<=bound_right:\n\t\t\t\t\tif bound_left in group_mjd:\n\t\t\t\t\t\tidx = mjd_list.index(bound_left)+ind_list[0]\n\t\t\t\t\t\tband_data.append([data[x][idx] for x in features])\n\t\t\t\t\telse:\n\t\t\t\t\t\tband_data.append([0 for x in features])\n\t\t\t\t\tbound_left +=1\n\t\t\t\tif i in train_id:\n\t\t\t\t\tX_train.append(band_data)\n\t\t\t\t\tY_train.append(data['type'][ind_list[0]])\n\t\t\t\telif i in test_id:\n\t\t\t\t\tX_test.append(band_data)\n\t\t\t\t\tY_test.append(data['type'][ind_list[0]])\n\n\n\n\treturn X_train, Y_train, X_test, Y_test\n\n\n\n\n\ndef format_group(data = list, train_id = list, test_id = list, features = list, sequence_length = None, remove_check_delta = 233, combine_check_delta = 0.7, group_size = None, min_size = 3, preprocess = 's', set_GPR = False):\n\t'''\n\tThis function groups the observations for one objects into groups.\n\tWe assume that the data are discontinuous and that a large-scale observation is made once a year.\n\tThen, we could regard one year's observations as a group.\n\tIn this way, we could figure out the year trend of the variable object.\n\tThe length of the sequence is the number of groups.\n\tThe vector in each timestep is the combination of a group's observation(one band's magnitudes differences)\n\tInputs:\n\t- data: the data from the preprocessed file\n\t- train_id: the id of train objects\n\t- test_id: the id of test objects\n\t- features: the targeted feature list. Other features will be emitted.\n\t- sequence_length: the number of data points in one sequence\n\t- remove_check_delta: used in remove_alone_mjd function. if the difference between two neighboring data points' mjd is larger than this value, the earlier one will be removed.\n\t- combine_check_delta: used in combine_narrow_mjd function. if the difference between two neighboring data points' mjd is smaller than this value, the later one will be removed.\n\t- group_size: the fixed size of one group after padding\n\t- min_size: the minimal size of the group before padding\n\t- preprocess: choose the method of preprocessing: \n\t\t\t\t  's': standardization\n\t\t\t\t  'n': normalization\n\t\t\t\t  'd': difference_data - set the difference between neighboring values as the preprocessed value\n\t- set_GPR: True means using Gaussian process regression to the magnitude data.\n\t\n\tReturns:\n\t- X_train/X_test: the train/test data with 3 dimensions: \n\t\tobjects, the number of observations in an objects, the vector in 1 observation\n\t- Y_train/Y_test: the label of train/test data. \n\n\t'''\n\n\tdata, max_delta, num_group, max_group, min_group, group_size, group_bound = group_observations(data = data, remove_check_delta = remove_check_delta, \n\t\tcombine_check_delta = combine_check_delta, sequence_length = sequence_length, group_size = group_size, min_size = min_size)\n\n\tdata = make_int_mjd(data)\n\n\tcombine = data.groupby(data['id']) \n\tid_list = []\n\tfor _id, group in combine:\n\t    id_list.append(_id)\n\n\tif set_GPR == True:\n\t\tdata = GPR(id_list, data, group_bound, features, withErr = False)\n\n\t# normalization/ differences\n\tif preprocess == 's':\n\t\tdata = standardization(id_list, data, features)\n\telif preprocess == 'n':\n\t\tdata = normalization(id_list, data, features)\n\telif preprocess == 'd':\n\t\tdata = difference_data(id_list, data, features, group_bound)\n\telse:\n\t\tprint('ERROR: wrong preprocess method input!')\n\n\n\tobj_data = []\n\tX_train = []\n\tY_train = []\n\tX_test = []\n\tY_test = []\n\n\n\n\tfor i in id_list:\n\t\tind_list = data[data.id == i].index.tolist()\n\t\tmjd_list = data['mjd'][ind_list].tolist()\n\t\tfeature = features[0]\n\t\tvalue_dict = {}\n\t\tfor f in features:\n\t\t\tvalue_dict[f] = data[f][ind_list].tolist()\n\t\t# value_list = data[feature][ind_list].tolist()\n\n\t\tn = 0\n\t\tobj_data = []\n\n\t\tfor bound in group_bound:\n\t\t\tbound_right = bound[1]\n\t\t\tbound_left = bound[0]\n\t\t\tgroup_mjd = [x for x in mjd_list if bound_left<=x<=bound_right]\n\n\t\t\tif len(group_mjd)==0:\n\t\t\t\tsequence_data = [0 for x in range(0,(group_size+1)*len(features))]\n\t\t\t\tobj_data.append(sequence_data)\n\t\t\telse:\n\t\t\t\tsequence_data = []\n\t\t\t\twhile bound_left<=bound_right:\n\t\t\t\t\tif bound_left in group_mjd:\n\t\t\t\t\t\tidx = mjd_list.index(bound_left)\n\t\t\t\t\t\tfor f in features:\n\t\t\t\t\t\t\tsequence_data.append(value_dict[f][idx])\n\t\t\t\t\telse:\n\t\t\t\t\t\tfor f in features:\n\t\t\t\t\t\t\tsequence_data.append(0)\n\t\t\t\t\t\n\t\t\t\t\tbound_left +=1\n\n\t\t\t\tobj_data.append(sequence_data)\n\n\n\t\tif i in train_id:\n\t\t\tX_train.append(obj_data)\n\t\t\tY_train.append(data['type'][ind_list[0]])\n\t\telif i in test_id:\n\t\t\tX_test.append(obj_data)\n\t\t\tY_test.append(data['type'][ind_list[0]])\n\n\t\n\n\treturn X_train, Y_train, X_test, Y_test\n\n\n\n\n\n\ndef format_simple(data = list, train_id = list, test_id = list, features = list, preprocess = 's'):\n\t'''\n\tThis is the simple way of generating the input data.\n\tThe length of the sequence is the biggest number of observations among all objects.\n\tThe vector in each timestep is the combination of targeted features. eg.v1 = [u,g,r]\n\tInputs:\n\t- data: the data from the preprocessed file\n\t- train_id: the id of train objects\n\t- test_id: the id of test objects\n\t- features: the targeted feature list. Other features will be emitted.\n\t- preprocess: choose the method of preprocessing: \n\t\t\t  's': standardization\n\t\t\t  'n': normalization\n\t\t\t  'd': difference_data - set the difference between neighboring values as the preprocessed value\n\tReturns: \n\t- X_train/X_test: the train/test data with 3 dimensions: \n\t\tobjects, the number of observations in an objects, the vector in 1 observation\n\t- Y_train/Y_test: the label of train/test data. \n\t- test_list: the final object id order in the test set\n\t- test_type: the final object type order in the test set\n\t'''\n\t\n\tcombine = data.groupby(data['id'])\n\n\tid_list = []\n\tfor _id, group in combine:\n\t    id_list.append(_id)\n\n\tif preprocess == 's':\n\t\tdata = standardization(id_list, data, features)\n\telif preprocess == 'n':\n\t\tdata = normalization(id_list, data, features)\n\telif preprocess == 'd':\n\t\tdata = difference_data(id_list, data, features, group_bound=[[np.min(data['mjd']),np.max(data['mjd'])]])\n\telse:\n\t\tprint('ERROR: wrong preprocess method input!')\n\n\tdata = data.fillna(0)\n\tprint(data[:10])\n\n\tobj_data = []\n\tX_train = []\n\tY_train = []\n\tX_test = []\n\tY_test = []\n\t\n\n\tlast_id = None\n\tlast_type = None\n\tn = 0\n\twhile n < len(data):\n\t\tobj_id = data['id'][n]\n\t\tif n == len(data)-1:\n\t\t\tobj_data.append([float(data[x][n]) for x in features])\n\t\tif obj_id != last_id or n == len(data)-1:\n\t\t\tif last_id in train_id:\n\t\t\t\tX_train.append(obj_data)\n\t\t\t\tif last_type != None:\n\t\t\t\t\tY_train.append(last_type)\n\t\t\telif last_id in test_id:\n\t\t\t\tX_test.append(obj_data)\t\n\t\t\t\tif last_type != None:\n\t\t\t\t\tY_test.append(last_type)\n\t\t\t\t\t\n\n\t\t\tobj_data = []\n\t\t\tlast_id = obj_id\n\t\t\t\n\t\t\tif 'type' in list(data.columns):\n\t\t\t\tlast_type = data['type'][n]\n\t\t\telse:\n\t\t\t\tlast_type = data['label'][n]\n\t\t\n\t\tobj_data.append([float(data[x][n]) for x in features])\n\t\tn = n + 1\n\n\n\treturn X_train, Y_train, X_test, Y_test\n\n\ndef load_data(path, test_fraction = 0.2, seed = None, features = list, set_format = 'group', preprocess = 's', group_size = None,\n\tmin_group =3, sequence_length = None, remove_check_delta = 233, combine_check_delta = 0.7, set_GPR = False): \n\t'''\n\tgenerate train/test set, divide the data into train and test sets.\n\tInputs:\n\t- path: the processed file address\n\t- test fraction: float, the fraction of test set among all data\n\t- seed: int, when the seed is the same, the random result is the same\n\t- features: list, the features which will be fed into the neural network\n\t- set_format: choose an input format: simple, group, season\n\t- preprocess: choose the method of preprocessing: \n\t\t\t\t  's': standardization\n\t\t\t\t  'n': normalization\n\t\t\t\t  'd': difference_data - set the difference between neighboring values as the preprocessed value\n\t- min_group: the minimal size of the group before padding\n\t- remove_check_delta: used in remove_alone_mjd function. if the difference between two neighboring data points' mjd is larger than this value, the earlier one will be removed.\n\t- combine_check_delta: used in combine_narrow_mjd function. if the difference between two neighboring data points' mjd is smaller than this value, the later one will be removed.\n\t- set_GPR: True means using Gaussian process regression to the magnitude data.\n\tReturns:\n\t- X_train/X_test: the train/test data with 3 dimensions: \n\t\tobjects, the number of observations in an objects, the vector in 1 observation\n\t- Y_train/Y_test: the label of train/test data. \n\t- length_train/test: the number of sequences in the train/test set\n\t- time_sequence: the number of data points in one sequence\n\t- input_dim: the dimension of input data(train/test sets)\n\t- num_classes: the number of classes/labels\n\n\t'''\n\n\tprint('Start constructing the input data...')\n\t\n\tlast_id = None\n\tids = []\n\tpadding = False\n\n\t# Used to generate a specified random number\n\tnp.random.seed(seed)\n\n\t# read the csv file\n\tdata = pd.read_csv(path)  \n\n\t# select columns\n\tif 'type' in list(data.columns):\n\t\tdata = data[['id','mjd','type']+features]\n\telse:\n\t\tdata = data[['id','mjd','label']+features]\n\n\t#remove NaN band values, replace NaN with 0\n\tdata = data.dropna(subset = features, axis = 'rows', how='all').reset_index()\n\tdata = data.fillna(0)\n\n\t# combine the data grouped by its id\n\tcombine = data.groupby(data['id']) \n\n\t# create the id list\n\tid_list = []\n\tfor _id, group in combine:\n\t\tid_list.append(_id)\n\n\t# shuffle the index\n\tnp.random.shuffle(id_list)\n\n\t# divide the id into test and train id lists\n\ttest_len = int(len(id_list)*test_fraction)\n\ttest_id = id_list[:test_len]\n\ttrain_id = id_list[test_len:]\n\ttrain_len = len(train_id)\n\n\t\n\tif set_format == 'group':\n\t\tX_train, Y_train, X_test, Y_test = format_group(data, train_id, test_id, features, min_size = min_group, sequence_length = sequence_length, remove_check_delta = remove_check_delta, combine_check_delta = combine_check_delta, group_size = group_size, preprocess = preprocess, set_GPR = set_GPR)\n\telif set_format == 'season':\n\t\tX_train, Y_train, X_test, Y_test = format_season(data,train_id, test_id, features, min_size = min_group, sequence_length = sequence_length, remove_check_delta = remove_check_delta, combine_check_delta = combine_check_delta, group_size = group_size, preprocess = preprocess, set_GPR = set_GPR)\n\telif set_format == 'simple':\n\t\tX_train, Y_train, X_test, Y_test = format_simple(data, train_id, test_id, features, preprocess = preprocess)\n\t\tpadding = True\n\telse:\n\t\tprint('ERROR: wrong format!')\n\n\n\n\tif padding == True:\n\t\t# obtain the max number of candidates among all objects, and pad the train and test sets to have the same number of candidates for each object\n\t\tif sequence_length == None:\n            maxlen = np.max([len(x) for x in X_train and X_test])\n        else:\n            maxlen = sequence_length\n\t\tX_train = ks.preprocessing.sequence.pad_sequences(X_train,maxlen=maxlen,padding='pre',value=0,dtype='Float32')\n\t\tX_test = ks.preprocessing.sequence.pad_sequences(X_test,maxlen=maxlen,padding='pre',value=0,dtype='Float32')\n\n\n\n\telse:\n\t\t# padding csv: convert to array\n\t\tX_train = np.array(X_train)\n\t\tX_test = np.array(X_test)\n\n\tnum_classes = np.unique(Y_train).shape[0]\n\n\tY_train = ks.utils.to_categorical(Y_train, num_classes = num_classes, dtype = 'Float32')\n\tY_test = ks.utils.to_categorical(Y_test, num_classes = num_classes, dtype = 'Float32')\n\n\tlength_train = X_train.shape[0]\n\tlength_test = X_test.shape[0]\n\n\ttime_sequence = X_train.shape[1]\n\tinput_dim = X_train.shape[2]\n\n\n\n\n\treturn (X_train, X_test, Y_train, Y_test),(length_train, length_test, time_sequence, input_dim, num_classes)\n\n\ndef cut_test_data(X_test, Y_test, cut_fraction, set_format, features, sequence_length == None):\n\t'''\n\tThis function is used for testing the accuracy/AUC with the increasing number of observations in one group\n\t'''\n\tnum_classes = np.unique(Y_test).shape[0]\n\tif cut_fraction == None or cut_fraction == 0 or set_format=='simple':\n\t\tif set_format != 'simple':\n\t\t\tX_test = np.array(X_test)\n\t\telse:\n\t\t\tif sequence_length==None:\n\t\t\t    maxlen = np.max([len(x) for x in X_test])\n\t\t\telse:\n\t\t\t\tmaxlen = sequence_length\n\t\t\tX_test = ks.preprocessing.sequence.pad_sequences(X_test,maxlen=maxlen,padding='pre',value=0)\n\t\tY_test = ks.utils.to_categorical(Y_test, num_classes = num_classes, dtype = 'float32')\n\n\telif float(cut_fraction)<1.0 and float(cut_fraction)>0:\n\t\tcut_fraction = float(cut_fraction)\n\t\tnew_test = []\n\t\tif set_format == 'season':\n\t\t\tseason_size = len(X_test[0])\n\t\t\tfeature_num = len(features)\n\t\t\tsave_num = int(season_size*(1-cut_fraction))\n\t\t\tfor season in X_test:\n\t\t\t\tsave = season[:save_num]\n\t\t\t\tremove = season[save_num:]\n\t\t\t\tzero_padding = [[0 for x in range(feature_num)] for y in range(len(remove))]\n\t\t\t\tnew_test.append(save+zero_padding)\n\n\t\telif set_format == 'group':\n\t\t\tgroup_size = int(len(X_test[0][0])/len(features))\n\t\t\tsave_num = int(group_size*(1-cut_fraction))\n\t\t\tremove_num = group_size - save_num\n\t\t\tfor obj in X_test:\n\t\t\t\tobj_sequence = []\n\t\t\t\tfor group in obj:\n\t\t\t\t\tn = 0\n\t\t\t\t\tcheck = 0\n\t\t\t\t\tgroup_obs = []\n\t\t\t\t\twhile n<len(features):\n\t\t\t\t\t\tgroup_obs +=group[check:check+save_num]+[0 for x in range(remove_num)]\n\t\t\t\t\t\tcheck += group_size\n\t\t\t\t\t\tn += 1\n\t\t\t\t\tobj_sequence.append(group_obs)\n\t\t\t\tnew_test.append(obj_sequence)\n\t\telse:\n\t\t\tprint('\\nCut function doesn\\'t support Simple input format.\\n')\n\n\t\tX_test = new_test\n\n\t\tX_test = np.array(X_test)\n\n\t\tY_test = ks.utils.to_categorical(Y_test, num_classes = num_classes, dtype = 'float32')\n\n\treturn X_test, Y_test\n\n\ndef load_test_data(path, seed = None, features = list, set_format = 'group', preprocess = 's', min_group =3, group_size = None, sequence_length = None, remove_check_delta = 233, combine_check_delta = 0.7, set_GPR = False, padding = False):\n\t'''\n\tThis function is used for the construction of test data. It will be used in the model prediction to test and compare the performance of different classifiers.\n\t'''\n\tprint('Start constructing the test data...')\n\tdata = pd.read_csv(path) \n\t\n\t\n\t# select columns\n\tif 'type' in list(data.columns):\n\t\tdata = data.loc[:,['id','mjd','type']+features]\n\telse:\n\t\tdata = data.loc[:,['id','mjd','label']+features]\n\n\t# print('test')\n\t#remove NaN band values, replace NaN with 0\n\t# data = data.dropna(subset = features, axis = 'rows', how='all').reset_index()\n\t# data = data.fillna(0)\n# \tprint(data[:10])\n\n\t\n\tcombine = data.groupby(data['id']) \n\tid_list = []\n\tfor _id, group in combine:\n\t\tid_list.append(_id)\n\n\tif set_format == 'group':\n\t\tremove1, remove2, X_test, Y_test = format_group(data, [], id_list, features, min_size = min_group, sequence_length = sequence_length, remove_check_delta = remove_check_delta, combine_check_delta = combine_check_delta, group_size = group_size, preprocess = preprocess, set_GPR = set_GPR)\n\telif set_format == 'season':\n\t\tremove1, remove2, X_test, Y_test = format_season(data, [], id_list, features, min_size = min_group, sequence_length = sequence_length, remove_check_delta = remove_check_delta, combine_check_delta = combine_check_delta, group_size = group_size, preprocess = preprocess, set_GPR = set_GPR)\n\telif set_format == 'simple':\n\t\tremove1, remove2, X_test, Y_test = format_simple(data,  [], id_list, features, preprocess = preprocess)\n\telse:\n\t\tprint('ERROR: wrong format!')\n\n\treturn X_test, Y_test\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# model #","metadata":{}},{"cell_type":"code","source":"from numpy.random import seed\n# seed(1)\nimport tensorflow as tf\nfrom tensorflow.keras import backend as K\n# tensorflow.random.set_seed(1)\nfrom tensorflow.keras.models import Sequential, Model\nfrom tensorflow.keras.layers import Dense, Dropout, Input\nfrom tensorflow.keras.layers import Embedding, Activation, Masking\nfrom tensorflow.keras.layers import LSTM, Bidirectional, GRU, SimpleRNN, TimeDistributed\nfrom tensorflow.keras.callbacks import EarlyStopping\nimport numpy as np\nimport time\nimport csv\nimport os\n\nprint(tf.config.list_physical_devices('GPU'))\n\n'''\nThis file is used for building and training the classifier and saving the results and models.\n'''\n\nclass LossHistory(tensorflow.keras.callbacks.Callback):\n    '''\n    This class is used for recording the loss, accuracy, AUC, f1_score value during training.\n    '''\n    def on_train_begin(self, logs={}):\n        self.epoch_loss = []\n        self.epoch_accuracy = []\n        self.epoch_AUC = []\n        self.epoch_f1_score = []\n\n        self.epoch_val_loss = []\n        self.epoch_val_accuracy = []\n        self.epoch_val_AUC = []\n        self.epoch_val_f1_score = []\n\n        self.batch_losses = []\n        self.batch_accuracy = []\n        self.batch_AUC = []\n        self.batch_f1_score = []\n      \n\n    def on_epoch_end(self, batch, logs={}):\n        self.epoch_loss.append(logs.get('loss'))\n        self.epoch_accuracy.append(logs.get('accuracy'))\n        self.epoch_AUC.append(logs.get('auc'))\n        self.epoch_f1_score.append(logs.get('f1_score'))\n\n        self.epoch_val_loss.append(logs.get('val_loss'))\n        self.epoch_val_accuracy.append(logs.get('val_accuracy'))\n        self.epoch_val_AUC.append(logs.get('val_auc'))\n        self.epoch_val_f1_score.append(logs.get('val_f1_score'))\n\n    def on_batch_end(self, batch,logs={}):\n        self.batch_losses.append(logs.get('loss'))\n        self.batch_accuracy.append(logs.get('accuracy'))\n        self.batch_AUC.append(logs.get('auc'))\n        self.batch_f1_score.append(logs.get('f1_score'))\n\n\ndef parser_config(config):\n    '''\n    parse the config file\n    '''\n#     config = {}\n#     f = open(path, 'r').readlines()\n#     n = 1\n#     while n<len(f):\n#         if ':' in f[n]:\n#             x = f[n][:-1].split(':')\n#             x[1] = x[1].strip()\n#             config[x[0]] = x[1]\n#             if x[0] == 'path':\n#                 if x[1] == '':\n#                     config[x[0]] = '../data/processed/balanced/final_v1.csv'\n#                 else:\n#                     config[x[0]] = x[1]\n#             elif x[0] == 'features':\n#                 if x[1]!= '':\n#                     config[x[0]] = x[1].split(',')\n#                 else:\n#                     print(\"WARNING: empty feature list, automatically set the feature as g\")\n#                     config[x[0]] = ['g']\n#             elif x[0] == 'format':\n#                 if x[1]!= '':\n#                     config[x[0]] = x[1]\n#                 else:\n#                     config[x[0]] = 'group'\n#             elif x[0] == 'processed':\n#                 if x[1]!= '':\n#                     if x[1][0] == 's':\n#                         config[x[0]] = 's'\n#                     elif x[1][0] == 'n':\n#                         config[x[0]] = 'n'\n#                     elif  x[1][0] == 'd':\n#                         config[x[0]] = 'd'\n#                     else:\n#                         print(\"WARNING: wrong preprocess input, automatically set as standardization\")\n#                         config[x[0]] = 's'\n#                 else:\n#                     print(\"WARNING: empty preprocess input, automatically set as standardization\")\n#                     config[x[0]] = 's'\n#             elif x[0] == 'hidden_layers':\n#                 config[x[0]] = [int(l) for l in x[1].strip('[').strip(']').split(',')]\n#             elif x[0] == 'metrics':\n#                 config[x[0]] = [l for l in x[1].split(',') if l != 'f1_score']\n#                 if 'f1_score' in x[1].split(','):\n#                     config[x[0]].append(f1_score)\n#             else:\n#                 if x[1].lower() == 'true':\n#                     config[x[0]] = True\n#                 if x[1].lower() == 'false':\n#                     config[x[0]] = False\n#                 if x[1] == ''or x[1] == 'None':\n#                     config[x[0]] = None\n#             if x[0] == 'seed' or x[0] == 'group_size' or x[0] == 'group_num':\n#                 if config[x[0]] != None:\n#                     config[x[0]] = int(config[x[0]])\n               \n#         n +=1\n#     print('configuration: ',config)\n    \n    # generate a file tree\n    path = config['save_path']\n    if os.path.exists(path) == False:\n        os.mkdir(os.getcwd()+'/'+path)\n    if os.path.exists(path+'/'+config['rnn_type'].lower()) == False: \n        os.mkdir(os.getcwd()+'/'+path+'/'+config['rnn_type'].lower())\n    path = path+'/'+config['rnn_type'].lower()\n    if os.path.exists(path+'/'+config['format']) == False:\n        os.mkdir(os.getcwd()+'/'+path+'/'+config['format'])\n    path = path+'/'+config['format']\n    if os.path.exists(path+'/'+str(config['features'])) == False:\n        os.mkdir(os.getcwd()+'/'+path+'/'+str(config['features']))\n    path = path + '/'+str(config['features'])\n    if os.path.exists(path+'/GP') == False and config['set_GPR'] == True:\n        os.mkdir(os.getcwd()+'/'+path+'/GP')\n    elif os.path.exists(path+'/non_GP') == False and config['set_GPR'] == False:\n        os.mkdir(os.getcwd()+'/'+path+'/non_GP')\n    if config['set_GPR'] == True:\n        path = path + '/GP'\n    else:\n        path = path + '/non_GP'\n    print(config['hidden_layers'])\n    if os.path.exists(path+'/'+str(config['hidden_layers'])) == False:\n        os.mkdir(os.getcwd()+'/'+path+'/'+str(config['hidden_layers']))\n    path = path + '/'+str(config['hidden_layers'])\n    if os.path.exists(path+'/'+config['processed']) == False:\n        os.mkdir(os.getcwd()+'/'+path+'/'+config['processed'])\n    path = path + '/'+config['processed']\n\n    print('save path: ',path)\n\n    return path\n\ndef binary_focal_loss(gamma=2, alpha=0.25):\n    '''\n    Binary form of focal loss.\n    \n    focal_loss(p_t) = -alpha_t * (1 - p_t)**gamma * log(p_t)\n        where p = sigmoid(x), p_t = p or 1 - p depending on if the label is 1 or 0, respectively.\n    References:\n        https://arxiv.org/pdf/1708.02002.pdf\n    Usage:\n     model.compile(loss=[binary_focal_loss(alpha=.25, gamma=2)], metrics=[\"accuracy\"], optimizer=adam)\n    '''\n    alpha = tf.constant(alpha, dtype=tf.float32)\n    gamma = tf.constant(gamma, dtype=tf.float32)\n\n    def binary_focal_loss_fixed(y_true, y_pred):\n        y_true = tf.cast(y_true, tf.float32)\n        alpha_t = y_true*alpha + (K.ones_like(y_true)-y_true)*(1-alpha)\n    \n        p_t = y_true*y_pred + (K.ones_like(y_true)-y_true)*(K.ones_like(y_true)-y_pred) + K.epsilon()\n        focal_loss = - alpha_t * K.pow((K.ones_like(y_true)-p_t),gamma) * K.log(p_t)\n        return K.mean(focal_loss)\n    return binary_focal_loss_fixed\n\n\n# def f1_score(y_true, y_pred):\n#     '''\n#     calculate F1 value\n#     '''\n\n#     TP = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))\n#     label_true = K.sum(K.round(K.clip(y_pred, 0, 1)))\n#     pred_true = K.sum(K.round(K.clip(y_true, 0, 1)))\n\n#     precision = TP/pred_true\n#     recall = TP/label_true\n#     f1_score = 2 * (precision * recall)/(precision + recall)\n\n#     return f1_score\n\n\n# def precision(y_true,y_pred): \n#     TP=tf.reduce_sum(y_true*tf.round(y_pred))\n#     TN=tf.reduce_sum((1-y_true)*(1-tf.round(y_pred)))\n#     FP=tf.reduce_sum((1-y_true)*tf.round(y_pred))\n#     FN=tf.reduce_sum(y_true*(1-tf.round(y_pred)))\n#     precision=TP/(TP+FP)\n#     return precision\n \n# def recall(y_true,y_pred): \n#     TP=tf.reduce_sum(y_true*tf.round(y_pred))\n#     TN=tf.reduce_sum((1-y_true)*(1-tf.round(y_pred)))\n#     FP=tf.reduce_sum((1-y_true)*tf.round(y_pred))\n#     FN=tf.reduce_sum(y_true*(1-tf.round(y_pred)))\n#     recall=TP/(TP+FN)\n#     return recall\n\n\ndef build(time_sequence = int, input_dim = int, num_class = 2, hidden_layer = [256,256,256], rnn_type = 'LSTM', dropout = 0.05,  activation = 'tanh', print_model = False, predict_trend = False):\n    '''\n    This is the function for building a classifer model.\n    Inputs:\n    - time_sequence: the length of the sequence\n    - input_dim: the number of dimensions for each vector\n    - num_class: the number of classes/labels\n    - hidden_layer: the number of layers in RNN architecture\n    - rnn_type: the type of RNN, for example, LSTM, GRU, SimpleRNN\n    - dropout: the fraction of objects dropped before being fed into the next layer to avoid overfitting\n    - activation: applies an activation function to an output\n    - print_model: True means the RNN architecture will be printed\n    - predict_trend: whether use the prediction function. ***This function is being developed***\n    Returns:\n    - model: the designed model with fixed layers\n    '''\n\n    # tensorflow.keras.initializers.Constant(value=0)\n\n    if rnn_type.lower() == 'lstm':\n        RNN = LSTM\n    elif rnn_type.lower() == 'gru':\n        RNN = GRU\n    else:\n        RNN = SimpleRNN\n\n    input_data = Input(shape=(time_sequence, input_dim))\n\n    x = Masking(mask_value=0.0000)(input_data)\n\n    n = 0\n    while n<len(hidden_layer):\n        x = RNN(units=hidden_layer[n], activation = activation, return_sequences=(n<len(hidden_layer)-1))(x)\n        x = Dropout(dropout)(x)\n        n +=1\n    y = Dense(units = num_class, activation = 'softmax')(x)\n\n\n    model = Model([input_data],[y])\n\n    if print_model == True:\n        print(model.summary())\n\n\n    return model\n\n\ndef train(model, path, config, batch_size = int, epochs = int, X_train = list, Y_train = list, X_test = list, Y_test = list, features = list,\n    optimizer = 'Adam', lr = None, decay = False, loss_function = 'binary_crossentropy', metrics = ['accuracy','AUC']):\n    '''\n    train the classifier model.\n    Inputs:\n    - model: the designed model generated from the build function\n    - batch_size: the number of sequences fed into the layer for each time\n    - epoches: the times for the input data being processed\n    - X_train/X_test: the train/test data with 3 dimensions: \n        objects, the number of observations in an objects, the vector in 1 observation\n    - Y_train/Y_test: the label of train/test data. \n    = features: the targeted feature list\n    - optimizer: the optimization method for the loss function\n    - lr: learning rate\n    - decay: whether the learning rate will decrease with the increasing nunmber of epochs\n    - loss_function: the function used to measure the degree to which the predicted value f(x) of the model is inconsistent with the true value Y\n    - metrics: the metrics during training for testing the performance of the classifier\n    Return:\n    - model: the trained model\n\n    '''\n\n    d_value = 0.0\n    if lr != None:\n        lr = float(lr)\n    else:\n        lr = 0.001\n    if decay == True:\n        d_value = lr/epochs\n    else:\n        d_value = 0\n\n\n    opt = 'Adam'\n    if optimizer.lower == 'adam':\n        opt = tf.keras.optimizers.Adam(lr=lr, beta_1=0.9, beta_2=0.999, epsilon=1e-08, decay=d_value)\n    elif optimizer.lower == 'sgd':\n        opt = tf.keras.optimizers.SGD(lr=lr, momentum=0.0, nesterov=False, name='SGD', decay=d_value)\n\n    model.compile(loss = loss_function, optimizer = opt, metrics = metrics)\n\n    history = LossHistory()\n\n    print(\"Start training...\")\n\n    start_time = time.time()\n\n    class_weight = {0: sum(x[1] for x in Y_train)/len(Y_train), 1: sum(x[0] for x in Y_train)/len(Y_train)}\n\n    earlystop = EarlyStopping(monitor = 'val_loss', patience = 3)\n    \n    model.fit(X_train, Y_train, batch_size=batch_size, epochs=epochs, validation_data=(X_test, Y_test), verbose=1, class_weight = class_weight, callbacks=[history,earlystop])\n\n    sum_time = time.time() - start_time\n\n    print(\"--- %s seconds ---\" % sum_time)\n\n    #plot graphs: loss, accuracy, AUC\n\n    plot_two_graph(history.epoch_val_loss, history.epoch_loss, path+'/loss_epoch.png', 'num of epochs', 'loss', 'validation_loss','training_loss')\n    plot_one_graph(history.batch_losses, path+'/loss_batch.png', 'num of batches', 'loss')\n\n    plot_two_graph(history.epoch_val_accuracy, history.epoch_accuracy, path+'/acc_epoch.png', 'num of epochs','accuracy', 'validation_acc', 'training_acc')\n    plot_one_graph(history.batch_accuracy, path+'/acc_batch.png', 'num of batches', 'accuracy')\n\n    plot_two_graph(history.epoch_val_AUC, history.epoch_AUC, path+'/AUC_epoch.png','num of epochs', 'AUC','validation_AUC','training_AUC')\n    plot_one_graph(history.batch_AUC, path+'/AUC_batch.png','num of batches', 'AUC')\n\n    scores = model.evaluate(X_test, Y_test, verbose=1)\n\n    # record all loss, acc, AUC information\n    with open(path+'/record.txt','w') as f:\n        f.write('loss_epoch_train: ')\n        f.write(str(history.epoch_loss))\n        f.write('\\nloss_epoch_valid: ')\n        f.write(str(history.epoch_val_loss))\n        f.write('\\nloss_batch_train: ')\n        f.write(str(history.batch_losses))\n        f.write('\\nacc_epoch_train: ')\n        f.write(str(history.epoch_accuracy))\n        f.write('\\nacc_epoch_valid: ')\n        f.write(str(history.epoch_val_accuracy))\n        f.write('\\nacc_batch_train: ')\n        f.write(str(history.batch_accuracy))\n        f.write('\\nauc_epoch_train: ')\n        f.write(str(history.epoch_AUC))\n        f.write('\\nauc_epoch_valid: ')\n        f.write(str(history.epoch_val_AUC))\n        f.write('\\nauc_batch_train: ')\n        f.write(str(history.batch_AUC))\n    f.close()\n\n    y_pred = model.predict(X_test)\n\n    actual_value = [x[1] for x in Y_test]\n    predict_value = [x[1] for x in y_pred]\n    auc = plot_ROC(actual_value, predict_value, path+'/val_')\n\n    y_pred = [np.argmax(x) for x in y_pred]\n    y_true = [np.argmax(x) for x in Y_test]\n\n    cm = confusion_matrix(y_true, y_pred)\n    cm_normalized = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n    plot_confusion_matrix(cm_normalized, path+'/val_confusion_matrix', title='confusion matrix',classes=['non-QSO', 'QSO'])\n\n    tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()\n    recall = tp/(tp+fn)\n    specificity = tn/(fp+tn)\n    precision = tp/(tp+fp)\n    f1_score = 2*precision*recall/(precision+recall)\n\n    key_list = ['features','format','processed','set_GPR','rnn_type','hidden_layers','dropout','batch_size','num_epochs','test_fraction','optimizer','learning_rate','decay']\n\n    with open(config['save_path']+'/train_results.csv','a+',newline='') as csvfile:\n        writer = csv.writer(csvfile)\n        with open(config['save_path']+'/train_results.csv','r',newline='') as f:\n            reader = csv.reader(f)\n            if not [row for row in reader]:\n                writer.writerow(key_list+['train_time']+model.metrics_names + ['val_recall','val_specificity','val_precision','val_f1_score']+ ['TN','FP','FN','TP'])\n            writer.writerow([config[x] for x in key_list] + [str(sum_time)] + [scores[i] for i in range(len(model.metrics_names))] +[recall, specificity, precision, f1_score]+ [tn, fp, fn, tp])\n\n\n\n\n    return model\n\ndef predict(model, X_test, Y_test, path, config):\n    '''\n    This function is used to test the classifier's performance with a confusion metrix graph and ROC graph.\n    '''\n\n    y_pred = model.predict(X_test)\n    actual_value = [x[1] for x in Y_test]\n    predict_value = [x[1] for x in y_pred]\n    auc = plot_ROC(actual_value, predict_value, path+'/test_')\n    y_pred = [np.argmax(x) for x in y_pred]\n    y_true = [np.argmax(x) for x in Y_test]\n    cm = confusion_matrix(y_true, y_pred)\n    cm_normalized = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n    plot_confusion_matrix(cm_normalized, path+'/test_confusion_matrix', title='confusion matrix',classes=['non-QSO', 'QSO'])\n\n    title = ','.join(config['features'])\n    key_list = ['features','format','processed','set_GPR','rnn_type','hidden_layers','dropout','batch_size','num_epochs','test_fraction','optimizer','learning_rate','decay']\n\n    tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()\n    recall = tp/(tp+fn)\n    specificity = tn/(fp+tn)\n    precision = tp/(tp+fp)\n    f1_score = 2*precision*recall/(precision+recall)\n\n    key_list = ['features','format','processed','set_GPR','rnn_type','hidden_layers','dropout','batch_size','num_epochs','test_fraction','optimizer','learning_rate','decay']\n\n    with open(config['save_path']+'/test_results.csv','a+',newline='') as csvfile:\n        writer = csv.writer(csvfile)\n        with open(config['save_path']+'/test_results.csv','r',newline='') as f:\n            reader = csv.reader(f)\n            if not [row for row in reader]:\n                writer.writerow(key_list+['AUC','recall','specificity','precision','f1_score']+ ['TN','FP','FN','TP'])\n            writer.writerow([config[x] for x in key_list]+[auc, recall, specificity, precision, f1_score]+ [tn, fp, fn, tp])\n\n\n   \n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### plots","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import confusion_matrix\nfrom sklearn.metrics import roc_curve, auc \nimport matplotlib.pyplot as plt\nfrom matplotlib.pyplot import MultipleLocator\nimport numpy as np\nimport sys\n\n\ndef plot_confusion_matrix(cm, savename, title='Confusion Matrix', classes = int):\n\n    plt.figure(figsize=(12, 8), dpi=100)\n    np.set_printoptions(precision=2)\n\n    ind_array = np.arange(len(classes))\n    x, y = np.meshgrid(ind_array, ind_array)\n    for x_val, y_val in zip(x.flatten(), y.flatten()):\n        c = cm[y_val][x_val]\n        if c > 0.001:\n            plt.text(x_val, y_val, \"%0.2f\" % (c,), color='red', fontsize=15, va='center', ha='center')\n    \n    plt.imshow(cm, interpolation='nearest', cmap='Blues')\n    plt.title(title,fontsize=20)\n    plt.colorbar()\n    xlocations = np.array(range(len(classes)))\n    plt.xticks(xlocations, classes, fontsize=18)\n    plt.yticks(xlocations, classes, fontsize=18)\n    plt.ylabel('Actual label', fontsize=20)\n    plt.xlabel('Predict label',fontsize=20)\n    \n    # offset the tick\n    tick_marks = np.array(range(len(classes))) + 0.5\n    plt.gca().set_xticks(tick_marks, minor=True)\n    plt.gca().set_yticks(tick_marks, minor=True)\n    plt.gca().xaxis.set_ticks_position('none')\n    plt.gca().yaxis.set_ticks_position('none')\n    plt.grid(True, which='minor', linestyle='-')\n    plt.gcf().subplots_adjust(bottom=0.15)\n    \n    # show confusion matrix\n    plt.savefig(savename+'.png', format='png')\n\n\n\ndef plot_ROC(actual, predict, path):\n    false_positive_rate, true_positive_rate, thresholds = roc_curve(actual, predict)\n    roc_auc = auc(false_positive_rate, true_positive_rate)\n    plt.figure(figsize=(10, 10), dpi=100)\n    plt.title('Receiver Operating Characteristic',fontsize=16)\n    plt.plot(false_positive_rate, true_positive_rate, 'b',label='ROC curve (area = %0.2f)' % roc_auc)\n    plt.legend(loc='lower right')\n    plt.plot([0,1],[0,1],'r--')\n    plt.xlim([-0.05,1.05])\n    plt.ylim([-0.05,1.05])\n    plt.xticks(fontsize=18)\n    plt.yticks(fontsize=18)\n    plt.ylabel('True Positive Rate',fontsize=20)\n    plt.xlabel('False Positive Rate',fontsize=20)\n    plt.savefig(path+'ROC.png', format = 'png')\n\n    with open(path+'roc_record.txt','w+') as f:\n        f.write('false_positive_rate: ')\n        f.write(str(list(false_positive_rate)))\n        f.write('\\ntrue_positive_rate: ')\n        f.write(str(list(true_positive_rate)))\n        f.write('\\nthresholds: ')\n        f.write(str(list(thresholds)))\n        f.write('\\nroc_auc: ')\n        f.write(str(roc_auc))\n    f.close()\n    return roc_auc\n\n\ndef plot_two_graph(value1, value2, savename, x_label, y_label, y1, y2):\n    '''\n    plot for epochs\n    '''\n    plt.figure(figsize=(12, 8), dpi=100)\n    plt.clf()\n    plt.plot([n for n in list(range(1,len(value1)+1))], value1, marker='.', label = y1)\n    plt.plot([n for n in list(range(1,len(value2)+1))], value2, marker='.', label = y2)\n    plt.legend(loc='upper left')\n    plt.xticks(fontsize=18)\n    plt.yticks(fontsize=18)\n    plt.xlabel(x_label,fontsize=20)\n    plt.ylabel(y_label,fontsize=20)\n    if 'epoch' in x_label and len(value1)<=50:\n        ax = plt.gca()\n        x_major_locator=MultipleLocator(1)\n        ax.xaxis.set_major_locator(x_major_locator)\n    plt.xlim(1, len(value1))\n    if 'loss' not in savename:\n        plt.ylim(0,1)\n    plt.savefig(savename, format = 'png')\n\ndef plot_one_graph(value, savename, x_label, y_label):\n    '''\n    plot for batches\n    '''\n    plt.figure(figsize=(12, 8), dpi=100)\n    plt.clf()\n    plt.plot([n for n in list(range(1,len(value)+1))], value, label = y_label)\n    plt.legend(loc='upper left')\n    plt.xlabel(x_label,fontsize=20)\n    plt.ylabel(y_label,fontsize=20)\n    plt.xticks(fontsize=18)\n    plt.yticks(fontsize=18)\n    plt.xlim(1, len(value))\n    if 'loss' not in savename:\n        plt.ylim(0,1)\n    plt.savefig(savename, format = 'png')   \n\n    ","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Start running codes\n### Configuration setting","metadata":{}},{"cell_type":"code","source":"config = {}\n\n#input config\nconfig['train_path'] = '../input/sdss-stripe82-quasar-targeted-dataset/train_set.csv'\nconfig['test_path']  = '../input/sdss-stripe82-quasar-targeted-dataset/test_set.csv'\nconfig['save_path'] = '20210510sdss'\nconfig['seed'] = 34\nconfig['features'] = ['g','r']\nconfig['format'] = 'simple'\nconfig['processed'] = 's'\nconfig['set_GPR'] = False\nconfig['group_size'] = 67\nconfig['group_num'] = None\nconfig['cut_fraction'] = None\n\n#network config\nconfig['rnn_type'] = 'LSTM'\nconfig['hidden_layers'] = [64,64,64,64]\nconfig['dropout'] = 0.35\nconfig['plot_model'] = True\n\n#training config\nconfig['batch_size'] = 512\nconfig['num_epochs'] = 50\nconfig['test_fraction'] = 0.2\nconfig['optimizer'] = 'adam'\nconfig['learning_rate'] = 0.0001\nconfig['decay'] = ''\nconfig['metrics'] = ['accuracy','AUC']\n","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training","metadata":{}},{"cell_type":"code","source":"path = parser_config(config)\n\n(X_train, X_valid, Y_train, Y_valid), (length_train, length_test, time_sequence, input_dim, num_classes) = load_data(path = config['train_path'], test_fraction = float(config['test_fraction']), seed = config['seed'], features = config['features'], set_format = config['format'], preprocess = config['processed'], group_size = config['group_size'], sequence_length = config['group_num'],set_GPR = config['set_GPR'])\n\nmodel = build(time_sequence = time_sequence, input_dim = input_dim, num_class = num_classes, hidden_layer = config['hidden_layers'], rnn_type = config['rnn_type'], dropout = float(config['dropout']),  print_model = config['plot_model'])\n\nmodel = train(model, path, config, batch_size = int(config['batch_size']), epochs = int(config['num_epochs']), X_train = X_train, Y_train = Y_train, X_test = X_valid, Y_test = Y_valid, features = config['features'], optimizer = config['optimizer'],  lr = config['learning_rate'], decay = config['decay'], loss_function = 'binary_crossentropy', metrics = config['metrics'])\n\nX_test, Y_test = load_test_data(config['test_path'], seed = config['seed'], features = config['features'], set_format = config['format'], preprocess = config['processed'], group_size = config['group_size'], sequence_length = config['group_num'], set_GPR = config['set_GPR']) \n\nX_test, Y_test= cut_test_data(X_test, Y_test, cut_fraction=config['cut_fraction'], set_format=config['format'], features=config['features'])\n\npredict(model, X_test, Y_test, path, config)\n\ntitle = ','.join(config['features'])    \n# model.save(path+'/model_'+title+'_ep_'+str(config['num_epochs'])+'_bs_'+str(config['batch_size'])+'.h5')\n# print(\"Save model to the disk.\") ","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}