text:
import os
import torch
from torch.utils.tensorboard import SummaryWriter
from torch import nn
from torch import distributions
from torch.utils.data.sampler import SubsetRandomSampler
from flows.NVP_flow.NVP_conformer import RealNVP
import schnetpack as spk
import numpy as np
import os
import sys
from GGNN.graphinvent import models as gnn_model
from torch.utils.data import DataLoader
from transBG_model import TransBG_model
import random
import time
from statistics import mean
from openff.toolkit.topology import Molecule
from openff.toolkit.topology import Topology
from openff.toolkit.typing.engines.smirnoff import ForceField
from rdkit.Chem.rdmolfiles import PDBWriter
from rdkit.Geometry import Point3D

from utils.conf_dataset import conformer_dataset, collate_conf
from utils.z_matrix import deconstruct_z_matrix, correct_conf_indexes
from utils.conf_energy import OpenMMEnergyWrapper, get_conformations_energy

class TransBG:
    def __init__(self, params):
        self.params = params

    def build_model(self):
        p = self.params

        #Setting a fixed seed for reproducibility
        torch.manual_seed(p.random_seed)
        random.seed(p.random_seed)

        print('Preparing data...', end = '', flush = True)
        sys.stdout.flush()

        #Split data into train and validation set
        indices = [i for i in list(range(p.dataset_size)) if i not in p.energy_val_indices]
        split_1 = int(np.floor(p.val_set_size * p.dataset_size))
        split_2 =  int(np.floor( (p.val_set_size + p.test_set_size) * p.dataset_size))
        np.random.seed(p.random_seed)
        np.random.shuffle(indices)
        train_indices, val_indices, test_indices = indices[split_2:], indices[:split_1], indices[split_1:split_2]
        # Create data samplers and loaders:
        train_sampler = SubsetRandomSampler(train_indices)
        val_sampler = SubsetRandomSampler(val_indices)
        test_sampler = SubsetRandomSampler(test_indices)

        self.dataset = conformer_dataset(p.sdf_path, p.xyz_path, pdb_path = p.pdb_path)
        self.train_loader = DataLoader(self.dataset, batch_size = p.batch_size, shuffle = False, collate_fn = collate_conf, sampler = train_sampler) #Set shuffle to False as is mutually exclussive with a sampler
        self.val_loader = DataLoader(self.dataset, batch_size = p.batch_size, shuffle = False, collate_fn = collate_conf, sampler = val_sampler)
        self.test_loader = DataLoader(self.dataset, batch_size = p.batch_size, shuffle = False, collate_fn = collate_conf, sampler = test_sampler)

        #Now the same for energy based learning

        #Energy learning set
        self.energy_train_indices = p.energy_train_indices
        self.energy_val_indices = p.energy_val_indices
        self.max_energies = torch.tensor(p.max_energies, dtype = torch.float32, device = p.device)

        # Create data samplers and loaders:
        energy_train_sampler = SubsetRandomSampler(self.energy_train_indices)
        energy_val_sampler = SubsetRandomSampler(self.energy_val_indices)
        self.energy_train_loader = DataLoader(self.dataset, batch_size = p.e_batch_size, shuffle = False, collate_fn = collate_conf, sampler = energy_train_sampler) #Set shuffle to False as is mutually exclussive with a sampler
        self.energy_val_loader = DataLoader(self.dataset, batch_size = p.e_batch_size, shuffle = False, collate_fn = collate_conf, sampler = energy_val_sampler)

        print('Done!', flush = True)

        print('Building the model...', end = '', flush = True)
        sys.stdout.flush()

        #Length of the vector codifying the information that is fed into the flow: context representation
        rep_size = p.mol_rep_size + p.atom_rep_size + p.schnet_rep_size + p.ref_frame_rep_size + p.symmetry_token_size
        #Generator functions for the scaling and translation network (simple feed forward neural nets):
        nets = lambda: nn.Sequential(nn.Linear(3 + rep_size, p.intermediate_rep_size), nn.LeakyReLU(), nn.Linear(p.intermediate_rep_size, p.intermediate_rep_size), nn.LeakyReLU(), nn.Linear(p.intermediate_rep_size, 3 + rep_size), nn.Tanh())
        nett = lambda: nn.Sequential(nn.Linear(3 + rep_size, p.intermediate_rep_size), nn.LeakyReLU(), nn.Linear(p.intermediate_rep_size, p.intermediate_rep_size), nn.LeakyReLU(), nn.Linear(p.intermediate_rep_size, 3 + rep_size))
        #masking (alternate the dimensions that are unchanged and the one which is transformed): 0 means transformed. The context representation is always masked.
        mask = np.ones((3, rep_size + 3))
        mask[0, rep_size + 3 - 3] = 0
        mask_d = mask.tolist()
        masks_d = torch.from_numpy(np.array(mask_d * p.num_coupling_layers).astype(np.float32)) # Mask for transforming only distances
        # Note: this can seem not very efficient since we forward just 0s multiplications when not transforming, but this is handled at the flow model.
        mask[1, rep_size + 3 - 2] = 0
        mask_da = mask.tolist()
        masks_da = torch.from_numpy(np.array(mask_da * p.num_coupling_layers).astype(np.float32)) # Mask for transforming distances and angles
        mask[2, rep_size + 3 - 1] = 0
        mask_dat = mask.tolist()
        mask_dat = torch.from_numpy(np.array(mask_dat * p.num_coupling_layers).astype(np.float32))
        # define base distribution
        base_distribution = distributions.MultivariateNormal(torch.zeros(3), torch.eye(3))

        # ggnn instance
        gnn = gnn_model.initialize_GGNN(p)

        # flow instance
        flow = RealNVP(nets, nett, mask_d = masks_d, mask_da = masks_da, mask_dat = mask_dat, base_dist = base_distribution, device = p.device).to(p.device)
        # schnet instance
        atomic_reps = spk.representation.SchNet(n_atom_basis = p.schnet_rep_size, n_filters = p.schnet_rep_size, n_interactions = p.n_interactions,
                                                cutoff = p.cutoff, n_gaussians = p.n_gaussians, normalize_filter = p.normalize_filter, coupled_interactions = p.coupled_interactions, # cutoff is given in Angstrong
                                                max_z = p.max_z + 1, trainable_gaussians = p.trainable_gaussians)  # max_z needs one unit more or raises errors
        schnet_summarizer = spk.atomistic.Atomwise(n_in = p.schnet_rep_size, n_out = p.schnet_rep_size)
        schnet = spk.atomistic.AtomisticModel(atomic_reps, schnet_summarizer).to(p.device)

        self.model = TransBG_model(flow, schnet, gnn, p,  torch.device(p.device))
        #self.best_model = ConformationsGenerator(flow, schnet, gnn, p,  torch.device(p.device))    #We do not generate this for memory constrains.
        self.l_min_val_loss = float('inf')
        self.e_min_val_loss = float('inf')
        self.ft_min_val_loss = float('inf')
        self.l_training_metrics = {'train_loss' : [], 'val_loss' : [] }
        self.e_training_metrics = {'train_loss' : [], 'val_loss' : [], 'train_average_energy' : [], 'val_average_energy' : [] }

        print("Done!", flush = True)

    def train_likelihood(self):
        p = self.params

        #Setting a fixed seed for reproducibility
        torch.manual_seed(p.random_seed)
        random.seed(p.random_seed)

        #Define optimizer
        self.l_optimizer = torch.optim.Adam(self.model.parameters(), lr = p.l_learning_rate, weight_decay = p.weight_decay)

        #Create  learning rate scheduler

        self.l_scheduler = torch.optim.lr_scheduler.OneCycleLR(optimizer=self.l_optimizer,
                                                               max_lr= p.l_max_rel_lr * p.l_learning_rate,
                                                               div_factor= 1. / p.l_max_rel_lr,
                                                               final_div_factor = 1. / p.l_min_rel_lr,
                                                               pct_start = 0.01,
                                                               total_steps = p.l_epochs
                                                               )

        #Create a folder for the best model
        if not os.path.isdir("models"):
            os.mkdir("models/")

        #Create a folder for the training metrics
        if not os.path.isdir("metrics"):
            os.mkdir("metrics/")

        for epoch in range(1, p.l_epochs + 1):

            t0 = time.time()
            loss_sum = 0
            discarded_batches = 0

            for i_batch, batch in enumerate(self.train_loader ,1):

                mol_batch = batch
                self.model.train()
                z, log_det_invJ = self.model.pushback(mol_batch)

                if not torch.any( torch.isnan(z) ):
                    z = z.squeeze( dim = 1 )
                    #nan_flag = [ i_row for i_row in range( z.size()[0] ) if not torch.any( torch.isnan(z[i_row])) or torch.any( torch.isnan(log_det_invJ[i_row])) ]
                    #z = z[nan_flag]
                    #log_det_invJ = log_det_invJ[nan_flag]

                    #Likelihood loss:
                    armonic_en = torch.tensor(0.5, dtype = torch.float32) * z * z
                    armonic_en = armonic_en.squeeze(dim = 1).sum(dim = 1)
                    loss = armonic_en - log_det_invJ

                    loss = loss[loss == loss] # Ignore nan values, in case they accidentally appear
                    loss = torch.mean( loss )

                    self.l_optimizer.zero_grad()
                    loss.backward()
                    self.l_optimizer.step()

                    train_loss = loss.item()
                    loss_sum += train_loss
                else:
                    discarded_batches += 1
                    print('Discarded batch due to the presence of nan values')

                print('.', end = '')
                sys.stdout.flush()

            #Evaluate:
            print( ' \n Evaluating: ')
            self.model.eval()
            val_loss_sum = 0
            discarded_batches = 0
            with torch.no_grad():
                for i_batch, batch in enumerate(self.val_loader ,1):

                    mol_batch = batch
                    z, log_det_invJ = self.model.pushback(mol_batch)

                    if not torch.any( torch.isnan(z) ):
                        z = z.squeeze( dim = 1 )
                        nan_flag = [ i_row for i_row in range( z.size()[0] ) if not torch.any( torch.isnan(z[i_row]) ) or torch.any( torch.isnan(log_det_invJ[i_row]) ) ]
                        z = z[nan_flag]
                        log_det_invJ = log_det_invJ[nan_flag]

                        #Likelihood loss:
                        val_armonic_en = torch.tensor(0.5, dtype = torch.float32) * z * z
                        val_armonic_en = val_armonic_en.squeeze(dim = 1).sum(dim = 1)
                        val_loss = torch.mean(val_armonic_en - log_det_invJ)

                        val_loss_sum += val_loss.item()

                    else:
                        discarded_batches += 1
                        print('Discarded batch due to the presence of nan values')

                    print('.', end = '')
                    sys.stdout.flush()

            #Update learning rate
            self.l_scheduler.step()

            #Print some statistics:
            t1 = time.time()
            #train_loss = loss_sum/len(self.train_loader) # We print and use the train loss of the last batch (only fair to compare)
            val_loss = val_loss_sum/( len(self.val_loader) - discarded_batches )
            self.l_training_metrics['train_loss'].append(train_loss)
            np.save('metrics/train_loss_' + p.model_name + '.npy', self.l_training_metrics['train_loss'])
            self.l_training_metrics['val_loss'].append(val_loss)
            np.save('metrics/val_loss_' + p.model_name + '.npy', self.l_training_metrics['val_loss'])
            print(f' \n Epoch: {epoch}: training loss: {train_loss:.4f}, validation_loss: {val_loss:.4f}, time = {t1-t0:.4f}')
            sys.stdout.flush()

            #Keep the best model
            if val_loss < self.l_min_val_loss:
                self.l_min_val_loss = val_loss
                #self.best_model.load_state_dict( self.model.state_dict() )
                torch.save( self.model.state_dict(), 'models/' + p.model_name )

    def finetune_likelihood(self, min_e = False):
        p = self.params

        #Setting a fixed seed for reproducibility
        torch.manual_seed(p.random_seed)
        random.seed(p.random_seed)

        #Define optimizer
        self.ft_optimizer = torch.optim.Adam(self.model.parameters(), lr = p.ft_learning_rate, weight_decay = p.weight_decay)

        #Create  learning rate scheduler

        self.ft_scheduler = torch.optim.lr_scheduler.OneCycleLR(optimizer=self.ft_optimizer,
                                                            max_lr= p.ft_max_rel_lr * p.ft_learning_rate,
                                                            div_factor= 1. / p.ft_max_rel_lr,
                                                            final_div_factor = 1. / p.ft_min_rel_lr,
                                                            pct_start = 0.01,
                                                            total_steps = p.ft_epochs
                                                            )

        #Create a folder for the best model
        if not os.path.isdir("models"):
            os.mkdir("models/")

        #Create a folder for the training metrics
        if not os.path.isdir("metrics"):
            os.mkdir("metrics/")

        for epoch in range(1, p.ft_epochs + 1):

            t0 = time.time()
            loss_sum = 0
            discarded_batches = 0

            for i_batch, batch in enumerate(self.energy_train_loader ,1):

                mol_batch = batch
                self.model.train()
                z, log_det_invJ = self.model.pushback(mol_batch, min_e)

                if not torch.any( torch.isnan(z) ):
                    z = z.squeeze( dim = 1 )
                    #nan_flag = [ i_row for i_row in range( z.size()[0] ) if not torch.any( torch.isnan(z[i_row])) or torch.any( torch.isnan(log_det_invJ[i_row])) ]
                    #z = z[nan_flag]
                    #log_det_invJ = log_det_invJ[nan_flag]

                    #Likelihood loss:
                    armonic_en = torch.tensor(0.5, dtype = torch.float32) * z * z
                    armonic_en = armonic_en.squeeze(dim = 1).sum(dim = 1)
                    loss = armonic_en - log_det_invJ

                    loss = loss[loss == loss] # Ignore nan values, in case they accidentally appear
                    loss = torch.mean( loss )

                    self.ft_optimizer.zero_grad()
                    loss.backward()
                    self.ft_optimizer.step()

                    train_loss = loss.item()
                    loss_sum += train_loss
                else:
                    discarded_batches += 1
                    print('Discarded batch due to the presence of nan values')

                print('.', end = '')
                sys.stdout.flush()

            #Evaluate:
            print( ' \n Evaluating: ')
            self.model.eval()
            val_loss_sum = 0
            discarded_batches = 0
            with torch.no_grad():
                for i_batch, batch in enumerate(self.energy_val_loader ,1):

                    mol_batch = batch
                    z, log_det_invJ = self.model.pushback(mol_batch, min_e)

                    if not torch.any( torch.isnan(z) ):
                        z = z.squeeze( dim = 1 )
                        nan_flag = [ i_row for i_row in range( z.size()[0] ) if not torch.any( torch.isnan(z[i_row]) ) or torch.any( torch.isnan(log_det_invJ[i_row]) ) ]
                        z = z[nan_flag]
                        log_det_invJ = log_det_invJ[nan_flag]

                        #Likelihood loss:
                        val_armonic_en = torch.tensor(0.5, dtype = torch.float32) * z * z
                        val_armonic_en = val_armonic_en.squeeze(dim = 1).sum(dim = 1)
                        val_loss = torch.mean(val_armonic_en - log_det_invJ)

                        val_loss_sum += val_loss.item()

                    else:
                        discarded_batches += 1
                        print('Discarded batch due to the presence of nan values')

                    print('.', end = '')
                    sys.stdout.flush()

            #Update learning rate
            self.ft_scheduler.step()

            #Print some statistics:
            t1 = time.time()
            train_loss = loss_sum/len(self.energy_train_loader)
            val_loss = val_loss_sum/( len(self.energy_val_loader) - discarded_batches )
            #self.l_training_metrics['train_loss'].append(train_loss)
            #np.save('metrics/train_loss_' + p.model_name + '.npy', self.l_training_metrics['train_loss'])
            #self.l_training_metrics['val_loss'].append(val_loss)
            #np.save('metrics/val_loss_' + p.model_name + '.npy', self.l_training_metrics['val_loss'])
            print(f' \n Epoch: {epoch}: training loss: {train_loss:.4f}, validation_loss: {val_loss:.4f}, time = {t1-t0:.4f}')
            sys.stdout.flush()

    def generate_conformations(self, mol_batch, n_conf = 1):
        '''

        '''
        conformations, log_det_J = self.model.pushforward(mol_batch, n_conf)
        return conformations, log_det_J

    def train_energy(self, min_e = False):
        '''
            The model needs to have been pretrained (i.e. using likelihood-based learning).
        '''

        p = self.params

        #Setting a fixed seed for reproducibility
        torch.manual_seed(p.random_seed)
        random.seed(p.random_seed)

        #Select the trainable models:
        if not p.train_GNN_ebl:
            for parameter in self.model.gnn.parameters():
                parameter.requires_grad = False

        if not p.train_schnet_ebl:
            for parameter in self.model.schnet.parameters():
                parameter.requires_grad = False

        #Define optimizer
        self.e_optimizer = torch.optim.Adam(self.model.parameters(), lr = p.e_learning_rate, weight_decay = p.weight_decay)

        #Create  learning rate scheduler

        self.e_scheduler = torch.optim.lr_scheduler.OneCycleLR(optimizer=self.e_optimizer,
                                                               max_lr= p.e_max_rel_lr * p.e_learning_rate,
                                                               div_factor= 1. / p.e_max_rel_lr,
                                                               final_div_factor = 1. / p.e_min_rel_lr,
                                                               pct_start = 0.01,
                                                               total_steps = p.e_epochs
                                                               )

        #Instance of the energy wrapper
        energy_wrapper = OpenMMEnergyWrapper()

        # Create the OpenMM systems
        forcefield = ForceField('openff-1.2.0.offxml')
        systems_dict = {}
        for i_mol in p.energy_train_indices + p.energy_val_indices:
            mol_dict = self.dataset[i_mol]
            rd_mol = mol_dict['rdkit_mol']
            partial_charges = mol_dict['partial_charges']

            off_mol = Molecule.from_rdkit(rd_mol)
            off_mol.partial_charges = partial_charges
            topology = Topology.from_molecules(off_mol)
            openmm_system = forcefield.create_openmm_system(topology, charge_from_molecules = [off_mol], allow_nonintegral_charges=True)
            systems_dict[str(mol_dict['qm9_idx'])] = openmm_system

        #Create a folder for the best model
        if not os.path.isdir("models"):
            os.mkdir("models/")

        #Create a folder for the training metrics
        if not os.path.isdir("metrics"):
            os.mkdir("metrics/")

        #Tensorboard writer instance: Writer will output to ./runs/ directory by default
        if not os.path.isdir("runs/"+self.params.model_name):
            os.mkdir("runs/"+self.params.model_name)
        e_writer = SummaryWriter(log_dir="runs/"+self.params.model_name, comment='_rep')

        for epoch in range(1, p.e_epochs + 1):

            t0 = time.time()
            train_loss_sum = 0
            train_average_energy_sum = 0
            l_train_loss_sum = 0
            e_train_loss_sum = 0
            constr_train_loss_sum = 0

            pm_e_train_loss = torch.zeros((len(self.energy_train_indices)), dtype=torch.float32).to(p.device)

            train_energy = torch.zeros((len(self.energy_train_indices)), dtype=torch.float32).to(p.device)
            train_conf_rate = [0.0 for _ in range( len(self.energy_train_indices) )]
            train_discarded_batches = 0

            for i_batch, batch in enumerate(self.energy_train_loader ,1):

                mol_batch = batch
                self.model.train()

                #Likelihood loss:
                z, log_det_invJ = self.model.pushback(mol_batch, min_e)
                z = z.squeeze( dim = 1 )
                nan_flag = [ i_row for i_row in range( z.size()[0] ) if not torch.any( torch.isnan(z[i_row]) ) or torch.any( torch.isnan(log_det_invJ[i_row]) ) ]
                z = z[nan_flag]
                log_det_invJ = log_det_invJ[nan_flag]
                armonic_en = torch.tensor(0.5, dtype = torch.float32) * z * z
                armonic_en = armonic_en.squeeze(dim = 1).sum(dim = 1)
                l_train_loss = armonic_en - log_det_invJ
                l_train_loss = l_train_loss[l_train_loss == l_train_loss] # Ignore nan values, in case they accidentally appear
                l_train_loss = torch.mean( l_train_loss )
                l_train_loss_sum += l_train_loss.item()

                #Energy loss:

                conformations, log_det_J = self.model.pushforward(mol_batch = mol_batch, n_conf = p.n_conformations)

                #nan protection
                for i_mol, mol_dict in enumerate(mol_batch):
                    nan_conf_indexes = torch.unique( torch.nonzero(torch.isnan(conformations[i_mol]))[:,0] ).tolist()
                    nan_det_indexes = torch.unique( torch.nonzero(torch.isnan(log_det_J[i_mol]))[:,0] ).tolist()
                    nan_conf_indexes = list(set(nan_conf_indexes) | set(nan_det_indexes))
                    nan_conf_indexes.sort()
                    ok_conf_indexes = list( set(list(range(p.n_conformations))) - set(nan_conf_indexes) )
                    if nan_conf_indexes:
                        nan_mol_index = mol_dict['dataset_idx']
                        print(f'Found nans in {len(nan_conf_indexes)} conformations of molecule {nan_mol_index}. Ignored for training.')
                        conformations[i_mol] = conformations[i_mol][ok_conf_indexes]
                        log_det_J[i_mol] = log_det_J[i_mol][ok_conf_indexes]

                mol_indexes = [ mol['dataset_idx'] for mol in mol_batch]
                mol_indexes_train_energy_set = [ self.energy_train_indices.index( mol['dataset_idx'] ) for mol in mol_batch ]
                mol_indexes_batch = []
                e_train_loss_backward = torch.zeros((len(mol_batch)), dtype=torch.float32).to(p.device)

                for i_mol, mol in enumerate(mol_batch):
                    mol_index_train_energy_set = self.energy_train_indices.index( mol['dataset_idx'] )
                    mol_train_energy = torch.zeros( (conformations[i_mol].size()[0] ), dtype=torch.float32).to(p.device)#len(conformations[i_mol])
                    mol_constr_loss = torch.zeros( (conformations[i_mol].size()[0] ), dtype=torch.float32).to(p.device)#
                    energy_comp_indexes = []
                    for i_conf, conformation in enumerate(conformations[i_mol]):
                        cartesian_coord, mol_constr_loss[i_conf] = deconstruct_z_matrix( conformation.clone(), mol['ref_atoms'].copy() ) #The Z-matrixes are clamped at deconstruct_z_matrix
                        try:
                            energy = energy_wrapper.apply(cartesian_coord, systems_dict[str(mol['qm9_idx'])], p.device)
                            if not torch.isnan(energy):
                                mol_train_energy[i_conf] = energy.clone()
                                energy_comp_indexes.append(i_conf)
                        except:
                            #print('+', flush = True, end = '')
                            mol_train_energy[i_conf] = self.max_energies[ mol_index_train_energy_set ] + 1.
                            #print('\n There was a problem computing energies with conformation {idx_conf} of molecule {idx_mol}. This conformation is skipped for training.'.format(idx_mol = mol['dataset_idx'], idx_conf = i_conf), flush=True)
                    if len(energy_comp_indexes) != conformations[i_mol].size()[0]:
                        idx_mol = mol['dataset_idx']
                        energy_comp_rate = float(len(energy_comp_indexes)) / float(conformations[i_mol].size()[0])
                        print( f'OpenMM was able to compute the energy for {energy_comp_rate*100:.4f} % of the conformations of molecule {idx_mol}.' )
                    constr_train_loss_sum += torch.mean(mol_constr_loss).item()
                    mol_train_energy_backward = mol_train_energy[ mol_train_energy < self.max_energies[ mol_index_train_energy_set ] ].clone()
                    mol_log_det_J_backward = log_det_J[i_mol][mol_train_energy < self.max_energies[mol_index_train_energy_set]].clone()
                    mol_constr_loss_backward = mol_constr_loss[mol_train_energy < self.max_energies[mol_index_train_energy_set]].clone()
                    mol_e_train_loss_backward = mol_train_energy_backward - mol_log_det_J_backward + mol_constr_loss_backward
                    if mol_e_train_loss_backward.numel():
                        e_train_loss_backward[i_mol] = torch.mean(mol_e_train_loss_backward)
                        mol_indexes_batch.append(i_mol)
                    else:
                        dataset_idx = mol['dataset_idx']
                        print(f'\n No learnable conformations generated for molecule {dataset_idx}.')
                        mol_indexes.remove(dataset_idx)
                        mol_indexes_train_energy_set.remove(mol_index_train_energy_set)

                    mol_train_energy = mol_train_energy[energy_comp_indexes]
                    mol_log_det_J = log_det_J[i_mol][energy_comp_indexes].clone()#.detach()
                    mol_e_train_loss = mol_train_energy - mol_log_det_J
                    pm_e_train_loss[ mol_index_train_energy_set ] = torch.mean(mol_e_train_loss)
                    train_energy[ mol_index_train_energy_set ] = torch.mean(mol_train_energy)

                    # Percentage of used conformations
                    mol_conf_rate = mol_train_energy_backward.size()[0]
                    if mol_conf_rate:
                        mol_conf_rate = float(mol_conf_rate)/p.n_conformations
                    else:
                        mol_conf_rate = 0.

                    train_conf_rate[ mol_index_train_energy_set ] = mol_conf_rate

                e_train_loss_mean =  torch.mean(pm_e_train_loss[mol_indexes_train_energy_set])  # No track gradients # Of the batch
                e_train_loss_backward_mean =  torch.mean(e_train_loss_backward[mol_indexes_batch])

                if e_train_loss_backward_mean.numel():
                    total_loss_backward = p.l_weight*l_train_loss + p.e_weight*e_train_loss_backward_mean
                    self.e_optimizer.zero_grad()
                    total_loss_backward.backward()
                    self.e_optimizer.step()

                else:
                    train_discarded_batches += 1
                    print('Discarded batch due to non-learnable configurations.')

                total_loss = p.l_weight*l_train_loss + p.e_weight*e_train_loss_mean

                train_loss = total_loss.item()
                train_loss_sum += train_loss

                train_average_energy = torch.mean(train_energy).item()
                train_average_energy_sum += train_average_energy

                e_train_loss_sum += e_train_loss_mean.item()

                print('.', end = '')
                sys.stdout.flush()

            #Evaluate:
            print( ' \n Evaluating: ')
            self.model.eval()
            val_loss_sum = 0
            val_average_energy_sum = 0
            l_val_loss_sum = 0
            e_val_loss_sum = 0
            constr_val_loss_sum = 0

            pm_e_val_loss = torch.zeros((len(self.energy_val_indices)), dtype=torch.float32).to(p.device)

            val_energy = torch.zeros((len(self.energy_val_indices)), dtype=torch.float32).to(p.device)
            val_conf_rate = [0.0 for _ in range( len(self.energy_val_indices) )]
            val_discarded_batches = 0

            with torch.no_grad():
                for i_batch, mol_batch in enumerate(self.energy_val_loader ,1):

                    #Likelihood loss:
                    z, log_det_invJ = self.model.pushback(mol_batch, min_e)
                    z = z.squeeze( dim = 1 )
                    nan_flag = [ i_row for i_row in range( z.size()[0] ) if not torch.any( torch.isnan(z[i_row]) ) or torch.any( torch.isnan(log_det_invJ[i_row]) ) ]
                    z = z[nan_flag]
                    log_det_invJ = log_det_invJ[nan_flag]
                    armonic_en = torch.tensor(0.5, dtype = torch.float32) * z * z
                    armonic_en = armonic_en.squeeze(dim = 1).sum(dim = 1)
                    l_val_loss = armonic_en - log_det_invJ
                    l_val_loss = l_val_loss[l_val_loss == l_val_loss] # Ignore nan values, in case they accidentally appear
                    l_val_loss = torch.mean( l_val_loss )
                    l_val_loss_sum += l_val_loss.item()

                    #Energy loss:

                    conformations, log_det_J = self.model.pushforward(mol_batch = mol_batch, n_conf = p.n_conformations)

                    #nan protection
                    for i_mol, mol_dict in enumerate(mol_batch):
                        nan_conf_indexes = torch.unique( torch.nonzero(torch.isnan(conformations[i_mol]))[:,0] ).tolist()
                        nan_det_indexes = torch.unique( torch.nonzero(torch.isnan(log_det_J[i_mol]))[:,0] ).tolist()
                        nan_conf_indexes = list(set(nan_conf_indexes) | set(nan_det_indexes))
                        nan_conf_indexes.sort()
                        ok_conf_indexes = list( set(list(range(p.n_conformations))) - set(nan_conf_indexes) )
                        if nan_conf_indexes:
                            nan_mol_index = mol_dict['dataset_idx']
                            print(f'Found nans in {len(nan_conf_indexes)} conformations of molecule {nan_mol_index}. Ignored for training.')
                            conformations[i_mol] = conformations[i_mol][ok_conf_indexes]
                            log_det_J[i_mol] = log_det_J[i_mol][ok_conf_indexes]

                    mol_indexes = [ mol['dataset_idx'] for mol in mol_batch]
                    mol_indexes_val_energy_set = [ self.energy_val_indices.index( mol['dataset_idx'] ) for mol in mol_batch ]
                    mol_indexes_batch = []
                    #e_val_loss_backward = torch.zeros((len(mol_batch)), dtype=torch.float32).to(p.device)

                    for i_mol, mol in enumerate(mol_batch):
                        mol_index_val_energy_set = self.energy_val_indices.index( mol['dataset_idx'] )
                        mol_val_energy = torch.zeros( conformations[i_mol].size()[0], dtype=torch.float32).to(p.device)
                        mol_constr_val_loss = torch.zeros( conformations[i_mol].size()[0], dtype=torch.float32).to(p.device)#p.n_conformations
                        energy_comp_indexes = []
                        for i_conf, conformation in enumerate(conformations[i_mol]):
                            cartesian_coord, mol_constr_val_loss[i_conf] = deconstruct_z_matrix( conformation.clone(), mol['ref_atoms'].copy() ) #The Z-matrixes are clamped at deconstruct_z_matrix
                            try:
                                energy = energy_wrapper.apply(cartesian_coord, systems_dict[str(mol['qm9_idx'])], p.device)
                                if not torch.isnan(energy):
                                    mol_val_energy[i_conf] = energy.clone()
                                    energy_comp_indexes.append(i_conf)
                            except:
                                pass
                                #print('+', flush = True, end = '')
                                #print('\n There was a problem computing energies with conformation {idx_conf} of molecule {idx_mol}. This conformation is skipped for evaluation.'.format(idx_mol = mol['dataset_idx'], idx_conf = i_conf), flush=True)

                        if len(energy_comp_indexes) != conformations[i_mol].size()[0]:
                            idx_mol = mol['dataset_idx']
                            energy_comp_rate = float( len(energy_comp_indexes) ) / float(conformations[i_mol].size()[0])
                            print( f'OpenMM was able to compute the energy for {energy_comp_rate*100:.4f} % of the conformations of molecule {idx_mol}.' )

                        mol_val_energy = mol_val_energy[energy_comp_indexes]
                        constr_val_loss_sum += torch.mean(mol_constr_val_loss).item()
                        mol_constr_val_loss = mol_constr_val_loss[energy_comp_indexes]
                        mol_log_det_J = log_det_J[i_mol][energy_comp_indexes].clone().detach()
                        mol_e_val_loss = mol_val_energy - mol_log_det_J + mol_constr_val_loss
                        pm_e_val_loss[ mol_index_val_energy_set ] = torch.mean(mol_e_val_loss)

                        val_energy[ mol_index_val_energy_set ] = torch.mean(mol_val_energy)

                    # Compute statistics
                    e_val_loss_mean =  torch.mean(pm_e_val_loss[mol_indexes_val_energy_set])
                    total_loss = p.l_weight*l_val_loss + p.e_weight*e_val_loss_mean

                    val_loss = total_loss.item()
                    val_loss_sum += val_loss

                    val_average_energy = torch.mean(val_energy).item()
                    val_average_energy_sum += val_average_energy

                    e_val_loss_sum += e_val_loss_mean.item()

                    print('.', end = '')
                    sys.stdout.flush()

            #Update learning rate
            self.e_scheduler.step()

            #Print some (a lot of :) )  statistics:
            t1 = time.time()
            train_loss = train_loss_sum/ ( len(self.energy_train_loader) -  train_discarded_batches )
            val_loss = val_loss_sum / ( len(self.energy_val_loader) -  val_discarded_batches )
            train_average_energy = train_average_energy_sum / len(self.energy_train_loader)
            val_average_energy = val_average_energy_sum / len(self.energy_val_loader)
            l_train_loss = l_train_loss_sum / len(self.energy_train_loader)
            e_train_loss = e_train_loss_sum / len(self.energy_train_loader)
            l_val_loss = l_val_loss_sum / len(self.energy_val_loader)
            e_val_loss = e_val_loss_sum / len(self.energy_val_loader)
            constr_train_loss = constr_train_loss_sum / len(self.energy_train_indices)
            constr_val_loss = constr_val_loss_sum / len(self.energy_val_indices)

            self.e_training_metrics['train_loss'].append(train_loss)
            np.save('metrics/e_train_loss_' + p.model_name + '.npy', self.e_training_metrics['train_loss'])
            self.e_training_metrics['val_loss'].append(val_loss)
            np.save('metrics/e_val_loss_' + p.model_name + '.npy', self.e_training_metrics['val_loss'])
            self.e_training_metrics['train_average_energy'].append(train_average_energy)
            np.save('metrics/e_train_average_energy_' + p.model_name + '.npy', self.e_training_metrics['train_average_energy'])
            self.e_training_metrics['val_average_energy'].append(val_average_energy)
            np.save('metrics/e_val_average_energy_' + p.model_name + '.npy', self.e_training_metrics['val_average_energy'])
            print(f' \n Epoch: {epoch}: training loss: {train_loss:.4f}, validation_loss: {val_loss:.4f}, training average energy: {train_average_energy:.4f}, validation average energy: {val_average_energy:.4f}, \n likelihood train loss: {l_train_loss:.4f}, likelihood validation loss: {l_val_loss:.4f}, energy train loss: {e_train_loss:.4f}, energy validation loss: {e_val_loss:.4f}, \n train conformations rate: {mean(train_conf_rate):4f}, training constraint loss: {constr_train_loss:.4f}, validation constraint loss: {constr_val_loss:.4f}, time = {t1-t0:.4f}')
            sys.stdout.flush()

            e_writer.add_scalars('Loss/Total loss', {'Total train loss':train_loss, 'Total val loss':val_loss}, epoch)
            e_writer.add_scalars('Loss/Likelihood loss', {'Likelihood train loss':l_train_loss, 'Likelihood val loss':l_val_loss}, epoch)
            e_writer.add_scalars('Loss/Energy loss', {'Energy train loss':e_train_loss, 'Energy val loss':e_val_loss}, epoch)

            # Per molecule statistics
            print('\n Metrics for molecules in the training set:')

            for mol_index_dataset, mol_index in enumerate(self.energy_train_indices):
                print(f' \n Molecule: {mol_index}: average energy: {train_energy[mol_index_dataset].item():.4f}, energy (train) loss: {pm_e_train_loss[mol_index_dataset].item():.4f} \n conformations rate: {train_conf_rate[mol_index_dataset]:4f}')

            print('\n Metrics for molecules in the validation set:')
            for mol_index_dataset, mol_index in enumerate(self.energy_val_indices):
                print(f' \n Molecule: {mol_index}: average energy: {val_energy[mol_index_dataset].item():.4f}, energy (validation) loss: {pm_e_val_loss[mol_index_dataset].item():.4f}')

            #Save metrics so that they can be plotted
            sys.stdout.flush()

            #Compute the energy histograms at this epoch

            if self.params.hist_evo:
                print('Computing training set histograms...', flush=True)
                for mol_index in self.params.energy_train_indices:
                    conformations, _ = self.generate_histogram(mol_index=mol_index, n_conf=self.params.hist_conf, step_conf= 5000, visualize = False, save_files = False, verbose = False, ignore_constraints = False)
                    e_writer.add_histogram('Train histograms/Molecule '+str(mol_index),np.array(conformations), epoch)
                print('Done!', flush=True)
                print('Computing validation set histograms...', flush=True)
                for mol_index in self.params.energy_val_indices:
                    conformations, _ = self.generate_histogram(mol_index=mol_index, n_conf=self.params.hist_conf, step_conf= 5000, visualize = False, save_files = False, verbose = False, ignore_constraints = False)
                    e_writer.add_histogram('Val histograms/Molecule '+str(mol_index),np.array(conformations), epoch)
                print('Done!', flush=True)

            #Keep the best model
            if val_loss < self.e_min_val_loss:
                self.e_min_val_loss = val_loss
                #self.best_model.load_state_dict( self.model.state_dict() )
                torch.save( self.model.state_dict(), 'models/e_tuned_' + p.model_name )

        e_writer.close()

    def generate_correct_conformations(self, mol_index, n_conf, step_conf):
        '''

        '''
        mol_dict = self.dataset[mol_index]
        mol_batch = [mol_dict]
        n_ok_conf = 0
        conformations = torch.empty(0, dtype=torch.float32)
        log_det_J = torch.empty(0, dtype=torch.float32)
        while n_ok_conf < n_conf:
            if n_conf - n_ok_conf > step_conf:
                now_conf = step_conf
            else:
                now_conf = n_conf - n_ok_conf
            with torch.no_grad():
                new_conformations, new_log_det_J = self.generate_conformations(mol_batch, n_conf = now_conf)
                new_conformations = new_conformations[0]
                new_log_det_J = new_log_det_J[0]
            ok_conf_indexes = correct_conf_indexes(new_conformations)
            new_conformations = new_conformations[ok_conf_indexes]
            new_log_det_J = new_log_det_J[ok_conf_indexes]
            conformations = torch.cat( (conformations, new_conformations.cpu()), dim=0)
            log_det_J = torch.cat( (log_det_J, new_log_det_J.cpu()), dim=0 )
            n_ok_conf += len(ok_conf_indexes)
        return conformations, log_det_J

    def generate_histogram(self, mol_index, n_conf = 100, step_conf = 50, visualize = False, file_name = None, save_files = True, verbose = True, ignore_constraints = False):
        '''

        '''
        if save_files:
            if not file_name:
                file_name =  'mol_'+ str(mol_index) + '_' + self.params.model_name
            energy_file_name = 'histograms_data/' + file_name + '_energies.npy'
            log_det_J_file_name = 'histograms_data/' + file_name + '_log_det_J.npy'
            if visualize:
                pdb_writer = PDBWriter('histograms_data/' + file_name + '.pdb', flavor = 4)

        mol_dict = self.dataset[mol_index]
        self.model.eval()
        with torch.no_grad():
            if verbose: print(f'Generating conformations of molecule {mol_index}...', flush=True)
            if ignore_constraints:
                conformations, log_det_J = self.generate_conformations([mol_dict], n_conf)
                conformations = conformations[0]
                log_det_J = conformations[0]
            else:
                conformations, log_det_J = self.generate_correct_conformations(mol_index, n_conf, step_conf)
            if verbose: print(f'Conformations of molecule {mol_index} generated, computing energies...', flush=True)

        forcefield = ForceField('openff-1.2.0.offxml')
        rdkit_mol = mol_dict['rdkit_mol']
        partial_charges = mol_dict['partial_charges']
        ref_atoms = mol_dict['ref_atoms']

        off_mol = Molecule.from_rdkit(rdkit_mol)
        off_mol.partial_charges = partial_charges
        topology = Topology.from_molecules(off_mol)
        openmm_system = forcefield.create_openmm_system(topology, charge_from_molecules = [off_mol], allow_nonintegral_charges=True)

        energies = []

        for i_conf, conformation in enumerate(conformations):
            cartesian_coord = deconstruct_z_matrix(conformation, ref_atoms, loss = False)
            try:
                energy = get_conformations_energy(cartesian_coord, openmm_system)
            except:
                energy = np.nan
                if verbose: print(f'Failed to compute the energy of coformation {i_conf}.')
            if visualize:
                #Set the rdkit molecule to the conformation coordinates
                conformer = rdkit_mol.GetConformer()
                for i_atom in range(rdkit_mol.GetNumAtoms()):
                    x = cartesian_coord[i_atom, 0].item()
                    y = cartesian_coord[i_atom, 1].item()
                    z = cartesian_coord[i_atom, 2].item()
                    conformer.SetAtomPosition(i_atom, Point3D(x,y,z) )
                pdb_writer.write(rdkit_mol)

            energies.append(energy)
        log_det_J = log_det_J.cpu().numpy()
        if save_files:
            np.save(energy_file_name, energies)
            np.save(log_det_J_file_name, log_det_J)

        if verbose: print('Done!', flush=True)

        return energies, log_det_J

max_stars_repo_name: olsson-group/transBG

max_issues_repo_name: olsson-group/transBG

max_forks_repo_name: olsson-group/transBG


text:
# !/usr/bin/env python

'''
ENPM 673 Spring 2019: Robot Perception
Project 3 GMM Bouy Detection

Author:
Ashwin Varghese Kuruttukulam(ashwinvk94@gmail.com)
Rachith Prakash
Graduate Students in Robotics,
University of Maryland, College Park
'''

import glob
from skimage.feature import hog
import numpy as np
import cv2
from sklearn import svm
from sklearn.metrics import classification_report,accuracy_score
import pickle

def showImage(image):
	cv2.namedWindow("Display frame",cv2.WINDOW_NORMAL)
	cv2.imshow("Display frame",image)
	cv2.waitKey(0)

winSize = (64,64)
blockSize = (16,16)
blockStride = (8,8)
cellSize = (8,8)
nbins = 16
derivAperture = 1
winSigma = -1.
histogramNormType = 0
L2HysThreshold = 2.0000000000000001e-01
gammaCorrection = 2
nlevels = 64
SignedGradients = True
hog = cv2.HOGDescriptor(winSize,blockSize,blockStride,cellSize,nbins,derivAperture,winSigma,histogramNormType,L2HysThreshold,gammaCorrection,nlevels,SignedGradients)

TRAINING_IMAGE_SIZE_X = 64
TRAINING_IMAGE_SIZE_Y = 64

trainingLabelShortened = [45, 21, 38, 35, 17, 1, 14, 19]
trainingFeaturesArr = []
trainingLabelsArr= []
nImageCounter = 0

# Getting images in training
pathToImages = 'Training'
foldersTraining = glob.glob(pathToImages+"/*")
# Removing readme file
try:
	foldersTraining.remove(pathToImages+'/Readme.txt')
except:
	print 'Readme file not present'
foldersTraining.sort()

# Getting images in training
pathToImages = 'Testing'
foldersTesting = glob.glob(pathToImages+"/*")
# Removing readme file
try:
	foldersTesting.remove(pathToImages+'/Readme.txt')
except:
	print 'Readme file not present'
foldersTesting.sort()
nImagesCutoff = 1800
nImageCounter = 0
# Iterate through images in each fodlers to compile on list of traffic sign images
for label,folderTraining in enumerate(foldersTraining):
	print nImageCounter
	if label not in trainingLabelShortened:
		continue
	folderTesting = foldersTesting[label]
	print 'label is:',label
	nImageCounter1 = 0
	PositiveImages = glob.glob(folderTraining+"/*.ppm") + glob.glob(folderTesting+"/*.ppm")
	print 'Training and testing folders are'
	print folderTraining
	print folderTesting
	print 'Total number of images in this label is:'
	print len(PositiveImages)
	for imagePath in PositiveImages:
		nImageCounter += 1
		nImageCounter1 += 1

		# load training image
		image = cv2.imread(imagePath)

		# convert to grayscale
		trainImage = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

		# resize
		trainImage = cv2.resize(trainImage, (TRAINING_IMAGE_SIZE_X, TRAINING_IMAGE_SIZE_Y), interpolation=cv2.INTER_AREA)
		# showImage(trainImage)

		# We have to tune these

		# fd = hog.compute(trainImage,winStride,padding,locations)
		fd = hog.compute(trainImage)
		fd = fd.T
		if nImageCounter==1:
			trainingFeaturesArr = fd
			trainingLabelsArr = np.array(1)
		else:
			trainingFeaturesArr = np.vstack((trainingFeaturesArr,fd))
			trainingLabelsArr = np.append(trainingLabelsArr,label)
	print 'Number of images is:',nImageCounter1
nPositiveImageCounter = nImageCounter
nImageCounter = 0
NegativeImages1 = glob.glob("negative_images_1/*.jpg")
NegativeImages2 = glob.glob("negative_images_2/*.jpg")
NegativeImages = NegativeImages1+NegativeImages2
for imagePath in NegativeImages:
	print nImageCounter
	nImageCounter += 1

	# load training image
	image = cv2.imread(imagePath)

	# convert to grayscale
	trainImage = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

	# resize
	trainImage = cv2.resize(trainImage, (TRAINING_IMAGE_SIZE_X, TRAINING_IMAGE_SIZE_Y), interpolation=cv2.INTER_AREA)
	# showImage(trainImage)

	# We have to tune these
	# fd = hog.compute(trainImage,winStride,padding,locations)
	fd = hog.compute(trainImage)
	fd = fd.T
	trainingFeaturesArr = np.vstack((trainingFeaturesArr, fd))
	trainingLabelsArr = np.append(trainingLabelsArr, label+1)
	if nImageCounter>=nImagesCutoff:
		break

print 'Number of positive', nPositiveImageCounter
print 'Number of negative_images', nImageCounter

#Initializing svm classifier
# model = svm.SVC(C=1.0, kernel='rbf', degree=3, gamma='auto_deprecated', coef0=0.0, shrinking=True, probability=False, tol=0.001, cache_size=200, class_weight='balanced')
model = svm.SVC(gamma=0.01,class_weight='balanced')
trainingLabelsArr = trainingLabelsArr.reshape(-1,1)
data_frame = np.hstack((trainingFeaturesArr,trainingLabelsArr))

np.random.shuffle(data_frame)

#Spliting into training and testing data
percentage = 80
partitionIndex = int(trainingLabelsArr.shape[0]*percentage/100)
print partitionIndex

x_train, x_test = data_frame[:partitionIndex,:-1],  data_frame[partitionIndex:,:-1]
y_train, y_test = data_frame[:partitionIndex,-1:].ravel() , data_frame[partitionIndex:,-1:].ravel()

print x_train.shape
print y_train.shape
print y_train
model.fit(x_train,y_train)

y_pred = model.predict(x_test)

# y_pred = model.predict(x_test)
print("Accuracy: "+str(accuracy_score(y_test, y_pred)))
print('\n')
print(classification_report(y_test, y_pred))

filename = 'nine_class_svm.sav'
pickle.dump(model, open(filename, 'wb'))

ext: py

max_stars_repo_name: RachithP/traffic-sign-recognition

max_issues_repo_name: RachithP/traffic-sign-recognition

max_forks_repo_name: RachithP/traffic-sign-recognition

reason: import numpy


text:
import numpy as np
import cross_sections as xs

xs1 = xs.IBeam(1, 1, 0.1, 0.1)

L = 10
p = 1
E = 29000

def constant(x, **kwargs):
	return 1

def linearup(s, **kwargs):
	return x

load = constant

def simpsons(f, a, b, n): #function, start, stop, intervals
	if n % 2 == 0:
		h = (b-a)/n
		k = 0.0
		x = a + h

		for i in range(1, int(n/2) + 1):
			k += 4*f(x)
			x += 2*h

		x = a + 2*h
		for i in range(1, n//2):
			k += 2*f(x)
			x += 2*h
		return (h/3)*(f(a) + f(b) + k)
	else:
		print('n must be even')

I0 = lambda x: p * load(x, L=L)
I1 = lambda x: p * (L-x) * load(x, L=L)
I2 = lambda x: p * (L-x)**2 * load(x, L=L)
I3 = lambda x: p * (L-x)**3 * load(x, L=L)

Int0 = simpsons(I0, 0, L, 100)
Int1 = simpsons(I0, 0, L, 100)
Int2 = simpsons(I0, 0, L, 100)
Int3 = -simpsons(I0, 0, L, 100)

z = np.array([Int0, Int1, Int2, Int3])

a = L
b = L/(E * xs1.Ix)
c = L**2/(2 * E * xs1.Ix)
d = L**3/(6 * E * xs1.Ix)

B = np.matrix([[1, 0, 0, 0, -1, 0, 0, 0],
							 [a, 1, 0, 0, 0, -1, 0, 0],
							 [c, b, 1, 0, 0, 0, -1, 0],
							 [-d, -c, -a, 1, 0, 0, 0, -1]])

fixed = [1, 1, 0, 0]
free = [0, 0, 1, 1]

BC = np.array(fixed + free)

C = B[:, BC==1]

s = np.linalg.solve(C, z)

max_stars_repo_path: Examples/In progress/CEE213 CP4.py

max_stars_repo_name: Komutsou/StructPy

max_issues_repo_path: Examples/In progress/CEE213 CP4.py

max_issues_repo_name: Komutsou/StructPy

max_forks_repo_path: Examples/In progress/CEE213 CP4.py

max_forks_repo_name: Komutsou/StructPy


text:
import numpy as np
from numpy.linalg import solve

x=[19,25,31,38,44]
y=[19.0,32.3,49.0,73.3,97.8]

phi1=[i*i for i in x]
phi0=[1]*5

phi00=5
phi11=sum([i**2 for i in phi1])
phi01=sum(phi1)
yphi0=sum(y)
yphi1=sum([y[i]*phi1[i] for i in range(0,5)])

print(phi00)
print(phi01)
print(phi11)
print(yphi0)
print(yphi1)

a=np.mat([[phi00,phi01],[phi01,phi11]])#系数矩阵
b=np.mat([yphi0,yphi1]).T    #常数项列矩阵
x,y=solve(a,b)        #方程组的解

print(x,"\n",y)

reason: import numpy,from numpy


text:
from sklearn.metrics import mean_squared_error
import math
import os
import gc
import numpy as np
from keras.layers import Dense, LSTM, Activation, Dropout, ActivityRegularization, TimeDistributed, AveragePooling1D, Flatten
from keras.models import Sequential, model_from_json

# fix random seed for reproducibility purposes
# np.random.seed(7)

def getModel(trainX, trainY, testX, testY, model_path, weights_path, epochs=100, batch_size=16, look_back=1):
	'''
	If the model and weights already exist, load them from the cache
	If the model and weights DO NOT exist, then create and fit the model, saving the weights and model
	to the specified model_path and weights_path
	:param trainX: The train data along the x axis
	:param trainY: The train data along the y axis
	:param testX: The test data along the x axis
	:param testY: The test data along the y axis
	:param model_path: The path to either read or save the model
	:param weights_path: The path to either read or save the weights
	:param epochs: The number of epochs used for training (default is 100)
	:param batch_size: The batch size (default is 32)
	:param look_back: The look back value (default 1, should always be the same as the  dimenstion of your trainX data)
	:return:  The model
	'''
	# if model path exists, load it
	if os.path.exists(model_path):
		json_file = open(model_path, 'r')
		loaded_model_json = json_file.read()
		json_file.close()
		model = model_from_json(loaded_model_json)
		print("Loaded model from {}".format(model_path))
		if os.path.exists(weights_path):
			model.load_weights(weights_path)
	else:
		model = Sequential()
		model.add(LSTM(128, activation='relu', input_shape=(1, look_back), return_sequences=True, implementation=1))
		model.add(LSTM(32, input_shape=(1, look_back), return_sequences=True, implementation=1))
		# model.add(LSTM(8, input_shape=(1, look_back), return_sequences=True, implementation=1, activation='tanh'))
		model.add(LSTM(8, input_shape=(1, look_back), return_sequences=False, implementation=1, activation='tanh'))
		model.add(Dense(1, kernel_initializer='normal', activation='tanh'))
		if not os.path.exists(model_path):
			# serialize model to JSON
			model_json = model.to_json()
			with open(model_path, "w") as json_file:
				json_file.write(model_json)
				print("Saved model to {}".format(model_path))

	# must compile the model before training it
	model.compile(loss='mean_squared_error', optimizer='adam', metrics=['mse'])
	# if the weights path exists, use the weights!
	if os.path.exists(weights_path):
		model.load_weights(weights_path)
	else:
		# if no weights exist, train the model
		model.fit(trainX, trainY, epochs=epochs, batch_size=batch_size, validation_data=(testX, testY), verbose=3)
		# serialize weights to HDF5
		model.save_weights(weights_path)
		print("Saved weights to {}".format(weights_path))

	print("Collecting garbage.")
	gc.collect()

	return model

def invert_predictions(trainPredict, testPredict, trainY, testY, scaler):
	trainPredict = scaler.inverse_transform(trainPredict)
	trainY = scaler.inverse_transform([trainY])
	testPredict = scaler.inverse_transform(testPredict)
	testY = scaler.inverse_transform([testY])
	return trainPredict, testPredict, trainY, testY

def scorePrediction(trainPredict, testPredict, trainY, testY):
	# calculate root mean squared error
	trainScore = math.sqrt(mean_squared_error(trainY[0], trainPredict[:, 0]))
	print('Train Score: %.2f RMSE' % (trainScore))
	testScore = math.sqrt(mean_squared_error(testY[0], testPredict[:, 0]))
	print('Test Score: %.2f RMSE' % (testScore))
	return trainPredict, testPredict

def predictFuture(model, initialPredictionData, num_days, scaler):
	predictionXData = initialPredictionData
	predictions = []
	for i in range(0, num_days):
		futurePredict = model.predict(predictionXData)
		predictions.append(scaler.inverse_transform(futurePredict)[0])
		# shift everything to the right
		predictionXData = np.roll(predictionXData, 1)
		# place future predict at beginning of array
		predictionXData[0][0][0] = futurePredict[0][0]

	return predictions

def getOptimalNumHiddenNeurons(training_data_size, num_input_neurons, num_output_neurons):
	return math.floor(training_data_size/(5*(num_input_neurons + num_output_neurons)))

def calculatePercentChange(future_price, current_price):
	if future_price > current_price:
		percent_change = current_price / future_price[0]
	elif future_price < current_price:
		percent_change = -(current_price / future_price[0])
	else:
		percent_change = 0
	return percent_change

def calculateMomentum(future_data, historic_data):
	return future_data - historic_data

def determineAction(momentum, threshold):
	if momentum <= threshold and momentum >= -threshold:
		return 'HOLD'
	elif momentum > 0:
		return 'BUY'
	elif momentum < 0:
		return 'SELL'

max_stars_repo_name: driemworks/agatha

max_forks_repo_name: driemworks/agatha


text:
# *****************************************************************************
# Copyright (c) 2020, Intel Corporation All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#     Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#
#     Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
# EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
# *****************************************************************************

import numba
import numpy as np
import pandas as pd

from sdc.tests.test_base import TestCase
from sdc.tests.test_utils import skip_numba_jit, skip_sdc_jit

DATA = [1.0, 2., 3., 4., 5.]
INDEX = [5, 4, 3, 2, 1]
NAME = "sname"

def series_apply_square_usecase(S):

    def square(x):
        return x ** 2

    return S.apply(square)

class TestSeries_apply(object):

    def test_series_apply(self):
        test_impl = series_apply_square_usecase
        hpat_func = self.jit(test_impl)

        S = pd.Series(DATA)
        pd.testing.assert_series_equal(hpat_func(S), test_impl(S))

    def test_series_apply_convert_type(self):
        def test_impl(S):
            def to_int(x):
                return int(x)
            return S.apply(to_int)
        hpat_func = self.jit(test_impl)

        S = pd.Series(DATA)
        pd.testing.assert_series_equal(hpat_func(S), test_impl(S))

    @skip_sdc_jit("Series.index values are different")
    def test_series_apply_index(self):
        test_impl = series_apply_square_usecase
        hpat_func = self.jit(test_impl)

        S = pd.Series(DATA, INDEX)
        pd.testing.assert_series_equal(hpat_func(S), test_impl(S))

    @skip_sdc_jit('Attribute "name" are different')
    def test_series_apply_name(self):
        test_impl = series_apply_square_usecase
        hpat_func = self.jit(test_impl)

        S = pd.Series(DATA, name=NAME)
        pd.testing.assert_series_equal(hpat_func(S), test_impl(S))

    def test_series_apply_lambda(self):
        def test_impl(S):
            return S.apply(lambda a: 2 * a)
        hpat_func = self.jit(test_impl)

        S = pd.Series(DATA)
        pd.testing.assert_series_equal(hpat_func(S), test_impl(S))

    @skip_sdc_jit("'Var' object has no attribute 'py_func'")
    def test_series_apply_np(self):
        def test_impl(S):
            return S.apply(np.log)
        hpat_func = self.jit(test_impl)

        S = pd.Series(DATA)
        pd.testing.assert_series_equal(hpat_func(S), test_impl(S))

    @skip_sdc_jit("'function' object has no attribute 'py_func'")
    def test_series_apply_register_jitable(self):
        @numba.extending.register_jitable
        def square(x):
            return x ** 2

        def test_impl(S):
            return S.apply(square)
        hpat_func = self.jit(test_impl)

        S = pd.Series(DATA)
        pd.testing.assert_series_equal(hpat_func(S), test_impl(S))

    @skip_numba_jit("'args' in apply is not supported")
    @skip_sdc_jit("'args' in apply is not supported")
    def test_series_apply_args(self):
        @numba.extending.register_jitable
        def subtract_custom_value(x, custom_value):
            return x - custom_value

        def test_impl(S):
            return S.apply(subtract_custom_value, args=(5,))
        hpat_func = self.jit(test_impl)

        S = pd.Series(DATA)
        pd.testing.assert_series_equal(hpat_func(S), test_impl(S))

    @skip_numba_jit("'kwds' in apply is not supported")
    @skip_sdc_jit("'kwds' in apply is not supported")
    def test_series_apply_kwds(self):
        @numba.extending.register_jitable
        def add_custom_values(x, **kwargs):
            for month in kwargs:
                x += kwargs[month]
            return x

        def test_impl(S):
            return S.apply(add_custom_values, june=30, july=20, august=25)
        hpat_func = self.jit(test_impl)

        S = pd.Series(DATA)
        pd.testing.assert_series_equal(hpat_func(S), test_impl(S))

class _Test(TestSeries_apply, TestCase):
    pass

max_stars_repo_name: densmirn/sdc

max_issues_repo_name: densmirn/sdc

max_forks_repo_name: densmirn/sdc

reason: import numpy,import numba


text:
import tools.processing as pre
import numpy as np
from enum import Enum

opt = Enum("Option", "MEAN RANDOM")

glove_path = "data/embeddings/glove.840B.300d.txt"
phonem_path = "data/embeddings/embedding_phonems.txt"
GLOVE_SIZE = 300
PHONEM_SIZE = 3

def get_glove(path_to_glove, wordlist):
    embedding_weights = {}
    count_all_words = 0

    with open(path_to_glove,'r') as f:
        for line in f:
            vals = line.split(' ')
            word = str(vals[0])
            if word in wordlist:
                count_all_words += 1
                coefs = np.asarray(vals[1:],dtype='float32')
                coefs /= np.linalg.norm(coefs)
                embedding_weights[word] = coefs
            if count_all_words== len(wordlist) -1:
                break
    return embedding_weights

def get_embedding(words, option=opt.MEAN):
    embedding = get_glove(glove_path, words)

    unknowns = list(set(words) - set(embedding.keys()))

    mean = np.mean( np.array(list(embedding.values())), axis=0 )

    for word in unknowns:
        embedding[word] = mean

    return embedding

def get_embedding_matrix(word2index_map, vocabulary_size):
    wordlist = list(word2index_map.keys())
    word2embedding_dict = get_embedding(wordlist)
    embedding_matrix = np.zeros((vocabulary_size, GLOVE_SIZE))

    for word, index in word2index_map.items():
        if not word == "PAD_TOKEN":
            word_embedding = word2embedding_dict[word]
            embedding_matrix[index, :] = word_embedding

    return embedding_matrix

def get_phonem(path_to_phonem, phonemlist):
    embedding_weights = {}
    count_all_words = 0

    with open(path_to_phonem, 'r') as f:
        for line in f:
            vals = line.split(' ')
            word = str(vals[0])
            if word in phonemlist:
                count_all_words += 1
                coefs = np.asarray(vals[1:],dtype='float32') # already normalized
                embedding_weights[word] = coefs
            if count_all_words == len(phonemlist) -1:
                break
    return embedding_weights

def get_phonem_embedding(phonems):
    embedding = get_phonem(phonem_path, phonems)

    unknowns = list(set(phonems) - set(embedding.keys()))

    mean = np.mean( np.array(list(embedding.values())), axis=0 )

    for phonem in unknowns:
        embedding[phonem] = mean

    return embedding

def get_phonem_embedding_matrix(phonem2index_map, vocabulary_size):
    phonemlist = list(phonem2index_map.keys())
    word2embedding_dict = get_phonem_embedding(phonemlist)
    embedding_matrix = np.zeros((vocabulary_size, PHONEM_SIZE))

    for phonem, index in phonem2index_map.items():
        if not phonem == "PAD_TOKEN":
            phonem_embedding = word2embedding_dict[phonem]
            embedding_matrix[index, :] = phonem_embedding

    return embedding_matrix

def get_closest_words(embedded_word, embedding_dict, limit=10):

    weights = np.array(list(embedding_dict.values()))

    cosine_dists = np.dot( weights, embedded_word)
    ff = np.argsort(cosine_dists)[::-1][:10]

    words = list(embedding_dict.keys())
    index2word_map = dict( (key, value) for (key, value) in enumerate(words) )

    words = []
    distance = []
    for f in ff:
        words.append(index2word_map[f])
        distance.append(cosine_dists[f])

    return list(zip(words, distance))

def get_closest_word(embedding_word, embedding_dict):

    words = get_closest_words(embedding_word, embedding_dict, limit=1)

    return words[0][0]

max_stars_repo_name: frankzl/deep-rap

max_issues_repo_name: frankzl/deep-rap

max_forks_repo_name: frankzl/deep-rap


text:
import numpy as np

def run_viterbi(initial_scores, transition_scores, final_scores, emission_scores):
    length = np.size(emission_scores, 0) # Length of the sequence.
    num_states = np.size(initial_scores) # Number of states.

    # Variables storing the Viterbi scores.
    viterbi_scores = np.zeros([length, num_states])

    # Variables storing the paths to backtrack.
    viterbi_paths = -np.ones([length, num_states], dtype=int)

    # Most likely sequence.
    best_path = -np.ones(length, dtype=int)

    # Initialization.
    viterbi_scores[0,:] = emission_scores[0,:] * initial_scores

    # Viterbi loop.
    for pos in xrange(1,length):
        for current_state in xrange(num_states):
            viterbi_scores[pos, current_state] = \
                np.max(viterbi_scores[pos-1, :] * transition_scores[pos-1, current_state, :])
            viterbi_scores[pos, current_state] *= emission_scores[pos, current_state]
            viterbi_paths[pos, current_state] = \
                np.argmax(viterbi_scores[pos-1, :] * transition_scores[pos-1, current_state, :])

    # Termination.
    best_score = np.max(viterbi_scores[length-1,:] * final_scores)
    best_path[length-1] = np.argmax(viterbi_scores[length-1,:] * final_scores)

    # Backtrack.
    for pos in xrange(length-2, -1, -1):
        best_path[pos] = viterbi_paths[pos+1, best_path[pos+1]]

    return best_path, best_score

######
# Computes the viterbi path for a given sequence of lenght.
# N - Lenght of sequence
# H - Number of hidden states
# Receives:
# Node potentials (N,H) vector
# Edge potentials (N-1,H,H)
######
#def viterbi(node_potentials,edge_potentials):
#    H,N = node_potentials.shape
#    max_marginals = np.zeros([H,N])
#
#    ## Fil backpointers with zero to signal an errors
#    backpointers = np.zeros([H,N],dtype="int")
#    backpointers[:,1:] = -1
#    max_marginals[:,0] = node_potentials[:,0]
#    for pos in xrange(1,N):
#        max_value = 0
#
#        for current_state in xrange(H):
#            max_prob = -1
#            max_state = -1
#
#            for prev_state in xrange(H):
#                viter_v = max_marginals[prev_state,pos-1]
#                trans_v = edge_potentials[prev_state,current_state,pos-1]
#                prob = viter_v*trans_v
#                if(prob > max_prob):
#                    max_prob = prob
#                    max_state = prev_state
#            max_marginals[current_state,pos] = max_prob*node_potentials[current_state,pos]
#            backpointers[current_state,pos] = max_state
#    ##Recover the viterbi state
#    viterbi_path = np.zeros([N,1],dtype="int")
#    best = np.argmax(max_marginals[:,N-1])
#    viterbi_path[N-1] = best
#    for pos in xrange(N-2,-1,-1):
##        print viterbi_path
##        print viterbi_path[pos+1]
##        print "back"
##        print backpointers
#        viterbi_path[pos] = backpointers[viterbi_path[pos+1],pos+1]
#    return viterbi_path,max_marginals

max_stars_repo_name: SimonSuster/lxmls-toolkit

max_issues_repo_name: daviddao/LxMLS-labs-solution

max_forks_repo_name: daviddao/LxMLS-labs-solution


text:
import torch
import matplotlib.pyplot as plt
import numpy as np
import train
import argparse
import onmt.Constants
import matplotlib.ticker as ticker

def main(res):
    for k,v in res.items():
        if v['nse']['score'] is None:
            continue
        print(' prediction ::: '+' '.join(v['nse']['pred']))
        utt_locs = v['nse']['utt_attn'].data.transpose(0,1)[:len(v['nse']['pred'])].transpose(0,1).tolist()
        fig = plt.figure()

        plt.ion()

        #ax1 = fig.add_subplot(111)
        plt.imshow(utt_locs)

        pred = v['nse']['pred']
        #ax1.set_xlim(0,len(pred))
        plt.xlabel(' '.join(v['nse']['base']))
        plt.ylabel(' t ')
        plt.xticks(range(len(pred)),pred)
        #ax1.set_xticks(range(len(pred)))
        #ax1.xaxis.set_major_locator(ticker.MultipleLocator(1))
        #ax1.set_xticklabels(pred)

        #ax2 = ax1.twiny()
        #ax2.plt.matshow(utt_locs)
        #ax2.set_xlim(0,len(pred))
        #ax2.set_xlabel(' '.join(v['nse']['base']))#' base reply ')

        #ticks = v['nse']['base']
        #print(' === len ticks ==' + str(len(ticks)))
        #print(' === len pred ==' + str(len(pred)))
        #add_nulls = len(pred)-len(ticks)

        #if  add_nulls > 0:
        #    ticks += ['_' for i in range(add_nulls)]

        #print(' === len ticks ==' + str(len(ticks)))

        #ax2.set_xticks(range(len(ticks)))
        #ax2.xaxis.set_major_locator(ticker.MultipleLocator(1))
        #ax2.set_xticklabels(ticks)

        plt.show()

        print( ' size prediction : ' +str(len(v['nse']['pred'])))
        print( ' size locations : ' +str(len(utt_locs[0])))
        print(' utt_locs len : %d  --- utt_locs[0] len : %d '%(len(utt_locs),len(utt_locs[0])))

        ch= input('\n --> :')
        if ch == 'q':
            break
        plt.close()

if __name__ == '__main__':
    #res_file = 'predictions/utt_context.pt'
    res_file = 'predictions/nse_keys_cs3.pt'
    main(torch.load(res_file))


text:
# Copyright (c) 2019, MD2K Center of Excellence
# - Nasir Ali <nasir.ali08@gmail.com>
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

import numpy as np
import pandas as pd
import json
from geopy.distance import great_circle
from pyspark.sql.functions import pandas_udf, PandasUDFType
from shapely.geometry.multipoint import MultiPoint
from sklearn.cluster import DBSCAN
from pyspark.sql.types import StructField, StructType, StringType, FloatType
from pyspark.sql.types import StructField, StructType, StringType, FloatType, TimestampType, IntegerType,DateType,TimestampType
from cerebralcortex.algorithms.ema.ema_incentive_features import get_ema_incentive_features
schema = StructType([
    StructField("timestamp", TimestampType()),
    StructField("localtime", TimestampType()),
    StructField("user", StringType()),
    StructField("version", IntegerType()),
    StructField("status", StringType()),
    StructField("ema_id", StringType()),
    StructField("schedule_timestamp", TimestampType()),
    StructField("operation", StringType())
])

@pandas_udf(schema, PandasUDFType.GROUPED_MAP)
def get_ema_log_features(user_data):
    all_vals = []
    for index, row in user_data.iterrows():
        ema = row["log"]
        if not isinstance(ema, dict):
            ema = json.loads(ema)

        operation = ema["operation"].lower()
        if operation !="condition":
            status = ema.get("status","")
            ema_id = ema["id"]
            schedule_timestamp = ema.get("logSchedule",{}).get("scheduleTimestamp")
            if schedule_timestamp:
                schedule_timestamp = pd.to_datetime(schedule_timestamp, unit='ms')

            all_vals.append([row["timestamp"],row["localtime"], row["user"],1,status,ema_id,schedule_timestamp,operation])

    return pd.DataFrame(all_vals,columns=['timestamp','localtime', 'user', 'version','status','ema_id','schedule_timestamp','operation'])

max_stars_repo_name: brinnaebent/CerebralCortex-Kernel

max_issues_repo_name: brinnaebent/CerebralCortex-Kernel

max_forks_repo_name: brinnaebent/CerebralCortex-Kernel