{"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":"code","source":"from IPython.display import clear_output\n!pip install paddle-quantum\nclear_output()","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:38:52.278866Z","iopub.execute_input":"2022-03-31T01:38:52.279809Z","iopub.status.idle":"2022-03-31T01:39:54.902920Z","shell.execute_reply.started":"2022-03-31T01:38:52.279687Z","shell.execute_reply":"2022-03-31T01:39:54.901957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import time\nimport matplotlib\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport paddle\nimport astropy\nfrom numpy import pi as PI\nfrom matplotlib import pyplot as plt\n\nfrom paddle import matmul, transpose\nfrom paddle_quantum.circuit import UAnsatz\n\nimport sklearn\nfrom sklearn import svm\nfrom sklearn.datasets import fetch_openml, make_moons, make_circles\nfrom sklearn.model_selection import train_test_split\n\nfrom IPython.display import clear_output\nfrom tqdm import tqdm\n\nimport time\nimport matplotlib\nfrom numpy import pi as PI\n\nimport multiprocessing\nimport warnings\nfrom itertools import chain\nsns.set_style('whitegrid')\nwarnings.simplefilter('ignore', FutureWarning)\nwarnings.simplefilter('ignore', RuntimeWarning)\nfrom cesium.time_series import TimeSeries\nimport cesium.featurize as featurize\nfrom gatspy.periodic import LombScargleMultiband, LombScargleMultibandFast\nimport pdb","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:41:06.889237Z","iopub.execute_input":"2022-03-31T01:41:06.889573Z","iopub.status.idle":"2022-03-31T01:41:10.784671Z","shell.execute_reply.started":"2022-03-31T01:41:06.889531Z","shell.execute_reply":"2022-03-31T01:41:10.783800Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data visualization","metadata":{}},{"cell_type":"code","source":"from astropy.table import Table\n\nfilename = '../input/training_set.csv'\ndata = Table.read(filename, format='csv')\nnobjects = len(data)\ndata","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:10:18.754907Z","iopub.execute_input":"2022-03-31T02:10:18.755195Z","iopub.status.idle":"2022-03-31T02:10:20.567132Z","shell.execute_reply.started":"2022-03-31T02:10:18.755166Z","shell.execute_reply":"2022-03-31T02:10:20.566129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_series = pd.read_csv('../input/training_set.csv')\ntrain_metadata = pd.read_csv('../input/training_set_metadata.csv')\n","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:41:26.434414Z","iopub.execute_input":"2022-03-31T01:41:26.434753Z","iopub.status.idle":"2022-03-31T01:41:27.388110Z","shell.execute_reply.started":"2022-03-31T01:41:26.434721Z","shell.execute_reply":"2022-03-31T01:41:27.387101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def view_target(n):\n    obj_id = np.random.choice(train_metadata.object_id[train_metadata.target == n].values)\n    obj_df = train_series[train_series.object_id == obj_id]\n    fig, axes = plt.subplots(6,1, figsize=(10, 10))\n    axes[0].set_title(f'Class_{n}')\n    for i, ax in enumerate(axes):\n        ax.scatter(obj_df.mjd[obj_df.passband == i].values, obj_df.flux[obj_df.passband == i].values, alpha=0.5)    \n        ax.scatter(obj_df.mjd[obj_df.passband == i].values, obj_df.flux_err[obj_df.passband == i].values, alpha=0.5)\n    return obj_df","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:25.999022Z","iopub.execute_input":"2022-03-31T01:46:26.001052Z","iopub.status.idle":"2022-03-31T01:46:26.012747Z","shell.execute_reply.started":"2022-03-31T01:46:26.000946Z","shell.execute_reply":"2022-03-31T01:46:26.011628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"view_target(42)\nview_target(90)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:27.651320Z","iopub.execute_input":"2022-03-31T01:46:27.652316Z","iopub.status.idle":"2022-03-31T01:46:30.438479Z","shell.execute_reply.started":"2022-03-31T01:46:27.652262Z","shell.execute_reply":"2022-03-31T01:46:30.437787Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"groups = train_series.groupby(['object_id', 'passband'])\ntimes = groups.apply(\n    lambda block: block['mjd'].values).reset_index().rename(columns={0: 'seq'})\nflux = groups.apply(\n    lambda block: block['flux'].values\n).reset_index().rename(columns={0: 'seq'})\nerr = groups.apply(\n    lambda block: block['flux_err'].values\n).reset_index().rename(columns={0: 'seq'})\ndet = groups.apply(\n    lambda block: block['detected'].astype(bool).values\n).reset_index().rename(columns={0: 'seq'})\ntimes_list = times.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()\nflux_list = flux.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()\nerr_list = err.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()\ndet_list = det.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:30.439767Z","iopub.execute_input":"2022-03-31T01:46:30.440545Z","iopub.status.idle":"2022-03-31T01:46:45.286379Z","shell.execute_reply.started":"2022-03-31T01:46:30.440508Z","shell.execute_reply":"2022-03-31T01:46:45.285417Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We cannot make much sense of the data with the large observation gap, unsynchronised passband observations and sparsity of data. But We know that some objects have periodic bebaviours, so we can attempt to fold them by period. Here we will examine which classes are more likely to be periodic, and how they typically look like.","metadata":{}},{"cell_type":"code","source":"def fit_multiband_freq(tup):\n    idx, group = tup\n    t, f, e, b = group['mjd'], group['flux'], group['flux_err'], group['passband']\n    model = LombScargleMultiband(fit_period=True)\n    model.optimizer.period_range = (0.1, int((group['mjd'].max() - group['mjd'].min()) / 2))\n    model.fit(t, f, e, b)\n    return model","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:45.288437Z","iopub.execute_input":"2022-03-31T01:46:45.288739Z","iopub.status.idle":"2022-03-31T01:46:45.295251Z","shell.execute_reply.started":"2022-03-31T01:46:45.288708Z","shell.execute_reply":"2022-03-31T01:46:45.294128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_freq_features(N, subsetting_pos=None):\n    if subsetting_pos is None:\n        subset_times_list = times_list\n        subset_flux_list = flux_list\n    else:\n        subset_times_list = [v for i, v in enumerate(times_list) \n                             if i in set(subsetting_pos)]\n        subset_flux_list = [v for i, v in enumerate(flux_list) \n                            if i in set(subsetting_pos)]\n    feats = featurize.featurize_time_series(times=subset_times_list[:N],\n                                            values=subset_flux_list[:N],\n                                            features_to_use=['skew',\n                                                            'percent_beyond_1_std',\n                                                            'percent_difference_flux_percentile'\n                                                            ],\n                                            scheduler=None)\n    subset = train_series[train_series['object_id'].isin(\n        train_metadata['object_id'].iloc[subsetting_pos].iloc[:N])]\n    models = list(map(fit_multiband_freq, subset.groupby('object_id')))\n    feats['object_pos'] = subsetting_pos[:N]\n    feats['freq1_freq'] = [model.best_period for model in models]\n    return feats, models","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:45.296644Z","iopub.execute_input":"2022-03-31T01:46:45.296875Z","iopub.status.idle":"2022-03-31T01:46:45.309973Z","shell.execute_reply.started":"2022-03-31T01:46:45.296838Z","shell.execute_reply":"2022-03-31T01:46:45.308978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"unique_classes = train_metadata['target'].unique()\nunique_classes","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:45.312707Z","iopub.execute_input":"2022-03-31T01:46:45.313378Z","iopub.status.idle":"2022-03-31T01:46:45.328785Z","shell.execute_reply.started":"2022-03-31T01:46:45.313329Z","shell.execute_reply":"2022-03-31T01:46:45.327980Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_class_feats(label, N=10):\n    class_pos = train_metadata[train_metadata['target'] == label].index\n    class_feats, class_models = get_freq_features(N, class_pos)\n    return class_feats, class_models","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:45.329919Z","iopub.execute_input":"2022-03-31T01:46:45.330140Z","iopub.status.idle":"2022-03-31T01:46:45.338926Z","shell.execute_reply.started":"2022-03-31T01:46:45.330113Z","shell.execute_reply":"2022-03-31T01:46:45.337906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_phase_curves(feats, models, use_median_freq=False, hide_undetected=True, N=10):\n    for i in range(N):\n        freq = feats.loc[i, 'freq1_freq'].median()\n        freq_min = feats.loc[i, 'freq1_freq'].min()\n        freq_std = feats.loc[i, 'freq1_freq'].std()\n        skew = feats.loc[i, 'skew'].mean()\n        object_pos = int(feats.loc[i, 'object_pos'][0])\n        f, ax = plt.subplots(1, 2, figsize=(14, 4))\n        sample = train_series[train_series['object_id'] ==\n                              train_metadata['object_id'].iloc[object_pos]].copy()\n        colors = ['red', 'orange', 'yellow', 'green', 'blue', 'purple']\n        score = models[i].score(models[i].best_period)\n        \n        ax[0].scatter(x=sample['mjd'], \n                   y=sample['flux'], \n                   c=[colors[b] for b in sample['passband']],\n                   s=8, alpha=0.8)\n        ax[0].vlines(sample['mjd'], \n                  sample['flux'] - sample['flux_err'],\n                  sample['flux'] + sample['flux_err'],\n                  colors=[colors[b] for b in sample['passband']],\n                  linewidth=1, alpha=0.8)\n        \n        sample['phase'] = (sample['mjd'] / models[i].best_period) % 1\n        ax[1].scatter(x=sample['phase'], \n                   y=sample['flux'], \n                   c=[colors[b] for b in sample['passband']],\n                   s=8, alpha=0.8)\n        ax[1].vlines(sample['phase'], \n                  sample['flux'] - sample['flux_err'],\n                  sample['flux'] + sample['flux_err'],\n                  colors=[colors[b] for b in sample['passband']],\n                  linewidth=1, alpha=0.8)\n        x_range = np.linspace(sample['mjd'].min(), sample['mjd'].max(), 1000)\n        for band in range(6):\n            y = models[i].predict(x_range, band)\n            xs = (x_range / models[i].best_period) % 1\n            ords = np.argsort(xs)\n            ax[1].plot(xs[ords], y[ords], c=colors[band], alpha=0.4)\n        \n        title = ax[0].get_title()\n        ax[0].set_title('time')\n        ax[1].set_title('phase')\n        f.suptitle(title + f'object: {sample[\"object_id\"].iloc[0]}, '\n                   f'class: {train_metadata[\"target\"].iloc[object_pos]}\\n'\n                   f'period: {models[i].best_period: .4}, '\n                   f'period score: {score: .4}, '\n                   f'mean skew: {skew:.4}', y=1.1)\n        plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:45.342433Z","iopub.execute_input":"2022-03-31T01:46:45.343191Z","iopub.status.idle":"2022-03-31T01:46:45.363372Z","shell.execute_reply.started":"2022-03-31T01:46:45.343139Z","shell.execute_reply":"2022-03-31T01:46:45.362685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"warnings.simplefilter('ignore', UserWarning)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:45.715794Z","iopub.execute_input":"2022-03-31T01:46:45.716088Z","iopub.status.idle":"2022-03-31T01:46:45.721293Z","shell.execute_reply.started":"2022-03-31T01:46:45.716057Z","shell.execute_reply":"2022-03-31T01:46:45.720223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Class 42**","metadata":{}},{"cell_type":"code","source":"%%capture capt\nfeats, models = get_class_feats(42)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:46:48.337287Z","iopub.execute_input":"2022-03-31T01:46:48.337660Z","iopub.status.idle":"2022-03-31T01:53:42.971339Z","shell.execute_reply.started":"2022-03-31T01:46:48.337613Z","shell.execute_reply":"2022-03-31T01:53:42.969884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_phase_curves(feats, models)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:53:42.975344Z","iopub.execute_input":"2022-03-31T01:53:42.976261Z","iopub.status.idle":"2022-03-31T01:53:48.474909Z","shell.execute_reply.started":"2022-03-31T01:53:42.976163Z","shell.execute_reply":"2022-03-31T01:53:48.474057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Class 42 appears to be mostly flat. It is extragalactical, so high uncertainty happens in some cases. It most likely is a class of objects with \"burst\" events, although the actual burst is not always detected in the observation windows. When the burst happens, the object's magnitude increases dramatically across all bands and gradually falls back to normal levels in a few months. If the burst is indeed detected, it will be characterised by relatively low frequency std between bands together with very high detected periods. Other features like skewness can also be used to identify this class of objects, as a burst usually results in high skew in the light curve.","metadata":{}},{"cell_type":"markdown","source":"**Class 90**","metadata":{}},{"cell_type":"code","source":"%%capture capt\nfeats, models = get_class_feats(90)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T01:53:48.476145Z","iopub.execute_input":"2022-03-31T01:53:48.476359Z","iopub.status.idle":"2022-03-31T02:00:43.958226Z","shell.execute_reply.started":"2022-03-31T01:53:48.476333Z","shell.execute_reply":"2022-03-31T02:00:43.957097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_phase_curves(feats, models)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:00:43.961417Z","iopub.execute_input":"2022-03-31T02:00:43.962042Z","iopub.status.idle":"2022-03-31T02:00:49.642640Z","shell.execute_reply.started":"2022-03-31T02:00:43.961974Z","shell.execute_reply":"2022-03-31T02:00:49.641736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Class 90 is kind of similar to class 42 in that it is very likely a group of objects with \"burst\" behaviours, or objects with partially unobsered bursts. Like class 42, it is characterised by sudden peaks that usually last for a few months, high skew and often long fitted periods (or 1-day periods) with low score. So far we are still unable to tell the difference between class 90 and class 42.","metadata":{}},{"cell_type":"markdown","source":"choose class 42 and 90","metadata":{}},{"cell_type":"markdown","source":"# Data Preprocessing","metadata":{}},{"cell_type":"markdown","source":"# Transform into a binary classification","metadata":{}},{"cell_type":"code","source":"from sklearn import metrics\nfrom sklearn import preprocessing\nfrom sklearn.preprocessing import MultiLabelBinarizer\nimport numpy as np\nimport pandas as pd\ntrain_series = pd.read_csv('../input/training_set.csv')\ntrain_metadata = pd.read_csv('../input/training_set_metadata.csv')\n#test_data = pd.read_csv('../input/test_set_batch1.csv')\nobj_id_1 = np.random.RandomState(500).choice(train_metadata.object_id[train_metadata.target == 42].values)\nobj_df_1 = train_series[train_series.object_id == obj_id_1]\nobj_df_1_train = obj_df_1[:100]\nx_test_1 = obj_df_1[101:121]\nobj_id_2 = np.random.RandomState(500).choice(train_metadata.object_id[train_metadata.target == 90].values)\nobj_df_2 = train_series[train_series.object_id == obj_id_2]\nobj_df_2_train = obj_df_2[:100]\n#obj_id_1\n#obj_df_1\nx_test_2 = obj_df_2[101:121]\n#m = np.array(obj_df_1)\n#m\nx_test_1[:5]","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:00:49.643863Z","iopub.execute_input":"2022-03-31T02:00:49.644107Z","iopub.status.idle":"2022-03-31T02:00:51.070062Z","shell.execute_reply.started":"2022-03-31T02:00:49.644079Z","shell.execute_reply":"2022-03-31T02:00:51.069133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"obj_df_2[:5]","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:00:51.071219Z","iopub.execute_input":"2022-03-31T02:00:51.071563Z","iopub.status.idle":"2022-03-31T02:00:51.083033Z","shell.execute_reply.started":"2022-03-31T02:00:51.071529Z","shell.execute_reply":"2022-03-31T02:00:51.082406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"mb_1 = MultiLabelBinarizer()\nx_train_1 = mb_1.fit_transform(obj_df_1)\nx_test_1 = mb_1.fit_transform(test_id_1)\n\nmb_2 = MultiLabelBinarizer()\nx_train_2 = mb_2.fit_transform(obj_df_2)\nx_test_2 = mb_1.fit_transform(test_id_2)","metadata":{}},{"cell_type":"markdown","source":"# **Time Series Transformations**\n\n\ntransform each time series into the same set of derived quantities. \nThese include: \n1. the number of measurements; \n1. the minimum, maximum, mean, median, standard deviation, and skew of flux;\n1. the minimum, maximum, mean, median, standard deviation, and skew of flux error; \n1. the sum of the ratio between flux and flux error;\n1. the skew of the ratio between flux and flux error;\n1. the sum of the flux times squared flux ratio;\n1. the skew of the flux times squared flux ratio;\n1. the mean time between measurements; \n1. the maximum time between measurements;\n1. spectroscopic redshifts for the host galaxy; \n1. photometric redshifts for the host galaxy;\n1. the position of each object in the sky; \n1. the first two Fourier coefficients for each band, \n1. as well as kurtosis and skewness. ","metadata":{}},{"cell_type":"code","source":"def agg_func(x):\n    d = {}\n    flux, dflux = x[\"flux\"], x[\"flux_err\"]\n    flux_mean = np.sum(flux*np.square(flux/dflux))/np.sum(np.square(flux/dflux))\n    d[\"flux_mean\"] = flux_mean\n    d[\"flux_std\"] = np.std(flux/flux_mean, ddof = 1)\n    d[\"flux_amp\"] = (np.max(flux) - np.min(flux))/flux_mean\n    d[\"flux_beyond\"] = np.sum(np.abs(flux - flux_mean) > np.std(flux, ddof = 1))/flux.shape[0]\n    d[\"flux_mad\"] = np.median(np.abs((flux - np.median(flux))/flux_mean))\n    d[\"flux_skew\"] = skew(flux)\n    colnames = [\"flux_mean\", \"flux_std\", \"flux_amp\", \"flux_mad\", \"flux_beyond\", \"flux_skew\"]\n    return pd.Series(d, index = colnames)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def _finalize(self):\n    '''Store individual passband fluxes as object attributes'''\n    # in this example, we'll use the weighted mean to normalize the features\n    weighted_mean = lambda flux, dflux: np.sum(flux*(flux/dflux)**2)/np.sum((flux/dflux)**2)\n\n    # define some functions to compute simple descriptive statistics\n    normalized_flux_std = lambda flux, wMeanFlux: np.std(flux/wMeanFlux, ddof = 1)\n    normalized_amplitude = lambda flux, wMeanFlux: (np.max(flux) - np.min(flux))/wMeanFlux\n    normalized_MAD = lambda flux, wMeanFlux: np.median(np.abs((flux - np.median(flux))/wMeanFlux))\n    beyond_1std = lambda flux, wMeanFlux: sum(np.abs(flux - wMeanFlux) > np.std(flux, ddof = 1))/len(flux)\n\n    for pb in self._passbands:\n        ind = self.DFlc['passband'] == pb\n        pbname = self._pbnames[pb]\n\n        if len(self.DFlc[ind]) == 0:\n            setattr(self, f'{pbname}Std', np.nan)\n            setattr(self, f'{pbname}Amp', np.nan)\n            setattr(self, f'{pbname}MAD', np.nan)\n            setattr(self, f'{pbname}Beyond', np.nan)\n            setattr(self, f'{pbname}Skew', np.nan)\n            continue\n\n        f  = self.DFlc['flux'][ind]\n        df = self.DFlc['flux_err'][ind]\n        m  = weighted_mean(f, df)\n\n        # we'll save the measurements in each passband to simplify access.\n        setattr(self, f'{pbname}Flux', f)\n        setattr(self, f'{pbname}FluxUnc', df)\n        setattr(self, f'{pbname}Mean', m)\n\n        # compute the features\n        std = normalized_flux_std(f, df)\n        amp = normalized_amplitude(f, m)\n        mad = normalized_MAD(f, m)\n        beyond = beyond_1std(f, m)\n        skew = spstat.skew(f) \n\n        # and save the features\n        setattr(self, f'{pbname}Std', std)\n        setattr(self, f'{pbname}Amp', amp)\n        setattr(self, f'{pbname}MAD', mad)\n        setattr(self, f'{pbname}Beyond', beyond)\n        setattr(self, f'{pbname}Skew', skew)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_features(self):\n    '''Return all the features for this object'''\n    variables = ['Std', 'Amp', 'MAD', 'Beyond', 'Skew']\n    feats = []\n    for i, pb in enumerate(self._passbands):\n        pbname = self._pbnames[pb]\n        feats += [getattr(self, f'{pbname}{x}', np.nan) for x in variables]\n    return feats","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_features_to_agg(df):\n    # CPMP using the following feature was really silliy :)\n    # df['mjd_diff'] = df['mjd_max'] - df['mjd_min']\n    # see https://www.kaggle.com/c/PLAsTiCC-2018/discussion/69696\n    \n    # The others may be useful\n    df['flux_diff'] = df['flux_max'] - df['flux_min']\n    df['flux_dif2'] = (df['flux_max'] - df['flux_min']) / df['flux_mean']\n    df['flux_w_mean'] = df['flux_by_flux_ratio_sq_sum'] / df['flux_ratio_sq_sum']\n    df['flux_dif3'] = (df['flux_max'] - df['flux_min']) / df['flux_w_mean']\n\n    # del df['mjd_max'], df['mjd_min']\n\n    return df","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Logscale transformation**","metadata":{}},{"cell_type":"markdown","source":"from sklearn.preprocessing import FunctionTransformer\n\ntransformer = FunctionTransformer(np.log10, validate=True)","metadata":{}},{"cell_type":"markdown","source":"**MinMaxScaler Normalization/scaling and outliers**","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import MinMaxScaler\ntransformer = FunctionTransformer(np.log10)\nLogscale_data = transformer.transform(obj_df_1)\nLogscale_data = np.nan_to_num(Logscale_data) \nprint(Logscale_data[:5])\n\"\"\"\"\"\"\n# X_std = (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0))\n# X_scaled = X_std * (max - min) + min\n\"\"\"\"\"\"\nscaler = MinMaxScaler(feature_range=(-(np.pi/2),(np.pi/2)))\n\n# transform data\n\nNormal_data = scaler.fit_transform(Logscale_data)       \nNormal_data = np.array(Normal_data)\nprint(Normal_data[:5])","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:12:05.476817Z","iopub.execute_input":"2022-03-31T02:12:05.477164Z","iopub.status.idle":"2022-03-31T02:12:05.489088Z","shell.execute_reply.started":"2022-03-31T02:12:05.477125Z","shell.execute_reply":"2022-03-31T02:12:05.488240Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.preprocessing import FunctionTransformer\nfrom sklearn.preprocessing import MinMaxScaler\ndef data_preprocessing(origin_data):\n    transformer = FunctionTransformer(np.log10)\n    Logscale_data = transformer.transform(origin_data)\n    Logscale_data = np.nan_to_num(Logscale_data) \n    \"\"\"\"\"\"\n    # X_std = (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0))\n    # X_scaled = X_std * (max - min) + min\n    \"\"\"\"\"\"\n    scaler = MinMaxScaler(feature_range=(-(np.pi/2),(np.pi/2)))\n\n    # transform data\n\n    Normal_data = scaler.fit_transform(Logscale_data)                                          \n    return Normal_data","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:00:51.100928Z","iopub.execute_input":"2022-03-31T02:00:51.101164Z","iopub.status.idle":"2022-03-31T02:00:51.109528Z","shell.execute_reply.started":"2022-03-31T02:00:51.101136Z","shell.execute_reply":"2022-03-31T02:00:51.108425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Quantum Circuit test","metadata":{}},{"cell_type":"markdown","source":"**iswapgate^0.5 Test**","metadata":{}},{"cell_type":"code","source":"#iswapgate^0.5\nn = 2\n# 初始化电路\ncircuit_t = UAnsatz(n)\ntheta1 = paddle.to_tensor(np.array([np.pi/2], np.float64))\ntheta2 = paddle.to_tensor(np.array([(7*np.pi)/4], np.float64))\ncircuit_t.sdg(0)\ncircuit_t.h(0)\ncircuit_t.sdg(0)\ncircuit_t.rz(theta1[0], 0)\ncircuit_t.cnot([0, 1])\ncircuit_t.sdg(0)\ncircuit_t.h(0)\ncircuit_t.sdg(0)\ncircuit_t.sdg(1)\ncircuit_t.h(1)\ncircuit_t.sdg(1)\ncircuit_t.rz(theta2[0], 0)\ncircuit_t.rz(theta2[0], 1)\ncircuit_t.sdg(0)\ncircuit_t.h(0)\ncircuit_t.sdg(0)\ncircuit_t.rz(theta1[0], 0)\ncircuit_t.cnot([0, 1])\ncircuit_t.sdg(0)\ncircuit_t.h(0)\ncircuit_t.sdg(0)\ncircuit_t.sdg(1)\ncircuit_t.h(1)\ncircuit_t.sdg(1)\nprint(circuit_t)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:00:51.110717Z","iopub.execute_input":"2022-03-31T02:00:51.111512Z","iopub.status.idle":"2022-03-31T02:00:51.135852Z","shell.execute_reply.started":"2022-03-31T02:00:51.111438Z","shell.execute_reply":"2022-03-31T02:00:51.135067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Circuit Encode Test**","metadata":{}},{"cell_type":"code","source":"# 量子比特的数量\nn = 2\n# 初始化电路\ncircuit = UAnsatz(n)\n# x 是经典信息\nx1 = paddle.to_tensor([-1.45, 3, 2, -0.05], 'float64')\nx2 = paddle.to_tensor([-1.45, 3, 2, -0.05], 'float64')\nx3 = paddle.to_tensor([-1.45, 3, 2, -0.05], 'float64')\n\nfor i in range(n):\n    # 加上一层 Hadamard 门\n    circuit.superposition_layer()\n    # 加上一层旋转门 Rz\n    for j in range(n):\n        circuit.rz(x1[j] ,j)\n    # 加上一层旋转门 Ry\n        circuit.ry(x2[j] ,j)\n    # 加上一层旋转门 Rz\n        circuit.rz(x3[j] ,j)\nprint(circuit)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:08:52.822289Z","iopub.execute_input":"2022-03-31T02:08:52.822729Z","iopub.status.idle":"2022-03-31T02:08:52.880356Z","shell.execute_reply.started":"2022-03-31T02:08:52.822678Z","shell.execute_reply":"2022-03-31T02:08:52.879351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 导入训练集和测试集","metadata":{}},{"cell_type":"code","source":"X_train = data_preprocessing(np.array(obj_df_1_train))\nX_train = np.reshape(X_train, 600)\ny_train = data_preprocessing(np.array(obj_df_2_train))\ny_train = np.reshape(y_train, 600)\nX_test = data_preprocessing(np.array(x_test_1))\nX_test = np.reshape(X_test, 120)\ny_test = data_preprocessing(np.array(x_test_2))\ny_test = np.reshape(y_test, 120)\n#X_train\nprint(X_train[:20])\nprint(y_train[:5])","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:03:04.762212Z","iopub.execute_input":"2022-03-31T02:03:04.762588Z","iopub.status.idle":"2022-03-31T02:03:04.777946Z","shell.execute_reply.started":"2022-03-31T02:03:04.762550Z","shell.execute_reply":"2022-03-31T02:03:04.776933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 初始化进度条\nbar_format_string = '{l_bar}{bar}|[{elapsed}<{remaining}, ' '{rate_fmt}{postfix}]'\npbar = tqdm(total=100, bar_format=bar_format_string)\npbar.close()\nclear_output()","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:03:15.002055Z","iopub.execute_input":"2022-03-31T02:03:15.002941Z","iopub.status.idle":"2022-03-31T02:03:15.017114Z","shell.execute_reply.started":"2022-03-31T02:03:15.002899Z","shell.execute_reply":"2022-03-31T02:03:15.016279Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Quantum Circuit**","metadata":{}},{"cell_type":"code","source":"#x = paddle.to_tensor([0,0,0,0,0,0], 'float64')\n#X_train = np.array(X_train)\n#y_train = np.array(y_train)\nx1 = paddle.to_tensor(np.reshape(X_train, 600), 'float64')\nx2 = paddle.to_tensor(np.reshape(y_train, 600), 'float64')\ntheta = np.array([np.pi], np.float64)\ntheta = paddle.to_tensor(theta)\n# 量子比特的数量等于经典信息的长度\nn = 6\nr = 20\ntheta1 = paddle.to_tensor(np.array([np.pi/2], np.float64))\ntheta2 = paddle.to_tensor(np.array([(7*np.pi)/4], np.float64))\n# 初始化电路\ncircuit = UAnsatz(n)\n#x = paddle.to_tensor([lambda x : x, for x in range(n)], 'float64')\nfor q in range(r):\n    # 加上一层 Hadamard 门\n    circuit.superposition_layer()\n    for i in range(n):\n        for j in range(n):\n        # 加上一层旋转门 Rz\n            circuit.rz(x1[j] ,j)\n        # 加上一层旋转门 Ry\n            circuit.ry(x1[j] ,j)\n        # 加上一层旋转门 Rz\n            circuit.rz(x1[j] ,j)\n        \"\"\"\"\"\"\n        #iswap^0.5\n        for j in range(n-1):\n            circuit.sdg(j)\n            circuit.h(j)\n            circuit.sdg(j)\n            circuit.rz(theta1[0], j)\n            circuit.cnot([j, j+1])\n            circuit.sdg(j)\n            circuit.h(j)\n            circuit.sdg(j)\n            circuit.sdg(j+1)\n            circuit.h(j+1)\n            circuit.sdg(j+1)\n            circuit.rz(theta2[0], j)\n            circuit.rz(theta2[0], j+1)\n            circuit.sdg(j)\n            circuit.h(j)\n            circuit.sdg(j)\n            circuit.rz(theta1[0], j)\n            circuit.cnot([j, j+1])\n            circuit.sdg(j)\n            circuit.h(j)\n            circuit.sdg(j)\n            circuit.sdg(j+1)\n            circuit.h(j+1)\n            circuit.sdg(j+1)\n\n        \"\"\"\"\"\"\n        for j in range(n):\n            circuit.ry(x1[j] ,j)\n            circuit.rz(x1[j] ,j)\n            circuit.h(j)\n    #invert\n    for i in range(n):\n        for j in range(n):\n            circuit.h(j)\n            circuit.rz(x2[j] ,j)\n            circuit.ry(x2[j] ,j)\n        #iswap^0.5\n        for j in range(n-1):\n            circuit.sdg(j)\n            circuit.h(j)\n            circuit.sdg(j)\n            circuit.rz(theta1[0], j)\n            circuit.cnot([j, j+1])\n            circuit.sdg(j)\n            circuit.h(j)\n            circuit.sdg(j)\n            circuit.sdg(j+1)\n            circuit.h(j+1)\n            circuit.sdg(j+1)\n            circuit.rz(theta2[0], j)\n            circuit.rz(theta2[0], j+1)\n            circuit.sdg(j)\n            circuit.h(j)\n            circuit.sdg(j)\n            circuit.rz(theta1[0], j)\n            circuit.cnot([j, j+1])\n            circuit.sdg(j)\n            circuit.h(j)\n            circuit.sdg(j)\n            circuit.sdg(j+1)\n            circuit.h(j+1)\n            circuit.sdg(j+1)\n        \"\"\"\"\"\"\n        # 加上一层旋转门 Rz\n        for j in range(n):\n        # 加上一层旋转门 Rz\n            circuit.rz(x2[j] ,j)\n        # 加上一层旋转门 Ry\n            circuit.ry(x2[j] ,j)\n        # 加上一层旋转门 Rz\n            circuit.rz(x2[j] ,j)\n            circuit.h(j)\n        \"\"\"\"\"\"\n#print(circuit)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:03:26.347410Z","iopub.execute_input":"2022-03-31T02:03:26.347794Z","iopub.status.idle":"2022-03-31T02:03:26.883217Z","shell.execute_reply.started":"2022-03-31T02:03:26.347756Z","shell.execute_reply":"2022-03-31T02:03:26.882314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fin_state = circuit.run_state_vector()\nprint([np.round(i, 5) for i in fin_state.numpy()])","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:03:30.510006Z","iopub.execute_input":"2022-03-31T02:03:30.510689Z","iopub.status.idle":"2022-03-31T02:04:10.731984Z","shell.execute_reply.started":"2022-03-31T02:03:30.510648Z","shell.execute_reply":"2022-03-31T02:04:10.730998Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 返回测量结果为 0...0 的概率\npro = (fin_state[0].conj() * fin_state[0]).real().numpy()[0]\npro","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:04:10.733711Z","iopub.execute_input":"2022-03-31T02:04:10.733947Z","iopub.status.idle":"2022-03-31T02:04:10.740922Z","shell.execute_reply.started":"2022-03-31T02:04:10.733917Z","shell.execute_reply":"2022-03-31T02:04:10.740250Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def q_kernel_estimator(x1, x2):\n    #x1 = paddle.to_tensor(np.reshape(X_train, 600), 'float64')\n    #x2 = paddle.to_tensor(np.reshape(y_train, 600), 'float64')\n    theta = np.array([np.pi], np.float64)\n    theta = paddle.to_tensor(theta)\n    # 量子比特的数量等于经典信息的长度\n    n = 6\n    r = 20\n    theta1 = paddle.to_tensor(np.array([np.pi/2], np.float64))\n    theta2 = paddle.to_tensor(np.array([(7*np.pi)/4], np.float64))\n    # 初始化电路\n    circuit = UAnsatz(n)\n    #x = paddle.to_tensor([lambda x : x, for x in range(n)], 'float64')\n    for q in range(r):\n        # 加上一层 Hadamard 门\n        circuit.superposition_layer()\n        for i in range(n):\n            for j in range(n):\n            # 加上一层旋转门 Rz\n                circuit.rz(x1[j] ,j)\n            # 加上一层旋转门 Ry\n                circuit.ry(x1[j] ,j)\n            # 加上一层旋转门 Rz\n                circuit.rz(x1[j] ,j)\n            \"\"\"\"\"\"\n            #iswap^0.5\n            for j in range(n-1):\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.rz(theta1[0], j)\n                circuit.cnot([j, j+1])\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.sdg(j+1)\n                circuit.h(j+1)\n                circuit.sdg(j+1)\n                circuit.rz(theta2[0], j)\n                circuit.rz(theta2[0], j+1)\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.rz(theta1[0], j)\n                circuit.cnot([j, j+1])\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.sdg(j+1)\n                circuit.h(j+1)\n                circuit.sdg(j+1)\n\n            \"\"\"\"\"\"\n            for j in range(n):\n                circuit.ry(x1[j] ,j)\n                circuit.rz(x1[j] ,j)\n                circuit.h(j)\n        #invert\n        for i in range(n):\n            for j in range(n):\n                circuit.h(j)\n                circuit.rz(x2[j] ,j)\n                circuit.ry(x2[j] ,j)\n            #iswap^0.5\n            for j in range(n-1):\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.rz(theta1[0], j)\n                circuit.cnot([j, j+1])\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.sdg(j+1)\n                circuit.h(j+1)\n                circuit.sdg(j+1)\n                circuit.rz(theta2[0], j)\n                circuit.rz(theta2[0], j+1)\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.rz(theta1[0], j)\n                circuit.cnot([j, j+1])\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.sdg(j+1)\n                circuit.h(j+1)\n                circuit.sdg(j+1)\n            \"\"\"\"\"\"\n            # 加上一层旋转门 Rz\n            for j in range(n):\n            # 加上一层旋转门 Rz\n                circuit.rz(x2[j] ,j)\n            # 加上一层旋转门 Ry\n                circuit.ry(x2[j] ,j)\n            # 加上一层旋转门 Rz\n                circuit.rz(x2[j] ,j)\n                circuit.h(j)\n            \"\"\"\"\"\"\n    # 用态矢量模式运行电路\n    fin_state = circuit.run_state_vector()\n    # 更新进度条\n    global pbar\n    global N\n    pbar.update(100/N)\n    # 返回测量结果为 0...0 的概率\n    return (fin_state[0].conj() * fin_state[0]).real().numpy()[0]","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:04:10.742192Z","iopub.execute_input":"2022-03-31T02:04:10.742591Z","iopub.status.idle":"2022-03-31T02:04:10.767282Z","shell.execute_reply.started":"2022-03-31T02:04:10.742552Z","shell.execute_reply":"2022-03-31T02:04:10.766583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 创建进度条，并设置所需要的量子核函数计算数量 N\ndef q_kernel_matrix(X1, X2):\n    return np.array([[q_kernel_estimator(x1, x2) for x2 in X2] for x1 in X1])\npbar = tqdm(total=100, \n            desc='训练 QKE-SVM 并分类中', \n            bar_format=bar_format_string)\nN = len(X_train) ** 2 + len(X_train) ** 2 + len(X_train) * len(X_test)\n\n# 创建一个具有量子核函数的支持向量机\nsvm_qke = svm.SVC(kernel=q_kernel_matrix)\n\n# 根据训练数据计算支持向量机的决策平面\nsvm_qke.fit(X_train, y_train)\n\n# 计算支持向量机分别对于训练数据和测试数据的分类预测值\npredict_svm_qke_train = svm_qke.predict(X_train)\npredict_svm_qke_test = svm_qke.predict(X_test)\n\n# 计算准确率\naccuracy_train = np.array(predict_svm_qke_train == y_train, dtype=int).sum()/len(y_train)\naccuracy_test = np.array(predict_svm_qke_test == y_test, dtype=int).sum()/len(y_test)","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:04:10.769047Z","iopub.execute_input":"2022-03-31T02:04:10.769513Z","iopub.status.idle":"2022-03-31T02:04:10.816065Z","shell.execute_reply.started":"2022-03-31T02:04:10.769440Z","shell.execute_reply":"2022-03-31T02:04:10.814801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Transformed","metadata":{}},{"cell_type":"code","source":"from sklearn import metrics\nfrom sklearn import preprocessing\nfrom sklearn.preprocessing import MultiLabelBinarizer\nimport numpy as np\nimport pandas as pd\ntrain_series = pd.read_csv('../input/training_set.csv')\ntrain_metadata = pd.read_csv('../input/training_set_metadata.csv')\nobj_id_1 = np.random.RandomState(500).choice(train_metadata.object_id[train_metadata.target == 42].values)\nobj_df_1 = train_series[train_series.object_id == obj_id_1]\nobj_df_1_train = obj_df_1[:100]\nx_test_1 = obj_df_1[101:121]\nobj_id_2 = np.random.RandomState(500).choice(train_metadata.object_id[train_metadata.target == 90].values)\nobj_df_2 = train_series[train_series.object_id == obj_id_2]\nobj_df_2_train = obj_df_2[:100]\n#obj_id_1\n#obj_df_1\nx_test_2 = obj_df_2[101:121]\n#m = np.array(obj_df_1)\n#m\nx_test_1[:5]","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:00:51.188357Z","iopub.status.idle":"2022-03-31T02:00:51.188923Z","shell.execute_reply.started":"2022-03-31T02:00:51.188741Z","shell.execute_reply":"2022-03-31T02:00:51.188762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_series = pd.read_csv('../input/training_set.csv')\ntrain_metadata = pd.read_csv('../input/training_set_metadata.csv')\nobj_id_1 = np.random.RandomState(500).choice(train_metadata.object_id[train_metadata.target == 42].values)\nobj_df_1 = train_series[train_series.object_id == obj_id_1]\nobj_df_1_train = obj_df_1[:100]\nx_test_1 = obj_df_1[101:121]\nobj_id_2 = np.random.RandomState(500).choice(train_metadata.object_id[train_metadata.target == 90].values)\nobj_df_2 = train_series[train_series.object_id == obj_id_2]\nobj_df_2_train = obj_df_2[:100]\n#obj_id_1\n#obj_df_1\nx_test_2 = obj_df_2[101:121]\n#m = np.array(obj_df_1)\n#m\nx_test_1[:5]","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:00:51.189981Z","iopub.status.idle":"2022-03-31T02:00:51.190362Z","shell.execute_reply.started":"2022-03-31T02:00:51.190179Z","shell.execute_reply":"2022-03-31T02:00:51.190204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"obj_df_1_train['flux_ratio_sq'] = np.power(obj_df_1_train['flux'] / obj_df_1_train['flux_err'], 2.0)\nobj_df_1_train['flux_by_flux_ratio_sq'] = obj_df_1_train['flux'] * obj_df_1_train['flux_ratio_sq']","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:05:33.495497Z","iopub.execute_input":"2022-03-31T02:05:33.495863Z","iopub.status.idle":"2022-03-31T02:05:33.506474Z","shell.execute_reply.started":"2022-03-31T02:05:33.495822Z","shell.execute_reply":"2022-03-31T02:05:33.505521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_features = obj_df_1_train.columns[1:]","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:05:35.987822Z","iopub.execute_input":"2022-03-31T02:05:35.988209Z","iopub.status.idle":"2022-03-31T02:05:35.993910Z","shell.execute_reply.started":"2022-03-31T02:05:35.988168Z","shell.execute_reply":"2022-03-31T02:05:35.992533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"groupObjects = obj_df_1_train.groupby('object_id')[data_features]\n\n#print(\"Add constant object features\")\nfeatures = train_metadata.drop(['target'], axis=1)\n\nprint(\"Add sum of mutable object features\")\nfeatures = pd.merge(features, groupObjects.agg('sum'), how='right', on='object_id', suffixes=['', '_sum'])\n\nprint(\"Add mean of mutable object features\")\nfeatures = pd.merge(features, groupObjects.agg('mean'), how='right', on='object_id', suffixes=['', '_mean'])\n\nprint(\"Add median of mutable features\")\nfeatures = pd.merge(features, groupObjects.agg('median'), how='right', on='object_id', suffixes=['', '_median'])\n\nprint(\"Add minimum of mutable features\")\nfeatures = pd.merge(features, groupObjects.agg('min'), how='right', on='object_id', suffixes=['', '_min'])\n\nprint(\"Add maximum of mutable features\")\nfeatures = pd.merge(features, groupObjects.agg('max'), how='right', on='object_id', suffixes=['', '_max'])\n\nprint(\"Add range of mutable features\")\nfeatures = pd.merge(features, groupObjects.agg(lambda x: max(x) - min(x)), how='right', on='object_id', suffixes=['', '_range'])\n\nprint(\"Add standard deviation of mutable features\")\nfeatures = pd.merge(features, groupObjects.agg('std'), how='right', on='object_id', suffixes=['', '_stddev'])\n\nprint(\"Add skew of mutable features\")\nfeatures = pd.merge(features, groupObjects.agg('skew'), how='right', on='object_id', suffixes=['', '_skew'])","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:05:37.574499Z","iopub.execute_input":"2022-03-31T02:05:37.574899Z","iopub.status.idle":"2022-03-31T02:05:37.668228Z","shell.execute_reply.started":"2022-03-31T02:05:37.574855Z","shell.execute_reply":"2022-03-31T02:05:37.666900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features = features.fillna(features.mean())","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:05:39.924784Z","iopub.execute_input":"2022-03-31T02:05:39.925149Z","iopub.status.idle":"2022-03-31T02:05:39.947233Z","shell.execute_reply.started":"2022-03-31T02:05:39.925116Z","shell.execute_reply":"2022-03-31T02:05:39.946016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features","metadata":{"execution":{"iopub.status.busy":"2022-03-31T02:05:41.473312Z","iopub.execute_input":"2022-03-31T02:05:41.473901Z","iopub.status.idle":"2022-03-31T02:05:41.510801Z","shell.execute_reply.started":"2022-03-31T02:05:41.473831Z","shell.execute_reply":"2022-03-31T02:05:41.509326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.preprocessing import FunctionTransformer\nfrom sklearn.preprocessing import MinMaxScaler\ndef data_preprocessing(origin_data):\n    transformer = FunctionTransformer(np.log10)\n    Logscale_data = transformer.transform(origin_data)\n    Logscale_data = np.nan_to_num(Logscale_data) \n    \"\"\"\"\"\"\n    # X_std = (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0))\n    # X_scaled = X_std * (max - min) + min\n    \"\"\"\"\"\"\n    scaler = MinMaxScaler(feature_range=(-(np.pi/2),(np.pi/2)))\n\n    # transform data\n\n    Normal_data = scaler.fit_transform(Logscale_data)                                          \n    return Normal_data","metadata":{"execution":{"iopub.status.busy":"2022-03-30T16:23:59.663882Z","iopub.execute_input":"2022-03-30T16:23:59.664206Z","iopub.status.idle":"2022-03-30T16:23:59.670601Z","shell.execute_reply.started":"2022-03-30T16:23:59.664159Z","shell.execute_reply":"2022-03-30T16:23:59.669807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train = data_preprocessing(np.array(features))\nX_train","metadata":{"execution":{"iopub.status.busy":"2022-03-30T16:24:01.393396Z","iopub.execute_input":"2022-03-30T16:24:01.393657Z","iopub.status.idle":"2022-03-30T16:24:01.402981Z","shell.execute_reply.started":"2022-03-30T16:24:01.39363Z","shell.execute_reply":"2022-03-30T16:24:01.401996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train = data_preprocessing(np.array(obj_df_1_train))\nX_train = np.reshape(X_train, 600)\ny_train = data_preprocessing(np.array(obj_df_2_train))\ny_train = np.reshape(y_train, 600)\nX_test = data_preprocessing(np.array(x_test_1))\nX_test = np.reshape(X_test, 120)\ny_test = data_preprocessing(np.array(x_test_2))\ny_test = np.reshape(y_test, 120)\n#X_train\nprint(X_train[:20])\nprint(y_train[:5])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def q_kernel_estimator(x1, x2):\n    #x1 = paddle.to_tensor(np.reshape(X_train, 600), 'float64')\n    #x2 = paddle.to_tensor(np.reshape(y_train, 600), 'float64')\n    theta = np.array([np.pi], np.float64)\n    theta = paddle.to_tensor(theta)\n    # 量子比特的数量等于经典信息的长度\n    n = 6\n    r = 20\n    theta1 = paddle.to_tensor(np.array([np.pi/2], np.float64))\n    theta2 = paddle.to_tensor(np.array([(7*np.pi)/4], np.float64))\n    # 初始化电路\n    circuit = UAnsatz(n)\n    #x = paddle.to_tensor([lambda x : x, for x in range(n)], 'float64')\n    for q in range(r):\n        # 加上一层 Hadamard 门\n        circuit.superposition_layer()\n        for i in range(n):\n            for j in range(n):\n            # 加上一层旋转门 Rz\n                circuit.rz(x1[j] ,j)\n            # 加上一层旋转门 Ry\n                circuit.ry(x1[j] ,j)\n            # 加上一层旋转门 Rz\n                circuit.rz(x1[j] ,j)\n            \"\"\"\"\"\"\n            #iswap^0.5\n            for j in range(n-1):\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.rz(theta1[0], j)\n                circuit.cnot([j, j+1])\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.sdg(j+1)\n                circuit.h(j+1)\n                circuit.sdg(j+1)\n                circuit.rz(theta2[0], j)\n                circuit.rz(theta2[0], j+1)\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.rz(theta1[0], j)\n                circuit.cnot([j, j+1])\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.sdg(j+1)\n                circuit.h(j+1)\n                circuit.sdg(j+1)\n\n            \"\"\"\"\"\"\n            for j in range(n):\n                circuit.ry(x1[j] ,j)\n                circuit.rz(x1[j] ,j)\n                circuit.h(j)\n        #invert\n        for i in range(n):\n            for j in range(n):\n                circuit.h(j)\n                circuit.rz(x2[j] ,j)\n                circuit.ry(x2[j] ,j)\n            #iswap^0.5\n            for j in range(n-1):\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.rz(theta1[0], j)\n                circuit.cnot([j, j+1])\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.sdg(j+1)\n                circuit.h(j+1)\n                circuit.sdg(j+1)\n                circuit.rz(theta2[0], j)\n                circuit.rz(theta2[0], j+1)\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.rz(theta1[0], j)\n                circuit.cnot([j, j+1])\n                circuit.sdg(j)\n                circuit.h(j)\n                circuit.sdg(j)\n                circuit.sdg(j+1)\n                circuit.h(j+1)\n                circuit.sdg(j+1)\n            \"\"\"\"\"\"\n            # 加上一层旋转门 Rz\n            for j in range(n):\n            # 加上一层旋转门 Rz\n                circuit.rz(x2[j] ,j)\n            # 加上一层旋转门 Ry\n                circuit.ry(x2[j] ,j)\n            # 加上一层旋转门 Rz\n                circuit.rz(x2[j] ,j)\n                circuit.h(j)\n            \"\"\"\"\"\"\n    # 用态矢量模式运行电路\n    fin_state = circuit.run_state_vector()\n    # 更新进度条\n    global pbar\n    global N\n    pbar.update(100/N)\n    # 返回测量结果为 0...0 的概率\n    return (fin_state[0].conj() * fin_state[0]).real().numpy()[0]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.svm import SVR\n# 创建进度条，并设置所需要的量子核函数计算数量 N\ndef q_kernel_matrix(X1, X2):\n    return np.array([[q_kernel_estimator(x1, x2) for x2 in X2] for x1 in X1])\npbar = tqdm(total=100, \n            desc='训练 QKE-SVM 并分类中', \n            bar_format=bar_format_string)\nN = len(X_train) ** 2 + len(X_train) ** 2 + len(X_train) * len(X_test)\n\n# 创建一个具有量子核函数的支持向量机\nsvm_qke = svm.SVC(kernel=q_kernel_matrix)\n\n# 根据训练数据计算支持向量机的决策平面\nsvm_qke.fit(X_train, y_train)\n\n# 计算支持向量机分别对于训练数据和测试数据的分类预测值\npredict_svm_qke_train = svm_qke.predict(X_train)\npredict_svm_qke_test = svm_qke.predict(X_test)\n\n# 计算准确率\naccuracy_train = np.array(predict_svm_qke_train == y_train, dtype=int).sum()/len(y_train)\naccuracy_test = np.array(predict_svm_qke_test == y_test, dtype=int).sum()/len(y_test)","metadata":{},"execution_count":null,"outputs":[]}]}