library(gamlss)
library(topsis)
library(forecast)
library(plotrix)
library(tsintermittent)
library(sfsmisc)

# For each scenario

o <- runif(1, 17, 140) # order cost in each scenario
s <- runif(1, 5, 50) # shortage cost in each scenario
h <- runif(1, 0, 0.68) # holding cost in each scenario
prop <- runif(1, 0.01, 0.99) # Proportion of zeros in the sample of each scenario
l <- 2 # lead time period

# ZINO Statistical Distribution

muNO <- as.numeric(
  gamlss(
    subset(k, k > 0) ~ 1,
    sigma.formula = ~1,
    nu.formula = ~1,
    family = NO(mu.link = "identity", sigma.link = "identity"),
    method = CG()
  )$mu.coefficients
)
sigmaNO <- as.numeric(
  gamlss(
    subset(k, k > 0) ~ 1,
    sigma.formula = ~1,
    nu.formula = ~1,
    family = NO(mu.link = "identity", sigma.link = "identity"),
    method = CG()
  )$sigma.coefficients
)
muNODLTsc <- as.numeric(
  gamlss(
    subset(DLT, DLT > 0) ~ 1,
    sigma.formula = ~1,
    nu.formula = ~1,
    family = NO(mu.link = "identity", sigma.link = "identity"),
    method = CG()
  )$mu.coefficients
)
sigmaNODLTsc <- as.numeric(
  gamlss(
    subset(DLT, DLT > 0) ~ 1,
    sigma.formula = ~1,
    nu.formula = ~1,
    family = NO(mu.link = "identity", sigma.link = "identity"),
    method = CG()
  )$sigma.coefficients
)

muNO # estimate mu of subset of DPUT>0
sigmaNO # estimate sigmaof subset of DPUT>0
muNOZ <- muNO * (1 - prop) # estimate mean DPUT ZINO
muNOZ
muNODLT <- muNODLTsc * (1 - nudlt) # estimate mean LTD with ZINOsum
muNODLT
sigmaNODLT <- sigmaNODLTsc * (1 - nudlt) # estimate standard deviation of LTD with ZINOsum
sigmaNODLT

# Optimization with iterative process
# Stop iterations when the variation of results is less than the threshold of 0.1
threshold <- 0.1
NNO <- sqrt(2 * o * muNOZ * 365 / h) # Initial quantity to order
rNORMAL <- qNO((1 - NNO * h / (muNOZ * 365 * s)) * (1 - nudlt), muNODLT, sigmaNODLT) # Initial reorder point
SrNORMAL <- sum(
  (
    subset(DLT, DLT > rNORMAL) - rNORMAL
  ) * dNO(
    subset(DLT, DLT > rNORMAL), muNODLT, sigmaNODLT
  ) / sum(
    dNO(subset(DLT, DLT > 0), muNODLT, sigmaNODLT)
  )
) # Initial expected shortage by cycle
repeat {
  prevNNO <- NNO
  prevRNORMAL <- rNORMAL
  prevSrNORMAL <- SrNORMAL
  NNO <- sqrt(2 * o * muNOZ * 365 / h) # quantity to order
  rNORMAL <- qNO((1 - NNO * h / (muNOZ * 365 * s)) * (1 - nudlt), muNODLT, sigmaNODLT) # reorder point
  SrNORMAL <- sum(
    (
      subset(DLT, DLT > rNORMAL) - rNORMAL
    ) * dNO(
      subset(DLT, DLT > rNORMAL), muNODLT, sigmaNODLT
    ) / sum(
      dNO(subset(DLT, DLT > 0), muNODLT, sigmaNODLT)
    )
  ) # expected shortage by cycle
  if (abs(NNO - prevNNO) < threshold) {
    break
  }
}
NNO
rNORMAL
SrNORMAL

# Total cost with ZINO statistical distribution
TCNO <- (NNO / 2 + rNORMAL - muNOZ * l) * h + muNOZ * 365 * o / NNO + SrNORMAL * muNOZ * 365 * s / NNO

# Fill-rate with ZINO statistical distribution
FRNO <- (
    sum(
      dNO(subset(k, k < NNO), muNO, sigmaNO) * subset(k, k < NNO)
    ) + sum(
      dNO(subset(k, k > NNO), muNO, sigmaNO) * NNO
    )
  ) / sum(
    dNO(k, muNO, sigmaNO) * k
  )

# KL divergence of ZINOsum statistical distribution in LTD

KLNOsumn <-
  sum(
    (
      (
        1 - nudlt
      ) * density(
        subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0))
      )$y / sum(
        density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y
      )
    ) * (
      ifelse(
        log(
          (
            (1 - nudlt) * density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y / sum(density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y)
          ) / (
            (1 - nudlt) * dNO(
              qNO(runif(length(subset(DLT, DLT > 0)), 0, 1), muNODLT, sigmaNODLT), muNODLT, sigmaNODLT
            ) / sum(
              dNO(qNO(runif(length(subset(DLT, DLT > 0)), 0, 1), muNODLT, sigmaNODLT), muNODLT, sigmaNODLT)
            )
          )
        ) == -Inf,
        0,
        log(
          (
            (1 - nudlt) * density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y / sum(density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y)
          ) / (
            (1 - nudlt) * dNO(
              qNO(runif(length(subset(DLT, DLT > 0)), 0, 1), muNODLT, sigmaNODLT), muNODLT, sigmaNODLT
            ) / sum(
              dNO(qNO(runif(length(subset(DLT, DLT > 0)), 0, 1), muNODLT, sigmaNODLT), muNODLT, sigmaNODLT)
            )
          )
        )
      )
    )
  )


####### SES Approach

muzeroSES <- mean(ets(c(remDLT[1, ], remDLT[2, ]), model = "AZZ", opt.crit = "sigma")$fitted)
sigmaSES <- sqrt(
  sum(
    (
      data.frame(
        c(remDLT[1, ], remDLT[2, ])
      )[1:200 - 1, ] - data.frame(
        ets(c(remDLT[1, ], remDLT[2, ]), model = "AZZ", opt.crit = "sigma")$fitted
      )[2:200, ]
    )^2
  ) / 199
)
muzeroSES # estimate mean DPUT SES
sigmaSES # estimate standard deviation DPUT SES
# Initial
NSES <- sqrt(2 * o * muzeroSES * 365 / h)
rSES <- qNO((1 - NSES * h / (muzeroSES * 365 * s)) * (1 - nudlt), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt))
SrSES <- sum(
  (
    subset(DLT, DLT > rSES) - rSES
  ) * dNO(
    subset(DLT, DLT > rSES), muzeroSES * l, sigmaSES * sqrt(l)
  ) / sum(
    dNO(subset(DLT, DLT > 0), muzeroSES * l, sigmaSES * sqrt(l))
  )
)
# Iterations
threshold <- 0.1
repeat {
  prevNSES <- NSES
  prevRSES <- rSES
  prevSrSES <- SrSES
  NSES <- sqrt(2 * muzeroSES * 365 * (o + s * SrSES) / h)
  rSES <- qNO((1 - NSES * h / (muzeroSES * 365 * s)) * (1 - nudlt), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt))
  SrSES <- sum(
    (
      subset(DLT, DLT > rSES) - rSES
    ) * dNO(
      subset(DLT, DLT > rSES), muzeroSES * l, sigmaSES * sqrt(l)
    ) / sum(
      dNO(subset(DLT, DLT > 0), muzeroSES * l, sigmaSES * sqrt(l))
    )
  )
  if (abs(NSES - prevNSES) < threshold) {
    break
  }
}

NSES
rSES
SrSES

# Stop iterations when the variation of results is less than the threshold of 0.1

# TC SES method

TCSES <- (NSES / 2 + rSES - muzeroSES * l) * h + muzeroSES * 365 * o / NSES + SrSES * muzeroSES * 365 * s / NSES

# Fill-rate SES method

FRSES <- matrix(
  (
    sum(
      dNO(subset(k, k < NSES), muzeroSES, sigmaSES) * subset(k, k < NSES)
    ) + sum(
      dNO(subset(k, k > NSES), muzeroSES, sigmaSES) * NSES
    )
  ) / sum(dNO(k, muzeroSES, sigmaSES) * k), 1, 1
)


# KL divergence of SESsum statistical distribution in LTD

KLSESsumn <-
  sum(
    (
      (
        1 - nudlt
      ) * density(
        subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0))
      )$y / sum(
        density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y
      )
    ) * (
      ifelse(
        log(
          (
            (1 - nudlt) * density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y / sum(density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y)
          ) / (
            (1 - nudlt) * dNO(
              qNO(
                runif(length(subset(DLT, DLT > 0)), 0, 1), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt)
              ), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt)
            ) / sum(
              dNO(
                qNO(
                  runif(length(subset(DLT, DLT > 0)), 0, 1), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt)
                ), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt)
              )
            )
          )
        ) == -Inf,
        0,
        log(
          (
            (1 - nudlt) * density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y / sum(density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y)
          ) / (
            (1 - nudlt) * dNO(
              qNO(runif(length(subset(DLT, DLT > 0)), 0, 1), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt)), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt)
            ) / sum(
              dNO(qNO(runif(length(subset(DLT, DLT > 0)), 0, 1), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt)), muzeroSES * l, sigmaSES * sqrt(l) * (1 - nudlt))
            )
          )
        )
      )
    )
  )


####### CROSTON's Approach

muCROST <- mean(data.frame(crost(c(remDLT[1, ], remDLT[2, ]), type = "sba")$frc.in)[125:200, ])
sigmaCROST <- sqrt(
  sum(
    (
      data.frame(
        c(remDLT[1, ], remDLT[2, ])
      )[124:200 - 1, ] - data.frame(
        crost(c(remDLT[1, ], remDLT[2, ]), type = "sba")$frc.in
      )[125:200, ]
    )^2
  ) / 75
)
muCROST # estimate mean DPUT CROST
sigmaCROST # estimate mean DPUT CROST

# Initial
NCROST <- sqrt(2 * o * muCROST * 365 / h)
rCROST <- qNO((1 - NCROST * h / (muCROST * 365 * s)) * (1 - nudlt), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt))
SrCROST <- sum(
  (
    subset(DLT, DLT > rCROST) - rCROST
  ) * dNO(
    subset(DLT, DLT > rCROST), muCROST * l, sigmaCROST * sqrt(l)
  ) / sum(
    dNO(subset(DLT, DLT > 0), muCROST * l, sigmaCROST * sqrt(l))
  )
)
# Iterations
threshold <- 0.1
repeat {
  prevNCROST <- NCROST
  prevRCROST <- rCROST
  prevSrCROST <- SrCROST
  NCROST <- sqrt(2 * muCROST * 365 * (o + s * SrCROST) / h)
  rCROST <- qNO((1 - NCROST * h / (muCROST * 365 * s)) * (1 - nudlt), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt))
  SrCROST <- sum(
    (
      subset(DLT, DLT > rCROST) - rCROST
    ) * dNO(
      subset(DLT, DLT > rCROST), muCROST * l, sigmaCROST * sqrt(l)
    ) / sum(
      dNO(subset(DLT, DLT > 0), muCROST * l, sigmaCROST * sqrt(l))
    )
  )
  if (abs(NCROST - prevNCROST) < threshold) {
    break
  }
}

NCROST
rCROST
SrCROST

# TC Chroston's method
TCCROST <- (NCROST / 2 + rCROST - muCROST * l) * h + muCROST * 365 * o / NCROST + SrCROST * muCROST * 365 * s / NCROST

# Fill-rate Chroston's method
FRCROST <- (
  sum(
    dNO(subset(k, k < NCROST), muCROST, sigmaCROST) * subset(k, k < NCROST)
  ) + sum(
    dNO(subset(k, k > NCROST), muCROST, sigmaCROST) * NCROST
  )
) / sum(
  dNO(k, muCROST, sigmaCROST) * k
)

# KL divergence of CROST sum statistical distribution in LTD

KLCROSTsumn <-
  sum(
    (
      (
        1 - nudlt
      ) * density(
        subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0))
      )$y / sum(
        density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y
      )
    ) * (
      ifelse(
        log(
          (
            (1 - nudlt) * density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y / sum(density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y)
          ) / (
            (1 - nudlt) * dNO(
              qNO(
                runif(length(subset(DLT, DLT > 0)), 0, 1), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt)
              ), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt)
            ) / sum(
              dNO(
                qNO(
                  runif(length(subset(DLT, DLT > 0)), 0, 1), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt)
                ), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt)
              )
            )
          )
        ) == -Inf,
        0,
        log(
          (
            (1 - nudlt) * density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y / sum(density(subset(DLT, DLT > 0), n = length(subset(DLT, DLT > 0)))$y)
          ) / (
            (1 - nudlt) * dNO(
              qNO(
                runif(length(subset(DLT, DLT > 0)), 0, 1), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt)
              ), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt)
            ) / sum(
              dNO(
                qNO(
                  runif(length(subset(DLT, DLT > 0)), 0, 1), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt)
                ), muCROST * l, sigmaCROST * sqrt(l) * (1 - nudlt)
              )
            )
          )
        )
      )
    ))


# TOPSIS
d <- array(c(TCNO, TCSES, TCCROST, FRNO, FRSES, FRCROST, SrNORMAL, SrSES, SrCROST, KLNOsumn, KLSESsumn, KLCROSTsumn), c(1, 3, 4))


top <- topsis(d[1, , ], c(1, 1, 1, 1), c("-", "+", "-", "-"))[, 3]

# EXPORT RESULTS

x <- data.frame(
  Proportion = prop, Order.C = o, Holding.C = h, Shotgage.C = s,
  TCzinormal = TCNO, TCses = TCSES, TCcroston = TCCROST,
  FRzinormal = FRNO, FRses = FRSES, FRcroston = FRCROST,
  Srzinormal = SrNORMAL, Srses = SrSES, Srcroston = SrCROST,
  Qzinormal = NNO, Qses = NSES, Qcroston = NCROST,
  rzinormal = rNORMAL, rses = rSES, rcroston = rCROST,
  KLzinormalDLT = KLNOsumn, KLsesDLT = KLSESsumn, KLcrostonDLT = KLCROSTsumn,
  TOPSIS = t(top)
)

write.table(x, file = "result.csv", sep = ",", col.names = NA, qmethod = "double")