SurfDock / model /models /mdn_score_model_v6.py
Upload folder using huggingface_hub (part 4)
96f168d verified about 2 months ago
31.4 kB
import math

from e3nn import o3
import torch
from torch import nn
from torch.nn import functional as F
from torch_cluster import radius, radius_graph
from torch_scatter import scatter, scatter_mean,scatter_add
import numpy as np
from e3nn.nn import BatchNorm
from torch_geometric.utils import  to_dense_batch
from utils import so3, torus
from datasets.process_mols import lig_feature_dims, rec_residue_feature_dims
from utils.mdn_utils import compute_euclidean_distances_matrix,compute_euclidean_distances_matrix_TopN
from utils.training_mdn import mdn_loss_fn,calculate_probablity
""""
Version 6: this version use surface node to replace rec node to cal mdn ,correspond to surface_score_model version3
"""
class AtomEncoder(torch.nn.Module):

    def __init__(self, emb_dim, feature_dims, sigma_embed_dim, lm_embedding_type= None):
        # first element of feature_dims tuple is a list with the lenght of each categorical feature and the second is the number of scalar features
        super(AtomEncoder, self).__init__()
        self.atom_embedding_list = torch.nn.ModuleList()
        self.num_categorical_features = len(feature_dims[0])
        self.num_scalar_features = feature_dims[1] #+ sigma_embed_dim
        self.lm_embedding_type = lm_embedding_type
        for i, dim in enumerate(feature_dims[0]):
            emb = torch.nn.Embedding(dim, emb_dim)
            torch.nn.init.xavier_uniform_(emb.weight.data)
            self.atom_embedding_list.append(emb)

        if self.num_scalar_features > 0:
            self.linear = torch.nn.Linear(self.num_scalar_features, emb_dim)
        if self.lm_embedding_type is not None:
            if self.lm_embedding_type == 'esm':
                self.lm_embedding_dim = 1280
            else: raise ValueError('LM Embedding type was not correctly determined. LM embedding type: ', self.lm_embedding_type)
            self.lm_embedding_layer = torch.nn.Linear(self.lm_embedding_dim + emb_dim, emb_dim)

    def forward(self, x):
        x_embedding = 0
        if self.lm_embedding_type is not None:
            assert x.shape[1] == self.num_categorical_features + self.num_scalar_features + self.lm_embedding_dim,f'{x.shape[1]}=={self.num_categorical_features} + {self.num_scalar_features} + {self.lm_embedding_dim}'
        else:
            assert x.shape[1] == self.num_categorical_features + self.num_scalar_features
        for i in range(self.num_categorical_features):
            x_embedding += self.atom_embedding_list[i](x[:, i].long())

        if self.num_scalar_features > 0:
            x_embedding += self.linear(x[:, self.num_categorical_features:self.num_categorical_features + self.num_scalar_features])
        if self.lm_embedding_type is not None:
            x_embedding = self.lm_embedding_layer(torch.cat([x_embedding, x[:, -self.lm_embedding_dim:]], axis=1))
        return x_embedding


class TensorProductConvLayer(torch.nn.Module):
    def __init__(self, in_irreps, sh_irreps, out_irreps, n_edge_features, residual=True, batch_norm=True, dropout=0.0,
                 hidden_features=None):
        super(TensorProductConvLayer, self).__init__()
        self.in_irreps = in_irreps
        self.out_irreps = out_irreps
        self.sh_irreps = sh_irreps
        self.residual = residual
        if hidden_features is None:
            hidden_features = n_edge_features

        self.tp = tp = o3.FullyConnectedTensorProduct(in_irreps, sh_irreps, out_irreps, shared_weights=False)

        self.fc = nn.Sequential(
            nn.Linear(n_edge_features, hidden_features),
            nn.ReLU(),
            nn.Dropout(dropout),
            nn.Linear(hidden_features, tp.weight_numel)
        )
        self.batch_norm = BatchNorm(out_irreps) if batch_norm else None

    def forward(self, node_attr, edge_index, edge_attr, edge_sh, out_nodes=None, reduce='mean'):

        edge_src, edge_dst = edge_index
        tp = self.tp(node_attr[edge_dst], edge_sh, self.fc(edge_attr))

        out_nodes = out_nodes or node_attr.shape[0]
        out = scatter(tp, edge_src, dim=0, dim_size=out_nodes, reduce=reduce)

        if self.residual:
            padded = F.pad(node_attr, (0, out.shape[-1] - node_attr.shape[-1]))
            out = out + padded

        if self.batch_norm:

            out = self.batch_norm(out)
        return out


class TensorProductScoreModelV6(torch.nn.Module):
    def __init__(self, args,t_to_sigma, device, timestep_emb_func, in_lig_edge_features=10, in_rec_edge_features = 5,sigma_embed_dim=32, sh_lmax=2,
                 ns=16, nv=4, num_conv_layers=2, lig_max_radius=5, rec_max_radius=30, cross_max_distance=250,
                 center_max_distance=30, distance_embed_dim=32, cross_distance_embed_dim=32, no_torsion=False,
                 scale_by_sigma=True, use_second_order_repr=False, batch_norm=True,
                 dynamic_max_cross=False, dropout=0.0, lm_embedding_type=None, mdn_mode=True,
                 mdn_dropout=0, mdn_no_batchnorm=False,n_gaussians = 20):
        super(TensorProductScoreModelV6, self).__init__()
        self.args = args
        self.t_to_sigma = t_to_sigma
        self.in_lig_edge_features = in_lig_edge_features
        self.sigma_embed_dim = sigma_embed_dim
        self.lig_max_radius = lig_max_radius
        self.rec_max_radius = rec_max_radius
        self.cross_max_distance = cross_max_distance
        self.dynamic_max_cross = dynamic_max_cross
        self.center_max_distance = center_max_distance
        self.distance_embed_dim = distance_embed_dim
        self.cross_distance_embed_dim = cross_distance_embed_dim
        self.sh_irreps = o3.Irreps.spherical_harmonics(lmax=sh_lmax)
        self.ns, self.nv = ns, nv
        self.scale_by_sigma = scale_by_sigma
        self.device = device
        self.no_torsion = no_torsion
        self.timestep_emb_func = timestep_emb_func
        self.mdn_mode = mdn_mode
        self.num_conv_layers = num_conv_layers

        self.lig_node_embedding = AtomEncoder(emb_dim=ns, feature_dims=lig_feature_dims, sigma_embed_dim=sigma_embed_dim)
        self.lig_edge_embedding = nn.Sequential(nn.Linear(in_lig_edge_features + distance_embed_dim, ns),nn.ReLU(), nn.Dropout(dropout),nn.Linear(ns, ns))

        self.rec_node_embedding = AtomEncoder(emb_dim=ns, feature_dims=rec_residue_feature_dims, sigma_embed_dim=sigma_embed_dim, lm_embedding_type=lm_embedding_type)
        self.rec_edge_embedding = nn.Sequential(nn.Linear(in_rec_edge_features + distance_embed_dim, ns), nn.ReLU(), nn.Dropout(dropout),nn.Linear(ns, ns))

        # self.cross_edge_embedding = nn.Sequential(nn.Linear(sigma_embed_dim , ns), nn.ReLU(), nn.Dropout(dropout),nn.Linear(ns, ns))
        self.surface_node_embedding = AtomEncoder(emb_dim=ns, feature_dims=[[],4], sigma_embed_dim=sigma_embed_dim)
        self.surface_edge_embedding = nn.Sequential(nn.Linear(3 + distance_embed_dim, ns), nn.ReLU(), nn.Dropout(dropout),nn.Linear(ns, ns))
        
        self.surface_rec_cross_edge_embedding = nn.Sequential(nn.Linear(cross_distance_embed_dim, ns), nn.ReLU(), nn.Dropout(dropout),nn.Linear(ns, ns))
        
        self.lig_distance_expansion = GaussianSmearing(0.0, lig_max_radius, distance_embed_dim)
        self.rec_distance_expansion = GaussianSmearing(0.0, rec_max_radius, distance_embed_dim)
        self.cross_distance_expansion = GaussianSmearing(0.0, cross_max_distance, cross_distance_embed_dim)
        self.surface_distance_expansion = GaussianSmearing(0.0, rec_max_radius, distance_embed_dim)

        if use_second_order_repr:
            irrep_seq = [
                f'{ns}x0e',
                f'{ns}x0e + {nv}x1o + {nv}x2e',
                f'{ns}x0e + {nv}x1o + {nv}x2e + {nv}x1e + {nv}x2o',
                f'{ns}x0e + {nv}x1o + {nv}x2e + {nv}x1e + {nv}x2o + {ns}x0o'
            ]
        else:
            irrep_seq = [
                f'{ns}x0e',
                f'{ns}x0e + {nv}x1o',
                f'{ns}x0e + {nv}x1o + {nv}x1e',
                f'{ns}x0e + {nv}x1o + {nv}x1e + {ns}x0o'
            ]
        lig_conv_layers= []
        # surface modules
        surface_conv_layers=[]
        residue_to_surface_conv_layers = []
        lig_conv_layers, rec_conv_layers = [], []
        for i in range(num_conv_layers):
            in_irreps = irrep_seq[min(i, len(irrep_seq) - 1)]
            out_irreps = irrep_seq[min(i + 1, len(irrep_seq) - 1)]
            parameters = {
                'in_irreps': in_irreps,
                'sh_irreps': self.sh_irreps,
                'out_irreps': out_irreps,
                'n_edge_features': 3 * ns,
                'hidden_features': 3 * ns,
                'residual': False,
                'batch_norm': batch_norm,
                'dropout': dropout
            }
            if i ==0:
                residue_to_surface_conv_layers.append(TensorProductConvLayer(** {
                'in_irreps': f'{ns}x0e + {nv}x1o + {nv}x1e + {ns}x0o',
                'sh_irreps': self.sh_irreps,
                'out_irreps':  in_irreps,
                'n_edge_features': 3 * ns,
                'hidden_features': 3 * ns,
                'residual': False,
                'batch_norm': batch_norm,
                'dropout': dropout
            }))
                rec_conv_layers.append(TensorProductConvLayer(** {
                'in_irreps': in_irreps,
                'sh_irreps': self.sh_irreps,
                'out_irreps': f'{ns}x0e + {nv}x1o + {nv}x1e + {ns}x0o',
                'n_edge_features': 3 * ns,
                'hidden_features': 3 * ns,
                'residual': False,
                'batch_norm': batch_norm,
                'dropout': dropout
            }))
            lig_layer = TensorProductConvLayer(**parameters)
            lig_conv_layers.append(lig_layer)


            if i != num_conv_layers - 1:
                surface_layer = TensorProductConvLayer(**parameters)
                surface_conv_layers.append(surface_layer)

        self.residue_to_surface_conv_layers = nn.ModuleList(residue_to_surface_conv_layers)
        self.surface_conv_layers = nn.ModuleList(surface_conv_layers)
        self.lig_conv_layers = nn.ModuleList(lig_conv_layers)

        self.rec_conv_layers = nn.ModuleList(rec_conv_layers)


        if self.mdn_mode:
            # only mdn model without atom type and bond type predict 
            if self.num_conv_layers >=2:
                mdn_hidden_dim = 4 * ns
            else:
                mdn_hidden_dim = 2 * ns
            self.MLP = nn.Sequential(nn.Linear(mdn_hidden_dim, mdn_hidden_dim), nn.BatchNorm1d(mdn_hidden_dim), nn.ELU(), nn.Dropout(p= mdn_dropout)) 
            self.z_pi = nn.Linear(mdn_hidden_dim, n_gaussians)
            self.z_sigma = nn.Linear(mdn_hidden_dim, n_gaussians)
            self.z_mu = nn.Linear(mdn_hidden_dim, n_gaussians)
            # ligand distance prediction
            if self.args.ligand_distance_prediction:
                self.ligand_MLP = nn.Sequential(nn.Linear(mdn_hidden_dim, mdn_hidden_dim), nn.BatchNorm1d(mdn_hidden_dim), nn.ELU(), nn.Dropout(p= mdn_dropout)) 
                self.ligand_z_pi = nn.Linear(mdn_hidden_dim, n_gaussians)
                self.ligand_z_sigma = nn.Linear(mdn_hidden_dim, n_gaussians)
                self.ligand_z_mu = nn.Linear(mdn_hidden_dim, n_gaussians)
            if self.args.atom_type_prediction:
                self.atom_types = nn.Sequential(nn.Linear(mdn_hidden_dim//2, 2*mdn_hidden_dim), nn.BatchNorm1d(2*mdn_hidden_dim), nn.ELU(), nn.Dropout(p= mdn_dropout),nn.Linear(2*mdn_hidden_dim, 119)) 
            if self.args.bond_type_prediction:
                self.bond_types = nn.Sequential(nn.Linear(mdn_hidden_dim, 2*mdn_hidden_dim), nn.BatchNorm1d(2*mdn_hidden_dim), nn.ELU(), nn.Dropout(p= mdn_dropout),nn.Linear(2*mdn_hidden_dim, 4))
            if self.args.residue_type_prediction:
                self.residue_types = nn.Sequential(nn.Linear(mdn_hidden_dim//2, 2*mdn_hidden_dim), nn.BatchNorm1d(2*mdn_hidden_dim), nn.ELU(), nn.Dropout(p= mdn_dropout),nn.Linear(2*mdn_hidden_dim, 38)) 

        else:
            # center of mass translation and rotation components
            self.center_distance_expansion = GaussianSmearing(0.0, center_max_distance, distance_embed_dim)
            self.center_edge_embedding = nn.Sequential(
                nn.Linear(distance_embed_dim + sigma_embed_dim, ns),
                nn.ReLU(),
                nn.Dropout(dropout),
                nn.Linear(ns, ns)
            )

            self.final_conv = TensorProductConvLayer(
                in_irreps=self.lig_conv_layers[-1].out_irreps,
                sh_irreps=self.sh_irreps,
                out_irreps=f'2x1o + 2x1e',
                n_edge_features=2 * ns,
                residual=False,
                dropout=dropout,
                batch_norm=batch_norm
            )
            self.tr_final_layer = nn.Sequential(nn.Linear(1 + sigma_embed_dim, ns),nn.Dropout(dropout), nn.ReLU(), nn.Linear(ns, 1))
            self.rot_final_layer = nn.Sequential(nn.Linear(1 + sigma_embed_dim, ns),nn.Dropout(dropout), nn.ReLU(), nn.Linear(ns, 1))

            if not no_torsion:
                # torsion angles components
                self.final_edge_embedding = nn.Sequential(
                    nn.Linear(distance_embed_dim, ns),
                    nn.ReLU(),
                    nn.Dropout(dropout),
                    nn.Linear(ns, ns)
                )
                self.final_tp_tor = o3.FullTensorProduct(self.sh_irreps, "2e")
                self.tor_bond_conv = TensorProductConvLayer(
                    in_irreps=self.lig_conv_layers[-1].out_irreps,
                    sh_irreps=self.final_tp_tor.irreps_out,
                    out_irreps=f'{ns}x0o + {ns}x0e',
                    n_edge_features=3 * ns,
                    residual=False,
                    dropout=dropout,
                    batch_norm=batch_norm
                )
                self.tor_final_layer = nn.Sequential(
                    nn.Linear(2 * ns, ns, bias=False),
                    nn.Tanh(),
                    nn.Dropout(dropout),
                    nn.Linear(ns, 1, bias=False)
                )
    def forward(self, data):

        # build ligand graph
        lig_node_attr, lig_edge_index, lig_edge_attr, lig_edge_sh = self.build_lig_conv_graph(data)
        lig_src, lig_dst = lig_edge_index
        lig_node_attr = self.lig_node_embedding(lig_node_attr)
        lig_edge_attr = self.lig_edge_embedding(lig_edge_attr)
        # build receptor graph
        rec_node_attr, rec_edge_index, rec_edge_attr, rec_edge_sh = self.build_rec_conv_graph(data)
        rec_src, rec_dst = rec_edge_index
        rec_node_attr = self.rec_node_embedding(data['receptor'].x)
        rec_edge_attr = self.rec_edge_embedding(rec_edge_attr)
        # build surface graph
        surface_node_attr,surface_edge_index, surface_edge_attr, surface_edge_sh = self.build_surface_conv_graph(data)
        surface_src, surface_dst = surface_edge_index
        surface_node_attr = self.surface_node_embedding(surface_node_attr)
        surface_edge_attr = self.surface_edge_embedding(surface_edge_attr)

        # update receptor embedding and then update embedding to surface
        rec_edge_attr_ = torch.cat([rec_edge_attr, rec_node_attr[rec_src, :self.ns], rec_node_attr[rec_dst, :self.ns]], -1)
        rec_intra_update = self.rec_conv_layers[0](rec_node_attr, rec_edge_index, rec_edge_attr_, rec_edge_sh)
        rec_node_attr = F.pad(rec_node_attr, (0, rec_intra_update.shape[-1] - rec_node_attr.shape[-1]))
        rec_node_attr = rec_node_attr + rec_intra_update


        # surface ,residue cross graph builld this info will use one shot
        surface_rec_cross_edge_index, surface_rec_cross_edge_attr, surface_rec_cross_edge_sh = self.build_surface_rec_cross_conv_graph(data)
        surface_rec_cross_rec, surface_rec_cross_surface = surface_rec_cross_edge_index
        surface_rec_cross_edge_attr = self.surface_rec_cross_edge_embedding( surface_rec_cross_edge_attr)

        residue_to_surface_edge_attr_ = torch.cat([surface_rec_cross_edge_attr, rec_node_attr[ surface_rec_cross_rec, :self.ns], surface_node_attr[surface_rec_cross_surface, :self.ns]], -1)
        # just one layer for feature update ,maybe can add more layers?
        surface_inter_residue_update = self.residue_to_surface_conv_layers[0](rec_node_attr, torch.flip(surface_rec_cross_edge_index,dims = [0]), residue_to_surface_edge_attr_, surface_rec_cross_edge_sh,
                                                              out_nodes=surface_node_attr.shape[0])
        surface_node_attr = F.pad(surface_node_attr, (0, surface_inter_residue_update.shape[-1] - surface_node_attr.shape[-1]))
        surface_node_attr = surface_node_attr + surface_inter_residue_update


        for l in range(len(self.lig_conv_layers)):
            # intra graph message passing
            lig_edge_attr_ = torch.cat([lig_edge_attr, lig_node_attr[lig_src, :self.ns], lig_node_attr[lig_dst, :self.ns]], -1)
            lig_intra_update = self.lig_conv_layers[l](lig_node_attr, lig_edge_index, lig_edge_attr_, lig_edge_sh)
            if l != len(self.lig_conv_layers) - 1:
                surface_edge_attr_ = torch.cat([surface_edge_attr, surface_node_attr[surface_src, :self.ns], surface_node_attr[surface_dst, :self.ns]], -1)
                surface_intra_update = self.surface_conv_layers[l](surface_node_attr, surface_edge_index, surface_edge_attr_, surface_edge_sh)
            # padding original features
            lig_node_attr = F.pad(lig_node_attr, (0, lig_intra_update.shape[-1] - lig_node_attr.shape[-1]))
            # update features with residual updates
            lig_node_attr = lig_node_attr + lig_intra_update #+ lig_inter_update

            if l != len(self.lig_conv_layers) - 1:
                surface_node_attr = F.pad(surface_node_attr, (0, surface_intra_update.shape[-1] - surface_node_attr.shape[-1]))
                surface_node_attr = surface_node_attr + surface_intra_update
            # print('after pad',lig_node_attr.shape,rec_node_attr.shape)
        if self.mdn_mode:
            scalar_lig_attr = torch.cat([lig_node_attr[:,:self.ns],lig_node_attr[:,-self.ns:] ], dim=1) if self.num_conv_layers >= 3 else lig_node_attr[:,:self.ns]
            scalar_rec_attr = torch.cat([surface_node_attr[:,:self.ns],surface_node_attr[:,-self.ns:] ], dim=1) if self.num_conv_layers >= 3 else surface_node_attr[:,:self.ns]
            h_l_x, l_mask = to_dense_batch(scalar_lig_attr, data['ligand'].batch, fill_value=0)
            h_t_x, t_mask = to_dense_batch(scalar_rec_attr, data['surface'].batch, fill_value=0)
            h_l_pos, _ = to_dense_batch(data['ligand'].pos, data['ligand'].batch, fill_value=0)
            h_t_pos, _ = to_dense_batch(data['surface'].pos, data['surface'].batch, fill_value=0)
            # aeesrtEncountered unequal batch-sizes
            assert h_l_x.size(0) == h_t_x.size(0), 'Encountered unequal batch-sizes'
            (B, N_l, C_out), N_t = h_l_x.size(), h_t_x.size(1)
            h_l_x = h_l_x.unsqueeze(-2)
            h_l_x = h_l_x.repeat(1, 1, N_t, 1) # [B, N_l, N_t, C_out]
            h_t_x = h_t_x.unsqueeze(-3)
            h_t_x = h_t_x.repeat(1, N_l, 1, 1) # [B, N_l, N_t, C_out]
            ########
            C = torch.cat((h_l_x, h_t_x), -1)
            self.C_mask = C_mask = l_mask.view(B, N_l, 1) & t_mask.view(B, 1, N_t)
            self.C = C = C[C_mask]
            C = self.MLP(C)
            # Get batch indexes for ligand-target combined features
            C_batch = torch.tensor(range(B)).unsqueeze(-1).unsqueeze(-1).to(self.device)
            C_batch = C_batch.repeat(1, N_l, N_t)[C_mask]#.to(self.device)
            # Outputs interactions predictions
            pi = F.softmax(self.z_pi(C), -1)
            sigma = F.elu(self.z_sigma(C))+1.1
            mu = F.elu(self.z_mu(C))+1
            dist = compute_euclidean_distances_matrix(h_l_pos, h_t_pos)[C_mask]
            mdn_loss_interaction = mdn_loss_fn(pi, sigma, mu, dist.unsqueeze(1).detach(),dist_threhold=self.args.mdn_dist_threshold_train if self.args.mdn_dist_threshold_train is not None else 7.0)

            if self.args.ligand_distance_prediction:
                # output for ligand distance predictions
                h_l_x, l_mask = to_dense_batch(scalar_lig_attr, data['ligand'].batch, fill_value=0)
                h_l = h_l_x.unsqueeze(-2)
                h_l = h_l.repeat(1, 1, N_l, 1) # [B, N_l, N_t, C_out]
                h_t_x = h_l_x.unsqueeze(-3)
                h_t_x = h_t_x.repeat(1, N_l, 1, 1) # [B, N_l, N_t, C_out]
                C_ligand = torch.cat((h_l, h_t_x), -1)
                C_mask_ligand = l_mask.view(B, N_l, 1) & l_mask.view(B, 1, N_l)

                sample = 1-torch.eye(N_l).to(self.device)
                sample = sample.unsqueeze(0)
                C_mask_ligand = C_mask_ligand*sample.repeat(B, 1, 1).bool()

                self.C_ligand = C_ligand = C_ligand[C_mask_ligand]
                C_ligand = self.ligand_MLP(C_ligand)

                pi_ligand = F.softmax(self.ligand_z_pi(C_ligand), -1)
                sigma_ligand = F.elu(self.ligand_z_sigma(C_ligand))+1.1
                mu_ligand = F.elu(self.ligand_z_mu(C_ligand))+1
                dist_ligand = compute_euclidean_distances_matrix(h_l_pos, h_l_pos)[C_mask_ligand]
                mdn_loss_ligand = mdn_loss_fn(pi_ligand, sigma_ligand, mu_ligand, dist_ligand.unsqueeze(1).detach(),dist_threhold=self.args.mdn_dist_threshold_train if self.args.mdn_dist_threshold_train is not None else 7.0)
            else:
                mdn_loss_ligand = mdn_loss_interaction*0.0
            if self.args.atom_type_prediction:
            ####### node type predictions
            # ligand atom types 119
                atom_types_label = data['ligand'].x[:,0]
                atom_types_pred = self.atom_types(scalar_lig_attr)
                # classification loss
                atom_types_loss = F.cross_entropy(atom_types_pred, atom_types_label)
            else:
                atom_types_loss = mdn_loss_interaction*0.0
            if self.args.bond_type_prediction:
                    # bond types 4
                bond_types_label = torch.argmax(data['ligand', 'lig_bond', 'ligand'].edge_attr[:,:4], dim=-1, keepdim=False)
                bond_types_pred = self.bond_types(torch.cat([scalar_lig_attr[data['ligand', 'ligand'].edge_index[0]], scalar_lig_attr[data['ligand', 'ligand'].edge_index[1]]], axis=1))
                        # # classification loss
                bond_types_loss = F.cross_entropy(bond_types_pred, bond_types_label)
                # print('bond_types_loss',bond_types_loss)
            else:
                bond_types_loss = mdn_loss_interaction*0.0
            if self.args.residue_type_prediction:
            # residue types 38
                residue_types_label = data['receptor'].x[:,0].long()
                residue_types_pred = self.residue_types(scalar_rec_attr)
                residue_types_loss = F.cross_entropy(residue_types_pred, residue_types_label)
            else:
                residue_types_loss = mdn_loss_interaction*0.0
            if self.training:
                if torch.isnan(mdn_loss_interaction  + mdn_loss_ligand  + atom_types_loss + bond_types_loss + residue_types_loss):
                    print(mdn_loss_interaction , mdn_loss_ligand , atom_types_loss , bond_types_loss , residue_types_loss)
                return mdn_loss_interaction , mdn_loss_ligand , atom_types_loss , bond_types_loss , residue_types_loss

            else:
                # probx = 0.0
                # for dist in dists_list:
                prob = calculate_probablity(pi, sigma, mu, dist.unsqueeze(1).detach(),dist_threhold=self.args.mdn_dist_threshold_test if self.args.mdn_dist_threshold_test is not None else 5.0)
                probx = scatter_add(prob,C_batch, dim=0, dim_size=B)
                return mdn_loss_interaction , mdn_loss_ligand , atom_types_loss , bond_types_loss , residue_types_loss,probx
    def build_lig_conv_graph(self, data):
   
        # compute edges
        radius_edges = radius_graph(data['ligand'].pos, self.lig_max_radius, data['ligand'].batch)
        edge_index = torch.cat([data['ligand', 'ligand'].edge_index, radius_edges], 1).long()
        edge_attr = torch.cat([
            data['ligand', 'ligand'].edge_attr,
            torch.zeros(radius_edges.shape[-1], self.in_lig_edge_features, device=data['ligand'].x.device)
        ], 0)
        # compute initial features

        edge_attr = edge_attr#torch.cat([edge_attr, edge_sigma_emb], 1)
        node_attr = data['ligand'].x#torch.cat([data['ligand'].x, data['ligand'].node_sigma_emb], 1)
        src, dst = edge_index
        edge_vec = data['ligand'].pos[dst.long()] - data['ligand'].pos[src.long()]
        # 
        edge_length_emb = self.lig_distance_expansion(edge_vec.norm(dim=-1))
        edge_attr = torch.cat([edge_attr, edge_length_emb], 1)
        edge_sh = o3.spherical_harmonics(self.sh_irreps, edge_vec, normalize=True, normalization='component')
        return node_attr, edge_index, edge_attr, edge_sh
    def build_rec_conv_graph(self, data):
        # builds the receptor initial node and edge embeddings

        node_attr = data['receptor'].x#torch.cat([data['receptor'].x, data['receptor'].node_sigma_emb], 1)
        # this assumes the edges were already created in preprocessing since protein's structure is fixed
        edge_index = data['receptor', 'receptor'].edge_index
        src, dst = edge_index
        edge_vec = data['receptor'].pos[dst.long()] - data['receptor'].pos[src.long()]
        edge_length_emb = self.rec_distance_expansion(edge_vec.norm(dim=-1))
        # edge_sigma_emb = data['receptor'].node_sigma_emb[edge_index[0].long()]

        edge_attr =torch.cat([data['receptor', 'rec_contact', 'receptor'].edge_attr, edge_length_emb], 1).float()
        edge_sh = o3.spherical_harmonics(self.sh_irreps, edge_vec, normalize=True, normalization='component')
        return node_attr, edge_index, edge_attr, edge_sh
    def build_cross_conv_graph(self, data, cross_distance_cutoff):
        # builds the cross edges between ligand and receptor
        if torch.is_tensor(cross_distance_cutoff):
            # different cutoff for every graph (depends on the diffusion time)
            edge_index = radius(data['receptor'].pos / cross_distance_cutoff[data['receptor'].batch],
                                data['ligand'].pos / cross_distance_cutoff[data['ligand'].batch], 1,
                                data['receptor'].batch, data['ligand'].batch, max_num_neighbors=30)
        else:
            edge_index = radius(data['receptor'].pos, data['ligand'].pos, cross_distance_cutoff,
                            data['receptor'].batch, data['ligand'].batch, max_num_neighbors=30)
        src, dst = edge_index
        edge_vec = data['receptor'].pos[dst.long()] - data['ligand'].pos[src.long()]
        edge_length_emb = self.cross_distance_expansion(edge_vec.norm(dim=-1))
        # edge_sigma_emb = data['ligand'].node_sigma_emb[src.long()]
        edge_attr = edge_length_emb
        edge_sh = o3.spherical_harmonics(self.sh_irreps, edge_vec, normalize=True, normalization='component')
        return edge_index, edge_attr, edge_sh
    def build_center_conv_graph(self, data):
        # builds the filter and edges for the convolution generating translational and rotational scores
        edge_index = torch.cat([data['ligand'].batch.unsqueeze(0), torch.arange(len(data['ligand'].batch)).to(data['ligand'].x.device).unsqueeze(0)], dim=0)

        center_pos, count = torch.zeros((data.num_graphs, 3)).to(data['ligand'].x.device), torch.zeros((data.num_graphs, 3)).to(data['ligand'].x.device)
        center_pos.index_add_(0, index=data['ligand'].batch, source=data['ligand'].pos)
        center_pos = center_pos / torch.bincount(data['ligand'].batch).unsqueeze(1)

        edge_vec = data['ligand'].pos[edge_index[1]] - center_pos[edge_index[0]]
        edge_attr = self.center_distance_expansion(edge_vec.norm(dim=-1))
        edge_sigma_emb = data['ligand'].node_sigma_emb[edge_index[1].long()]
        edge_attr = torch.cat([edge_attr, edge_sigma_emb], 1)
        edge_sh = o3.spherical_harmonics(self.sh_irreps, edge_vec, normalize=True, normalization='component')
        return edge_index, edge_attr, edge_sh
    
    def build_bond_conv_graph(self, data):
        # builds the graph for the convolution between the center of the rotatable bonds and the neighbouring nodes
        bonds = data['ligand', 'ligand'].edge_index[:, data['ligand'].edge_mask].long()
        bond_pos = (data['ligand'].pos[bonds[0]] + data['ligand'].pos[bonds[1]]) / 2
        bond_batch = data['ligand'].batch[bonds[0]]
        edge_index = radius(data['ligand'].pos, bond_pos, self.lig_max_radius, batch_x=data['ligand'].batch, batch_y=bond_batch)
        edge_vec = data['ligand'].pos[edge_index[1]] - bond_pos[edge_index[0]]
        edge_attr = self.lig_distance_expansion(edge_vec.norm(dim=-1))
        edge_attr = self.final_edge_embedding(edge_attr)
        edge_sh = o3.spherical_harmonics(self.sh_irreps, edge_vec, normalize=True, normalization='component')
        return bonds, edge_index, edge_attr, edge_sh
    def build_surface_conv_graph(self, data):
        node_attr = torch.nan_to_num(data['surface'].x)
        # this assumes the edges were already created in preprocessing since protein's structure is fixed
        edge_index = data['surface','surface_edge','surface'].edge_index
        src, dst = edge_index
        edge_vec = data['surface'].pos[dst.long()] - data['surface'].pos[src.long()]
        edge_length_emb = self.surface_distance_expansion(edge_vec.norm(dim=-1))
        # edge_sigma_emb = data['surface'].node_sigma_emb[edge_index[0].long()]
        edge_attr = torch.cat([data['surface','surface_edge','surface'].edge_attr, edge_length_emb], 1).float()
        edge_sh = o3.spherical_harmonics(self.sh_irreps, edge_vec, normalize=True, normalization='component')
        return node_attr, edge_index, edge_attr, edge_sh
    
    def build_surface_rec_cross_conv_graph(self, data, cross_distance_cutoff = 15):
        edge_index = radius(data['surface'].pos, data['receptor'].pos, cross_distance_cutoff,
                        data['surface'].batch, data['receptor'].batch, max_num_neighbors=30)
        src, dst = edge_index
        edge_vec = data['surface'].pos[dst.long()] - data['receptor'].pos[src.long()]

        edge_length_emb = self.cross_distance_expansion(edge_vec.norm(dim=-1))
        # edge_sigma_emb = data['receptor'].node_sigma_emb[src.long()]
        edge_attr = edge_length_emb#torch.cat([edge_sigma_emb, edge_length_emb], 1)
        edge_sh = o3.spherical_harmonics(self.sh_irreps, edge_vec, normalize=True, normalization='component')
        return edge_index, edge_attr, edge_sh
class GaussianSmearing(torch.nn.Module):
    # used to embed the edge distances
    def __init__(self, start=0.0, stop=5.0, num_gaussians=50):
        super().__init__()
        offset = torch.linspace(start, stop, num_gaussians)
        self.coeff = -0.5 / (offset[1] - offset[0]).item() ** 2
        self.register_buffer('offset', offset)

    def forward(self, dist):
        dist = dist.view(-1, 1) - self.offset.view(1, -1)
        return torch.exp(self.coeff * torch.pow(dist, 2))