orient-bias-filter: Orientation Bias Filter

# Evaluate whether asymptote squashing or double-exponential transform
# can deal with singularities at endpoints during integration


# apply linear change of variable y = c x + d
# map from [-1, 1] to [a, b]
linear_transform <- function(x, a=0, b=1) {
	c <- 0.5 * (b - a);
	d <- 0.5 * (a + b);
	c * x + d	
}

linear_atransform <- function(y, a=0, b=1) {
	c <- 0.5 * (b - a);
	d <- 0.5 * (a + b);
	(y - d) / c
}

# apply change of variable y = tanh( pi sinh(x) / 2 )
# map from (-Inf, Inf) to (-1, 1)
dexp_transform <- function(x) {
	tanh( pi * sinh(x) / 2)
}

dexp_atransform <- function(y) {
	asinh(atanh(y) * 2 / pi)
}

integrate_trapezoid <- function(f, a, b, n=64, eps=1e-6) {
	h <- (b - a - 2*eps) / n;
	x <- seq(a + eps, b - eps - h, length=n);

	h * sum(f(x))
}

integrate_midpoint <- function(f, a, b, n=64) {
	h <- (b - a) / n;
	x <- seq(a, b - h, length=n) + h/2;
	
	h * sum(f(x))
}

integrate_dexp <- function(f, a, b, n=64) {
	c <- 0.5 * (b - a);
	d <- 0.5 * (a + b);

	h <- (3 - -3) / n;
	k <- seq(-n/2, n/2, by=1);
	t <- seq(-3, 3, length=n + 1);
	x <- dexp_transform(t);

	ws <- pi/2 * cosh(k*h) / cosh(pi/2 * sinh(k*h))^2;
	f(c * dexp_transform(x) + d) * ws
	integral <- h * sum( f(c * x + d) * ws );

	c * integral
}

####

e <- 0.001;
alpha <- 0.1;
beta <- 10;

# null model: theta = 0

# all bases are reference, so that maximum likelihood of phi = 0

f <- function(phi) {
	#(1 - e)*(1 - phi) * phi^(alpha-1) * (1 - phi)^(beta-1)
	(1 - e)*(1 - phi) * dbeta(phi, alpha, beta)
}

f_reg <- function(phi) {
	#(1 - e)*(1 - phi) * phi^(alpha-1) * (1 - phi)^(beta-1)  -  (1 - e)*phi^(alpha -1)
	(1 - e)*(1 - phi) * dbeta(phi, alpha, beta)  -  (1 - e)*phi^(alpha -1)
}

phis <- seq(0, 0.01, by=1e-4);

plot(phis, f(phis), type="l")
f(1e-10)
f(1e-100)
f(0)

plot(phis, f_reg(phis), type="l")
f_reg(1e-10)
f_reg(1e-100)
f_reg(0)

phis <- seq(0, 1, by=1e-2);
xs <- seq(-1, 1, by=2e-2);

phis <- seq(0, 0.01, by=1e-4);
plot(phis, f(phis))
plot(phis, f(map_domain(dexp_transform(phis))))

phis <- seq(0, 1, by=1e-2);
plot(phis, f(phis))
plot(phis, f(map_domain(dexp_transform(phis))))
f(map_domain(dexp_transform(0)))


integrate(function(x) f(map_domain(dexp_transform(x))), -3, 3)

f3 <- function(x) {
	exp(-x / 5.0) * (2.0 + sin(2.0 * x))
}

integrate(f3, 0, 10)


plot(phis, f(dexp_transform(map_domain(xs))))


integrate(f, 0, 1)$value
integrate(f_reg, 0, 1)$value + (1 - e)/alpha



plot(curve(f3(x), xlim=c(0, 10)))

integrate(f3, 0, 10)
integrate_trapezoid(f3, 0, 10)
integrate_trapezoid(f3, 0, 10, eps=0)
integrate_midpoint(f3, 0, 10)
integrate_dexp(f3, 0, 10)

plot(curve(f(x), xlim=c(0, 1)))
plot(curve(f(linear_transform(dexp_transform(x))), xlim=c(-4, 4)))

plot(curve(f(x), xlim=c(0, 1e-6)))
plot(curve(f_reg(x), xlim=c(0, 1e-6)))

integrate(f, 0, 1)$value
integrate_trapezoid(f, 0, 1, eps=1e-10, n=1024)
integrate_midpoint(f, 0, 1, n=1024)
integrate_dexp(f, 0, 1)

# asymptote squashing helps the standard method but hurts the other
integrate(f_reg, 0, 1)$value + (1 - e)/alpha
integrate_trapezoid(f_reg, 0, 1) + (1 - e)/alpha
integrate_midpoint(f_reg, 0, 1, n=1024) + (1 - e)/alpha
integrate_dexp(f_reg, 0, 1) + (1 - e)/alpha

library(devtools)
load_all("..")

xs <- centered_grid(exp(-15))
midpoint_integration(xs, f_reg) + (1 - e)/alpha

midpoint_integration(xs, f)


# accuracy ranking:
# 1. standard integrate (Gauss quadrature) with asymptote squashing
# 2. double exponential transform integrate

# integrate with midpoint is woefully *inaccurate* with asymptote at x = 0
# because it is sensitive to where the grid is placed