# Some R codes related to the paper 
#
# M. Chae, R. Martin, and S. G. Walker (2019). On an algorithm for solving Fredholm equations 
# of the first kind. Statistics and Computing, volume 29, pages 645-654; arXiv:1709.02695


myint <- function(fx, x, m=2) {

  h <- x[2] - x[1]
  n <- length(x)
  if(!(m %in% c(0, 1, 2))) stop("Only Riemann, Trapezoid, or Simpson is implemented!")
  if(m == 0) intf <- h * sum(fx[-n])
  else if(m == 1) intf <- (h / 2) * sum(fx[-1] + fx[-n])
  else {

    if(n %% 2 != 0) stop("Simpson's rule requires grid of even length!")
    intf <- (h / 3) * (sum(fx) + sum(fx[2:(n-1)]) + 2 * sum(fx[seq(2, n-1, by=2)]))

  }
  return(intf)

}

newpr <- function(x, fx, u, K, p0, M) {

  if(length(u) %% 2 != 0) stop("u-grid length should be of even length!")
  if(length(x) %% 2 != 0) stop("x-grid length should be of even length!")
  tK <- t(K)
  p <- p0
  Kp <- t(tK * p)
  fp <- apply(Kp, 1, myint, x=u)
  D <- myint(fx * log(fx / fp), x=x)
  #D.diff <- numeric(M)
  for(m in 1:M) {

    p <- apply(Kp * fx / fp, 2, myint, x=x)
    Kp <- t(tK * p)
    fp <- apply(Kp, 1, myint, x=u)
    Dm <- myint(fx * log(fx / fp), x=x)
    #D.diff[m] <- D[length(D)] - Dm
    D <- c(D, Dm)

  }
  return(list(u=u, pu=p, x=x, fpx=fp, D=D))

}


## Galaxy data illustration

M <- 25
library(MASS)
data(galaxies)
Y <- galaxies / 1000
d <- density(Y, n=2**9)
y <- x <- d$x
fy <- d$y
p0 <- dt(x, df=1)
K <- dnorm(outer(y, x, "-"))
o <- newpr(y, fy, x, K, p0, M)
plot(x, o$pu, type="l", xlab=expression(theta), ylab=expression(p(theta)))
ylim <- c(0, max(c(o$fpx, fy)))
hist(Y, freq=FALSE, breaks=10, xlab="x", ylab="f(x)", col="gray", border="white", main="", ylim=ylim)
lines(y, o$fpx, lty=2)
lines(y, fy)
plot(o$D, type="l", xlab="m", ylab=expression(D(f,f[m])))