# adopted from
# https://github.com/openai/improved-diffusion/blob/main/improved_diffusion/gaussian_diffusion.py
# and
# https://github.com/lucidrains/denoising-diffusion-pytorch/blob/7706bdfc6f527f58d33f84b7b522e61e6e3164b3/denoising_diffusion_pytorch/denoising_diffusion_pytorch.py
# and
# https://github.com/openai/guided-diffusion/blob/0ba878e517b276c45d1195eb29f6f5f72659a05b/guided_diffusion/nn.py
#
# thanks!


import math
import torch
import torch.nn as nn
import numpy as np
from einops import repeat

from ldm.util import instantiate_from_config

class MPCExp(nn.Module):
    def __init__(self, iters=8):
        super().__init__()
        self.iters = iters

    def forward(self, x):
        b0 = (x < - 14.0).to(torch.float32)
        result = (1 + (x / 2**self.iters)) * (1-b0)
        for _ in range(self.iters):
            result = result.square()
        return result
    
class MPCReciprocal(nn.Module):
    def __init__(self, exp_iters=8, nr_iters=10):
        super().__init__()
        self.nr_iters = nr_iters
        self.exp = MPCExp(iters=exp_iters)

    def forward(self, x):
        result = 3 * self.exp(1 - 2 * x) + 0.003
        for _ in range(self.nr_iters):
            result += result - result.square() * x
        return result

class MPCInvsqrt(nn.Module):
    def __init__(self, iters=4, exp_iters=8):
        super().__init__()
        self.iters = iters
        self.exp = MPCExp(iters=exp_iters)

    def forward(self, x):
        y = self.exp(-((x/2)+0.2))*2.2 + 0.2
        y = y - x / 1024
        for _ in range(self.iters):
            y = y * (3 - x * y.square()) / 2
        return y

class MPCExpLoose(nn.Module):
    def __init__(self, iters=8):
        super().__init__()
        self.iters = iters
        coef = [0.00488676, 0.06755453, 0.36066246, 0.91766503, 0.98805643]
        self.coef = coef[::-1]  # reverse order

    def forward(self, x):
        b0 = (x > - 4.8).to(x.dtype)
        x2 = x.square()
        seg1 = (self.coef[4] * x2 + self.coef[3] * x + self.coef[2]) * x2 + self.coef[1] * x + self.coef[0]
        ret = seg1 * b0
        return ret

class MPCSoftmax(nn.Module):
    def __init__(self, exp_iters=8, reci_exp_iters=8, nr_iters=10, dim=-1, layer_idx=28):
        super().__init__()
        self.dim = dim
        self.exp = MPCExp()
        self.reciprocal = MPCReciprocal(exp_iters=reci_exp_iters, nr_iters=nr_iters)
        self.mpc_mode = True

    def forward(self, x):
        if self.mpc_mode:
            delta = 1 / (math.sqrt(x.size(self.dim)) * 2**3)
            x_max = x.max(dim=self.dim, keepdim=True)[0]
            x = x - x_max
            x_exp = self.exp(x) #.div(math.sqrt(x.size(self.dim)))
            x_sum = x_exp.sum(dim=self.dim, keepdim=True)  * delta
            x_sum = self.reciprocal(x_sum) * delta
            x_exp = x_exp * x_sum
            return x_exp
        else:
            return torch.nn.functional.softmax(x, dim=self.dim)

class MPCSiLU(nn.Module):
    def __init__(self):
        super().__init__()
        self.a_coeffs = [-0.3067541139982155,  -0.0819767021525476, -0.0055465625580307]
        self.b_coeffs = [0.0085064025895951, 0.5, 0.2281430841728270, -0.011113046708173, 0.0002743776353465]
        self.mpc_mode = True
    
    def forward(self, x):
        if self.mpc_mode:
            b0 = (x < -8.0).to(torch.float32)
            b1 = (x < -4.0).to(torch.float32)
            b2 = (x > 4.0).to(torch.float32)
            b3 = (1-b1) * (1-b2)
            b4 = (1-b0) * b1  # x in [-8.0, -4.0)

            x2 = x * x
            seg1 = self.a_coeffs[2] * x2 + self.a_coeffs[1] * x + self.a_coeffs[0]
            seg2 = ((self.b_coeffs[4] * x2 + self.b_coeffs[3]) * x2 + self.b_coeffs[2]) * x2 + self.b_coeffs[1] * x + self.b_coeffs[0]
            ret = b2 * x + b4 * seg1 + b3 * seg2
            return ret
        else:
            return x * torch.sigmoid(x)

class MPCSiLULoose(nn.Module):
    def __init__(self):
        super().__init__()
        coef = [-5.28080234e-06,  3.97939016e-04, -1.18558684e-02,  2.28920817e-01,  5.00011828e-01,  8.82801754e-03]
        self.coef = coef[::-1]
        self.mpc_mode = True
    
    def forward(self, x):
        if self.mpc_mode:
            b0 = (x > -5.4).to(x.dtype)
            b1 = (x > 4.0).to(x.dtype)
            b2 = b0 * (1-b1)

            x2 = x.square()
            x4 = x2.square()
            seg1 = (self.coef[5] * x4 + self.coef[4] * x2 + self.coef[3]) * x4 + self.coef[2] * x2 + self.coef[1] * x + self.coef[0]
            ret = b1 * x + b2 * seg1
            return ret
        else:
            return x * torch.sigmoid(x)

class MPCGeLU(nn.Module): # poly 6
    def __init__(self):
        super().__init__()
        self.a_coeffs = [-0.5054031199708174, -0.4222658115198386, -0.1180761295118195, -0.0110341340306157] # 0, 1, 2, 3
        self.b_coeffs = [0.0085263215410380, 0.5,  0.3603292692789629, -0.037688200365904, 0.0018067462606141] # 0, 1, 2, 4, 6
        self.mpc_mode = True

    def forward(self, x):
        if self.mpc_mode:
            b0 = (x < -5.0).to(torch.float32)
            b1 = (x < -1.97).to(torch.float32)
            b2 = (x > 3.0).to(torch.float32)
            b3 = (1-b1) * (1-b2)
            b4 = (1-b0) * b1

            x2 = x * x
            seg1 = (self.a_coeffs[3] * x + self.a_coeffs[2]) * x2 + self.a_coeffs[1] * x + self.a_coeffs[0]
            seg2 = ((self.b_coeffs[4] * x2 + self.b_coeffs[3]) * x2 + self.b_coeffs[2]) * x2 + self.b_coeffs[1] * x + self.b_coeffs[0]
            ret = b2 * x + b4 * seg1 + b3 * seg2
            return ret
        else:
            return nn.functional.gelu(x)
        
class MPCGeLULoose(nn.Module): # seg 1 # deg 4 (4, 2, 1, 0), c*c=2, p*c=1
    def __init__(self):
        super().__init__()
        coef = [-0.01625818,  0.30074739,  0.49693733,  0.03222096]
        self.coef = coef[::-1]  # reverse order
        self.mpc_mode = True

    def forward(self, x):
        if self.mpc_mode:
            b0 = (x > -2.1).to(x.dtype)
            b1 = (x > 3.0).to(x.dtype)
            b2 = b0 * (1-b1)

            x2 = x.square()
            seg = (self.coef[3] * x2 + self.coef[2]) * x2 + self.coef[1] * x + self.coef[0]
            ret = x * b1 + seg * b2 
            return ret
        else:
            return nn.functional.gelu(x)


class MPCGroupNorm(nn.GroupNorm):
    __constants__ = ['num_groups', 'num_channels', 'eps', 'affine']
    num_groups: int
    num_channels: int
    eps: float
    affine: bool

    def __init__(self, num_groups: int, num_channels: int, eps: float = 1e-5, affine: bool = True,
                 device=None, dtype=None) -> None:
        factory_kwargs = {'device': device, 'dtype': dtype}
        super().__init__(num_groups, num_channels)
        if num_channels % num_groups != 0:
            raise ValueError('num_channels must be divisible by num_groups')

        self.num_groups = num_groups
        self.num_channels = num_channels
        self.eps = eps
        self.affine = affine

        self.invsqrt = MPCInvsqrt(iters=6)
        self.mpc_mode = True

    def forward(self, x):
        if self.mpc_mode:
            orig_shape = x.shape
            x = x.reshape(orig_shape[0], self.num_groups, -1)

            denom = min(np.sqrt(x.size(-1) / (122880)), 1)

            x = x.mul(denom)
            mean = x.mean(dim=-1, keepdim=True)
            var = (x - mean).square().mean(dim=-1, keepdim=True)

            x = (x - mean) * self.invsqrt(var + self.eps)
            x = x.reshape(orig_shape)

            weight = self.weight.reshape(1,-1,1,1).to(torch.float32)
            bias = self.bias.reshape(1,-1,1,1).to(torch.float32)
            
            x = x * weight + bias

            return x
    
        else:
            return nn.functional.group_norm(x, self.num_groups, self.weight, self.bias, self.eps)

class MPCGroupNorm32(MPCGroupNorm):
    def forward(self, x):
        if self.mpc_mode:
            return super().forward(x)
        else:
            return super().forward(x.float()).type(x.dtype)


def make_beta_schedule(schedule, n_timestep, linear_start=1e-4, linear_end=2e-2, cosine_s=8e-3):
    if schedule == "linear":
        betas = (
                torch.linspace(linear_start ** 0.5, linear_end ** 0.5, n_timestep, dtype=torch.float64) ** 2
        )

    elif schedule == "cosine":
        timesteps = (
                torch.arange(n_timestep + 1, dtype=torch.float64) / n_timestep + cosine_s
        )
        alphas = timesteps / (1 + cosine_s) * np.pi / 2
        alphas = torch.cos(alphas).pow(2)
        alphas = alphas / alphas[0]
        betas = 1 - alphas[1:] / alphas[:-1]
        betas = np.clip(betas, a_min=0, a_max=0.999)

    elif schedule == "sqrt_linear":
        betas = torch.linspace(linear_start, linear_end, n_timestep, dtype=torch.float64)
    elif schedule == "sqrt":
        betas = torch.linspace(linear_start, linear_end, n_timestep, dtype=torch.float64) ** 0.5
    else:
        raise ValueError(f"schedule '{schedule}' unknown.")
    return betas.numpy()


def make_ddim_timesteps(ddim_discr_method, num_ddim_timesteps, num_ddpm_timesteps, verbose=True):
    if ddim_discr_method == 'uniform':
        c = num_ddpm_timesteps // num_ddim_timesteps
        ddim_timesteps = np.asarray(list(range(0, num_ddpm_timesteps, c)))
    elif ddim_discr_method == 'quad':
        ddim_timesteps = ((np.linspace(0, np.sqrt(num_ddpm_timesteps * .8), num_ddim_timesteps)) ** 2).astype(int)
    else:
        raise NotImplementedError(f'There is no ddim discretization method called "{ddim_discr_method}"')

    # add one to get the final alpha values right (the ones from first scale to data during sampling)
    steps_out = ddim_timesteps + 1
    if verbose:
        print(f'Selected timesteps for ddim sampler: {steps_out}')
    return steps_out


def make_ddim_sampling_parameters(alphacums, ddim_timesteps, eta, verbose=True):
    # select alphas for computing the variance schedule
    alphas = alphacums[ddim_timesteps]
    alphas_prev = np.asarray([alphacums[0]] + alphacums[ddim_timesteps[:-1]].tolist())

    # according the the formula provided in https://arxiv.org/abs/2010.02502
    sigmas = eta * np.sqrt((1 - alphas_prev) / (1 - alphas) * (1 - alphas / alphas_prev))
    sigmas_orig = np.sqrt((1 - alphas_prev) / (1 - alphas) * (1 - alphas / alphas_prev))
    if verbose:
        print(f'Selected alphas for ddim sampler: a_t: {alphas}; a_(t-1): {alphas_prev}')
        print(f'For the chosen value of eta, which is {eta}, '
              f'this results in the following sigma_t schedule for ddim sampler {sigmas}')
    return sigmas, alphas, alphas_prev, sigmas_orig


def betas_for_alpha_bar(num_diffusion_timesteps, alpha_bar, max_beta=0.999):
    """
    Create a beta schedule that discretizes the given alpha_t_bar function,
    which defines the cumulative product of (1-beta) over time from t = [0,1].
    :param num_diffusion_timesteps: the number of betas to produce.
    :param alpha_bar: a lambda that takes an argument t from 0 to 1 and
                      produces the cumulative product of (1-beta) up to that
                      part of the diffusion process.
    :param max_beta: the maximum beta to use; use values lower than 1 to
                     prevent singularities.
    """
    betas = []
    for i in range(num_diffusion_timesteps):
        t1 = i / num_diffusion_timesteps
        t2 = (i + 1) / num_diffusion_timesteps
        betas.append(min(1 - alpha_bar(t2) / alpha_bar(t1), max_beta))
    return np.array(betas)


def extract_into_tensor(a, t, x_shape):
    b, *_ = t.shape
    out = a.gather(-1, t)
    return out.reshape(b, *((1,) * (len(x_shape) - 1)))


def checkpoint(func, inputs, params, flag):
    """
    Evaluate a function without caching intermediate activations, allowing for
    reduced memory at the expense of extra compute in the backward pass.
    :param func: the function to evaluate.
    :param inputs: the argument sequence to pass to `func`.
    :param params: a sequence of parameters `func` depends on but does not
                   explicitly take as arguments.
    :param flag: if False, disable gradient checkpointing.
    """
    if flag:
        args = tuple(inputs) + tuple(params)
        return CheckpointFunction.apply(func, len(inputs), *args)
    else:
        return func(*inputs)


class CheckpointFunction(torch.autograd.Function):
    @staticmethod
    def forward(ctx, run_function, length, *args):
        ctx.run_function = run_function
        ctx.input_tensors = list(args[:length])
        ctx.input_params = list(args[length:])

        with torch.no_grad():
            output_tensors = ctx.run_function(*ctx.input_tensors)
        return output_tensors

    @staticmethod
    def backward(ctx, *output_grads):
        ctx.input_tensors = [x.detach().requires_grad_(True) for x in ctx.input_tensors]
        with torch.enable_grad():
            # Fixes a bug where the first op in run_function modifies the
            # Tensor storage in place, which is not allowed for detach()'d
            # Tensors.
            shallow_copies = [x.view_as(x) for x in ctx.input_tensors]
            output_tensors = ctx.run_function(*shallow_copies)
        input_grads = torch.autograd.grad(
            output_tensors,
            ctx.input_tensors + ctx.input_params,
            output_grads,
            allow_unused=True,
        )
        del ctx.input_tensors
        del ctx.input_params
        del output_tensors
        return (None, None) + input_grads


def timestep_embedding(timesteps, dim, max_period=10000, repeat_only=False):
    """
    Create sinusoidal timestep embeddings.
    :param timesteps: a 1-D Tensor of N indices, one per batch element.
                      These may be fractional.
    :param dim: the dimension of the output.
    :param max_period: controls the minimum frequency of the embeddings.
    :return: an [N x dim] Tensor of positional embeddings.
    """
    if not repeat_only:
        half = dim // 2
        freqs = torch.exp(
            -math.log(max_period) * torch.arange(start=0, end=half, dtype=torch.float32) / half
        ).to(device=timesteps.device)
        args = timesteps[:, None].float() * freqs[None]
        embedding = torch.cat([torch.cos(args), torch.sin(args)], dim=-1)
        if dim % 2:
            embedding = torch.cat([embedding, torch.zeros_like(embedding[:, :1])], dim=-1)
    else:
        embedding = repeat(timesteps, 'b -> b d', d=dim)
    return embedding


def zero_module(module):
    """
    Zero out the parameters of a module and return it.
    """
    for p in module.parameters():
        p.detach().zero_()
    return module


def scale_module(module, scale):
    """
    Scale the parameters of a module and return it.
    """
    for p in module.parameters():
        p.detach().mul_(scale)
    return module


def mean_flat(tensor):
    """
    Take the mean over all non-batch dimensions.
    """
    return tensor.mean(dim=list(range(1, len(tensor.shape))))


def he_normalization(channels):
    """
    Make a standard normalization layer.
    :param channels: number of input channels.
    :return: an nn.Module for normalization.
    """
    return MPCGroupNorm32(32, channels)


def he_conv_nd(dims, *args, **kwargs):
    """
    Create a 1D, 2D, or 3D convolution module.
    """
    if dims == 1:
        return nn.Conv1d(*args, **kwargs)
    elif dims == 2:
        return nn.Conv2d(*args, **kwargs)
    # elif dims == 3:
    #     return nn.Conv3d(*args, **kwargs)
    raise ValueError(f"unsupported dimensions: {dims}")


def he_linear(*args, **kwargs):
    """
    Create a linear module.
    """
    return nn.Linear(*args, **kwargs)


def he_avg_pool_nd(dims, *args, **kwargs):
    """
    Create a 1D, 2D, or 3D average pooling module.
    """
    if dims == 1:
        return nn.AvgPool1d(*args, **kwargs)
    elif dims == 2:
        return nn.AvgPool2d(*args, **kwargs)
    raise ValueError(f"unsupported dimensions: {dims}")


class HybridConditioner(nn.Module):

    def __init__(self, c_concat_config, c_crossattn_config):
        super().__init__()
        self.concat_conditioner = instantiate_from_config(c_concat_config)
        self.crossattn_conditioner = instantiate_from_config(c_crossattn_config)

    def forward(self, c_concat, c_crossattn):
        c_concat = self.concat_conditioner(c_concat)
        c_crossattn = self.crossattn_conditioner(c_crossattn)
        return {'c_concat': [c_concat], 'c_crossattn': [c_crossattn]}


def noise_like(shape, device, repeat=False):
    repeat_noise = lambda: torch.randn((1, *shape[1:]), device=device).repeat(shape[0], *((1,) * (len(shape) - 1)))
    noise = lambda: torch.randn(shape, device=device)
    return repeat_noise() if repeat else noise()