import pickle
import glob
from tqdm import tqdm
import torch
import torch.nn as nn
import torch.nn.functional as F
import numpy as np
import torch
from torch.utils.data import Dataset, DataLoader
from tqdm import tqdm
from scipy.interpolate import splev
import logging

path = '/gpfs/mvoronina/MT/Project/Seq/sequences_7959_sequences.pt'
data = torch.load(path)

labels = torch.load('/gpfs/mvoronina/MT/Project/Seq/sequences_labels.pt')
labels_part = labels[-459:]


def bs(x, df=None, knots=None, degree=3, intercept=False):
    """
    df : int
        The number of degrees of freedom to use for this spline. The
        return value will have this many columns. You must specify at least
        one of `df` and `knots`.
    knots : list(float)
        The interior knots of the spline. If unspecified, then equally
        spaced quantiles of the input data are used. You must specify at least
        one of `df` and `knots`.
    degree : int
        The degree of the piecewise polynomial. Default is 3 for cubic splines.
    intercept : bool
        If `True`, the resulting spline basis will span the intercept term
        (i.e. the constant function). If `False` (the default) then this
        will not be the case, which is useful for avoiding overspecification
        in models that include multiple spline terms and/or an intercept term.

    """

    order = degree + 1
    inner_knots = []
    if df is not None and knots is None:
        n_inner_knots = df - order + (1 - intercept)
        if n_inner_knots < 0:
            n_inner_knots = 0
            print("df was too small; have used %d"
                  % (order - (1 - intercept)))

        if n_inner_knots > 0:
            inner_knots = np.percentile(
                x, 100 * np.linspace(0, 1, n_inner_knots + 2)[1:-1])

    elif knots is not None:
        inner_knots = knots

    all_knots = np.concatenate(
        ([np.min(x), np.max(x)] * order, inner_knots))

    all_knots.sort()

    n_basis = len(all_knots) - (degree + 1)
    basis = np.empty((x.shape[0], n_basis), dtype=float)

    for i in range(n_basis):
        coefs = np.zeros((n_basis,))
        coefs[i] = 1
        basis[:, i] = splev(x, (all_knots, coefs, degree))

    if not intercept:
        basis = basis[:, 1:]
    return basis

class BSplineTransformation(nn.Module):

    def __init__(self, degrees_of_freedom, log=False, scaled=False):
        super(BSplineTransformation, self).__init__()
        self._spline_tr = None
        self._log = log
        self._scaled = scaled
        self._df = degrees_of_freedom

    def forward(self, input):
        if self._spline_tr is None:
            spatial_dim = input.size()[-1]
            self._spline_tr = spline_factory(spatial_dim, self._df, log=self._log)
            if self._scaled:
                self._spline_tr = self._spline_tr / spatial_dim
            if input.is_cuda:
                self._spline_tr = self._spline_tr.cuda()
        
        return  torch.matmul(input, self._spline_tr)


def spline_factory(n, df, log=False):
    if log:
        dist = np.array(np.arange(n) - n/2.0)
        dist = np.log(np.abs(dist) + 1) * ( 2*(dist>0)-1)
        n_knots = df - 4
        knots = np.linspace(np.min(dist),np.max(dist),n_knots+2)[1:-1]
        return torch.from_numpy(bs(
            dist, knots=knots, intercept=True)).float()
    else:
        dist = np.arange(n)
        return torch.from_numpy(bs(
            dist, df=df, intercept=True)).float()    
    
class Sei(nn.Module):
    def __init__(self, n_classes=3):
        super(Sei, self).__init__()

        self.lconv1 = nn.Sequential(
            nn.Conv1d(4, 860, kernel_size=9, padding=4),
            nn.Conv1d(860, 860, kernel_size=9, padding=4))

        self.conv1 = nn.Sequential(
            nn.Conv1d(860, 860, kernel_size=9, padding=4),
            nn.ReLU(inplace=True),
            nn.Conv1d(860, 860, kernel_size=9, padding=4),
            nn.ReLU(inplace=True))

        self.lconv2 = nn.Sequential(
            nn.MaxPool1d(kernel_size=4, stride=4),
            nn.Dropout(p=0.2),
            nn.Conv1d(860, 480, kernel_size=9, padding=4),
            nn.Conv1d(480, 480, kernel_size=9, padding=4))

        self.conv2 = nn.Sequential(
            nn.Dropout(p=0.2),
            nn.Conv1d(480, 480, kernel_size=9, padding=4),
            nn.ReLU(inplace=True),
            nn.Conv1d(480, 480, kernel_size=9, padding=4),
            nn.ReLU(inplace=True))

        self.lconv3 = nn.Sequential(
            nn.MaxPool1d(kernel_size=4, stride=4),
            nn.Dropout(p=0.2),
            nn.Conv1d(480, 128, kernel_size=9, padding=4),
            nn.Conv1d(128, 128, kernel_size=9, padding=4))

        self.conv3 = nn.Sequential(
            nn.Dropout(p=0.2),
            nn.Conv1d(128, 128, kernel_size=9, padding=4),
            nn.ReLU(inplace=True),
            nn.Conv1d(128, 128, kernel_size=9, padding=4),
            nn.ReLU(inplace=True))

        self.dconv1 = nn.Sequential(
            nn.Dropout(p=0.10),
            nn.Conv1d(128, 128, kernel_size=5, dilation=2, padding=4),
            nn.ReLU(inplace=True))
        self.dconv2 = nn.Sequential(
            nn.Dropout(p=0.10),
            nn.Conv1d(128, 128, kernel_size=5, dilation=4, padding=8),
            nn.ReLU(inplace=True))
        self.dconv3 = nn.Sequential(
            nn.Dropout(p=0.10),
            nn.Conv1d(128, 128, kernel_size=5, dilation=8, padding=16),
            nn.ReLU(inplace=True))
        self.dconv4 = nn.Sequential(
            nn.Dropout(p=0.10),
            nn.Conv1d(128, 128, kernel_size=5, dilation=16, padding=32),
            nn.ReLU(inplace=True))
        self.dconv5 = nn.Sequential(
            nn.Dropout(p=0.10),
            nn.Conv1d(128, 128, kernel_size=5, dilation=25, padding=50),
            nn.ReLU(inplace=True))

        self._spline_df = int(128/8)        
        self.spline_tr = nn.Sequential(
            nn.Dropout(p=0.5),
            BSplineTransformation(self._spline_df, scaled=False))

        self.classifiers = nn.ModuleList([nn.Sequential(
            nn.Linear(128 * self._spline_df, 3),
            nn.ReLU(inplace=True),
            nn.Linear(3, 3),
            nn.Softmax(dim=1)) for _ in range(12)])

    def forward(self, x):
        """Forward propagation of a batch.
        """
        lout1 = self.lconv1(x)
        out1 = self.conv1(lout1)

        lout2 = self.lconv2(out1 + lout1)
        out2 = self.conv2(lout2)

        lout3 = self.lconv3(out2 + lout2)
        out3 = self.conv3(lout3)

        dconv_out1 = self.dconv1(out3 + lout3)
        cat_out1 = out3 + dconv_out1
        dconv_out2 = self.dconv2(cat_out1)
        cat_out2 = cat_out1 + dconv_out2
        dconv_out3 = self.dconv3(cat_out2)
        cat_out3 = cat_out2 + dconv_out3
        dconv_out4 = self.dconv4(cat_out3)
        cat_out4 = cat_out3 + dconv_out4
        dconv_out5 = self.dconv5(cat_out4)
        out = cat_out4 + dconv_out5
        
        spline_out = self.spline_tr(out)
        reshape_out = spline_out.view(spline_out.size(0), 128 * self._spline_df)
        
        predictions = []
        for classifier in self.classifiers:
            prediction = classifier(reshape_out)
            predictions.append(prediction)
        
        return torch.stack(predictions, dim=1)


def get_optimizer(lr):
    """
    The optimizer and the parameters with which to initialize the optimizer.
    At a later time, we initialize the optimizer by also passing in the model
    parameters (`model.parameters()`). We cannot initialize the optimizer
    until the model has been initialized.
    """
    return (torch.optim.SGD,
            {"lr": lr, "weight_decay": 1e-7, "momentum": 0.9})


def one_hot_encode_sequence(sequence):
    mapping = {'A': [1, 0, 0, 0], 'T': [0, 1, 0, 0], 'G': [0, 0, 1, 0], 'C': [0, 0, 0, 1]}
    encoded_sequence = []
    for base in sequence:
        encoded_sequence.append(mapping[base])
    return torch.tensor(encoded_sequence, dtype=torch.float).reshape((1, 4, len(encoded_sequence)))




class SequenceDataset(torch.utils.data.Dataset):
    
    def __init__(self, data_file, label_file, logger, save_file_prefix=None):
        self.sequences = []
        self.labels = []
        with open(label_file, 'r') as f:
            next(f)
            
            for line in f:
                values = line.strip().split(",")
                encoded_labels = []
                for value in values[1:]:
                    encoded_labels.append(one_hot_dict[value])
                self.labels.append(torch.from_numpy(np.array(encoded_labels)).float())
                
        try:
            with open(data_file, 'r') as f:
                i = 0
                for line in f:
                    self.sequences.append(one_hot_encode_sequence(line.rstrip()))
                    i+=1
                    logger.info(f"{i} line proceeded")


            if save_file_prefix is not None:
                torch.save(self.sequences, f'{save_file_prefix}_sequences.pt')
                torch.save(self.labels, f'{save_file_prefix}_labels.pt')
        except:
            logger.error("An error occurred:", exc_info=True)
        
    def __getitem__(self, idx):
        sequence = self.sequences[idx]
        label = self.labels[idx]
        return sequence, label

    def __len__(self):
        return len(self.labels)
    
    
one_hot_dict = {"": [1, 0, 0], '0.0': [0, 1, 0], '1.0': [0, 0, 1]}

logger = logging.getLogger("./SEI_run")
logger.setLevel(logging.INFO)

# create the logging file handler
fh = logging.FileHandler("./SEI_run")

formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
fh.setFormatter(formatter)

# add handler to logger object
logger.addHandler(fh)

model = Sei()
optimizer_class, optimizer_params = get_optimizer(lr=0.01)
optimizer = optimizer_class(model.parameters(), **optimizer_params)

epochs=100

model.train()
logger.info(f'Start training: n_epochs={epochs}')
try:
    for epoch in range(epochs):
        running_loss = 0.0
        for i in tqdm(range(len(data))):
            sequence, label = data[i], labels_part[i]
            optimizer.zero_grad()
            output = model.forward(sequence)
            loss = F.binary_cross_entropy_with_logits(output.squeeze(0), label)
            loss.backward()
            optimizer.step()
            running_loss += loss.item()
            logger.info(f"{i} sequence, {epoch} epoch: loss = {loss.item()}")

        logger.info('Epoch %d loss: %.3f' % (epoch + 1, running_loss / len(dataset)))
        torch.save(model.state_dict(), f'./model_epoch_{epoch+1}.pt')
                        
except:
    logger.error("An error occurred:", exc_info=True)               
                        

Faster the code. Write optimal version