# Demo: Comprehensive Test of Quantile Regression Constraints
# Author: Alexandre Abbes
#
# This demo runs a complete test of all constraint types:
# - Partial monotonicity (increasing on first intervals only)
# - Full monotonicity (decreasing everywhere)
# - Convexity constraint
# - Unconstrained fit
#
# The results are displayed in a 2x2 panel for comparison.


oldpar <- par(mfrow = c(2,2))


library(BsplineQuantReg)

cat("==============================================\n")
cat("Demo: Comprehensive Quantile Regression Tests\n")
cat("==============================================\n\n")

if (!exists("degree")){degree=3}


# Generate synthetic data
#set.seed(42)
n_points <- 30
xtab <- seq(0, 1, length.out = n_points + 1)  # 51 points from 0 to 1

# Data: increasing trend with oscillation + noise
# y = 2x + 0.5*sin(6 pi x) + noise
ytab <- 2 * xtab + 0.5 * sin(6 * pi * xtab) + 0.2 * rnorm(n_points + 1)

x_eval <- seq(0, 1, length.out = 200)
y_true <-  2 * x_eval + 0.5 * sin(6 * pi * x_eval)


# Number of intervals
kn <- 12
cat(sprintf("Number of observations: %d\n", length(xtab)))
cat(sprintf("Number of intervals (kn): %d\n", kn))
cat(sprintf("Number of basis functions: %d\n\n", kn + 3))

# Define knots (quantile-based)
knots <- quantile(xtab, probs = seq(0, 1, length.out = kn + 1))

# Define partial monotonicity: increasing only on first 7 intervals
n_constrained <- 7
monot_partial <- c(rep(1, n_constrained), rep(0, kn - n_constrained))

cat("Constraint configurations:\n")
cat("-------------------------\n")
cat("1. Unconstrained:        no constraints\n")
cat("2. Partial increasing:   increasing on intervals 1-7 only\n")
cat("3. Full decreasing:      decreasing everywhere\n")
cat("4. Convexity:            second derivative >= 0 everywhere\n\n")

cat("Fitting models...\n")

# Fit 1: Unconstrained (tau = 0.9, 0.5, 0.1)
cat("  - Unconstrained (tau = 0.9, 0.5, 0.1)")
res_uncon1 <- quantile_spline(xtab, ytab, knots, tau = 0.9,
                                    monot = 0, convcons = 0,degree=degree )
res_uncon2 <- quantile_spline(xtab, ytab, knots, tau = 0.1,
                                    monot = 0, convcons = 0,degree=degree)
res_uncon3 <- quantile_spline(xtab, ytab, knots, tau = 0.5,
                                    monot = 0, convcons = 0,degree=degree)

cat(" done\n")

# Fit 2: Partial increasing (tau = 0.5)
cat("  - Partial increasing (tau = 0.5, 0.1, intervals 1-7)...")
res_croissant1 <- quantile_spline(xtab, ytab, knots, tau = 0.5,
                                        monot = monot_partial, convcons = 0,
                                  degree=degree)
res_croissant2 <- quantile_spline(xtab, ytab, knots, tau = 0.1,
                                         monot = monot_partial,degree=degree, convcons = 0)
cat(" done\n")

# Fit 3: Full decreasing (tau = 0.5)
cat("  - Full decreasing (tau = 0.5)...")
res_decroissant <- quantile_spline(xtab, ytab, knots, tau = 0.5,
                                          monot = -1, convcons = 0,
                                   degree=degree,
                                          solver = "HIGHS")
cat(" done\n")

# Fit 4: Convexity constraint (tau = 0.5)
cat("  - Convexity constraint (tau = 0.5)...")
res_convexe <- quantile_spline(xtab, ytab, knots, tau = 0.5,
                                      monot = 0, convcons = 1,
                               degree=degree,
                                      solver = "HIGHS")
cat(" done\n\n")

# Evaluate all fits on a fine grid

y_uncon1 <- spline_eval(res_uncon1, x_eval)
y_uncon2 <- spline_eval(res_uncon2, x_eval)
y_uncon3 <- spline_eval(res_uncon3, x_eval)
y_croiss1 <- spline_eval(res_croissant1, x_eval)
y_croiss2 <- spline_eval(res_croissant2, x_eval)
y_decroiss <- spline_eval(res_decroissant, x_eval)
y_convexe <- spline_eval(res_convexe, x_eval)

# Compute derivatives for verification
compute_derivative <- function(y, x) {
  diff(y) / diff(x)
}

deriv_croiss <- compute_derivative(y_croiss1, x_eval)
deriv_decroiss <- compute_derivative(y_decroiss, x_eval)

cat("Verification of constraints:\n")
cat("----------------------------\n")
cat(sprintf("  Partial increasing: min derivative = %.4f (should be >= 0 on constrained region)\n",
            min(deriv_croiss[1:140])))  # first 140 points ~ first 7 intervals
cat(sprintf("  Full decreasing:    max derivative = %.4f (should be <= 0)\n",
            max(deriv_decroiss)))
cat("\n")

# Create plots

# Ajouter un espace en haut de la page pour le titre
par(mfrow = c(2, 2), mar = c(4, 4, 4, 2), oma = c(0, 0, 2, 0))

# Plot 1: Partial increasing constraint
plot(xtab, ytab, pch = 16,  col = "black",
     xlab = "x", ylab = "y",
     main = "Partial Increasing Constraint\n(first 7 intervals only)\n tau=0.5, 0.1")
lines(x_eval, y_croiss1, col = "blue", lwd = 2)
lines(x_eval, y_croiss2, col = "blue", lwd = 2)
abline(v = knots, col = "blue", lty = 2, lwd = 0.2)
legend("topleft",c("true","fitted"),col = c("black", "blue"),lwd = 1.5)
lines(x_eval, y_true, col = "black", lwd = 1, )
# Highlight the constrained region
abline(v = knots[n_constrained + 1], col = "red", lty = 2, lwd = 2)
text(knots[n_constrained + 1] + 0.1, max(ytab) - 1.7,
     "Constrained \n region", col = "red", cex = 0.7, srt = 90)
grid()

# Plot 2: Full decreasing constraint
plot(xtab, ytab, pch = 16,  col = "black",
     xlab = "x", ylab = "y",
     main = "Full Decreasing Constraint\n(everywhere)")
lines(x_eval, y_decroiss, col = "darkgreen", lwd = 2)
abline(v = knots, col = "blue", lty = 2, lwd = 0.5)
legend("topleft",c("true","fitted"),col = c("black", "darkgreen"),lwd = 2)
lines(x_eval, y_true, col = "black", lwd = 1, )
grid()

# Plot 3: Convexity constraint
plot(xtab, ytab, pch = 16,  col = "black",
     xlab = "x", ylab = "y",
     main = "Convexity Constraint\n(second derivative >= 0)")
lines(x_eval, y_convexe, col = "purple", lwd = 2)
legend("topleft",c("true","fitted"),lwd = 2,col=c("black","purple"))
lines(x_eval, y_true, col = "black", lwd = 1, )
abline(v = knots, col = "blue", lty = 2, lwd = 0.5)
grid()

# Plot 4: Unconstrained (tau = 0.9)
plot(xtab, ytab, pch = 16,  col = "black",
     xlab = "x", ylab = "y",
     main = "Unconstrained\n(tau = 0.9, , 0.5, 0.1)")
legend("topleft",c("true","fitted"),col=c("black","red"),lwd = 2)
lines(x_eval, y_true, col = "black", lwd = 1, )
lines(x_eval, y_uncon1, col = "red", lwd = 2)
lines(x_eval, y_uncon2, col = "red", lwd = 2)
lines(x_eval, y_uncon3, col = "red", lwd = 2)
abline(v = knots, col = "blue", lty = 2, lwd = 0.5)
grid()

#unlink("/home/abbes/R/x86_64-pc-linux-gnu-library/4.6/BsplineQuantReg", recursive = TRUE)
# Ajouter le titre global
mtext(paste("Comprehensive Quantile Regression Tests - Degree", degree),
      side = 3, line = 0.5, outer = TRUE, cex = 1.2, font = 2)


cat("==============================================\n")
cat("Demo completed. Check the plots for comparison.\n")
cat("==============================================\n")

par(oldpar)
