## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, fig.width = 7, fig.height = 5)
library(BsplineQuantReg)

## ----polymul------------------------------------------------------------------
# (1 + x) * (1 + x) = 1 + 2x + x^2
p1 <- c(1, 1)  # 1 + x
p2 <- c(1, 1)  # 1 + x
product <- polymul(p1, p2)
print(product)  # c(1, 2, 1) → 1 + 2x + x^2

# (x^2 + 2x + 1) * (x - 1) = x^3 + x^2 - x - 1
p3 <- c(1, 2, 1)  # x^2 + 2x + 1
p4 <- c(1, -1)    # x - 1
polymul(p3, p4)

## ----polyadd------------------------------------------------------------------
# (1 + x) + (1 - x) = 2
polyadd(c(1, 1), c(1, -1))

# (x^2 + 1) + (x + 1) = x^2 + x + 2
polyadd(c(1, 0, 1), c(1, 1))

## ----poly_eval----------------------------------------------------------------
# P(x) = 1 + x + x^2
p <- c(1, 1, 1)  # x^2 + x + 1
poly_eval(p, c(0, 1, 2))  # returns c(1, 3, 7)

# Evaluate at many points
x <- seq(-2, 2, length.out = 10)
y <- poly_eval(c(1, 0, -1), x)  # 1 - x^2
plot(x, y, type = "l", main = "1 - x^2")

## ----polyderiv----------------------------------------------------------------
# P(x) = x^2 → P'(x) = 2x
p <- c(1, 0, 0)  # x^2
polyderiv(p, der = 1)  # c(2, 0) → 2x

# Second derivative: P''(x) = 2
polyderiv(p, der = 2)  # c(2)

# Higher order derivatives
p <- c(1, 2, 3, 4)  # 4x^3 + 3x^2 + 2x + 1
polyderiv(p, der = 2)  # 24x + 6

## ----reduce_pol---------------------------------------------------------------
# Remove leading zeros
reduce_pol(c(0, 0, 1, 2, 1))  # c(1, 2, 1) → x^2 + 2x + 1
reduce_pol(c(0, 0, 0, 5))     # c(5) → constant polynomial

## ----change_polynomial_base_taylor--------------------------------------------
# Convert (x-1)^2 to expansion around 0
# (x-1)^2 = x^2 - 2x + 1
coeffs_a <- c(1, 0, 0)  # (x-1)^2 in basis centered at a=1
change_polynomial_base_taylor(coeffs_a, a = 1, b = 0)
# Returns c(1, -2, 1) → x^2 - 2x + 1

## ----makpp--------------------------------------------------------------------
# Create a piecewise polynomial with two intervals
# Interval 1: x^2 on [0, 1]
# Interval 2: 2x - 1 on [1, 2]
coeff <- matrix(c(
  1, 0, 0,   # x^2
  0, 2, -1   # 2x - 1
), nrow = 2, byrow = TRUE)

knots <- c(0, 1, 2)

pp <- makpp(coeff, knots)
print(pp)

# Evaluate the PP form
x <- seq(0, 2, length.out = 100)
y <- evalpp(pp, x)

plot(x, y, type = "l", main = "Piecewise Polynomial")
abline(v = knots, col = "red", lty = 2)

## ----non_callable_pp----------------------------------------------------------
# A non-callable PP object is a list with components:
# - coeff: matrix of polynomial coefficients
# - knot: knot positions
# - degree: polynomial degree

pp <- makpp(coeff, knots, callable = FALSE)
class(pp)  # "non_callable_pp"

# Access components
pp$coeff
pp$knot
pp$degree

## ----callable_pp--------------------------------------------------------------
# Create a callable PP
pp_call <- makpp(coeff, knots, callable = TRUE)
class(pp_call)  # "callable_pp" "function"

# Evaluate directly
x <- seq(0, 2, length.out = 10)
y <- pp_call(x)

# Print shows information
print(pp_call)

# Can be used in plots
plot(pp_call, xlim = c(0, 2))

## -----------------------------------------------------------------------------
# P(x) = 3x^3 - 2x^2 + x - 5
p <- c(3, -2, 1, -5)
show_poly(p)
# "3*x^3 - 2*x^2 + 1*x - 5"

## ----local_basis--------------------------------------------------------------
# P(x) = 3(x-2)^3 - 2(x-2)^2 + (x-2) - 5
p <- c(3, -2, 1, -5)
show_poly(p, b = 2)
# "3*(x-2)^3 - 2*(x-2)^2 + 1*(x-2) - 5"

## ----convert_bases------------------------------------------------------------
# Coefficients are in basis (x-2)^k, display in canonical basis
p <- c(3, -2, 1, -5)  # In basis (x-2)^k
show_poly(p, a = 2, b = 0)
# This expands to: 3x^3 - 18x^2 + 37x - 27

## ----centers------------------------------------------------------------------
# Center at positive value
p <- c(1, 0, 2)  # (x-3)^2 + 2
show_poly(p, b = 3)
# "(x-3)^2 + 2"

# Center at negative value
p <- c(1, 0, 2)  # (x+2)^2 + 2
show_poly(p, b = -2)
# "(x+2)^2 + 2"

## ----digits-------------------------------------------------------------------
p <- c(1/3, -2/7, 1/5)
show_poly(p, digits = 2)
# "0.33*x^2 - 0.29*x + 0.2"

show_poly(p, digits = 4)
# "0.3333*x^2 - 0.2857*x + 0.2"

## ----negative-----------------------------------------------------------------
# All coefficients negative
p <- c(-3, -2, -1)  # -3x^2 - 2x - 1
show_poly(p)
# "-3*x^2 - 2*x - 1"

# Mixed coefficients
p <- c(3, -2, 1, -5)  # 3x^3 - 2x^2 + x - 5
show_poly(p)
# "3*x^3 - 2*x^2 + 1*x - 5"

## ----bspline_example----------------------------------------------------------
# Create a B-spline and display its polynomial form
sn <- c(0, 0, 0, 0, 0.3, 0.6, 1, 1, 1, 1)
basis <- Bspline_base(sn, degree = 3)
basis$coeff <- c(1, -0.5, 2, 0.5, -1, 0.75)

# Convert to PP
pp <- Bsplinetopp(basis, callable = FALSE)

# Display each piece
cat("Piece 1: ", show_poly(pp$coeff[1, ], b = 0), "\n")
cat("Piece 2: ", show_poly(pp$coeff[2, ], b = 0), "\n")
cat("Piece 3: ", show_poly(pp$coeff[3, ], b = 0), "\n")

## ----show_pp_integration------------------------------------------------------
# Create a PP with two intervals
knot <- c(0, 0.5, 1)
coeff <- matrix(c(1, 2, 0.5, 0, 1, -1), nrow = 2, ncol = 3, byrow = TRUE)
pp <- makpp(coeff, knot)

# Display in local basis
show_pp(pp, local = TRUE, digits = 3)

## ----show_pp------------------------------------------------------------------
# Create a PP with two intervals
knot <- c(0, 0.5, 1)
coeff <- matrix(c(1, 2, 0.5, 0, 1, -1), nrow = 2, ncol = 3, byrow = TRUE)
pp <- makpp(coeff, knot)

# Display in canonical basis (1, x, x², ...)
show_pp(pp, local = FALSE)
# [0.000, 0.500] 0.5x^2 + 2x + 1
# [0.500, 1.000] -1x^2 + 1x + 0

# Display in local basis ((x-a)^i)
show_pp(pp, local = TRUE)
# [0.000, 0.500] 0.5(x-0)^2 + 2(x-0) + 1
# [0.500, 1.000] -1(x-0.5)^2 + 1(x-0.5) + 0

## -----------------------------------------------------------------------------

# Create a B-spline and display its polynomial form
sn <- c(0, 0, 0, 0, 0.3, 0.6, 1, 1, 1, 1)
basis <- Bspline_base(sn, degree = 3)
basis$coeff <- c(1, -0.5, 2, 0.5, -1, 0.75)

show_pp(basis, local = TRUE, verbose = TRUE)
# PP Information:
# Degree: 3
# Knots: 0, 0.3, 0.6, 1
# Number of intervals: 3
# 
# [0.000, 0.300] 1(x-0)^3 + 2.5(x-0)^2 + 1.5(x-0) + 1
# [0.300, 0.600] -4.63(x-0.3)^3 + ... 
# [0.600, 1.000] ...

## ----get_parameters-----------------------------------------------------------
# For callable PP objects, use get_parameters()
params <- get_parameters(pp_call)
print(params$degree)
print(params$knot)
print(params$coeff)

# For non-callable PP, access directly
pp=makpp(pp,callable=FALSE,verbose=TRUE)
pp$degree
pp$knot
pp$coeff

## ----extract_interval---------------------------------------------------------
# Get polynomial for interval i
pp=makpp(pp,callable=FALSE)
interval_poly <- pp$coeff[1, ]  # First interval
print(interval_poly)

# Evaluate on a subinterval
x_sub <- seq(0, 1, length.out = 20)
y_sub <- poly_eval(interval_poly, x_sub)

## ----combine_pp---------------------------------------------------------------
# Create two PP objects and combine them
# This is useful for constructing complex piecewise functions
coeff1 <- matrix(c(1, 0), nrow = 1)
coeff2 <- matrix(c(0, 1), nrow = 1)

pp1 <- makpp(coeff1, c(0, 1))
pp2 <- makpp(coeff2, c(1, 2))

# Combine coefficients and knots
combined_coeff <- rbind(coeff1, coeff2)
combined_knots <- c(0, 1, 2)
pp_combined <- makpp(combined_coeff, combined_knots)

# Evaluate combined function
x <- seq(0, 2, length.out = 100)
y_combined <- evalpp(pp_combined, x)
plot(x, y_combined, type = "l")
abline(v = c(0, 1, 2), col = "red", lty = 2)

## ----visualize_pp-------------------------------------------------------------
# Create a PP with multiple pieces
knots <- seq(0, 1, length.out = 5)
degree <- 2
n_pieces <- length(knots) - 1

# Generate random coefficients
coeff <- matrix(rnorm(n_pieces * (degree + 1)), nrow = n_pieces)

pp <- makpp(coeff, knots, callable = TRUE)

# Evaluate and plot
x <- seq(0, 1, length.out = 200)
y <- pp(x)

plot(x, y, type = "l", lwd = 2,
     main = "Random Piecewise Quadratic")
abline(v = knots, col = "red", lty = 2)
grid()

## ----pp_vs_bspline------------------------------------------------------------
# Create a B-spline and convert to PP form
sn <- c(0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1)
basis <- Bspline_base(sn, degree = 3)
basis$coeff <- runif(basis$n_splines)

# Convert to PP
pp_from_bspline <- Bsplinetopp(basis)

# Both represent the same function
x <- seq(0, 1, length.out = 100)
y_bspline <- spline_eval(basis, x)
y_pp <- evalpp(pp_from_bspline, x)

# They should match
max(abs(y_bspline - y_pp))