{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install tensorflow-gpu -U\n# !pip install tensorflow-gpu==2.0a0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import tensorflow as tf\nprint(tf.__version__)\n#is it using the gpu?\ntf.test.is_gpu_available(\n    cuda_only=False,\n    min_cuda_compute_capability=None\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# !pip install keras -U","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install megnet ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip show keras","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"conda install -c openbabel openbabel","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# import sys\n# !conda install --yes --prefix {sys.prefix} -c openbabel openbabel==2.4.1\n","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"import sys\n!conda install --yes --prefix {sys.prefix} -c rdkit rdkit==2018.03.3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sys.prefix","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# # only reading 10% of data for debug\n# train = pd.read_csv('../input/train.csv')\n# test = pd.read_csv('../input/champs-scalar-coupling/test.csv')\n# len(train['molecule_name'].unique())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test['id'].count()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# structure = pd.read_csv('../input/structures.csv')\n# len(structure['molecule_name'].unique())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# from keras_radam import RAdam","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train = pd.read_csv('../input/champs-scalar-coupling/train.csv')\nmax_value =204.88 \nmin_value = -36.2186\n\n\n\nprint(max_value,min_value)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"graph code"},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"Abstract classes and utility operations for building graph representations and\ndata loaders (known as Sequence objects in Keras).\n\nMost users will not need to interact with this module.\"\"\"\nimport json\n\n\nfrom keras.callbacks import ModelCheckpoint\nfrom keras.optimizers import Adam\nfrom keras.layers import Dense, Input, Concatenate, Add, Embedding, Dropout\nfrom megnet.layers import MEGNetLayer, Set2Set\nfrom megnet.activations import softplus2\nfrom keras.regularizers import l2\nfrom keras.backend import int_shape\nfrom megnet.callbacks import ModelCheckpointMAE, ManualStop, ReduceLRUponNan\nfrom megnet.data.graph import GraphBatchDistanceConvert, GraphBatchGenerator, GaussianDistance\nfrom megnet.data.crystal import CrystalGraph\nfrom megnet.utils.preprocessing import DummyScaler\nimport numpy as np\nimport os\nfrom warnings import warn\nfrom monty.serialization import dumpfn, loadfn\nfrom abc import ABCMeta, abstractmethod\nfrom operator import itemgetter\nfrom megnet.utils.general import expand_1st, to_list\nfrom monty.json import MSONable\nfrom megnet.data import local_env\nfrom inspect import signature\nfrom pymatgen.analysis.local_env import NearNeighbors\nfrom keras.utils import Sequence\n\nfrom keras.models import Model\nimport pickle\nfrom sklearn.model_selection import train_test_split\n\nimport pandas as pd\n\nimport six\n# #from pybel import *\nimport openbabel as ob\nimport pybel\nfrom multiprocessing import Pool\nfrom tqdm import *\nfrom pathlib import Path\nimport sys\nimport itertools\nfrom typing import List\nfrom functools import partial\nfrom collections import deque\nfrom pymatgen import Molecule, Element\nfrom pymatgen.io.babel import BabelMolAdaptor\nfrom megnet.data.qm9 import ring_to_vector\nfrom megnet.utils.general import fast_label_binarize\n\n\nclass StructureGraph(MSONable):\n    \"\"\"\n    This is a base class for converting converting structure into graphs or model inputs\n    Methods to be implemented are follows:\n        1. convert(self, structure)\n            This is to convert a structure into a graph dictionary\n        2. get_input(self, structure)\n            This method convert a structure directly to a model input\n        3. get_flat_data(self, graphs, targets)\n            This method process graphs and targets pairs and output model input list.\n\n    \"\"\"\n\n    # TODO (wardlt): Consider making \"num_*_features\" funcs to simplify making a MEGNet model\n\n    def __init__(self,\n                 nn_strategy=None,\n                 atom_converter=None,\n                 bond_converter=None,\n                 **kwargs):\n\n        if isinstance(nn_strategy, str):\n            strategy = local_env.get(nn_strategy)\n            parameters = signature(strategy).parameters\n            param_dict = {i: j.default for i, j in parameters.items()}\n            for i, j in kwargs.items():\n                if i in param_dict:\n                    setattr(self, i, j)\n                    param_dict.update({i: j})\n            self.nn_strategy = strategy(**param_dict)\n        elif isinstance(nn_strategy, NearNeighbors):\n            self.nn_strategy = nn_strategy\n        elif nn_strategy is None:\n            self.nn_strategy = None\n        else:\n            raise RuntimeError(\"Strategy not valid\")\n\n        self.atom_converter = atom_converter\n        self.bond_converter = bond_converter\n        if self.atom_converter is None:\n            self.atom_converter = self._get_dummy_converter()\n        if self.bond_converter is None:\n            self.bond_converter = self._get_dummy_converter()\n\n    def convert(self, structure, state_attributes=None):\n        \"\"\"\n        Take a pymatgen structure and convert it to a index-type graph representation\n        The graph will have node, distance, index1, index2, where node is a vector of Z number\n        of atoms in the structure, index1 and index2 mark the atom indices forming the bond and separated by\n        distance.\n        For state attributes, you can set structure.state = [[xx, xx]] beforehand or the algorithm would\n        take default [[0, 0]]\n\n        Args:\n            state_attributes: (list) state attributes\n            structure: (pymatgen structure)\n            (dictionary)\n        \"\"\"\n        state_attributes = state_attributes or [[0, 0]]\n        index1 = []\n        index2 = []\n        bonds = []\n        if self.nn_strategy is None:\n            raise RuntimeError(\"NearNeighbor strategy is not provided!\")\n        for n, neighbors in enumerate(self.nn_strategy.get_all_nn_info(structure)):\n            index1.extend([n] * len(neighbors))\n            for neighbor in neighbors:\n                index2.append(neighbor['site_index'])\n                bonds.append(neighbor['weight'])\n\n        atoms = [i.specie.Z for i in structure]\n\n        if np.size(np.unique(index1)) < len(atoms):\n            raise RuntimeError(\"Isolated atoms found in the structure\")\n        else:\n            return {'atom': np.array(atoms, dtype='int32').tolist(),\n                    'bond': bonds,\n                    'state': state_attributes,\n                    'index1': index1,\n                    'index2': index2\n                    }\n\n    def convert(self, structure, attributes, state_attributes=None):\n        \"\"\"\n        Take a pymatgen structure and convert it to a index-type graph representation\n        The graph will have node, distance, index1, index2, where node is a vector of Z number\n        of atoms in the structure, index1 and index2 mark the atom indices forming the bond and separated by\n        distance.\n        For state attributes, you can set structure.state = [[xx, xx]] beforehand or the algorithm would\n        take default [[0, 0]]\n\n        Args:\n            state_attributes: (list) state attributes\n            structure: (pymatgen structure)\n            (dictionary)\n        \"\"\"\n        state_attributes = state_attributes or [[0, 0]]\n        index1 = []\n        index2 = []\n        bonds = []\n        if self.nn_strategy is None:\n            raise RuntimeError(\"NearNeighbor strategy is not provided!\")\n        for n, neighbors in enumerate(self.nn_strategy.get_all_nn_info(structure)):\n            index1.extend([n] * len(neighbors))\n            for neighbor in neighbors:\n                index2.append(neighbor['site_index'])\n                bonds.append(neighbor['weight'])\n\n        atoms = [i.specie.Z for i in structure]\n\n        if np.size(np.unique(index1)) < len(atoms):\n            raise RuntimeError(\"Isolated atoms found in the structure\")\n        else:\n            return {'atom': np.array(atoms, dtype='int32').tolist(),\n                    'bond': bonds,\n                    'state': state_attributes,\n                    'index1': index1,\n                    'index2': index2\n                    }\n\n    def __call__(self, structure):\n        return self.convert(structure)\n\n    def get_input(self, structure):\n        \"\"\"\n        Turns a structure into model input\n        \"\"\"\n        graph = self.convert(structure)\n        return self.graph_to_input(graph)\n\n    def graph_to_input(self, graph):\n        \"\"\"\n        Turns a graph into model input\n\n        Args:\n            (dict): Dictionary description of the graph\n        Return:\n            ([np.ndarray]): Inputs in the form needed by MEGNet\n        \"\"\"\n        gnode = [0] * len(graph['atom'])\n        gbond = [0] * len(graph['index1'])\n\n        return [expand_1st(self.atom_converter.convert(graph['atom'])),\n                expand_1st(self.bond_converter.convert(graph['bond'])),\n                expand_1st(np.array(graph['state'])),\n                expand_1st(np.array(graph['index1'])),\n                expand_1st(np.array(graph['index2'])),\n                expand_1st(np.array(gnode)),\n                expand_1st(np.array(gbond))]\n\n    def get_flat_data(self, graphs, targets=None):\n        \"\"\"\n        Expand the graph dictionary to form a list of features and targets tensors.\n        This is useful when the model is trained on assembled graphs on the fly.\n\n        Args:\n            graphs: (list of dictionary) list of graph dictionary for each structure\n            targets: (list of float or list) Optional: corresponding target\n                values for each structure\n\n        Returns:\n            tuple(node_features, edges_features, global_values, index1, index2, targets)\n        \"\"\"\n\n        output = []  # Will be a list of arrays\n\n        # Convert the graphs to matrices\n        for feature in ['atom', 'bond', 'state', 'index1', 'index2']:\n            output.append([np.array(x[feature]) for x in graphs])\n\n        # If needed, add the targets\n        if targets is not None:\n            output.append([to_list(t) for t in targets])\n\n        return tuple(output)\n\n    def get_flat_scarlar_data(self, graphs, targets=None):\n        \"\"\"\n        Expand the graph dictionary to form a list of features and targets tensors.\n        This is useful when the model is trained on assembled graphs on the fly.\n\n        Args:\n            graphs: (list of dictionary) list of graph dictionary for each structure\n            targets: (list of float or list) Optional: corresponding target\n                values for each structure\n\n        Returns:\n            tuple(node_features, edges_features, global_values, index1, index2, targets)\n        \"\"\"\n\n        output = []  # Will be a list of arrays\n\n        # Convert the graphs to matrices\n        for feature in ['atom', 'bond', 'state', 'index1', 'index2']:\n            output.append([np.array(x[feature]) for x in graphs])\n\n        #  gnode = [0] * len(graph['atom'])\n        # gbond = [0] * len(graph['index1'])\n        # elf.graph_converter.graph_to_input(graph)\n        # # If needed, add the targets\n        # count=128\n        # if targets is not None:\n        #     for bond in output[1]:\n        #         c = len(bond)\n        #         if(c>count):\n        #             count = c\n        # print(\"bond max count:\"+str(c))\n        # current_words = list(current_words + [0] * (10 - len(current_words)))\n        # output.append([list(to_list(t['scalar_coupling_constant'])+[0]*(count-len(to_list(t['scalar_coupling_constant'])))) for t in targets])\n\n        for feature in ['output']:\n            ff = np.array([np.array(x[feature]) for x in targets])\n            print('mmmmm')\n            print(ff.shape)\n\n            # scal_coupling_constant = ff.reshape(-1, 1)\n            # print(scal_coupling_constant.shape)\n\n            output.append(ff)\n            print(ff)\n\n        return tuple(output)\n\n    @staticmethod\n    def _get_dummy_converter():\n        return DummyConverter()\n\n    def as_dict(self):\n        all_dict = super().as_dict()\n        if 'nn_strategy' in all_dict:\n            nn_strategy = all_dict.pop('nn_strategy')\n            all_dict.update({'nn_strategy': local_env.serialize(nn_strategy)})\n        return all_dict\n\n    @classmethod\n    def from_dict(cls, d):\n        if 'nn_strategy' in d:\n            nn_strategy = d.pop('nn_strategy')\n            nn_strategy_obj = local_env.deserialize(nn_strategy)\n            d.update({'nn_strategy': nn_strategy_obj})\n            return super().from_dict(d)\n        return super().from_dict(d)\n\n\nclass DistanceConverter(MSONable):\n    \"\"\"\n    Base class for distance conversion. The class needs to have a convert method.\n    \"\"\"\n\n    def convert(self, d):\n        raise NotImplementedError\n\n\nclass DummyConverter(DistanceConverter):\n    \"\"\"\n    Dummy converter as a placeholder\n    \"\"\"\n\n    def convert(self, d):\n        return d\n\n\nclass GaussianDistance(DistanceConverter):\n    \"\"\"\n    Expand distance with Gaussian basis sit at centers and with width 0.5.\n\n    Args:\n        centers: (np.array)\n        width: (float)\n    \"\"\"\n\n    def __init__(self, centers=np.linspace(0, 5, 100), width=0.5):\n        self.centers = centers\n        self.width = width\n\n    def convert(self, d):\n        \"\"\"\n        expand distance vector d with given parameters\n\n        Args:\n            d: (1d array) distance array\n\n        Returns\n            (matrix) N*M matrix with N the length of d and M the length of centers\n        \"\"\"\n        d = np.array(d)\n        return np.exp(-(d[:, None] - self.centers[None, :]) ** 2 / self.width ** 2)\n\n\nclass MoorseLongRange(DistanceConverter):\n    \"\"\"\n    This is an attempt to implement a Moorse/long range interactomic potential like\n    distance expansion. The distance is expanded with this basis at different equilibrium\n    bond distance, r_eq. It is still a work in progress. Do not use if you do not know\n    much about the parameters\n    ref: https://en.wikipedia.org/wiki/Morse/Long-range_potential#Function\n\n    Args:\n        d_e: (float) dissociate energy\n        r_ref: (float) reference bond length\n        r_eq: (list) equilibrium bond length\n        p: (int) exponential term in the original equation, see ref\n        q: (int) exponential term in the original equaiton, see ref\n        cm: (list) long range coefficients in u_LR = \\Sigma_i_N (cm_i / r^i)\n        betas: (list) parameters determining the transition between long range and short range\n    \"\"\"\n\n    def __init__(self, d_e=1, r_ref=2, r_eq=[1, 2, 3],\n                 p=2, q=2, cm=[1, 2, 3, 4],\n                 betas=[0.1, 0.2, 0.3, 0.4]):\n        self.d_e = d_e\n        self.r_ref = r_ref\n        self.r_eq = np.array(r_eq)\n        self.p = p\n        self.q = q\n        self.cm = np.array(cm).ravel()\n        self.n_cm = len(self.cm)\n        self.betas = np.array(betas).ravel()\n\n    def convert(self, d):\n        return self.d_e * (1 - self.u(d)[:, None] / self.u(self.r_eq)[None, :] *\n                           np.exp(-self.beta(d) * self.y(d[:, None], self.r_eq[None, :], self.p))) ** 2\n\n    def u(self, r):\n        m_i = np.arange(1, self.n_cm + 1)\n        if np.size(r) == 1:\n            return np.sum(self.cm / r ** m_i)\n        return np.sum(self.cm[None, :] / r[:, None] ** m_i[None, :], axis=1).ravel()\n\n    @staticmethod\n    def y(r, r_ref, p):\n        return (r ** p - r_ref ** p) / (r ** p + r_ref ** p)\n\n    def beta(self, r):\n        y_p_ref = self.y(r, self.r_ref, self.p)\n        y_q_ref = self.y(r, self.r_ref, self.q)\n        return self.beta_inf[None, :] * y_p_ref[:, None] + (1 - y_p_ref[:, None]) * \\\n               np.sum(self.betas[None, :] * y_q_ref[:, None] ** np.arange(0, len(self.betas))[None, :], axis=1).ravel()[\n               :, None]\n\n    @property\n    def beta_inf(self):\n        return np.log(2 * self.d_e / self.u(self.r_eq))\n\n\nclass BaseGraphBatchGenerator(Sequence):\n    \"\"\"Base class for classes that generate batches of training data for MEGNet.\n    Based on the Sequence class, which is the data loader equivalent for Keras.\n\n    Implementations of this base class must implement the :meth:`_generate_inputs`,\n    which generates the lists of graph descriptions for a batch.\n\n    The :meth:`process_atom_features` function and related functions are used to modify\n    the features for each atom, bond, and global features when creating a batch.\n    \"\"\"\n\n    def __init__(self, dataset_size, targets, batch_size=128, shuffle=True):\n        \"\"\"\n        Args:\n            dataset_size (int): Number of entries in dataset\n            targets (ndarray): Feature to be predicted for each network\n            batch_size (int): Maximum batch size\n            shuffle (bool): Whether to shuffle the data after each step\n        \"\"\"\n        if targets is not None:\n            self.targets = np.array(targets)\n        else:\n            self.targets = None\n        self.batch_size = batch_size\n        self.total_n = dataset_size\n        self.is_shuffle = shuffle\n        self.max_step = int(np.ceil(self.total_n / batch_size))\n        self.mol_index = np.arange(self.total_n)\n        if self.is_shuffle:\n            self.mol_index = np.random.permutation(self.mol_index)\n\n    def __len__(self):\n        return self.max_step\n\n    def _combine_graph_data(self, feature_list_temp, connection_list_temp, global_list_temp,\n                            index1_temp, index2_temp):\n        \"\"\"Compile the matrices describing each graph into single matrices for the entire graph\n\n        Beyond concatenating the graph descriptions, this operation updates the indices of each\n        node to be sequential across all graphs so they are not duplicated between graphs\n\n        Args:\n            feature_list_temp ([ndarray]): List of features for each node\n            connection_list_temp ([ndarray]): List of features for each connection\n            global_list_temp ([ndarray]): List of global state for each graph\n            index1_temp ([ndarray]): List of indices for the start of each bond\n            index2_temp ([ndarray]): List of indices for the end of each bond\n        Returns:\n            (tuple): Input arrays describing the entire batch of networks:\n                - ndarray: Features for each node\n                - ndarray: Features for each connection\n                - ndarray: Global state for each graph\n                - ndarray: Indices for the start of each bond\n                - ndarray: Indices for the end of each bond\n                - ndarray: Index of graph associated with each node\n                - ndarray: Index of graph associated with each connection\n        \"\"\"\n        # get atom's structure id\n        gnode = []\n        for i, j in enumerate(feature_list_temp):\n            gnode += [i] * len(j)\n        # get bond features from a batch of structures\n        # get bond's structure id\n        gbond = []\n        for i, j in enumerate(connection_list_temp):\n            gbond += [i] * len(j)\n\n        # assemble atom features together\n        feature_list_temp = np.concatenate(feature_list_temp, axis=0)\n        feature_list_temp = self.process_atom_feature(feature_list_temp)\n\n        # assemble bond feature together\n        connection_list_temp = np.concatenate(connection_list_temp, axis=0)\n        connection_list_temp = self.process_bond_feature(connection_list_temp)\n\n        # assemble state feature together\n        global_list_temp = np.concatenate(global_list_temp, axis=0)\n        global_list_temp = self.process_state_feature(global_list_temp)\n\n        # assemble bond indices\n        index1 = []\n        index2 = []\n        offset_ind = 0\n        for ind1, ind2 in zip(index1_temp, index2_temp):\n            index1 += [i + offset_ind for i in ind1]\n            index2 += [i + offset_ind for i in ind2]\n            offset_ind += (max(ind1) + 1)\n\n        # Compile the inputs in needed order\n        inputs = [expand_1st(feature_list_temp),\n                  expand_1st(connection_list_temp),\n                  expand_1st(global_list_temp),\n                  expand_1st(index1),\n                  expand_1st(index2),\n                  expand_1st(gnode),\n                  expand_1st(gbond)]\n        return inputs\n\n    def on_epoch_end(self):\n        if self.is_shuffle:\n            self.mol_index = np.random.permutation(self.mol_index)\n\n    def process_atom_feature(self, x):\n        return x\n\n    def process_bond_feature(self, x):\n        return x\n\n    def process_state_feature(self, x):\n        return x\n\n    def __getitem__(self, index):\n        # Get the indices for this batch\n        batch_index = self.mol_index[index * self.batch_size:(index + 1) * self.batch_size]\n\n        # Get the inputs for each batch\n        inputs = self._generate_inputs(batch_index)\n\n        # Make the graph data\n        inputs = self._combine_graph_data(*inputs)\n\n        # Return the batch\n        if self.targets is None:\n            return inputs\n        else:\n            # get targets\n            it = itemgetter(*batch_index)\n            target_temp = itemgetter_list(self.targets, batch_index)\n            target_temp = np.atleast_2d(target_temp)\n\n            return inputs, expand_1st(target_temp)\n\n    @abstractmethod\n    def _generate_inputs(self, batch_index):\n        \"\"\"Get the graph descriptions for each batch\n\n        Args:\n             batch_index ([int]): List of indices for training batch\n        Returns:\n            (tuple): Input arrays describing each network:\n                - [ndarray]: List of features for each node\n                - [ndarray]: List of features for each connection\n                - [ndarray]: List of global state for each graph\n                - [ndarray]: List of indices for the start of each bond\n                - [ndarray]: List of indices for the end of each bond\n        \"\"\"\n        pass\n\n\nclass GraphBatchGenerator(BaseGraphBatchGenerator):\n    \"\"\"\n    A generator class that assembles several structures (indicated by\n    batch_size) and form (x, y) pairs for model training.\n\n    Args:\n        atom_features: (list of np.array) list of atom feature matrix,\n        bond_features: (list of np.array) list of bond features matrix\n        state_features: (list of np.array) list of [1, G] state features,\n            where G is the global state feature dimension\n        index1_list: (list of integer) list of (M, ) one side atomic index of the bond,\n        M is different for different structures\n        index2_list: (list of integer) list of (M, ) the other side atomic\n            index of the bond, M is different for different structures,\n            but it has to be the same as the corresponding index1.\n        targets: (numpy array), N*1, where N is the number of structures\n        batch_size: (int) number of samples in a batch\n    \"\"\"\n\n    def __init__(self,\n                 atom_features,\n                 bond_features,\n                 state_features,\n                 index1_list,\n                 index2_list,\n                 targets=None,\n                 batch_size=128,\n                 is_shuffle=True):\n        super().__init__(len(atom_features), targets, batch_size, is_shuffle)\n        self.atom_features = atom_features\n        self.bond_features = bond_features\n        self.state_features = state_features\n        self.index1_list = index1_list\n        self.index2_list = index2_list\n\n    def _generate_inputs(self, batch_index):\n        \"\"\"Get the graph descriptions for each batch\n\n        Args:\n             batch_index ([int]): List of indices for training batch\n        Returns:\n            (tuple): Input arrays describe each network:\n                - [ndarray]: List of features for each nodes\n                - [ndarray]: List of features for each connection\n                - [ndarray]: List of global state for each graph\n                - [ndarray]: List of indices for the start of each bond\n                - [ndarray]: List of indices for the end of each bond\n        \"\"\"\n\n        # Get the features and connectivity lists for this batch\n        it = itemgetter(*batch_index)\n        feature_list_temp = itemgetter_list(self.atom_features, batch_index)\n        connection_list_temp = itemgetter_list(self.bond_features, batch_index)\n        global_list_temp = itemgetter_list(self.state_features, batch_index)\n        index1_temp = itemgetter_list(self.index1_list, batch_index)\n        index2_temp = itemgetter_list(self.index2_list, batch_index)\n\n        return feature_list_temp, connection_list_temp, global_list_temp, \\\n               index1_temp, index2_temp\n\n\nclass GraphBatchDistanceConvert(GraphBatchGenerator):\n    \"\"\"\n    Generate batch of structures with bond distance being expanded using a Expansor\n\n    Args:\n        atom_features: (list of np.array) list of atom feature matrix,\n        bond_features: (list of np.array) list of bond features matrix\n        state_features: (list of np.array) list of [1, G] state features, where G is the global state feature dimension\n        index1_list: (list of integer) list of (M, ) one side atomic index of the bond, M is different for differentstructures\n        index2_list: (list of integer) list of (M, ) the other side atomic index of the bond, M is different for different\n            structures, but it has to be the same as the correponding index1.\n        targets: (numpy array), N*1, where N is the number of structures\n        batch_size: (int) number of samples in a batch\n        is_shuffle: (bool) whether to shuffle the structure, default to True\n        distance_converter: (bool) converter for processing the distances\n\n    \"\"\"\n\n    def __init__(self,\n                 atom_features,\n                 bond_features,\n                 state_features,\n                 index1_list,\n                 index2_list,\n                 targets=None,\n                 batch_size=128,\n                 is_shuffle=True,\n                 distance_converter=None):\n        super().__init__(atom_features=atom_features,\n                         bond_features=bond_features,\n                         state_features=state_features,\n                         index1_list=index1_list,\n                         index2_list=index2_list,\n                         targets=targets,\n                         batch_size=batch_size,\n                         is_shuffle=is_shuffle)\n        self.distance_converter = distance_converter\n\n    def process_bond_feature(self, x):\n        return self.distance_converter.convert(x)\n\n\ndef itemgetter_list(l, indices):\n    \"\"\"\n    Get indices of l and return a tuple\n\n    Args:\n        l:  (list)\n        indices: (list) indices\n\n    Returns:\n        (tuple)\n    \"\"\"\n    it = itemgetter(*indices)\n    if np.size(indices) == 1:\n        return it(l),\n    else:\n        return it(l)\n\n\n\"\"\"\nTools for creating graph inputs from molecule data\n\"\"\"\n\n# from megnet.data.graph import (StructureGraph, GaussianDistance,\n#                                BaseGraphBatchGenerator, GraphBatchGenerator)\n#\ntry:\n    import pybel\nexcept ImportError:\n    pybel = None\n\ntry:\n    from rdkit import Chem\nexcept ImportError:\n    Chem = None\n\n__date__ = '12/01/2018'\n\n# List of features to use by default for each atom\n_ATOM_FEATURES = ['element', 'chirality', 'formal_charge', 'ring_sizes',\n                  'hybridization', 'donor', 'acceptor', 'aromatic']\n\n# List of features to use by default for each bond\n_BOND_FEATURES = ['bond_type', 'same_ring', 'spatial_distance', 'graph_distance', 'coupling_type']\n# _BOND_FEATURES = ['bond_type', 'same_ring', 'spatial_distance', 'graph_distance']\n\n# List of elements in library to use by default\n_ELEMENTS = ['H', 'C', 'N', 'O', 'F']\n\n\nclass SimpleMolGraph(StructureGraph):\n    \"\"\"\n    Default using all atom pairs as bonds. The distance between atoms are used\n    as bond features. By default the distance is expanded using a Gaussian\n    expansion with centers at np.linspace(0, 4, 20) and width of 0.5\n    \"\"\"\n\n    def __init__(self,\n                 nn_strategy='AllAtomPairs',\n                 atom_converter=None,\n                 bond_converter=None\n                 ):\n        if bond_converter is None:\n            bond_converter = GaussianDistance(np.linspace(0, 4, 20), 0.5)\n        super().__init__(nn_strategy=nn_strategy, atom_converter=atom_converter,\n                         bond_converter=bond_converter)\n\n\nclass MolecularGraph(StructureGraph):\n    \"\"\"Class for generating the graph inputs from a molecule\n\n    Computes many different features for the atoms and bonds in a molecule, and prepares them\n    in a form compatible with MEGNet models. The :meth:`convert` method takes a OpenBabel molecule\n    and, besides computing features, also encodes them in a form compatible with machine learning.\n    Namely, the `convert` method one-hot encodes categorical variables and concatenates\n    the atomic features\n\n    ## Atomic Features\n\n    This class can compute the following features for each atom\n\n    - `atomic_num`: The atomic number\n    - `element`: (categorical) Element identity. (Unlike `atomic_num`, element is one-hot-encoded)\n    - `chirality`: (categorical) R, S, or not a Chiral center (one-hot encoded).\n    - `formal_charge`: Formal charge of the atom\n    - `ring_sizes`: For rings with 9 or fewer atoms, how many unique rings\n    of each size include this atom\n    - `hybridization`: (categorical) Hybridization of atom: sp, sp2, sp3, sq.\n    planer, trig, octahedral, or hydrogen\n    - `donor`: (boolean) Whether the atom is a hydrogen bond donor\n    - `acceptor`: (boolean) Whether the atom is a hydrogen bond acceptor\n    - `aromatic`: (boolean) Whether the atom is part of an aromatic system\n\n    ## Atom Pair Features\n\n    The class also computes features for each pair of atoms\n\n    - `bond_type`: (categorical) Whether the pair are unbonded, or in a single, double, triple, or aromatic bond\n    - `same_ring`: (boolean) Whether the atoms are in the same aromatic ring\n    - `graph_distance`: Distance of shortest path between atoms on the bonding graph\n    - `spatial_distance`: Euclidean distance between the atoms. By default, this distance is expanded into\n        a vector of 20 different values computed using the `GaussianDistance` converter\n\n    \"\"\"\n\n    def __init__(self, atom_features=None, bond_features=None, distance_converter=None,\n                 known_elements=None):\n        \"\"\"\n        Args:\n            atom_features ([str]): List of atom features to compute\n            bond_features ([str]): List of bond features to compute\n            distance_converter (DistanceCovertor): Tool used to expand distances\n                from a single scalar vector to an array of values\n            known_elements ([str]): List of elements expected to be in dataset. Used only if the\n                feature `element` is used to describe each atom\n        \"\"\"\n\n        # Check if openbabel and RDKit are installed\n        if Chem is None or pybel is None:\n            raise RuntimeError('RDKit and openbabel must be installed')\n\n        super().__init__()\n        if bond_features is None:\n            bond_features = _BOND_FEATURES\n        if atom_features is None:\n            atom_features = _ATOM_FEATURES\n        if distance_converter is None:\n            distance_converter = GaussianDistance(np.linspace(0, 4, 20), 0.5)\n        if known_elements is None:\n            known_elements = _ELEMENTS\n\n        # Check if all feature names are valid\n        if any(i not in _ATOM_FEATURES for i in atom_features):\n            bad_features = set(atom_features).difference(_ATOM_FEATURES)\n            raise ValueError('Unrecognized atom features: {}'.format(', '.join(bad_features)))\n        self.atom_features = atom_features\n        if any(i not in _BOND_FEATURES for i in bond_features):\n            bad_features = set(bond_features).difference(_BOND_FEATURES)\n            raise ValueError('Unrecognized bond features: {}'.format(', '.join(bad_features)))\n        self.bond_features = bond_features\n        self.known_elements = known_elements\n        self.distance_converter = distance_converter\n\n    def convert(self, mol, state_attributes=None, full_pair_matrix=True):\n        \"\"\"\n        Compute the representation for a molecule\n\n        Args：\n            mol (pybel.Molecule): Molecule to generate features for\n            state_attributes (list): State attributes. Uses average mass and number of bonds per atom as default\n            full_pair_matrix (bool): Whether to generate info for all atom pairs, not just bonded ones\n        Returns:\n            (dict): Dictionary of features\n        \"\"\"\n\n        # Get the features features for all atoms and bonds\n        atom_features = []\n        atom_pairs = []\n        for idx, atom in enumerate(mol.atoms):\n            f = self.get_atom_feature(mol, atom)\n            atom_features.append(f)\n        atom_features = sorted(atom_features, key=lambda x: x[\"coordid\"])\n        num_atoms = mol.OBMol.NumAtoms()\n        for i, j in itertools.combinations(range(0, num_atoms), 2):\n            bond_feature = self.get_pair_feature(mol, i, j, full_pair_matrix)\n            if bond_feature:\n                atom_pairs.append(bond_feature)\n            else:\n                continue\n\n        # Compute the graph distance, if desired\n        if 'graph_distance' in self.bond_features:\n            graph_dist = self._dijkstra_distance(atom_pairs)\n            for i in atom_pairs:\n                i.update({'graph_distance': graph_dist[i['a_idx'], i['b_idx']]})\n\n        # Generate the state attributes (that describe the whole network)\n        state_attributes = state_attributes or [\n            [mol.molwt / num_atoms,\n             len([i for i in atom_pairs if i['bond_type'] > 0]) / num_atoms]\n        ]\n\n        # Get the atom features in the order they are requested by the user as a 2D array\n        atoms = []\n        for atom in atom_features:\n            atoms.append(self._create_atom_feature_vector(atom))\n\n        # Get the bond features in the order request by the user\n        bonds = []\n        index1_temp = []\n        index2_temp = []\n        for bond in atom_pairs:\n            # Store the index of each bond\n            index1_temp.append(bond.pop('a_idx'))\n            index2_temp.append(bond.pop('b_idx'))\n\n            # Get the desired bond features\n            bonds.append(self._create_pair_feature_vector(bond))\n\n        # Given the bonds (i,j), make it so (i,j) == (j, i)\n        index1 = index1_temp + index2_temp\n        index2 = index2_temp + index1_temp\n        bonds = bonds + bonds\n\n        # Sort the arrays by the beginning index\n        sorted_arg = np.argsort(index1)\n        index1 = np.array(index1)[sorted_arg].tolist()\n        index2 = np.array(index2)[sorted_arg].tolist()\n        bonds = np.array(bonds)[sorted_arg].tolist()\n\n        return {'atom': atoms,\n                'bond': bonds,\n                'state': state_attributes,\n                'index1': index1,\n                'index2': index2}\n\n    def convert_coupling1(self, mol, coupling_attributes, energy_dipole, charges_magnetic, target_convert=True,\n                          state_attributes=None, full_pair_matrix=True):\n        \"\"\"\n        Compute the representation for a molecule\n\n        Args：\n            mol (pybel.Molecule): Molecule to generate features for\n            state_attributes (list): State attributes. Uses average mass and number of bonds per atom as default\n            full_pair_matrix (bool): Whether to generate info for all atom pairs, not just bonded ones\n        Returns:\n            (dict): Dictionary of features\n        \"\"\"\n\n        # Get the features features for all atoms and bonds\n        atom_features = []\n        atom_pairs = []\n        for idx, atom in enumerate(mol.atoms):\n            get_data = charges_magnetic[charges_magnetic['atom_index'] == idx]\n            f = self.get_atom_feature(mol, get_data, atom)\n            atom_features.append(f)\n        atom_features = sorted(atom_features, key=lambda x: x[\"coordid\"])\n        num_atoms = mol.OBMol.NumAtoms()\n        for i, j in itertools.combinations(range(0, num_atoms), 2):\n            bond_feature = self.get_pair_feature_coupling(mol, coupling_attributes, i, j, full_pair_matrix,\n                                                          target_convert)\n            if bond_feature:\n                atom_pairs.append(bond_feature)\n            else:\n                continue\n\n        # Compute the graph distance, if desired\n        if 'graph_distance' in self.bond_features:\n            graph_dist = self._dijkstra_distance(atom_pairs)\n            for i in atom_pairs:\n                i.update({'graph_distance': graph_dist[i['a_idx'], i['b_idx']]})\n\n        # Generate the state attributes (that describe the whole network)\n        state_attributes = state_attributes or [\n            [mol.molwt / num_atoms,\n             len([i for i in atom_pairs if i['bond_type'] > 0]) / num_atoms, energy_dipole['dipole_moments'][0],\n             energy_dipole['potential_energy'][0]]\n        ]\n\n        # Get the atom features in the order they are requested by the user as a 2D array\n        atoms = []\n        for atom in atom_features:\n            atoms.append(self._create_atom_feature_vector(atom))\n\n        # Get the bond features in the order request by the user\n        bonds = []\n        index1_temp = []\n        index2_temp = []\n        scal_coupling_constant = []\n        bonds_pair_count = []\n        for bond in atom_pairs:\n            # Store the index of each bond\n            index1_temp.append(bond.pop('a_idx'))\n            index2_temp.append(bond.pop('b_idx'))\n            # train out data\n            if bond.get('bond_index_pair'):\n                bonds_pair_count.append(bond.pop('bond_index_pair'))\n\n            scal_coupling_constant.append(bond.pop('scalar_coupling_constant'))\n            # Get the desired bond features\n            bonds.append(self._create_pair_feature_vector(bond))\n\n        # Given the bonds (i,j), make it so (i,j) == (j, i)\n        index1 = index1_temp + index2_temp\n        index2 = index2_temp + index1_temp\n        bonds = bonds + bonds\n        # bonds_train_out = bonds_train_out+bonds_train_out\n        scal_coupling_constant = scal_coupling_constant + scal_coupling_constant\n\n        # Sort the arrays by the beginning index\n        sorted_arg = np.argsort(index1)\n        index1 = np.array(index1)[sorted_arg].tolist()\n        index2 = np.array(index2)[sorted_arg].tolist()\n\n        bonds = np.array(bonds)[sorted_arg].tolist()\n        # bonds_train_out = np.array(bonds_train_out)[sorted_arg].tolist()\n        scal_coupling_constant = np.array(scal_coupling_constant)[sorted_arg].tolist()\n\n        scal_coupling_constant = list(scal_coupling_constant + [0] * (756 - len(scal_coupling_constant)))\n        # print(\"vvvvv\")\n        # print(scal_coupling_constant.shape)\n        # print(scal_coupling_constant)\n        # scal_coupling_constant = scal_coupling_constant.reshape(-1,1)\n        # print(scal_coupling_constant.shape)\n\n        gnode = [0] * len(atoms)\n        gbond = [0] * len(index1)\n        if target_convert == False:\n            return ({'atom': atoms,\n                     'bond': bonds,\n                     'state': state_attributes,\n                     'index1': index1,\n                     'index2': index2\n\n                     }, {'index1': index1,\n                         'index2': index2,\n                         'coupling_attributes': coupling_attributes\n                         })\n\n        return ({'atom': atoms,\n                 'bond': bonds,\n                 'state': state_attributes,\n                 'index1': index1,\n                 'index2': index2,\n\n                 }, {'scal_coupling_constant': scal_coupling_constant})\n\n    def convert_coupling(self, mol, coupling_attributes, energy_dipole, charges_magnetic, target_convert=True,\n                         state_attributes=None, full_pair_matrix=True):\n        \"\"\"\n        Compute the representation for a molecule\n\n        Args：\n            mol (pybel.Molecule): Molecule to generate features for\n            state_attributes (list): State attributes. Uses average mass and number of bonds per atom as default\n            full_pair_matrix (bool): Whether to generate info for all atom pairs, not just bonded ones\n        Returns:\n            (dict): Dictionary of features\n        \"\"\"\n\n        # Get the features features for all atoms and bonds\n        atom_features = []\n        atom_pairs = []\n        for idx, atom in enumerate(mol.atoms):\n            get_data = charges_magnetic[charges_magnetic['atom_index'] == idx]\n            f = self.get_atom_feature(mol, get_data, atom)\n            atom_features.append(f)\n        atom_features = sorted(atom_features, key=lambda x: x[\"coordid\"])\n        num_atoms = mol.OBMol.NumAtoms()\n        # get the bonds\n        # for index,row in coupling_attributes.iterrows():\n        #     atom_index_0 = row[\"atom_index_0\"]\n        #     atom_index_1 = row['atom_index_1']\n        #     type = row['type']\n        #     if target_convert:\n        #         scalar_coupling_constant = row['scalar_coupling_constant']\n        #     else:\n        #         scalar_coupling_constant=0\n        for i, j in itertools.combinations(range(0, num_atoms), 2):\n            bond_feature = self.get_pair_feature_coupling(mol, coupling_attributes, i, j, full_pair_matrix,\n                                                          target_convert)\n            if bond_feature:\n                atom_pairs.append(bond_feature)\n            else:\n                continue\n\n        # Compute the graph distance, if desired\n        if 'graph_distance' in self.bond_features:\n            graph_dist = self._dijkstra_distance(atom_pairs)\n            for i in atom_pairs:\n                i.update({'graph_distance': graph_dist[i['a_idx'], i['b_idx']]})\n\n        # Generate the state attributes (that describe the whole network)\n        state_attributes = state_attributes or [\n            [mol.molwt / num_atoms,\n             len([i for i in atom_pairs if i['bond_type'] > 0]) / num_atoms\n           ]\n        ]\n       ##mol_out\n        mol_out=[]\n\n        mol_out.append( energy_dipole['dipole_moments'].tolist()[0])\n        mol_out.append(energy_dipole['potential_energy'].tolist()[0])\n\n        # Get the atom features in the order they are requested by the user as a 2D array\n        atoms = []\n        atmos_out=[]\n\n        for atom in atom_features:\n\n            atmos_out.append(atom.pop('sigma_iso'))\n            atmos_out.append(atom.pop('omega'))\n            atmos_out.append(atom.pop('kappa'))\n            atmos_out.append(atom.pop('mulliken_charge'))\n\n            atoms.append(self._create_atom_feature_vector(atom))\n\n        # Get the bond features in the order request by the user\n        bonds = []\n        index1_temp = []\n        index2_temp = []\n        scal_coupling_constant = []\n        bonds_pair_count = []\n        for bond in atom_pairs:\n            # Store the index of each bond\n            index1_temp.append(bond.pop('a_idx'))\n            index2_temp.append(bond.pop('b_idx'))\n            # train out data\n\n            scal_coupling_constant.append(bond.pop('scalar_coupling_constant'))\n            # Get the desired bond features\n            bonds.append(self._create_pair_feature_vector(bond))\n\n        # Given the bonds (i,j), make it so (i,j) == (j, i)\n        index1 = index1_temp + index2_temp\n        index2 = index2_temp + index1_temp\n        bonds = bonds + bonds\n        # bonds_train_out = bonds_train_out+bonds_train_out\n        scal_coupling_constant = scal_coupling_constant + scal_coupling_constant\n\n        # Sort the arrays by the beginning index\n        sorted_arg = np.argsort(index1)\n        index1 = np.array(index1)[sorted_arg].tolist()\n        index2 = np.array(index2)[sorted_arg].tolist()\n\n        bonds = np.array(bonds)[sorted_arg].tolist()\n        # bonds_train_out = np.array(bonds_train_out)[sorted_arg].tolist()\n        scal_coupling_constant = np.array(scal_coupling_constant)[sorted_arg].tolist()\n\n        scal_coupling_constant = list(scal_coupling_constant + [0] * (820 - len(scal_coupling_constant)))\n\n        # atmos_out = list(atmos_out + [0] * (112 - len(atmos_out)))\n\n        # output = np.concatenate([mol_out,atmos_out,scal_coupling_constant],axis=0)\n\n        # print(scal_coupling_constant.shape)\n        # print(scal_coupling_constant)\n        # scal_coupling_constant = scal_coupling_constant.reshape(-1,1)\n        # print(scal_coupling_constant.shape)\n\n        if target_convert == False:\n            return ({'atom': atoms,\n                     'bond': bonds,\n                     'state': state_attributes,\n                     'index1': index1,\n                     'index2': index2\n\n                     }, {'index1': index1,\n                         'index2': index2,\n                         'coupling_attributes': coupling_attributes\n                         })\n\n        return ({'atom': atoms,\n                 'bond': bonds,\n                 'state': state_attributes,\n                 'index1': index1,\n                 'index2': index2\n\n                 }, {'scal_coupling_constant': scal_coupling_constant})\n\n    def _create_pair_feature_vector(self, bond: dict) -> List[float]:\n        \"\"\"Generate the feature vector from the bond feature dictionary\n\n        Handles the binarization of categorical variables, and performing the distance conversion\n\n        Args:\n            bond (dict): Features for a certain pair of atoms\n        Returns:\n            ([float]) Values converted to a vector\n            \"\"\"\n        bond_temp = []\n        for i in self.bond_features:\n            # Some features require conversion (e.g., binarization)\n            if i in bond:\n                if i == \"bond_type\":\n                    bond_temp.extend(fast_label_binarize(bond[i], [0, 1, 2, 3, 4]))\n                elif i == \"same_ring\":\n                    bond_temp.append(int(bond[i]))\n                elif i == \"coupling_type\":\n                    bond_temp.extend(fast_label_binarize(bond[i], [0, 1, 2, 3, 4, 5, 6, 7]))\n                elif i == \"scalar_coupling_constant\":\n                    bond_temp.append(float(bond[i]))\n\n                elif i == \"spatial_distance\":\n                    expanded = self.distance_converter.convert([bond[i]])[0]\n                    if isinstance(expanded, np.ndarray):\n                        # If we use a distance expansion\n                        bond_temp.extend(expanded.tolist())\n                    else:\n                        # If not\n                        bond_temp.append(expanded)\n                else:\n                    bond_temp.append(bond[i])\n        return bond_temp\n\n    def _create_atom_feature_vector(self, atom: dict) -> List[int]:\n        \"\"\"Generate the feature vector from the atomic feature dictionary\n\n        Handles the binarization of categorical variables, and transforming the ring_sizes to a list\n\n        Args:\n            atom (dict): Dictionary of atomic features\n        Returns:\n            ([int]): Atomic feature vector\n        \"\"\"\n        atom_temp = []\n        for i in self.atom_features:\n            if i == 'chirality':\n                atom_temp.extend(fast_label_binarize(atom[i], [0, 1, 2]))\n            elif i == 'element':\n                atom_temp.extend(fast_label_binarize(atom[i], self.known_elements))\n            elif i in ['aromatic', 'donor', 'acceptor']:\n                atom_temp.append(int(atom[i]))\n            elif i == 'hybridization':\n                atom_temp.extend(fast_label_binarize(atom[i], [1, 2, 3, 4, 5, 6]))\n            elif i == 'ring_sizes':\n                atom_temp.extend(ring_to_vector(atom[i]))\n            else:  # It is a scalar\n                atom_temp.append(atom[i])\n        return atom_temp\n\n    def _dijkstra_distance(self, pairs):\n        \"\"\"\n        Compute the graph distance between each pair of atoms,\n        using the network defined by the bonded atoms.\n\n        Args:\n            pairs ([dict]): List of bond information\n        Returns:\n            ([int]) Distance for each pair of bonds\n        \"\"\"\n        bonds = []\n        for p in pairs:\n            if p['bond_type'] > 0:\n                bonds.append([p['a_idx'], p['b_idx']])\n\n        return dijkstra_distance(bonds)\n\n    def get_atom_feature(self, mol, charges_magnetic, atom):\n        \"\"\"\n        Generate all features of a particular atom\n\n        Args:\n            mol (pybel.Molecule): Molecule being evaluated\n            atom (pybel.Atom): Specific atom being evaluated\n        Return:\n            (dict): All features for that atom\n        \"\"\"\n\n        # Get the link to the OpenBabel representation of the atom\n        obatom = atom.OBAtom\n        atom_idx = atom.idx - 1  # (pybel atoms indices start from 1)\n\n        # Get the element\n        element = Element.from_Z(obatom.GetAtomicNum()).symbol\n\n        # Get the fast-to-compute properties\n        output = {\"element\": element,\n                  \"atomic_num\": obatom.GetAtomicNum(),\n                  \"formal_charge\": obatom.GetFormalCharge(),\n                  \"hybridization\": 6 if element == 'H' else obatom.GetHyb(),\n                  \"acceptor\": obatom.IsHbondAcceptor(),\n                  \"donor\": obatom.IsHbondDonorH() if atom.type == 'H' else obatom.IsHbondDonor(),\n                  \"aromatic\": obatom.IsAromatic(),\n                  \"coordid\": atom.coordidx,\n                  \"sigma_iso\": charges_magnetic['sigma_iso'].tolist()[0],\n                  \"omega\": charges_magnetic['omega'].tolist()[0],\n                  \"kappa\": charges_magnetic['kappa'].tolist()[0],\n                  \"mulliken_charge\": charges_magnetic['mulliken_charge'].tolist()[0]\n\n                  }\n\n        # Get the chirality, if desired\n        if 'chirality' in self.atom_features:\n            # Determine whether the molecule has chiral centers\n            chiral_cc = self._get_chiral_centers(mol)\n            if atom_idx not in chiral_cc:\n                output['chirality'] = 0\n            else:\n                # 1 --> 'R', 2 --> 'S'\n                output['chirality'] = 1 if chiral_cc[atom_idx] == 'R' else 2\n\n        # Find the rings, if desired\n        if 'ring_sizes' in self.atom_features:\n            rings = mol.OBMol.GetSSSR()  # OpenBabel caches ring computation internally, no need to cache ourselves\n            output['ring_sizes'] = [r.Size() for r in rings if r.IsInRing(atom.idx)]\n\n        return output\n\n    def create_bond_feature(self, mol, bid, eid):\n        \"\"\"\n        Create information for a bond for a pair of atoms that are not actually bonded\n\n        Args:\n            mol (pybel.Molecule): Molecule being featurized\n            bid (int): Index of atom beginning of the bond\n            eid (int): Index of atom at the end of the bond\n        \"\"\"\n        a1 = mol.OBMol.GetAtom(bid + 1)\n        a2 = mol.OBMol.GetAtom(eid + 1)\n        same_ring = mol.OBMol.AreInSameRing(a1, a2)\n        return {\"a_idx\": bid,\n                \"b_idx\": eid,\n                \"bond_type\": 0,\n                \"same_ring\": True if same_ring else False,\n                \"spatial_distance\": a1.GetDistance(a2)}\n\n    def create_bond_feature_coupling(self, mol, bid, eid, coupling_type, scalar_coupling_constant, bond_index_pair):\n        \"\"\"\n        Create information for a bond for a pair of atoms that are not actually bonded\n\n        Args:\n            mol (pybel.Molecule): Molecule being featurized\n            bid (int): Index of atom beginning of the bond\n            eid (int): Index of atom at the end of the bond\n        \"\"\"\n        a1 = mol.OBMol.GetAtom(bid + 1)\n        a2 = mol.OBMol.GetAtom(eid + 1)\n        same_ring = mol.OBMol.AreInSameRing(a1, a2)\n        return {\"a_idx\": bid,\n                \"b_idx\": eid,\n                \"bond_type\": 0,\n                \"same_ring\": True if same_ring else False,\n                \"spatial_distance\": a1.GetDistance(a2),\n                \"coupling_type\": coupling_type,\n                \"scalar_coupling_constant\": scalar_coupling_constant,\n                \"bond_index_pair\": bond_index_pair\n                }\n\n    def get_pair_feature(self, mol, bid, eid, full_pair_matrix):\n        \"\"\"\n        Get the features for a certain bond\n\n        Args:\n            mol (pybel.Molecule): Molecule being featurized\n            bid (int): Index of atom beginning of the bond\n            eid (int): Index of atom at the end of the bond\n            full_pair_matrix (bool): Whether to compute the matrix for every atom - even those that\n                are not actually bonded\n        \"\"\"\n        # Find the bonded pair of atoms\n        bond = mol.OBMol.GetBond(bid + 1, eid + 1)\n        if not bond:  # If the bond is ordered in the other direction\n            bond = mol.OBMol.GetBond(eid + 1, bid + 1)\n\n        # If the atoms are not bonded\n        if not bond:\n            if full_pair_matrix:\n                return self.create_bond_feature(mol, bid, eid)\n            else:\n                return None\n\n        # Compute bond features\n        a1 = mol.OBMol.GetAtom(bid + 1)\n        a2 = mol.OBMol.GetAtom(eid + 1)\n        same_ring = mol.OBMol.AreInSameRing(a1, a2)\n        return {\"a_idx\": bid,\n                \"b_idx\": eid,\n                \"bond_type\": 4 if bond.IsAromatic() else bond.GetBondOrder(),\n                \"same_ring\": True if same_ring else False,\n                \"spatial_distance\": a1.GetDistance(a2)}\n\n    def get_pair_feature_coupling1(self, mol, type, bid, eid, full_pair_matrix, target_convert=False,\n                                   scalar_coupling_constant=0):\n        \"\"\"\n        Get the features for a certain bond\n\n        Args:\n            mol (pybel.Molecule): Molecule being featurized\n            bid (int): Index of atom beginning of the bond\n            eid (int): Index of atom at the end of the bond\n            full_pair_matrix (bool): Whether to compute the matrix for every atom - even those that\n                are not actually bonded\n        \"\"\"\n        # print(coupling_attr.columns)\n\n        # coupling_arr0 = coupling_attr[(coupling_attr['atom_index_0']==bid) & (coupling_attr['atom_index_1']==eid)]\n        # coupling_arr1 = coupling_attr[(coupling_attr['atom_index_0'] == eid) & (coupling_attr['atom_index_1'] == bid)]\n        #\n        # # print(bid,eid)\n        #\n        # coupling_type=-1\n        # scalar_coupling_constant=0\n        # bond_index_pair=[]\n        # if(len(coupling_arr0)>0):\n        #     coupling_type = coupling_arr0['type'].tolist()[0]\n        #     if(target_convert):\n        #         scalar_coupling_constant = coupling_arr0['scalar_coupling_constant'].tolist()[0]\n        #     bond_index_pair.append((bid,eid))\n        #     # print(\"ddddd\")\n        #     # print(bond_index_pair)\n        # if (len(coupling_arr1) > 0):\n        #     coupling_type = coupling_arr1['type'].tolist()[0]\n        #     if (target_convert):\n        #         scalar_coupling_constant = coupling_arr1['scalar_coupling_constant'].tolist()[0]\n        #     bond_index_pair.append((eid, bid))\n        # print(\"fffff\")\n        # print(bond_index_pair)\n\n        # print(coupling_arr)\n        ##get coupling attr\n\n        # Find the bonded pair of atoms\n        bond = mol.OBMol.GetBond(bid + 1, eid + 1)\n        if not bond:  # If the bond is ordered in the other direction\n            bond = mol.OBMol.GetBond(eid + 1, bid + 1)\n\n        # If the atoms are not bonded\n        if not bond:\n            if full_pair_matrix:\n                return self.create_bond_feature_coupling(mol, bid, eid, type, scalar_coupling_constant)\n            else:\n                return None\n\n        # Compute bond features\n        a1 = mol.OBMol.GetAtom(bid + 1)\n        a2 = mol.OBMol.GetAtom(eid + 1)\n        same_ring = mol.OBMol.AreInSameRing(a1, a2)\n\n        if (target_convert == False):\n            scalar_coupling_constant = 0\n        # else:\n        #     scalar_coupling_constant = coupling_arr1['scalar_coupling_constant'].tolist()[0]\n        return {\"a_idx\": bid,\n                \"b_idx\": eid,\n                \"bond_type\": 4 if bond.IsAromatic() else bond.GetBondOrder(),\n                \"same_ring\": True if same_ring else False,\n                \"spatial_distance\": a1.GetDistance(a2),\n                \"coupling_type\": type,\n                \"scalar_coupling_constant\": scalar_coupling_constant,\n                }\n\n    def get_pair_feature_coupling(self, mol, coupling_attr, bid, eid, full_pair_matrix, target_convert=False):\n        \"\"\"\n        Get the features for a certain bond\n\n        Args:\n            mol (pybel.Molecule): Molecule being featurized\n            bid (int): Index of atom beginning of the bond\n            eid (int): Index of atom at the end of the bond\n            full_pair_matrix (bool): Whether to compute the matrix for every atom - even those that\n                are not actually bonded\n        \"\"\"\n        # print(coupling_attr.columns)\n        coupling_arr0 = coupling_attr[(coupling_attr['atom_index_0'] == bid) & (coupling_attr['atom_index_1'] == eid)]\n        coupling_arr1 = coupling_attr[(coupling_attr['atom_index_0'] == eid) & (coupling_attr['atom_index_1'] == bid)]\n\n        # print(bid,eid)\n\n        coupling_type = -1\n        scalar_coupling_constant = 0\n        bond_index_pair = []\n        if (len(coupling_arr0) > 0):\n            coupling_type = coupling_arr0['type'].tolist()[0]\n            if (target_convert):\n                scalar_coupling_constant = coupling_arr0['scalar_coupling_constant'].tolist()[0]\n            bond_index_pair.append((bid, eid))\n            # print(\"ddddd\")\n            # print(bond_index_pair)\n        if (len(coupling_arr1) > 0):\n            coupling_type = coupling_arr1['type'].tolist()[0]\n            if (target_convert):\n                scalar_coupling_constant = coupling_arr1['scalar_coupling_constant'].tolist()[0]\n            bond_index_pair.append((eid, bid))\n            # print(\"fffff\")\n            # print(bond_index_pair)\n\n        # print(coupling_arr)\n        ##get coupling attr\n\n        # Find the bonded pair of atoms\n        bond = mol.OBMol.GetBond(bid + 1, eid + 1)\n        if not bond:  # If the bond is ordered in the other direction\n            bond = mol.OBMol.GetBond(eid + 1, bid + 1)\n\n        # If the atoms are not bonded\n        if not bond:\n            if full_pair_matrix:\n                return self.create_bond_feature_coupling(mol, bid, eid, coupling_type, scalar_coupling_constant,\n                                                         bond_index_pair)\n            else:\n                return None\n\n        # Compute bond features\n        a1 = mol.OBMol.GetAtom(bid + 1)\n        a2 = mol.OBMol.GetAtom(eid + 1)\n        same_ring = mol.OBMol.AreInSameRing(a1, a2)\n\n        if (target_convert == False):\n            scalar_coupling_constant = 0\n        # else:\n        #     scalar_coupling_constant = coupling_arr1['scalar_coupling_constant'].tolist()[0]\n        return {\"a_idx\": bid,\n                \"b_idx\": eid,\n                \"bond_type\": 4 if bond.IsAromatic() else bond.GetBondOrder(),\n                \"same_ring\": True if same_ring else False,\n                \"spatial_distance\": a1.GetDistance(a2),\n                \"coupling_type\": coupling_type,\n                \"scalar_coupling_constant\": scalar_coupling_constant,\n                \"bond_index_pair\": bond_index_pair\n                }\n\n    def _get_rdk_mol(self, mol, format='smiles'):\n        \"\"\"\n        Return: RDKit Mol (w/o H)\n        \"\"\"\n        if format == 'pdb':\n            return Chem.rdmolfiles.MolFromPDBBlock(mol.write(\"pdb\"))\n        elif format == 'smiles':\n            return Chem.rdmolfiles.MolFromSmiles(mol.write(\"smiles\"))\n\n    def _get_chiral_centers(self, mol):\n        \"\"\"\n        Use RDKit to find the chiral centers with CIP(R/S) label\n\n        This provides the absolute stereochemistry.  The chiral label obtained\n        from pybabel and rdkit.mol.getchiraltag is relative positions of the bonds as provided\n\n        Args:\n            mol (Molecule): Molecule to asses\n        Return:\n            (dict): Keys are the atom index and values are the CIP label\n        \"\"\"\n        mol_rdk = self._get_rdk_mol(mol, 'smiles')\n        if mol_rdk is None:\n            # Conversion to RDKit has failed\n            return {}\n        else:\n            chiral_cc = Chem.FindMolChiralCenters(mol_rdk)\n            return dict(chiral_cc)\n\n\ndef dijkstra_distance(bonds):\n    \"\"\"\n    Compute the graph distance based on the dijkstra algorithm\n\n    Args:\n        bonds: (list of list), for example [[0, 1], [1, 2]] means two bonds formed by atom 0, 1 and atom 1, 2\n\n    Returns:\n        full graph distance matrix\n    \"\"\"\n    nb_atom = max(itertools.chain(*bonds)) + 1\n    graph_dist = np.ones((nb_atom, nb_atom), dtype=np.int32) * np.infty\n    for bond in bonds:\n        graph_dist[bond[0], bond[1]] = 1\n        graph_dist[bond[1], bond[0]] = 1\n\n    queue = deque()  # Queue used in all loops\n    visited = set()  # Used in all loops\n    for i in range(nb_atom):\n        graph_dist[i, i] = 0\n        visited.clear()\n        queue.append(i)\n        while queue:\n            s = queue.pop()\n            visited.add(s)\n\n            for k in np.where(graph_dist[s, :] == 1)[0]:\n                if k not in visited:\n                    queue.append(k)\n                    graph_dist[i, k] = min(graph_dist[i, k],\n                                           graph_dist[i, s] + 1)\n                    graph_dist[k, i] = graph_dist[i, k]\n    return graph_dist\n\n\ndef mol_from_smiles(smiles):\n    mol = pybel.readstring(format='smi', string=smiles)\n    mol.make3D()\n    return mol\n\n\ndef mol_from_pymatgen(mol):\n    \"\"\"\n    Args:\n        mol(Molecule)\n    \"\"\"\n    mol = pybel.Molecule(BabelMolAdaptor(mol).openbabel_mol)\n    mol.make3D()\n    return mol\n\n\ndef mol_from_file(file_path, file_format='xyz'):\n    \"\"\"\n    Args:\n        file_path(str)\n        file_format(str): allow formats that open babel supports\n    \"\"\"\n    mol = [r for r in pybel.readfile(format=file_format,\n                                     filename=file_path)][0]\n    return mol\n\n\ndef _convert_mol(mol, molecule_format, converter):\n    \"\"\"Convert a molecule from string to its graph features\n\n    Utility function used in the graph generator.\n\n    The parse and convert operations are both in this function due to Pybel objects\n    not being serializable. By not using the Pybel representation of each molecule\n    as an input to this function, we can use multiprocessing to parallelize conversion\n    over molecules as strings can be passed as pickle objects to the worker threads but\n    but Pybel objects cannot.\n\n    Args:\n        mol (str): String representation of a molecule\n        molecule_format (str): Format of the string representation\n        converter (MolecularGraph): Tool used to generate graph representation\n    Returns:\n        (dict): Graph representation of the molecule\n    \"\"\"\n\n    # Convert molecule into pybel format\n    if molecule_format == 'smiles':\n        mol = mol_from_smiles(mol)  # Used to generate 3D coordinates/H atoms\n    else:\n        mol = pybel.readstring(molecule_format, mol)\n\n    return converter.convert(mol)\n\n\nclass MolecularGraphBatchGenerator(BaseGraphBatchGenerator):\n    \"\"\"Generator that creates batches of molecular data by computing graph properties on demand\n\n    If your dataset is small enough that the descriptions of the whole dataset fit in memory,\n    we recommend using :class:`megnet.data.graph.GraphBatchGenerator` instead to avoid\n    the computational cost of dynamically computing graphs.\"\"\"\n\n    def __init__(self, mols, targets=None, converter=None, molecule_format='xyz',\n                 batch_size=128, shuffle=True, n_jobs=1):\n        \"\"\"\n        Args:\n            mols ([str]): List of the string reprensetations of each molecule\n            targets ([ndarray]): Properties of each molecule to be predicted\n            converter (MolecularGraph): Converter used to generate graph features\n            molecule_format (str): Format of each of the string representations in `mols`\n            batch_size (int): Target size for each batch\n            shuffle (bool): Whether to shuffle the training data after each epoch\n            n_jobs (int): Number of worker threads (None to use all threads).\n        \"\"\"\n\n        super().__init__(len(mols), targets, batch_size, shuffle)\n        self.mols = np.array(mols)\n        if converter is None:\n            converter = MolecularGraph()\n        self.converter = converter\n        self.molecule_format = molecule_format\n        self.n_jobs = n_jobs\n\n        def mute():\n            sys.stdout = open(os.devnull, 'w')\n            sys.stderr = open(os.devnull, 'w')\n\n        self.pool = Pool(self.n_jobs, initializer=mute) if self.n_jobs != 1 else None\n\n    def __del__(self):\n        if self.pool is not None:\n            self.pool.close()  # Kill thread pool if generator is deleted\n\n    def _generate_inputs(self, batch_index):\n        # Get the molecules for this batch\n        mols = self.mols[batch_index]\n\n        # Generate the graphs\n        graphs = self._generate_graphs(mols)\n\n        # Return them as flattened into array format\n        return self.converter.get_flat_data(graphs)\n\n    def _generate_graphs(self, mols):\n        \"\"\"Generate graphs for a certain collection of molecules\n\n        Args:\n            mols ([string]): Molecules to process\n        Returns:\n            ([dict]): Graphs for all of the molecules\n        \"\"\"\n        if self.pool is None:\n            graphs = [_convert_mol(m, self.molecule_format, self.converter) for m in mols]\n        else:\n            func = partial(_convert_mol, molecule_format=self.molecule_format,\n                           converter=self.converter)\n            graphs = self.pool.map(func, mols)\n        return graphs\n\n    def create_cached_generator(self) -> GraphBatchGenerator:\n        \"\"\"Generates features for all of the molecules and stores them in memory\n\n        Returns:\n            (GraphBatchGenerator) Graph genereator that relies on having the graphs in memory\n        \"\"\"\n\n        # Make all the graphs\n        graphs = self._generate_graphs(self.mols)\n\n        # Turn them into a fat array\n        inputs = self.converter.get_flat_data(graphs, self.targets)\n\n        return GraphBatchGenerator(*inputs, is_shuffle=self.is_shuffle,\n                                   batch_size=self.batch_size)\n\n\nclass GraphModel:\n    \"\"\"\n    Composition of keras model and converter class for transfering structure\n    object to input tensors. We add methods to train the model from\n    (structures, targets) pairs\n\n    Args:\n        model: (keras model)\n        graph_converter: (object) a object that turns a structure to a graph,\n            check `megnet.data.crystal`\n        target_scaler: (object) a scaler object for converting targets, check\n            `megnet.utils.preprocessing`\n        metadata: (dict) An optional dict of metadata associated with the model.\n            Recommended to incorporate some basic information such as units,\n            MAE performance, etc.\n\n    \"\"\"\n\n    def __init__(self,\n                 model,\n                 graph_converter,\n                 target_scaler=DummyScaler(),\n                 metadata=None,\n                 **kwargs):\n        self.model = model\n        self.graph_converter = graph_converter\n        self.target_scaler = target_scaler\n        self.metadata = metadata or {}\n\n    def __getattr__(self, p):\n        return getattr(self.model, p)\n\n    def train(self,\n              train_structures,\n              train_targets,\n              validation_structures=None,\n              validation_targets=None,\n              epochs=1000,\n              batch_size=128,\n              verbose=1,\n              callbacks=None,\n              scrub_failed_structures=False,\n              prev_model=None,\n              lr_scaling_factor=0.5,\n              patience=500,\n              **kwargs):\n        \"\"\"\n        Args:\n            train_structures: (list) list of pymatgen structures\n            train_targets: (list) list of target values\n            validation_structures: (list) list of pymatgen structures as validation\n            validation_targets: (list) list of validation targets\n            epochs: (int) number of epochs\n            batch_size: (int) training batch size\n            verbose: (int) keras fit verbose, 0 no progress bar, 1 only at the epoch end and 2 every batch\n            callbacks: (list) megnet or keras callback functions for training\n            scrub_failed_structures: (bool) whether to scrub structures with failed graph computation\n            prev_model: (str) file name for previously saved model\n            lr_scaling_factor: (float, less than 1) scale the learning rate down when nan loss encountered\n            patience: (int) patience for early stopping\n            **kwargs:\n        \"\"\"\n        train_graphs, train_targets = self.get_all_graphs_targets(train_structures, train_targets,\n                                                                  scrub_failed_structures=scrub_failed_structures)\n        if validation_structures is not None:\n            val_graphs, validation_targets = self.get_all_graphs_targets(\n                validation_structures, validation_targets, scrub_failed_structures=scrub_failed_structures)\n        else:\n            val_graphs = None\n\n        self.train_from_graphs(train_graphs,\n                               train_targets,\n                               validation_graphs=val_graphs,\n                               validation_targets=validation_targets,\n                               epochs=epochs,\n                               batch_size=batch_size,\n                               verbose=verbose,\n                               callbacks=callbacks,\n                               prev_model=prev_model,\n                               lr_scaling_factor=lr_scaling_factor,\n                               patience=patience,\n                               **kwargs\n                               )\n\n    def train_coupling(self,\n                       mols,\n                       validation_structures=None,\n                       validation_targets=None,\n                       epochs=1000,\n                       batch_size=128,\n                       verbose=1,\n                       callbacks=None,\n                       scrub_failed_structures=False,\n                       prev_model=None,\n                       lr_scaling_factor=0.5,\n                       patience=500,\n                       **kwargs):\n        \"\"\"\n        Args:\n            train_structures: (list) list of pymatgen structures\n            train_targets: (list) list of target values\n            validation_structures: (list) list of pymatgen structures as validation\n            validation_targets: (list) list of validation targets\n            epochs: (int) number of epochs\n            batch_size: (int) training batch size\n            verbose: (int) keras fit verbose, 0 no progress bar, 1 only at the epoch end and 2 every batch\n            callbacks: (list) megnet or keras callback functions for training\n            scrub_failed_structures: (bool) whether to scrub structures with failed graph computation\n            prev_model: (str) file name for previously saved model\n            lr_scaling_factor: (float, less than 1) scale the learning rate down when nan loss encountered\n            patience: (int) patience for early stopping\n            **kwargs:\n        \"\"\"\n\n#         train_graphs, train_targets = self.get_all_graphs_targets_coupling(mols,\n#                                                                           scrub_failed_structures=scrub_failed_structures)\n\n        data_graphs=[]\n        data_targets=[]\n\n        with open('../input/shiyan-graph/train_graphs', 'rb') as fp:\n            data_graphs = pickle.load(fp)\n\n        with open('../input/shiyan-graph/train_targets', 'rb') as fp:\n            data_targets = pickle.load(fp)\n\n#         with open('../input/shiyan-graph-1/train_graphs', 'rb') as fp:\n#             data_graphs = pickle.load(fp)\n\n#         with open('../input/shiyan-graph-1/train_targets', 'rb') as fp:\n#             data_targets = pickle.load(fp)\n        \n#         max_value =204.88 \n#         min_value = -36.2186\n\n#         for target in  data_targets:\n#             target['output'] =  (np.array(target['output'])-min_value)/(max_value-min_value) \n#         data_targets[0]\n# #         print(len(data_targets))\n#         print(data_targets[0])\n        train_graphs, val_graphs, train_targets, validation_targets = train_test_split(data_graphs, data_targets,\n                                                                                       test_size=0.20, random_state=42)\n\n#         if validation_structures is not None:\n#             val_graphs, validation_targets = self.get_all_graphs_targets_coupling(\n#                 mols, scrub_failed_structures=scrub_failed_structures)\n#         else:\n#             val_graphs = None\n\n        self.train_from_graphs_coupling(train_graphs,\n                                        train_targets,\n                                        validation_graphs=val_graphs,\n                                        validation_targets=validation_targets,\n                                        epochs=epochs,\n                                        batch_size=batch_size,\n                                        verbose=verbose,\n                                        callbacks=callbacks,\n                                        prev_model=prev_model,\n                                        lr_scaling_factor=lr_scaling_factor,\n                                        patience=patience,\n                                        **kwargs\n                                        )\n\n    def train_from_graphs_coupling(self,\n                                   train_graphs,\n                                   train_targets,\n                                   validation_graphs=None,\n                                   validation_targets=None,\n                                   epochs=1000,\n                                   batch_size=1,\n                                   verbose=2,\n                                   callbacks=None,\n                                   prev_model=None,\n                                   lr_scaling_factor=0.5,\n                                   patience=500,\n                                   **kwargs\n                                   ):\n\n        # load from saved model\n        if prev_model:\n            self.load_weights(prev_model)\n        is_classification = False\n        monitor = 'val_acc' if is_classification else 'val_loss'\n        mode = 'max' if is_classification else 'min'\n        dirname = kwargs.pop('dirname', 'callback')\n        if not os.path.isdir(dirname):\n            os.makedirs(dirname)\n        if callbacks is None:\n            filepath = \"weights.best.hdf5\"\n            checkpoint = ModelCheckpoint(filepath, monitor='val_loss', verbose=1, save_best_only=True,\n                                         mode='min')\n#             stop = keras.callbacks.EarlyStopping( monitor='val_loss', patience=0, verbose=0, mode='auto')\n\n            callbacks= [checkpoint,stop_early,lrate]\n        train_nb_atoms = [len(i['atom']) for i in train_graphs]\n        train_targets = [self.target_scaler.transform(i, j) for i, j in zip(train_targets, train_nb_atoms)]\n\n        if validation_graphs is not None:\n            filepath = os.path.join(dirname, '%s_{epoch:05d}_{%s:.6f}.hdf5' % (monitor, monitor))\n            val_nb_atoms = [len(i['atom']) for i in validation_graphs]\n            validation_targets = [self.target_scaler.transform(i, j) for i, j in zip(validation_targets, val_nb_atoms)]\n            val_inputs = self.graph_converter.get_flat_scarlar_data(validation_graphs, validation_targets)\n\n            val_generator = self._create_generator(*val_inputs,\n                                                   batch_size=batch_size)\n            steps_per_val = int(np.ceil(len(validation_graphs) / batch_size))\n#             callbacks.extend([ReduceLRUponNan(filepath=filepath,\n#                                               monitor=monitor,\n#                                               mode=mode,\n#                                               factor=lr_scaling_factor,\n#                                               patience=patience,\n#                                               )])\n#             callbacks.extend([ModelCheckpointMAE(filepath=filepath,\n#                                                  monitor=monitor,\n#                                                  mode=mode,\n#                                                  save_best_only=True,\n#                                                  save_weights_only=False,\n#                                                  val_gen=val_generator,\n#                                                  steps_per_val=steps_per_val,\n#                                                  target_scaler=self.target_scaler)])\n        else:\n            val_generator = None\n            steps_per_val = None\n        # print(np.array(train_graphs[2]).shape)\n        train_inputs = self.graph_converter.get_flat_scarlar_data(train_graphs, train_targets)\n\n        # print(np.array(train_inputs[1]).shape)\n        # print(np.array(train_inputs[5]).shape)\n        # check dimension match\n        # print(np.array(train_graphs).shape)\n        self.check_dimension(train_graphs[0])\n        # ########################################new dynamics code####################\n        #         target_ph = tf.placeholder(tf.float32, shape=(None, 1))\n        #         target_data = tf.data.Dataset.from_tensor_slices(target_ph)\n        #         target_data = target_data.batch(batch_size)\n        #         target_iter = target_data.make_initializable_iterator()\n        #         target = target_iter.get_next()\n        #         learning_rate = 1e-3  # Learning rate\n        #\n        #         optimizer = Adam(lr=learning_rate)\n        #         self.compile(optimizer=optimizer, loss='mse', target_tensors=target)\n        #\n        #         # Training setup\n        #         sess = K.get_session()\n        #         print(batch_size)\n        #\n        #         # input = [train_inputs]\n        #         batches_train = batch_iterator([train_inputs[0],train_inputs[1],train_inputs[2],train_inputs[3],train_inputs[4],train_inputs[5],train_inputs[6],train_inputs[7]], 1, epochs)\n        #         loss = 0\n        #         batch_index = 0\n        #         batches_in_epoch = np.ceil(len(train_graphs[0]) / batch_size)\n        #\n        #         # Training loop\n        #         for b in batches_train:\n        #             # batch = Batch(b[0], b[1],b[2],b[3],b[4],b[5],b[6])\n        #             y_ = b[7]\n        #             print(y_[0])\n        #             print(y_[0].shape)\n        #             sess.run(target_iter.initializer, feed_dict={target_ph: y_})\n        #             loss += self.train_on_batch([b[0], b[1],b[2],b[3],b[4],b[5],b[6]], None)\n        #\n        #             batch_index += 1\n        #             if batch_index == batches_in_epoch:\n        #                 print(batch_index)\n        #                 print('Loss: {}'.format(loss / batches_in_epoch))\n        #                 loss = 0\n        #                 # batch_index = 0\n        #         print(batch_index)\n\n        train_generator = self._create_generator(*train_inputs, batch_size=batch_size)\n        steps_per_train = int(np.ceil(len(train_graphs) / batch_size))\n        self.fit_generator(train_generator, steps_per_epoch=steps_per_train,\n                           validation_data=val_generator, validation_steps=steps_per_val,\n                           epochs=epochs, initial_epoch = 0, verbose=verbose, callbacks=callbacks, **kwargs)\n\n    def train_from_graphs(self,\n                          train_graphs,\n                          train_targets,\n                          validation_graphs=None,\n                          validation_targets=None,\n                          epochs=1000,\n                          batch_size=128,\n                          verbose=1,\n                          callbacks=None,\n                          prev_model=None,\n                          lr_scaling_factor=0.5,\n                          patience=500,\n                          **kwargs\n                          ):\n\n        # load from saved model\n        if prev_model:\n            self.load_weights(prev_model)\n        is_classification = 'entropy' in self.model.loss\n        monitor = 'val_acc' if is_classification else 'val_mae'\n        mode = 'max' if is_classification else 'min'\n        dirname = kwargs.pop('dirname', 'callback')\n        if not os.path.isdir(dirname):\n            os.makedirs(dirname)\n        if callbacks is None:\n            # with this call back you can stop the model training by `touch STOP`\n            keras.callbacks.EarlyStopping(monitor=monitor, patience=0, verbose=0, mode='auto')\n\n            callbacks = [ManualStop(),EarlyStopping]\n        train_nb_atoms = [len(i['atom']) for i in train_graphs]\n        train_targets = [self.target_scaler.transform(i, j) for i, j in zip(train_targets, train_nb_atoms)]\n\n        if validation_graphs is not None:\n            filepath = os.path.join(dirname, '%s_{epoch:05d}_{%s:.6f}.hdf5' % (monitor, monitor))\n            val_nb_atoms = [len(i['atom']) for i in validation_graphs]\n            validation_targets = [self.target_scaler.transform(i, j) for i, j in\n                                  zip(validation_targets, val_nb_atoms)]\n            val_inputs = self.graph_converter.get_flat_data(validation_graphs, validation_targets)\n\n            val_generator = self._create_generator(*val_inputs,\n                                                   batch_size=batch_size)\n            steps_per_val = int(np.ceil(len(validation_graphs) / batch_size))\n            callbacks.extend([ReduceLRUponNan(filepath=filepath,\n                                              monitor=monitor,\n                                              mode=mode,\n                                              factor=lr_scaling_factor,\n                                              patience=patience,\n                                              )])\n            callbacks.extend([ModelCheckpointMAE(filepath=filepath,\n                                                 monitor=monitor,\n                                                 mode=mode,\n                                                 save_best_only=True,\n                                                 save_weights_only=False,\n                                                 val_gen=val_generator,\n                                                 steps_per_val=steps_per_val,\n                                                 target_scaler=self.target_scaler)])\n        else:\n            val_generator = None\n            steps_per_val = None\n\n        train_inputs = self.graph_converter.get_flat_data(train_graphs, train_targets)\n        # check dimension match\n        self.check_dimension(train_graphs[0])\n        train_generator = self._create_generator(*train_inputs, batch_size=batch_size)\n        steps_per_train = int(np.ceil(len(train_graphs) / batch_size))\n        self.fit_generator(train_generator, steps_per_epoch=steps_per_train,\n                           validation_data=val_generator, validation_steps=steps_per_val,\n                           epochs=epochs, verbose=verbose, callbacks=callbacks, **kwargs)\n\n    def check_dimension(self, graph):\n        \"\"\"\n        Check the model dimension against the graph converter dimension\n        Args:\n            graph: structure graph\n\n        Returns:\n\n        \"\"\"\n        test_inp = self.graph_converter.graph_to_input(graph)\n        input_shapes = [i.shape for i in test_inp]\n\n        model_input_shapes = [int_shape(i) for i in self.model.inputs]\n\n        def _check_match(real_shape, tensor_shape):\n            if len(real_shape) != len(tensor_shape):\n                return False\n            matched = True\n            for i, j in zip(real_shape, tensor_shape):\n                if j is None:\n                    continue\n                else:\n                    if i == j:\n                        continue\n                    else:\n                        matched = False\n            return matched\n\n        for i, j, k in zip(['atom features', 'bond features', 'state features'],\n                           input_shapes[:3], model_input_shapes[:3]):\n            matched = _check_match(j, k)\n            if not matched:\n                raise ValueError(\"The data dimension for %s is %s and does not match model \"\n                                 \"required shape of %s\" % (i, str(j), str(k)))\n\n    def get_all_graphs_targets(self, structures, targets, scrub_failed_structures=False):\n        \"\"\"\n        Compute the graphs from structures and spit out (graphs, targets) with options to\n        automatically remove structures with failed graph computations\n\n        Args:\n            structures: (list) pymatgen structure list\n            targets: (list) target property list\n            scrub_failed_structures: (bool) whether to scrub those failed structures\n\n        Returns:\n            graphs, targets\n\n        \"\"\"\n        graphs_valid = []\n        targets_valid = []\n\n        for i, (s, t) in enumerate(zip(structures, targets)):\n            try:\n                graph = self.graph_converter.convert(s, t)\n                graphs_valid.append(graph)\n                targets_valid.append(t)\n            except Exception as e:\n                if scrub_failed_structures:\n                    warn(\"structure with index %d failed the graph computations\" % i,\n                         UserWarning)\n                    continue\n                else:\n                    raise RuntimeError(str(e))\n        return graphs_valid, targets_valid\n\n    def get_all_graphs_targets_coupling(self, mols, scrub_failed_structures=False, target_convert=True):\n        \"\"\"\n        Compute the graphs from structures and spit out (graphs, targets) with options to\n        automatically remove structures with failed graph computations\n\n        Args:\n            structures: (list) pymatgen structure list\n            targets: (list) target property list\n            scrub_failed_structures: (bool) whether to scrub those failed structures\n\n        Returns:\n            graphs, targets\n\n        \"\"\"\n        graphs_valid = []\n        targets_valid = []\n        # target = np.zeros((len(mols), 135), dtype=np.float32)\n\n        with tqdm(total=len(mols)) as pbar:\n            for i, (s, t, energy_dipole, charges_magnetic) in enumerate(mols):\n                try:\n                    mol = read_ob_molecule(s)\n                    mol = pybel.Molecule(mol)\n                    graph, train = self.graph_converter.convert_coupling(mol, t, energy_dipole, charges_magnetic,\n                                                                         target_convert)\n                    # self.check_dimension(graph)\n                    # graph = self.graph_converter.graph_to_input(graph)\n                    graphs_valid.append(graph)\n\n                    #\n                    # # Get the coupling constant data in the order request by the user\n                    # train_coupling = np.array(t).tolist()\n                    # # train_coupling= t\n                    # # bonds = []\n                    # index1_temp = []\n                    # index2_temp = []\n                    # for train in train_coupling:\n                    #     # Store the index of each bond\n                    #     index1_temp.append(train[2])\n                    #     index2_temp.append(train[3])\n                    #\n                    #\n                    # # Given the train (i,j), make it so (i,j) == (j, i)\n                    # index1 = index1_temp + index2_temp\n                    # index2 = index2_temp + index1_temp\n                    # train_coupling = train_coupling + train_coupling\n                    #\n                    # # Sort the arrays by the beginning index\n                    # sorted_arg = np.argsort(index1)\n                    # index1 = np.array(index1)[sorted_arg].tolist()\n                    # index2 = np.array(index2)[sorted_arg].tolist()\n                    # train = np.array(train_coupling)[sorted_arg].tolist()\n                    # train_dict = {\n                    # 'train': train,\n                    # 'index1': index1,\n                    # 'index2': index2\n                    # }\n                    # # target_graph = self.graph_converter.convert(s, t,True)\n                    targets_valid.append(train)\n                except Exception as e:\n                    if scrub_failed_structures:\n                        warn(\"structure with index %d failed the graph computations\" % i,\n                             UserWarning)\n                        continue\n                    else:\n                        raise RuntimeError(str(e))\n                pbar.update()\n        return graphs_valid, targets_valid\n\n    def predict_structure(self, structure):\n        \"\"\"\n        Predict property from structure\n\n        Args:\n            structure: pymatgen structure or molecule\n\n        Returns:\n            predicted target value\n        \"\"\"\n        graph = self.graph_converter.convert(structure)\n        return self.predict_graph(graph)\n\n    def predict_graph(self, graph):\n        \"\"\"\n        Predict property from graph\n\n        Args:\n            graph: a graph dictionary, see megnet.data.graph\n\n        Returns:\n            predicted target value\n\n        \"\"\"\n        inp = self.graph_converter.graph_to_input(graph)\n        return self.target_scaler.inverse_transform(self.predict(inp).ravel(), len(graph['atom']))\n\n    def _create_generator(self, *args, **kwargs):\n        if hasattr(self.graph_converter, 'bond_converter'):\n            kwargs.update({'distance_converter': self.graph_converter.bond_converter})\n            return GraphBatchDistanceConvert(*args, **kwargs)\n        else:\n            return GraphBatchGenerator(*args, **kwargs)\n\n    def save_model(self, filename):\n        \"\"\"\n        Save the model to a keras model hdf5 and a json config for additional\n        converters\n\n        Args:\n            filename: (str) output file name\n\n        Returns:\n            None\n        \"\"\"\n        self.model.save(filename)\n        dumpfn(\n            {\n                'graph_converter': self.graph_converter,\n                'target_scaler': self.target_scaler,\n                'metadata': self.metadata\n            },\n            filename + '.json'\n        )\n\n    @classmethod\n    def from_file(cls, filename):\n        \"\"\"\n        Class method to load model from\n            filename for keras model\n            filename.json for additional converters\n\n        Args:\n            filename: (str) model file name\n\n        Returns\n            GraphModel\n        \"\"\"\n        configs = loadfn(filename + '.json')\n        from keras.models import load_model\n        from megnet.layers import _CUSTOM_OBJECTS\n        model = load_model(filename, custom_objects=_CUSTOM_OBJECTS)\n        configs.update({'model': model})\n        return GraphModel(**configs)\n\n    @classmethod\n    def from_url(cls, url):\n        \"\"\"\n        Download and load a model from a URL. E.g.\n        https://github.com/materialsvirtuallab/megnet/blob/master/mvl_models/mp-2019.4.1/formation_energy.hdf5\n\n        Args:\n            url: (str) url link of the model\n\n        Returns:\n            GraphModel\n        \"\"\"\n        import urllib.request\n        fname = url.split(\"/\")[-1]\n        urllib.request.urlretrieve(url, fname)\n        urllib.request.urlretrieve(url + \".json\", fname + \".json\")\n        return cls.from_file(fname)\n\n\nclass MEGNetModel(GraphModel):\n    \"\"\"\n    Construct a graph network model with or without explicit atom features\n    if n_feature is specified then a general graph model is assumed,\n    otherwise a crystal graph model with z number as atom feature is assumed.\n\n    Args:\n        nfeat_edge: (int) number of bond features\n        nfeat_global: (int) number of state features\n        nfeat_node: (int) number of atom features\n        nblocks: (int) number of MEGNetLayer blocks\n        lr: (float) learning rate\n        n1: (int) number of hidden units in layer 1 in MEGNetLayer\n        n2: (int) number of hidden units in layer 2 in MEGNetLayer\n        n3: (int) number of hidden units in layer 3 in MEGNetLayer\n        nvocal: (int) number of total element\n        embedding_dim: (int) number of embedding dimension\n        nbvocal: (int) number of bond types if bond attributes are types\n        bond_embedding_dim: (int) number of bond embedding dimension\n        ngvocal: (int) number of global types if global attributes are types\n        global_embedding_dim: (int) number of global embedding dimension\n        npass: (int) number of recurrent steps in Set2Set layer\n        ntarget: (int) number of output targets\n        act: (object) activation function\n        l2_coef: (float or None) l2 regularization parameter\n        is_classification: (bool) whether it is a classification task\n        loss: (object or str) loss function\n        metrics: (list or dict) List or dictionary of Keras metrics to be evaluated by the model during training and testing\n        dropout: (float) dropout rate\n        graph_converter: (object) object that exposes a \"convert\" method for structure to graph conversion\n        target_scaler: (object) object that exposes a \"transform\" and \"inverse_transform\" methods for transforming the target values\n        optimizer_kwargs (dict): extra keywords for optimizer, for example clipnorm and clipvalue\n    \"\"\"\n\n    def __init__(self,\n                 nfeat_edge=None,\n                 nfeat_global=None,\n                 nfeat_node=None,\n                 nblocks=3,\n                 lr=1e-3,\n                 n1=64,\n                 n2=32,\n                 n3=16,\n                 nvocal=95,\n                 embedding_dim=16,\n                 nbvocal=None,\n                 bond_embedding_dim=None,\n                 ngvocal=None,\n                 global_embedding_dim=None,\n                 npass=3,\n                 ntarget=1,\n                 act=softplus2,\n                 is_classification=False,\n                 loss=\"mae\",\n                 metrics=None,\n                 l2_coef=None,\n                 dropout=None,\n                 graph_converter=None,\n                 target_scaler=DummyScaler(),\n                 optimizer_kwargs=None,\n                 dropout_on_predict=False\n                 ):\n\n        # Build the MEG Model\n        model = make_megnet_model(nfeat_edge=nfeat_edge,\n                                  nfeat_global=nfeat_global,\n                                  nfeat_node=nfeat_node,\n                                  nblocks=nblocks,\n                                  n1=n1,\n                                  n2=n2,\n                                  n3=n3,\n                                  nvocal=nvocal,\n                                  embedding_dim=embedding_dim,\n                                  nbvocal=nbvocal,\n                                  bond_embedding_dim=bond_embedding_dim,\n                                  ngvocal=ngvocal,\n                                  global_embedding_dim=global_embedding_dim,\n                                  npass=npass,\n                                  ntarget=ntarget,\n                                  act=act,\n                                  is_classification=is_classification,\n                                  l2_coef=l2_coef,\n                                  dropout=dropout,\n                                  dropout_on_predict=dropout_on_predict)\n\n        # Compile the model with the optimizer\n        # loss = 'binary_crossentropy' if is_classification else loss\n        # loss = 'binary_crossentropy' if is_classification else loss\n        # loss = log_mae+mse\n        opt_params = {'lr': lr}\n        if optimizer_kwargs is not None:\n            opt_params.update(optimizer_kwargs)\n        # model.compile(Adam(**opt_params), loss, metrics=metrics)\n        model.compile(Adam(**opt_params), log_mae, metrics=[mse])\n#         model.compile(opt, log_mae, metrics=[mse])\n        if graph_converter is None:\n            graph_converter = CrystalGraph(cutoff=4, bond_converter=GaussianDistance(np.linspace(0, 5, 100), 0.5))\n\n        super().__init__(model=model, target_scaler=target_scaler, graph_converter=graph_converter)\n\n\ndef make_megnet_model(nfeat_edge=None, nfeat_global=None, nfeat_node=None, nblocks=3,\n                      n1=64, n2=32, n3=16, nvocal=95, embedding_dim=16, nbvocal=None,\n                      bond_embedding_dim=None, ngvocal=None, global_embedding_dim=None,\n                      npass=3, ntarget=1, act=softplus2, is_classification=False,\n                      l2_coef=None, dropout=None, dropout_on_predict=False):\n    \"\"\"Make a MEGNet Model\n\n    Args:\n        nfeat_edge: (int) number of bond features\n        nfeat_global: (int) number of state features\n        nfeat_node: (int) number of atom features\n        nblocks: (int) number of MEGNetLayer blocks\n        n1: (int) number of hidden units in layer 1 in MEGNetLayer\n        n2: (int) number of hidden units in layer 2 in MEGNetLayer\n        n3: (int) number of hidden units in layer 3 in MEGNetLayer\n        nvocal: (int) number of total element\n        embedding_dim: (int) number of embedding dimension\n        nbvocal: (int) number of bond types if bond attributes are types\n        bond_embedding_dim: (int) number of bond embedding dimension\n        ngvocal: (int) number of global types if global attributes are types\n        global_embedding_dim: (int) number of global embedding dimension\n        npass: (int) number of recurrent steps in Set2Set layer\n        ntarget: (int) number of output targets\n        act: (object) activation function\n        l2_coef: (float or None) l2 regularization parameter\n        is_classification: (bool) whether it is a classification task\n        dropout: (float) dropout rate\n        dropout_on_predict (bool): Whether to use dropout during prediction and training\n    Returns:\n        (Model) Keras model, ready to run\n    \"\"\"\n\n    # Get the setting for the training kwarg of Dropout\n    dropout_training = True if dropout_on_predict else None\n\n    # Create the input blocks\n    int32 = 'int32'\n    if nfeat_node is None:\n        x1 = Input(shape=(None,), dtype=int32)  # only z as feature\n        x1_ = Embedding(nvocal, embedding_dim)(x1)\n    else:\n        x1 = Input(shape=(None, nfeat_node))\n        x1_ = x1\n    if nfeat_edge is None:\n        x2 = Input(shape=(None,), dtype=int32)\n        x2_ = Embedding(nbvocal, bond_embedding_dim)(x2)\n    else:\n        x2 = Input(shape=(None, nfeat_edge))\n        x2_ = x2\n    if nfeat_global is None:\n        x3 = Input(shape=(None,), dtype=int32)\n        x3_ = Embedding(ngvocal, global_embedding_dim)(x3)\n    else:\n        x3 = Input(shape=(None, nfeat_global))\n        x3_ = x3\n    x4 = Input(shape=(None,), dtype=int32)\n    x5 = Input(shape=(None,), dtype=int32)\n    x6 = Input(shape=(None,), dtype=int32)\n    x7 = Input(shape=(None,), dtype=int32)\n    if l2_coef is not None:\n        reg = l2(l2_coef)\n    else:\n        reg = None\n\n    # two feedforward layers\n    def ff(x, n_hiddens=[n1, n2]):\n        out = x\n        for i in n_hiddens:\n            out = Dense(i, activation=act, kernel_regularizer=reg)(out)\n        return out\n\n    # a block corresponds to two feedforward layers + one MEGNetLayer layer\n    # Note the first block does not contain the feedforward layer since\n    # it will be explicitly added before the block\n    def one_block(a, b, c, has_ff=True):\n        if has_ff:\n            x1_ = ff(a)\n            x2_ = ff(b)\n            x3_ = ff(c)\n        else:\n            x1_ = a\n            x2_ = b\n            x3_ = c\n        out = MEGNetLayer(\n            [n1, n1, n2], [n1, n1, n2], [n1, n1, n2],\n            pool_method='mean', activation=act, kernel_regularizer=reg)(\n            [x1_, x2_, x3_, x4, x5, x6, x7])\n\n        x1_temp = out[0]\n        x2_temp = out[1]\n        x3_temp = out[2]\n        if dropout:\n            x1_temp = Dropout(dropout)(x1_temp, training=dropout_training)\n            x2_temp = Dropout(dropout)(x2_temp, training=dropout_training)\n            x3_temp = Dropout(dropout)(x3_temp, training=dropout_training)\n        return x1_temp, x2_temp, x3_temp\n\n    x1_ = ff(x1_)\n    x2_ = ff(x2_)\n    x3_ = ff(x3_)\n    for i in range(nblocks):\n        if i == 0:\n            has_ff = False\n        else:\n            has_ff = True\n        x1_1 = x1_\n        x2_1 = x2_\n        x3_1 = x3_\n        x1_1, x2_1, x3_1 = one_block(x1_1, x2_1, x3_1, has_ff)\n        # skip connection\n        x1_ = Add()([x1_, x1_1])\n        x2_ = Add()([x2_, x2_1])\n        x3_ = Add()([x3_, x3_1])\n    # set2set for both the atom and bond\n    node_vec = Set2Set(T=npass, n_hidden=n3, kernel_regularizer=reg)([x1_, x6])\n    edge_vec = Set2Set(T=npass, n_hidden=n3, kernel_regularizer=reg)([x2_, x7])\n    # concatenate atom, bond, and global\n    final_vec = Concatenate(axis=-1)([node_vec, edge_vec, x3_])\n    # final_vec = edge_vec\n    print(final_vec.shape)\n    if dropout:\n        final_vec = Dropout(dropout)(final_vec, training=dropout_training)\n    # final dense layers\n    final_vec = Dense(n2, activation=act, kernel_regularizer=reg)(final_vec)\n    final_vec = Dense(n3, activation=act, kernel_regularizer=reg)(final_vec)\n    if is_classification:\n        final_act = 'sigmoid'\n    else:\n        final_act = None\n    out = Dense(820, activation=final_act)(final_vec)\n    # out = final_vec\n    model = Model(inputs=[x1, x2, x3, x4, x5, x6, x7], outputs=out)\n\n    return model\n\n\ndef generate_arrays_from_file(model,batch_size=2):\n    graph_path = ['../input/all-graph-12/input_data0','../input/all-graph-12/input-data1','../input/all-graph-2/input-data2','../input/graph-data-3/input-data3']\n#     target_path = ['../input/shiyan-graph/train_targets','../input/shiyan-graph-1/train_targets','../input/shiyangraph2/train_targets','../input/shiyangraph3/train_targets']\n    while 1:\n\n        for graph_p in graph_path:\n            \n            with open(graph_p, 'rb') as fp:\n                input_data = pickle.load(fp)\n                \n            cnt = 0\n            feature_list_temp=[]\n            connection_list_tempp = []\n            global_list_temp =[]\n            index1_temp= []\n            index2_temp = []\n            temp_target = []\n            X=[]\n            for atom_features,bond_features,state_features,index1_list,index2_list,target in zip(*input_data):\n                # create Numpy arrays of input dat\n                # and labels, from each line in the file\n                feature_list_temp.append(atom_features)\n                connection_list_tempp.append(bond_features)\n                global_list_temp.append(state_features)\n                index1_temp.append(index1_list)\n                index2_temp.append(index2_list)\n                temp_target.append(target)\n                \n                \n                gnode = [0] * len(atom_features)\n                gbond = [0] * len(index1_list)\n\n       \n                \n                \n                cnt += 1\n\n                if cnt==batch_size:\n                    X.append(feature_list_temp)\n                    X.append(connection_list_tempp)\n                    X.append(global_list_temp)\n                    X.append(index1_temp)\n                    X.append(index2_temp)\n                  \n                    yield ([np.array(feature_list_temp),np.array(connection_list_tempp),np.array(global_list_temp),np.array(index1_temp),np.array(index2_temp),np.array(gnode),np.array(gbond)],np.array(temp_target))\n                    cnt = 0\n                    X=[]\n                    feature_list_temp=[] \n                    connection_list_tempp = []\n                    global_list_temp =[]\n\n                    index1_temp= []\n                    index2_temp = []\n                    temp_target = []\n                    \n# from megnet.data.molecule import MolecularGraph\n\n\npd.set_option('display.width', None)  # 设置字符显示宽度\npd.set_option('display.max_rows', None)  # 设置显示最大行\n# only reading 10% of data for debug\n# train = pd.read_csv('../../../../resources/train.csv')\n# test = pd.read_csv('../../../../resources/test.csv')\n\n##\n## Build molecules from files.xyz\n##\n\nobConversion = ob.OBConversion()\n#def read_ob_molecule(molecule_name, datadir=\"../input/champs-scalar-coupling/structures\"):\ndef read_ob_molecule(molecule_name, filename=None):\n    mol = ob.OBMol()\n    path = f\"../input/champs-scalar-coupling/structures/{molecule_name}.xyz\"\n    if not obConversion.ReadFile(mol, path):\n        raise FileNotFoundError(f\"Could not read molecule {filename}\")\n    return mol\n\n# print(read_ob_molecule('dsgdb9nsd_000001'))\n# print(train.columns)\n# print(train[['molecule_name', 'atom_index_0', 'atom_index_1', 'type',\n#        'scalar_coupling_constant']].head(10))\nn_cpu = 4\nmols={}\n\n\n    # molecule_names = np.concatenate([train.molecule_name.unique(), test.molecule_name.unique()])\n# molecule_names = train.molecule_name.unique()\n#     # xyzfiles = [Path(f'../../../../resources/structures/{f}.xyz') for f in molecule_names]\n# n = len(molecule_names)\n# i=0\n# with tqdm(total=n) as pbar:\n#     for name in molecule_names:\n#         # mol = read_ob_molecule(name)\n#             # mol1 = pybel.readfile('xyz', filename).\n#             # print(pybel.Molecule(mol))\n# #             print(train[train['molecule_name']==res[0]])\n# #             print(train[train['molecule_name'] == name][['molecule_name', 'atom_index_0', 'atom_index_1', 'type',\n# #                                                          'scalar_coupling_constant']])\n# #         mols[name] =(pybel.Molecule(mol), train[train['molecule_name'] == name])\n# #         outMDL = obConversion.WriteString(mol)\n#         mols[name] = (name, train[train['molecule_name'] == name])\n#         if(i==1):\n#              break\n#         i=i+1\n#         pbar.update()\n#\n#\n# # print(list(mols.values())[0])\n# mols = list(mols.values())\n# import pickle\n#\n# with open('train_outfile', 'wb') as fp:\n#     pickle.dump(mols, fp)\n\n# with open ('../input/new-train-data-set/train_outfile', 'rb') as fp:\n#     amols = pickle.load(fp)\n# mols=amols[:20]\n# print(mols[:2])\n# test_mols = amols[25:30]\n\n# with open ('train_graphs', 'rb') as fp:\n#     mols = pickle.load(fp)\n\n# for i in mols:\n#     print(i)\n# jsObj = json.dumps(mols)\n#\n# fileObject = open('jsonFile.json', 'w')\n# fileObject.write(jsObj)\n# fileObject.close()\n\n\n\n# targets = [[0,0,0,-2.8],[0,0,0,2.4]]\ndef mse(orig, preds):\n    # Mask values for which no scalar coupling exists\n    # orig=orig[:114]\n    # preds=preds[:114]\n    mask = tf.where(tf.equal(orig, 0), orig, tf.ones_like(orig))\n\n    nums = tf.boolean_mask(orig, mask)\n    preds = tf.boolean_mask(preds, mask)\n\n    reconstruction_error = tf.reduce_mean(tf.square(tf.subtract(nums, preds)))\n\n    return reconstruction_error\n\ndef log_mae(orig, preds):\n    # orig = orig[114:]\n    # preds= preds[114:]\n    # Mask values for which no scalar coupling exists\n    mask = tf.where(tf.equal(orig, 0), orig, tf.ones_like(orig))\n\n    nums = tf.boolean_mask(orig, mask)\n    preds = tf.boolean_mask(preds, mask)\n\n    reconstruction_error = tf.math.log(tf.reduce_mean(tf.abs(tf.subtract(nums, preds))))\n\n    return reconstruction_error\n\n\n\n# Define some callbacks, the initial learning rate and the optimizer\n# learning_rate = 0.001\ndef step_decay(epoch):\n    initial_lrate = 0.01\n    drop = 0.1\n    epochs_drop = 20.0\n    lrate = initial_lrate * np.power(drop,\n           np.floor((epoch)/epochs_drop))\n    tf.print(\"Learning rate: \", lrate)\n    return lrate\n\n# from keras.callbacks import LearningRateScheduler\n# lrate = LearningRateScheduler(step_decay)\n\n\nimport keras.backend as K\nfrom keras.callbacks import LearningRateScheduler\nfrom keras.callbacks import ReduceLROnPlateau\ndef scheduler(epoch):\n    # 每隔100个epoch，学习率减小为原来的1/10\n    if epoch % 100 == 0 and epoch != 0:\n        lr = K.get_value(model.optimizer.lr)\n        K.set_value(model.optimizer.lr, lr * 0.1)\n        print(\"lr changed to {}\".format(lr * 0.1))\n    return K.get_value(model.optimizer.lr)\n  \n# lrate = LearningRateScheduler(step_decay)\n\nstop_early = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience = 20, restore_best_weights=True)\n\nlrate  =  ReduceLROnPlateau(monitor='val_loss', factor=0.1,patience=10, min_lr=0.000001, verbose = 1)\n\n# opt = tf.optimizers.Adam(learning_rate=learning_rate)\n\n\nmodel = MEGNetModel(35, 2, 27, nblocks=1, lr=1e-2,\n                    n1=4, n2=32, n3=16, npass=1, ntarget=1,graph_converter=MolecularGraph())\n# model = MEGNetModel.from_file(\"./train_model_1\")\n# mpnn.compile(opt, log_mae, metrics = [mae, log_mse])\n\n\n\n# model.summary()\n#\n# model.load_weights('../input/shiyan-data-graph/val_mae_00009_0.119766.hdf5')\n\nmodel.load_weights('../input/tune-model/weights.best.hdf5')\nmodel.train_coupling(mols,epochs=200, verbose=1,batch_size=1)\n# model.save(\"train_model_1\")\n\n# ###test data begin#######\n#\n# test_mols={}\n# #     molecule_names = np.concatenate([train.molecule_name.unique(), test.molecule_name.unique()])\n# molecule_names = test.molecule_name.unique()\n#     # xyzfiles = [Path(f'../../../../resources/structures/{f}.xyz') for f in molecule_names]\n# n = len(molecule_names)\n# i=0\n# with tqdm(total=n) as pbar:\n#     for name in molecule_names:\n#         # mol = read_ob_molecule(name)\n#             # mol1 = pybel.readfile('xyz', filename).\n#             # print(pybel.Molecule(mol))\n# #             print(train[train['molecule_name']==res[0]])\n#\n#         test_mols[name] =(name, test[test['molecule_name'] == name])\n#         if(i==10):\n#             break\n#         i=i+1\n#         pbar.update()\n#\n#\n# # print(list(mols.values())[0])\n# test_mols = list(test_mols.values())\n#\n# test_graphs,  test_datas= model.get_all_graphs_targets_coupling(test_mols,target_convert=False)\n#\n\n\n# with open('../input/test-dat/test_graphs', 'rb') as fp:\n#     test_graphs = pickle.load(fp)\n\n# with open('../input/test-data/test_datas', 'rb') as fp:\n#     test_datas = pickle.load(fp)\n# for graph,test_data in zip(test_graphs,test_datas):\n#     # print(graph)\n#     # print(test_data)\n#     tt = model.predict_graph(graph)\n#     length = len(test_data['index1'])\n#     couple_constant_data = tt[:length]\n#     print(couple_constant_data)\n#     couple_constant_data= couple_constant_data*(max_value-min_value)+min_value\n#     print(couple_constant_data)\n#     df1 = pd.DataFrame({'index1':test_data['index1'],\n#                     'index2':test_data['index2'],\n#                     'couple_constant_data':couple_constant_data\n#                     })\n#     df2=test_data['coupling_attributes']\n#     df2.rename(columns={'atom_index_0':'index1', 'atom_index_1':'index2'}, inplace = True)\n#     print(df1)\n#     print(df2)\n#     data=[]\n#     for index,row in df2.iterrows():\n\n#         values1 = df1[(df1['index1'] == row['index1']) &(df1['index2'] == row['index2'])]['couple_constant_data'].tolist()[0]\n#         values2 = df1[(df1['index1'] == row['index2']) & (df1['index2'] == row['index1'])]['couple_constant_data'].tolist()[0]\n\n#         # coupling_arr0 = coupling_attr[(coupling_attr['atom_index_0'] == bid) & (coupling_attr['atom_index_1'] == eid)]\n#         # one = values1['couple_constant_data'].tolist()[0]\n#         # two = values2['couple_constant_data'].tolist()[0]\n\n#         data.append((values1+values2)/2)\n\n#     print(\"################\")\n#     df2['couple_data_value'] = data\n\n#     print(df2)\n#     # df= pd.merge(df2,df1)\n\n#     # print(df)\n#     df2.to_csv('test_result.csv',index = False,header=0,mode='a')\n# #     print('fininsed####################')\n# # #     # print(df[['molecule_name','index1', 'index2', 'couple_constant_data',\n# #     #    'type', 'scalar_coupling_constant']])\n# print('fininsed####################')\n\n# predicted_atom = [model.predict_structure(x) for x in structures]\n# print(predicted_atom)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import keras\nprint(tensorflow.__version__)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# import pickle\n\n# # with open('../input/shiyan-graph/train_graphs', 'rb') as fp:\n# #     data_graph = pickle.load(fp)\n\n# # with open('../input/shiyan-graph/train_targets', 'rb') as fp:\n# #     data_target = pickle.load(fp)\n\n# with open('../input/shiyan-graph-1/train_graphs', 'rb') as fp:\n#     data_graphs = pickle.load(fp)\n\n# with open('../input/shiyan-graph-1/train_targets', 'rb') as fp:\n#     data_targets = pickle.load(fp)\n    \n# # data_graph=None\n# # data_targets=None\n# # #         print(len(data_targets))\n# print(data_targets[0])\n# for target in  data_targets:\n#     target =  (target-min_value)/(max_value -min_value) \n# data_targets[:2]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip show keras","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls ../input/","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}