# USL 3-parameter model 
# Created by NJG on Thu Sep 21 14:28:21 2017

# SGI Origin data 
df.gcap <- read.table(text="N	X_N
                    8 	130
                    12 	170
                    16 	190
                    20 	210
                    24 	220
                    28 	245
                    32 	255
                    48 	280
                    64 	310",
                    header=TRUE
)


# Sensitivity produces -beta unless bounds are applied to USL parameters
usl.fit <- nls(X_N ~ gamma * N / (1 + alpha * (N - 1) + beta * N * (N - 1)), 
           df.gcap, start=list(gamma=10, alpha=0.01, beta=0.0001), 
           algorithm="port",
           lower=c(1,0,0), upper=c(100,0.05,0.0002))

summary(usl.fit)

A <- coef(usl.fit)['alpha']
B <- coef(usl.fit)['beta']
G <- coef(usl.fit)['gamma']

Nmax <- sqrt((1 - A) / B)
Xmax <- G * Nmax /  (1 + A * (Nmax - 1) + B * Nmax * (Nmax - 1))
Nopt <- abs(1 / A)
Xroof <-G * Nopt 


plot(df.gcap$N, df.gcap$X_N, 
     xlim=c(0,200), ylim=c(0,450), 
     xlab="Users (N)", ylab="Thruput X(N)")
abline(v=Nopt, col="gray")
abline(h=Xroof, col="gray")
abline(a=0, b=Xroof/Nopt, col="gray")
abline(v=Nmax, col="red")
abline(h=Xmax, col="red")
curve(G * x / (1 + A * (x - 1) + B * x * (x - 1)), 
      from=0, to=200, add=TRUE, col="blue", lwd=2)
title("USL Scalability Analysis - GCAP 2017")
text(100, 200, "NJG on Thu Sep 21 2017", col="lightgray",cex=0.75)

legend("bottom", 
       legend=eval(parse(text=sprintf("expression(alpha == %.4f, beta == %.6f, gamma == %.4f, Nmax==%.2f, Xmax==%.2f,Xroof==%.2f)", A, B, G, Nmax, Xmax, Xroof))), ncol=2, inset=0.05, cex=0.75)