Upload graph_decoder/diffusion_model.py with huggingface_hub

graph_decoder/diffusion_model.py

@@ -0,0 +1,422 @@
+# Copyright 2024 the Llamole team.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import os 
+import yaml
+import json
+
+import torch
+import torch.nn as nn
+import torch.nn.functional as F
+
+from . import diffusion_utils as utils
+from .molecule_utils import graph_to_smiles, check_valid
+from .transformer import Transformer
+from .visualize_utils import MolecularVisualization
+
+class GraphDiT(nn.Module):
+    def __init__(
+        self,
+        model_config_path,
+        data_info_path,
+        model_dtype,
+    ):
+        super().__init__()
+
+        dm_cfg, data_info = utils.load_config(model_config_path, data_info_path)
+
+        input_dims = data_info.input_dims
+        output_dims = data_info.output_dims
+        nodes_dist = data_info.nodes_dist
+        active_index = data_info.active_index
+
+        self.model_config = dm_cfg
+        self.data_info = data_info
+        self.T = dm_cfg.diffusion_steps
+        self.Xdim = input_dims["X"]
+        self.Edim = input_dims["E"]
+        self.ydim = input_dims["y"]
+        self.Xdim_output = output_dims["X"]
+        self.Edim_output = output_dims["E"]
+        self.ydim_output = output_dims["y"]
+        self.node_dist = nodes_dist
+        self.active_index = active_index
+        self.max_n_nodes = data_info.max_n_nodes
+        self.atom_decoder = data_info.atom_decoder
+        self.hidden_size = dm_cfg.hidden_size
+        self.mol_visualizer = MolecularVisualization(self.atom_decoder)
+
+        self.denoiser = Transformer(
+            max_n_nodes=self.max_n_nodes,
+            hidden_size=dm_cfg.hidden_size,
+            depth=dm_cfg.depth,
+            num_heads=dm_cfg.num_heads,
+            mlp_ratio=dm_cfg.mlp_ratio,
+            drop_condition=dm_cfg.drop_condition,
+            Xdim=self.Xdim,
+            Edim=self.Edim,
+            ydim=self.ydim,
+        )
+        self.model_dtype = model_dtype
+        # self.device = next(self.denoiser.parameters()).device
+
+        # model_params = torch.load(model_params_path, map_location='cpu')
+        # self.denoiser.load_state_dict(model_params)
+        
+        self.noise_schedule = utils.PredefinedNoiseScheduleDiscrete(
+            dm_cfg.diffusion_noise_schedule, timesteps=dm_cfg.diffusion_steps
+        )
+        x_marginals = data_info.node_types.to(self.model_dtype) / torch.sum(
+            data_info.node_types.to(self.model_dtype)
+        )
+        e_marginals = data_info.edge_types.to(self.model_dtype) / torch.sum(
+            data_info.edge_types.to(self.model_dtype)
+        )
+        x_marginals = x_marginals / x_marginals.sum()
+        e_marginals = e_marginals / e_marginals.sum()
+
+        xe_conditions = data_info.transition_E.to(self.model_dtype)
+        xe_conditions = xe_conditions[self.active_index][:, self.active_index]
+
+        xe_conditions = xe_conditions.sum(dim=1)
+        ex_conditions = xe_conditions.t()
+        xe_conditions = xe_conditions / xe_conditions.sum(dim=-1, keepdim=True)
+        ex_conditions = ex_conditions / ex_conditions.sum(dim=-1, keepdim=True)
+
+        self.transition_model = utils.MarginalTransition(
+            x_marginals=x_marginals,
+            e_marginals=e_marginals,
+            xe_conditions=xe_conditions,
+            ex_conditions=ex_conditions,
+            y_classes=self.ydim_output,
+            n_nodes=self.max_n_nodes,
+        )
+        self.limit_dist = utils.PlaceHolder(X=x_marginals, E=e_marginals, y=None)
+
+    # def to(self, *args, **kwargs):
+    #     self = super().to(*args, **kwargs)
+    #     self.model_dtype = next(self.denoiser.parameters()).dtype
+    #     return self
+
+    def init_model(self, model_dir, verbose=False):
+        model_file = os.path.join(model_dir, 'model.pt')
+        if os.path.exists(model_file):
+            self.denoiser.load_state_dict(torch.load(model_file, map_location='cpu', weights_only=True))
+        else:
+            raise FileNotFoundError(f"Model file not found: {model_file}")
+        
+        if verbose:
+            print('GraphDiT Denoiser Model initialized.')
+            print('Denoiser model:\n', self.denoiser)
+
+    def save_pretrained(self, output_dir):
+        if not os.path.exists(output_dir):
+            os.makedirs(output_dir)
+        
+        # Save model
+        model_path = os.path.join(output_dir, 'model.pt')
+        torch.save(self.denoiser.state_dict(), model_path)
+        
+        # Save model config
+        config_path = os.path.join(output_dir, 'model_config.yaml')
+        with open(config_path, 'w') as f:
+            yaml.dump(vars(self.model_config), f)
+        
+        # Save data info
+        data_info_path = os.path.join(output_dir, 'data.meta.json')
+        data_info_dict = {
+            "active_atoms": self.data_info.active_atoms,
+            "max_node": self.data_info.max_n_nodes,
+            "n_atoms_per_mol_dist": self.data_info.n_nodes.tolist(),
+            "bond_type_dist": self.data_info.edge_types.tolist(),
+            "transition_E": self.data_info.transition_E.tolist(),
+            "atom_type_dist": self.data_info.node_types.tolist(),
+            "valencies": self.data_info.valency_distribution.tolist()
+        }
+        with open(data_info_path, 'w') as f:
+            json.dump(data_info_dict, f, indent=2)
+        
+        print('GraphDiT Model and configurations saved to:', output_dir)
+
+    def disable_grads(self):
+        self.denoiser.disable_grads()
+    
+    def forward(
+        self, x, edge_index, edge_attr, graph_batch, properties, no_label_index
+    ):
+        raise ValueError('Not Implement')
+
+    def _forward(self, noisy_data, unconditioned=False):
+        noisy_x, noisy_e, properties = (
+            noisy_data["X_t"].to(self.model_dtype),
+            noisy_data["E_t"].to(self.model_dtype),
+            noisy_data["y_t"].to(self.model_dtype).clone(),
+        )
+        node_mask, timestep = (
+            noisy_data["node_mask"],
+            noisy_data["t"],
+        )
+        
+        pred = self.denoiser(
+            noisy_x,
+            noisy_e,
+            node_mask,
+            properties,
+            timestep,
+            unconditioned=unconditioned,
+        )
+        return pred
+
+    def apply_noise(self, X, E, y, node_mask):
+        """Sample noise and apply it to the data."""
+
+        # Sample a timestep t.
+        # When evaluating, the loss for t=0 is computed separately
+        lowest_t = 0 if self.training else 1
+        t_int = torch.randint(
+            lowest_t, self.T + 1, size=(X.size(0), 1), device=X.device
+        ).to(
+            self.model_dtype
+        )  # (bs, 1)
+        s_int = t_int - 1
+
+        t_float = t_int / self.T
+        s_float = s_int / self.T
+
+        # beta_t and alpha_s_bar are used for denoising/loss computation
+        beta_t = self.noise_schedule(t_normalized=t_float)  # (bs, 1)
+        alpha_s_bar = self.noise_schedule.get_alpha_bar(t_normalized=s_float)  # (bs, 1)
+        alpha_t_bar = self.noise_schedule.get_alpha_bar(t_normalized=t_float)  # (bs, 1)
+
+        Qtb = self.transition_model.get_Qt_bar(
+            alpha_t_bar, X.device
+        )  # (bs, dx_in, dx_out), (bs, de_in, de_out)
+
+        bs, n, d = X.shape
+        X_all = torch.cat([X, E.reshape(bs, n, -1)], dim=-1)
+        prob_all = X_all @ Qtb.X
+        probX = prob_all[:, :, : self.Xdim_output]
+        probE = prob_all[:, :, self.Xdim_output :].reshape(bs, n, n, -1)
+
+        sampled_t = utils.sample_discrete_features(
+            probX=probX, probE=probE, node_mask=node_mask
+        )
+
+        X_t = F.one_hot(sampled_t.X, num_classes=self.Xdim_output)
+        E_t = F.one_hot(sampled_t.E, num_classes=self.Edim_output)
+        assert (X.shape == X_t.shape) and (E.shape == E_t.shape)
+
+        y_t = y
+        z_t = utils.PlaceHolder(X=X_t, E=E_t, y=y_t).type_as(X_t).mask(node_mask)
+
+        noisy_data = {
+            "t_int": t_int,
+            "t": t_float,
+            "beta_t": beta_t,
+            "alpha_s_bar": alpha_s_bar,
+            "alpha_t_bar": alpha_t_bar,
+            "X_t": z_t.X,
+            "E_t": z_t.E,
+            "y_t": z_t.y,
+            "node_mask": node_mask,
+        }
+        return noisy_data
+
+    @torch.no_grad()
+    def generate(
+        self,
+        properties,
+        device,
+        guide_scale=1.,
+        num_nodes=None,
+        number_chain_steps=50,
+    ):
+        properties = [float('nan') if x is None else x for x in properties]
+        properties = torch.tensor(properties, dtype=torch.float).reshape(1, -1).to(device)
+        batch_size = properties.size(0)
+        assert batch_size == 1
+        # print('self.denoiser.dtype', self.model_dtype)
+        if num_nodes is None:
+            num_nodes = self.node_dist.sample_n(batch_size, device)
+        else:
+            num_nodes = torch.LongTensor([num_nodes]).to(device)
+
+        arange = (
+            torch.arange(self.max_n_nodes, device=device)
+            .unsqueeze(0)
+            .expand(batch_size, -1)
+        )
+        node_mask = arange < num_nodes.unsqueeze(1)
+
+        z_T = utils.sample_discrete_feature_noise(
+            limit_dist=self.limit_dist, node_mask=node_mask
+        )
+        X, E = z_T.X, z_T.E
+
+        assert (E == torch.transpose(E, 1, 2)).all()
+
+        if number_chain_steps > 0:
+            chain_X_size = torch.Size((number_chain_steps, X.size(1)))
+            chain_E_size = torch.Size((number_chain_steps, E.size(1), E.size(2)))
+            chain_X = torch.zeros(chain_X_size)
+            chain_E = torch.zeros(chain_E_size)
+
+        # Iteratively sample p(z_s | z_t) for t = 1, ..., T, with s = t - 1.
+        y = properties
+        for s_int in reversed(range(0, self.T)):
+            s_array = s_int * torch.ones((batch_size, 1)).type_as(y)
+            t_array = s_array + 1
+            s_norm = s_array / self.T
+            t_norm = t_array / self.T
+
+            # Sample z_s
+            sampled_s, discrete_sampled_s = self.sample_p_zs_given_zt(
+                s_norm, t_norm, X, E, y, node_mask, guide_scale, device
+            )
+            X, E, y = sampled_s.X, sampled_s.E, sampled_s.y
+            
+            if number_chain_steps > 0:
+                # Save the first keep_chain graphs
+                write_index = (s_int * number_chain_steps) // self.T
+                chain_X[write_index] = discrete_sampled_s.X[:1]
+                chain_E[write_index] = discrete_sampled_s.E[:1]
+
+        # Sample
+        sampled_s = sampled_s.mask(node_mask, collapse=True)
+        X, E, y = sampled_s.X, sampled_s.E, sampled_s.y
+
+        molecule_list = []
+        n = num_nodes[0]
+        atom_types = X[0, :n].cpu()
+        edge_types = E[0, :n, :n].cpu()
+        molecule_list.append([atom_types, edge_types])
+        smiles = graph_to_smiles(molecule_list, self.atom_decoder)[0]
+
+        # Visualize Chains
+        if number_chain_steps > 0:
+            final_X_chain = X[:1]
+            final_E_chain = E[:1]
+
+            chain_X[0] = final_X_chain                  # Overwrite last frame with the resulting X, E
+            chain_E[0] = final_E_chain
+
+            chain_X = utils.reverse_tensor(chain_X)
+            chain_E = utils.reverse_tensor(chain_E)
+
+            # Repeat last frame to see final sample better
+            chain_X = torch.cat([chain_X, chain_X[-1:].repeat(10, 1)], dim=0)
+            chain_E = torch.cat([chain_E, chain_E[-1:].repeat(10, 1, 1)], dim=0)
+            mol_img_list = self.mol_visualizer.visualize_chain(chain_X.numpy(), chain_E.numpy())
+        else:
+            mol_img_list = []
+
+        return smiles, mol_img_list
+
+    def check_valid(self, smiles):
+        return check_valid(smiles)
+    
+    def sample_p_zs_given_zt(
+        self, s, t, X_t, E_t, properties, node_mask, guide_scale, device
+    ):
+        """Samples from zs ~ p(zs | zt). Only used during sampling.
+        if last_step, return the graph prediction as well"""
+        bs, n, _ = X_t.shape
+        beta_t = self.noise_schedule(t_normalized=t)  # (bs, 1)
+        alpha_s_bar = self.noise_schedule.get_alpha_bar(t_normalized=s)
+        alpha_t_bar = self.noise_schedule.get_alpha_bar(t_normalized=t)
+
+        # Neural net predictions
+        noisy_data = {
+            "X_t": X_t,
+            "E_t": E_t,
+            "y_t": properties,
+            "t": t,
+            "node_mask": node_mask,
+        }
+
+        def get_prob(noisy_data, unconditioned=False):
+            pred = self._forward(noisy_data, unconditioned=unconditioned)
+
+            # Normalize predictions
+            pred_X = F.softmax(pred.X, dim=-1)  # bs, n, d0
+            pred_E = F.softmax(pred.E, dim=-1)  # bs, n, n, d0
+
+            # Retrieve transitions matrix
+            Qtb = self.transition_model.get_Qt_bar(alpha_t_bar, device)
+            Qsb = self.transition_model.get_Qt_bar(alpha_s_bar, device)
+            Qt = self.transition_model.get_Qt(beta_t, device)
+
+            Xt_all = torch.cat([X_t, E_t.reshape(bs, n, -1)], dim=-1)
+            predX_all = torch.cat([pred_X, pred_E.reshape(bs, n, -1)], dim=-1)
+
+            unnormalized_probX_all = utils.reverse_diffusion(
+                predX_0=predX_all, X_t=Xt_all, Qt=Qt.X, Qsb=Qsb.X, Qtb=Qtb.X
+            )
+
+            unnormalized_prob_X = unnormalized_probX_all[:, :, : self.Xdim_output]
+            unnormalized_prob_E = unnormalized_probX_all[
+                :, :, self.Xdim_output :
+            ].reshape(bs, n * n, -1)
+
+            unnormalized_prob_X[torch.sum(unnormalized_prob_X, dim=-1) == 0] = 1e-5
+            unnormalized_prob_E[torch.sum(unnormalized_prob_E, dim=-1) == 0] = 1e-5
+
+            prob_X = unnormalized_prob_X / torch.sum(
+                unnormalized_prob_X, dim=-1, keepdim=True
+            )  # bs, n, d_t-1
+            prob_E = unnormalized_prob_E / torch.sum(
+                unnormalized_prob_E, dim=-1, keepdim=True
+            )  # bs, n, d_t-1
+            prob_E = prob_E.reshape(bs, n, n, pred_E.shape[-1])
+
+            return prob_X, prob_E
+
+        prob_X, prob_E = get_prob(noisy_data)
+
+        ### Guidance
+        if guide_scale != 1:
+            uncon_prob_X, uncon_prob_E = get_prob(
+                noisy_data, unconditioned=True
+            )
+            prob_X = (
+                uncon_prob_X
+                * (prob_X / uncon_prob_X.clamp_min(1e-5)) ** guide_scale
+            )
+            prob_E = (
+                uncon_prob_E
+                * (prob_E / uncon_prob_E.clamp_min(1e-5)) ** guide_scale
+            )
+            prob_X = prob_X / prob_X.sum(dim=-1, keepdim=True).clamp_min(1e-5)
+            prob_E = prob_E / prob_E.sum(dim=-1, keepdim=True).clamp_min(1e-5)
+
+        # assert ((prob_X.sum(dim=-1) - 1).abs() < 1e-3).all()
+        # assert ((prob_E.sum(dim=-1) - 1).abs() < 1e-3).all()
+
+        sampled_s = utils.sample_discrete_features(
+            prob_X, prob_E, node_mask=node_mask, step=s[0, 0].item()
+        )
+
+        X_s = F.one_hot(sampled_s.X, num_classes=self.Xdim_output).to(self.model_dtype)
+        E_s = F.one_hot(sampled_s.E, num_classes=self.Edim_output).to(self.model_dtype)
+
+        assert (E_s == torch.transpose(E_s, 1, 2)).all()
+        assert (X_t.shape == X_s.shape) and (E_t.shape == E_s.shape)
+
+        out_one_hot = utils.PlaceHolder(X=X_s, E=E_s, y=properties)
+        out_discrete = utils.PlaceHolder(X=X_s, E=E_s, y=properties)
+
+        return out_one_hot.mask(node_mask).type_as(properties), out_discrete.mask(
+            node_mask, collapse=True
+        ).type_as(properties)
+