diff --git a/planetiler-benchmarks/src/main/java/com/onthegomap/planetiler/benchmarks/BenchmarkInterpolator.java b/planetiler-benchmarks/src/main/java/com/onthegomap/planetiler/benchmarks/BenchmarkInterpolator.java new file mode 100644 index 0000000000..94f113cac1 --- /dev/null +++ b/planetiler-benchmarks/src/main/java/com/onthegomap/planetiler/benchmarks/BenchmarkInterpolator.java @@ -0,0 +1,55 @@ +package com.onthegomap.planetiler.benchmarks; + +import com.onthegomap.planetiler.util.Scales; +import java.util.function.DoubleUnaryOperator; +import java.util.function.Supplier; + +public class BenchmarkInterpolator { + private static double sum = 0; + + public static void main(String[] args) { + long times = 10_000_000; + benchmarkInterpolator("linear", times, Scales::linear); + benchmarkInterpolator("sqrt", times, Scales::sqrt); + benchmarkInterpolator("pow2", times, () -> Scales.power(2)); + benchmarkInterpolator("pow10", times, () -> Scales.power(1)); + benchmarkInterpolator("log", times, () -> Scales.log(10)); + benchmarkInterpolator("log2", times, () -> Scales.log(2)); + System.err.println(sum); + } + + private static void benchmarkInterpolator(String name, long times, Supplier get) { + benchmarkAndInverted(name + "_2", 1, 2, times, () -> get.get().put(1, 1d).put(2, 2d)); + benchmarkAndInverted(name + "_3", 1, 2, times, () -> get.get().put(1, 1d).put(1.5, 2d).put(2, 3d)); + } + + private static void benchmarkAndInverted(String name, double start, double end, long steps, + Supplier build) { + benchmark(name + "_f", start, end, steps, build::get); + benchmark(name + "_i", start, end, steps, () -> build.get().invert()); + } + + private static void benchmark(String name, double start, double end, long steps, + Supplier build) { + double delta = (end - start) / steps; + double x = start; + double result = 0; + for (long i = 0; i < steps / 10_000; i++) { + result += build.get().applyAsDouble(x += delta); + } + x = start; + long a = System.currentTimeMillis(); + for (long i = 0; i < steps; i++) { + result += build.get().applyAsDouble(x += delta); + } + x = start; + var preBuilt = build.get(); + long b = System.currentTimeMillis(); + for (long i = 0; i < steps; i++) { + result += preBuilt.applyAsDouble(x += delta); + } + long c = System.currentTimeMillis(); + sum += result; + System.err.println(name + "\t" + (b - a) + "\t" + (c - b)); + } +} diff --git a/planetiler-core/src/main/java/com/onthegomap/planetiler/util/IInterpolator.java b/planetiler-core/src/main/java/com/onthegomap/planetiler/util/IInterpolator.java new file mode 100644 index 0000000000..870d7c9c7c --- /dev/null +++ b/planetiler-core/src/main/java/com/onthegomap/planetiler/util/IInterpolator.java @@ -0,0 +1,24 @@ +package com.onthegomap.planetiler.util; + +import java.util.function.BiFunction; +import java.util.function.DoubleFunction; +import java.util.function.DoubleUnaryOperator; + +public interface IInterpolator, V> extends DoubleFunction { + T self(); + + T put(double x, V y); + + + interface ValueInterpolator extends BiFunction> {} + interface Continuous & Continuous> + extends IInterpolator, DoubleUnaryOperator { + default DoubleUnaryOperator invert() { + return Interpolator.invertIt(this.self()); + } + + default T put(double x, double y) { + return put(x, Double.valueOf(y)); + } + } +} diff --git a/planetiler-core/src/main/java/com/onthegomap/planetiler/util/Interpolator.java b/planetiler-core/src/main/java/com/onthegomap/planetiler/util/Interpolator.java new file mode 100644 index 0000000000..4714638e8c --- /dev/null +++ b/planetiler-core/src/main/java/com/onthegomap/planetiler/util/Interpolator.java @@ -0,0 +1,272 @@ +package com.onthegomap.planetiler.util; + +import com.carrotsearch.hppc.DoubleArrayList; +import java.util.ArrayList; +import java.util.List; +import java.util.function.DoubleFunction; +import java.util.function.DoubleUnaryOperator; +import org.apache.commons.lang3.ArrayUtils; + +public class Interpolator, V> implements IInterpolator { + private static final ValueInterpolator INTERPOLATE_NUMERIC = + (a, b) -> t -> a * (1 - t) + b * t; + private static final DoubleUnaryOperator IDENTITY = x -> x; + private final ValueInterpolator valueInterpolator; + + DoubleUnaryOperator transform = IDENTITY; + DoubleUnaryOperator reverseTransform = IDENTITY; + boolean clamp = false; + private V defaultValue; + final DoubleArrayList domain = new DoubleArrayList(); + final List range = new ArrayList<>(); + private DoubleFunction fn; + double minKey = Double.POSITIVE_INFINITY; + double maxKey = Double.NEGATIVE_INFINITY; + + protected Interpolator(ValueInterpolator valueInterpolator) { + this.valueInterpolator = valueInterpolator; + } + + T setTransforms(DoubleUnaryOperator forward, DoubleUnaryOperator reverse) { + fn = null; + this.transform = forward; + this.reverseTransform = reverse; + return self(); + } + + @Override + public V apply(double operand) { + if (clamp) { + operand = Math.clamp(operand, minKey, maxKey); + } + if (Double.isNaN(operand)) { + return defaultValue; + } + if (fn == null) { + fn = rescale(); + } + return fn.apply(transform.applyAsDouble(operand)); + } + + public double applyAsDouble(double operand) { + return apply(operand) instanceof Number n ? n.doubleValue() : Double.NaN; + } + + private DoubleFunction rescale() { + if (domain.size() > 2) { + int j = Math.min(domain.size(), range.size()) - 1; + DoubleUnaryOperator[] d = new DoubleUnaryOperator[j]; + DoubleFunction[] r = new DoubleFunction[j]; + int i = -1; + + double[] domainItems = new double[domain.size()]; + for (int k = 0; k < domainItems.length; k++) { + domainItems[k] = transform.applyAsDouble(domain.get(k)); + } + List rangeItems = domainItems[j] < domainItems[0] ? range.reversed() : range; + + // Reverse descending domains. + if (domainItems[j] < domainItems[0]) { + ArrayUtils.reverse(domainItems); + } + + while (++i < j) { + d[i] = normalize(domainItems[i], domainItems[i + 1]); + r[i] = valueInterpolator.apply(rangeItems.get(i), rangeItems.get(i + 1)); + } + + return x -> { + int ii = bisect(domainItems, x, 1, j) - 1; + return r[ii].apply(d[ii].applyAsDouble(x)); + }; + } else { + double d0 = transform.applyAsDouble(domain.get(0)), d1 = transform.applyAsDouble(domain.get(1)); + V r0 = range.get(0), r1 = range.get(1); + boolean reverse = d1 < d0; + final double dlo = reverse ? d1 : d0; + final double dhi = reverse ? d0 : d1; + final V rlo = reverse ? r1 : r0; + final V rhi = reverse ? r0 : r1; + DoubleUnaryOperator normalize = normalize(dlo, dhi); + DoubleFunction interpolate = valueInterpolator.apply(rlo, rhi); + return x -> interpolate.apply(normalize.applyAsDouble(x)); + } + } + + private static int bisect(double[] a, double x, int lo, int hi) { + if (lo < hi) { + do { + int mid = (lo + hi) >>> 1; + if (a[mid] <= x) + lo = mid + 1; + else + hi = mid; + } while (lo < hi); + } + return lo; + } + + private static DoubleUnaryOperator interpolate(double a, double b) { + return t -> a * (1 - t) + b * t; + } + + private static DoubleUnaryOperator normalize(double a, double b) { + double delta = b - a; + return delta == 0 ? x -> 0.5 : Double.isNaN(delta) ? x -> Double.NaN : x -> (x - a) / delta; + } + + @SuppressWarnings("unchecked") + public T self() { + return (T) this; + } + + public T clamp(boolean clamp) { + this.clamp = clamp; + return self(); + } + + public T defaultValue(V value) { + this.defaultValue = value; + return self(); + } + + public T put(double stop, V value) { + fn = null; + minKey = Math.min(stop, minKey); + maxKey = Math.max(stop, maxKey); + domain.add(stop); + range.add(value); + return self(); + } + + // TODO + // private static class Inverted extends Interpolator {} + + public static > DoubleUnaryOperator invertIt( + Interpolator interpolator) { + var result = linear(); + int j = Math.min(interpolator.domain.size(), interpolator.range.size()); + for (int i = 0; i < j; i++) { + result.put(interpolator.range.get(i).doubleValue(), + interpolator.transform.applyAsDouble(interpolator.domain.get(i))); + } + DoubleUnaryOperator retVal = + interpolator.reverseTransform == IDENTITY ? result::applyAsDouble : + x -> interpolator.reverseTransform.applyAsDouble(result.apply(x)); + return interpolator.clamp ? retVal.andThen(x -> Math.clamp(x, interpolator.minKey, interpolator.maxKey)) : retVal; + } + + public static class Power, V> extends Interpolator { + + public Power(ValueInterpolator valueInterpolator) { + super(valueInterpolator); + } + + public T exponent(double exponent) { + double inverse = 1d / exponent; + return setTransforms(x -> Math.pow(x, exponent), y -> Math.pow(y, inverse)); + } + + public static class Numeric extends Power implements IInterpolator.Continuous { + + public Numeric(ValueInterpolator valueInterpolator) { + super(valueInterpolator); + } + } + public static class Other extends Power, V> { + + public Other(ValueInterpolator valueInterpolator) { + super(valueInterpolator); + } + } + } + + public static class Log, V> extends Interpolator { + + public Log(ValueInterpolator valueInterpolator) { + super(valueInterpolator); + } + + private static DoubleUnaryOperator log(double base) { + double logBase = Math.log(base); + return x -> Math.log(x) / logBase; + } + + public T base(double base) { + DoubleUnaryOperator forward = + base == 10d ? Math::log10 : + base == Math.E ? Math::log : + log(base); + DoubleUnaryOperator reverse = + base == Math.E ? Math::exp : + x -> Math.pow(base, x); + return setTransforms(forward, reverse); + } + + public static class Numeric extends Log implements IInterpolator.Continuous { + + public Numeric(ValueInterpolator valueInterpolator) { + super(valueInterpolator); + } + } + public static class Other extends Log, V> { + + public Other(ValueInterpolator valueInterpolator) { + super(valueInterpolator); + } + } + } + + public static class Linear, V> extends Interpolator { + + public Linear(ValueInterpolator valueInterpolator) { + super(valueInterpolator); + } + + public static class Numeric extends Linear + implements IInterpolator.Continuous { + + public Numeric(ValueInterpolator valueInterpolator) { + super(valueInterpolator); + } + } + public static class Other extends Linear, V> { + + public Other(ValueInterpolator valueInterpolator) { + super(valueInterpolator); + } + } + } + + public static Linear.Numeric linear(ValueInterpolator valueInterpolator) { + return new Linear.Numeric(valueInterpolator); + } + + public static Linear.Numeric linear() { + return new Linear.Numeric(INTERPOLATE_NUMERIC); + } + + public static Log.Numeric log(ValueInterpolator valueInterpolator) { + return new Log.Numeric(valueInterpolator).base(10); + } + + public static Log.Numeric log() { + return log(INTERPOLATE_NUMERIC).base(10); + } + + + // TODO specialized double implementation without boxing (interpolator, values?) ? + // TODO set special args before returning reusable thing? forget self type + + public static Power.Numeric npower(ValueInterpolator valueInterpolator) { + return new Power.Numeric(valueInterpolator).exponent(1); + } + + public static Power.Numeric power() { + return npower(INTERPOLATE_NUMERIC).exponent(1); + } + + public static Power.Numeric sqrt() { + return power().exponent(0.5); + } +} diff --git a/planetiler-core/src/main/java/com/onthegomap/planetiler/util/Scales.java b/planetiler-core/src/main/java/com/onthegomap/planetiler/util/Scales.java new file mode 100644 index 0000000000..5eab80e51c --- /dev/null +++ b/planetiler-core/src/main/java/com/onthegomap/planetiler/util/Scales.java @@ -0,0 +1,372 @@ +package com.onthegomap.planetiler.util; + +import com.carrotsearch.hppc.DoubleArrayList; +import java.util.ArrayList; +import java.util.List; +import java.util.function.BiFunction; +import java.util.function.DoubleFunction; +import java.util.function.DoubleUnaryOperator; +import org.apache.commons.lang3.ArrayUtils; + +public class Scales { + private static final Interpolator INTERPOLATE_NUMERIC = + (a, b) -> t -> a * (1 - t) + b * t; + + public static ThresholdScale quantize(int min, int max, V minValue, V... values) { + ThresholdScale result = threshold(minValue); + int n = values.length; + for (int i = 0; i < n; i++) { + result.putAbove(((i + 1d) * max - (i - n) * min) / (n + 1), values[i]); + } + return result; + } + + public static ThresholdScale quantize(int min, int max, List values) { + ThresholdScale result = threshold(values.getFirst()); + int n = values.size() - 1; + int i = -1; + while (++i < n) { + result.putAbove(((i + 1d) * max - (i - n) * min) / (n + 1), values.get(i + 1)); + } + return result; + } + + + public interface Interpolator extends BiFunction> {} + + private interface Scale, V> extends DoubleFunction { + + T defaultValue(V defaultValue); + + @SuppressWarnings("unchecked") + default T self() { + return (T) this; + } + } + + private interface Continuous, V> extends Scale { + + T put(double x, V y); + + T clamp(boolean clamp); + } + + private interface Threshold, V> extends Scale { + + T putAbove(double x, V y); + + double invertMin(V x); + + double invertMax(V x); + + Extent invertExtent(V x); + } + + public record Extent(double min, double max) {} + + public static class ThresholdScale implements Threshold, V> { + + private V defaultValue = null; + private final List domain = new ArrayList<>(); + private final List range = new ArrayList<>(); + private double[] domainArray = null; + + private ThresholdScale(V minValue) { + range.add(minValue); + } + + @Override + public ThresholdScale putAbove(double x, V y) { + domainArray = null; + domain.add(x); + range.add(y); + return self(); + } + + @Override + public double invertMin(V x) { + int idx = range.indexOf(x); + return idx < 0 ? Double.NaN : idx == 0 ? Double.NEGATIVE_INFINITY : domain.get(idx - 1); + } + + @Override + public double invertMax(V x) { + int idx = range.indexOf(x); + return idx < 0 ? Double.NaN : idx == range.size() - 1 ? Double.POSITIVE_INFINITY : domain.get(idx); + } + + @Override + public Extent invertExtent(V x) { + int idx = range.indexOf(x); + return new Extent( + idx < 0 ? Double.NaN : idx == 0 ? Double.NEGATIVE_INFINITY : domain.get(idx - 1), + idx < 0 ? Double.NaN : idx == range.size() - 1 ? Double.POSITIVE_INFINITY : domain.get(idx) + ); + } + + @Override + public ThresholdScale defaultValue(V defaultValue) { + this.defaultValue = defaultValue; + return self(); + } + + @Override + public V apply(double x) { + if (domainArray == null) { + domainArray = new double[domain.size()]; + int i = 0; + for (Double d : domain) { + domainArray[i++] = d; + } + } + if (Double.isNaN(x)) { + return defaultValue; + } + int idx = bisect(domainArray, x, 0, Math.min(domainArray.length, range.size() - 1)); + return range.get(idx); + } + } + + private interface ContinuousDoubleScale> + extends Continuous, DoubleUnaryOperator { + + default T put(double x, double y) { + return put(x, Double.valueOf(y)); + } + + @Override + default double applyAsDouble(double operand) { + return apply(operand); + } + + DoubleUnaryOperator invert(); + } + + + private static class BaseContinuous, V> implements Continuous { + static final DoubleUnaryOperator IDENTITY = x -> x; + + private final Interpolator valueInterpolator; + final DoubleUnaryOperator transform; + boolean clamp = false; + private V defaultValue; + final DoubleArrayList domain = new DoubleArrayList(); + final List range = new ArrayList<>(); + private DoubleFunction fn; + double minKey = Double.POSITIVE_INFINITY; + double maxKey = Double.NEGATIVE_INFINITY; + + private BaseContinuous(Interpolator valueInterpolator, DoubleUnaryOperator transform) { + this.valueInterpolator = valueInterpolator; + this.transform = transform; + } + + @Override + public V apply(double operand) { + if (clamp) { + operand = Math.clamp(operand, minKey, maxKey); + } + if (Double.isNaN(operand)) { + return defaultValue; + } + if (fn == null) { + fn = rescale(); + } + return fn.apply(transform.applyAsDouble(operand)); + } + + private DoubleFunction rescale() { + if (domain.size() > 2) { + int j = Math.min(domain.size(), range.size()) - 1; + DoubleUnaryOperator[] d = new DoubleUnaryOperator[j]; + DoubleFunction[] r = new DoubleFunction[j]; + int i = -1; + + double[] domainItems = new double[domain.size()]; + for (int k = 0; k < domainItems.length; k++) { + domainItems[k] = transform.applyAsDouble(domain.get(k)); + } + List rangeItems = domainItems[j] < domainItems[0] ? range.reversed() : range; + + // Reverse descending domains. + if (domainItems[j] < domainItems[0]) { + ArrayUtils.reverse(domainItems); + } + + while (++i < j) { + d[i] = normalize(domainItems[i], domainItems[i + 1]); + r[i] = valueInterpolator.apply(rangeItems.get(i), rangeItems.get(i + 1)); + } + + return x -> { + int ii = bisect(domainItems, x, 1, j) - 1; + return r[ii].apply(d[ii].applyAsDouble(x)); + }; + } else { + double d0 = transform.applyAsDouble(domain.get(0)), d1 = transform.applyAsDouble(domain.get(1)); + V r0 = range.get(0), r1 = range.get(1); + boolean reverse = d1 < d0; + final double dlo = reverse ? d1 : d0; + final double dhi = reverse ? d0 : d1; + final V rlo = reverse ? r1 : r0; + final V rhi = reverse ? r0 : r1; + DoubleFunction interpolate = valueInterpolator.apply(rlo, rhi); + return x -> { + double value = (x - dlo) / (dhi - dlo); + return interpolate.apply(value); + }; + } + } + + private static DoubleUnaryOperator normalize(double a, double b) { + double delta = b - a; + return delta == 0 ? x -> 0.5 : Double.isNaN(delta) ? x -> Double.NaN : x -> (x - a) / delta; + } + + @Override + public T put(double x, V y) { + fn = null; + minKey = Math.min(x, minKey); + maxKey = Math.max(x, maxKey); + domain.add(x); + range.add(y); + return self(); + } + + @Override + public T clamp(boolean clamp) { + this.clamp = clamp; + return self(); + } + + @Override + public T defaultValue(V defaultValue) { + this.defaultValue = defaultValue; + return self(); + } + } + + + private static int bisect(double[] a, double x, int lo, int hi) { + if (lo < hi) { + do { + int mid = (lo + hi) >>> 1; + if (a[mid] <= x) + lo = mid + 1; + else + hi = mid; + } while (lo < hi); + } + return lo; + } + + public static class DoubleContinuous extends BaseContinuous + implements ContinuousDoubleScale { + private final DoubleUnaryOperator reverseTransform; + + private DoubleContinuous(Interpolator valueInterpolator, DoubleUnaryOperator transform, + DoubleUnaryOperator reverseTransform) { + super(valueInterpolator, transform); + this.reverseTransform = reverseTransform; + } + + @Override + public DoubleUnaryOperator invert() { + var result = linear(); + int j = Math.min(domain.size(), range.size()); + for (int i = 0; i < j; i++) { + result.put(range.get(i).doubleValue(), + transform.applyAsDouble(domain.get(i))); + } + DoubleUnaryOperator retVal = + reverseTransform == IDENTITY ? result : x -> reverseTransform.applyAsDouble(result.apply(x)); + return clamp ? retVal.andThen(x -> Math.clamp(x, minKey, maxKey)) : retVal; + } + } + + public static class OtherContinuous extends BaseContinuous, V> { + + private OtherContinuous(Interpolator valueInterpolator, DoubleUnaryOperator transform) { + super(valueInterpolator, transform); + } + } + + public static DoubleContinuous linear() { + return new DoubleContinuous(INTERPOLATE_NUMERIC, BaseContinuous.IDENTITY, BaseContinuous.IDENTITY); + } + + public static DoubleContinuous linearDouble(Interpolator interp) { + return new DoubleContinuous(interp, BaseContinuous.IDENTITY, BaseContinuous.IDENTITY); + } + + public static OtherContinuous linear(Interpolator valueInterpolator) { + return new OtherContinuous<>(valueInterpolator, BaseContinuous.IDENTITY); + } + + + public static DoubleContinuous power(double exponent) { + double inverse = 1d / exponent; + return new DoubleContinuous(INTERPOLATE_NUMERIC, + x -> Math.pow(x, exponent), + y -> Math.pow(y, inverse) + ); + } + + public static OtherContinuous power(double exponent, Interpolator valueInterpolator) { + return new OtherContinuous<>(valueInterpolator, x -> Math.pow(x, exponent)); + } + + public static DoubleContinuous sqrt() { + return power(0.5); + } + + public static OtherContinuous sqrt(Interpolator valueInterpolator) { + return power(0.5, valueInterpolator); + } + + + private static DoubleUnaryOperator logBase(double base) { + double logBase = Math.log(base); + return x -> Math.log(x) / logBase; + } + + private static DoubleUnaryOperator expFn(double base) { + return base == Math.E ? Math::exp : + x -> Math.pow(base, x); + } + + private static DoubleUnaryOperator logFn(double base) { + return base == 10d ? Math::log10 : + base == Math.E ? Math::log : + logBase(base); + } + + public static DoubleContinuous log(double base) { + return new DoubleContinuous(INTERPOLATE_NUMERIC, logFn(base), expFn(base)); + } + + public static OtherContinuous log(double base, Interpolator valueInterpolator) { + return new OtherContinuous<>(valueInterpolator, logFn(base)); + } + + public static DoubleContinuous exponential(double base) { + return new DoubleContinuous(INTERPOLATE_NUMERIC, expFn(base), logFn(base)); + } + + public static OtherContinuous exponential(double base, Interpolator valueInterpolator) { + return new OtherContinuous<>(valueInterpolator, expFn(base)); + } + + public static OtherContinuous bezier(double x1, double y1, double x2, double y2) { + var bezier = new UnitBezier(x1, y1, x2, y2); + return new OtherContinuous<>((a, b) -> { + double diff = b - a; + return t -> bezier.solve(t) * diff + a; + }, BaseContinuous.IDENTITY); + } + + + public static ThresholdScale threshold(V minValue) { + return new ThresholdScale<>(minValue); + } +} diff --git a/planetiler-core/src/main/java/com/onthegomap/planetiler/util/UnitBezier.java b/planetiler-core/src/main/java/com/onthegomap/planetiler/util/UnitBezier.java new file mode 100644 index 0000000000..1425573c1d --- /dev/null +++ b/planetiler-core/src/main/java/com/onthegomap/planetiler/util/UnitBezier.java @@ -0,0 +1,94 @@ +package com.onthegomap.planetiler.util; + +public class UnitBezier { + private final double ax; + private final double bx; + private final double cx; + + private final double ay; + private final double by; + private final double cy; + + public UnitBezier(double p1x, double p1y, double p2x, double p2y) { + // Calculate the polynomial coefficients, implicit first and last control points are (0,0) and (1,1). + cx = 3.0 * p1x; + bx = 3.0 * (p2x - p1x) - cx; + ax = 1.0 - cx - bx; + + cy = 3.0 * p1y; + by = 3.0 * (p2y - p1y) - cy; + ay = 1.0 - cy - by; + } + + double sampleCurveX(double t) { + // `ax t^3 + bx t^2 + cx t' expanded using Horner's rule. + return ((ax * t + bx) * t + cx) * t; + } + + double sampleCurveY(double t) { + return ((ay * t + by) * t + cy) * t; + } + + double sampleCurveDerivativeX(double t) { + return (3.0 * ax * t + 2.0 * bx) * t + cx; + } + + // Given an x value, find a parametric value it came from. + double solveCurveX(double x, double epsilon) { + double t0; + double t1; + double t2; + double x2; + double d2; + int i; + + if (x < 0) + return 0; + if (x > 1) + return 1; + + // First try a few iterations of Newton's method -- normally very fast. + for (t2 = x, i = 0; i < 8; i++) { + x2 = sampleCurveX(t2) - x; + if (Math.abs(x2) < epsilon) + return t2; + d2 = sampleCurveDerivativeX(t2); + if (Math.abs(d2) < 1e-6) + break; + t2 = t2 - x2 / d2; + } + + // Fall back to the bisection method for reliability. + t0 = 0.0; + t1 = 1.0; + t2 = x; + + if (t2 < t0) + return t0; + if (t2 > t1) + return t1; + + while (t0 < t1) { + x2 = sampleCurveX(t2); + if (Math.abs(x2 - x) < epsilon) + return t2; + if (x > x2) + t0 = t2; + else + t1 = t2; + t2 = (t1 - t0) * .5 + t0; + } + + // Failure. + return t2; + } + + double solve(double x, double epsilon) { + return sampleCurveY(solveCurveX(x, epsilon)); + } + + + double solve(double x) { + return solve(x, 1e-6); + } +} diff --git a/planetiler-core/src/test/java/com/onthegomap/planetiler/util/InterpolatorTest.java b/planetiler-core/src/test/java/com/onthegomap/planetiler/util/InterpolatorTest.java new file mode 100644 index 0000000000..91c106dbe3 --- /dev/null +++ b/planetiler-core/src/test/java/com/onthegomap/planetiler/util/InterpolatorTest.java @@ -0,0 +1,236 @@ +package com.onthegomap.planetiler.util; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +import java.util.List; +import org.junit.jupiter.api.Test; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.ValueSource; + +class InterpolatorTest { + @ParameterizedTest + @ValueSource(booleans = {false, true}) + void testLinear(boolean clamp) { + var interpolator = Scales.linear() + .put(0, 10) + .put(1, 20) + .clamp(clamp); + assertTransform(interpolator, 0, 10); + assertTransform(interpolator, 1, 20); + assertTransform(interpolator, 0.5, 15); + // clamp + assertClose(clamp ? 20 : 30, interpolator.applyAsDouble(2)); + assertClose(clamp ? 1 : 2, interpolator.invert().applyAsDouble(30)); + assertClose(clamp ? 10 : 0, interpolator.applyAsDouble(-1)); + assertClose(clamp ? 0 : -1, interpolator.invert().applyAsDouble(0)); + } + + @ParameterizedTest + @ValueSource(booleans = {false, true}) + void testLog(boolean clamp) { + var interpolator = Scales.log(10) + .put(10, 1d) + .put(1_000, 2d) + .clamp(clamp); + assertTransform(interpolator, 10, 1); + assertTransform(interpolator, 100, 1.5); + assertTransform(interpolator, 1_000, 2); + // clamp + assertClose(clamp ? 1 : 0.5, interpolator.applyAsDouble(1)); + assertClose(clamp ? 10 : 1, interpolator.invert().applyAsDouble(0.5)); + assertClose(clamp ? 2 : 2.5, interpolator.applyAsDouble(10_000)); + assertClose(clamp ? 1_000 : 10_000, interpolator.invert().applyAsDouble(2.5)); + } + + @ParameterizedTest + @ValueSource(booleans = {false, true}) + void testLog2(boolean clamp) { + var interpolator = Scales.log(2) + .put(2, 1d) + .put(8, 2d) + .clamp(clamp); + assertTransform(interpolator, 2, 1); + assertTransform(interpolator, 4, 1.5); + assertTransform(interpolator, 8, 2); + // clamp + assertClose(clamp ? 1 : 0.5, interpolator.applyAsDouble(1)); + assertClose(clamp ? 2 : 1, interpolator.invert().applyAsDouble(0.5)); + assertClose(clamp ? 2 : 2.5, interpolator.applyAsDouble(16)); + assertClose(clamp ? 8 : 16, interpolator.invert().applyAsDouble(2.5)); + } + + @ParameterizedTest + @ValueSource(booleans = {false, true}) + void testPower(boolean clamp) { + var interpolator = Scales.power(2) + .put(Math.sqrt(2), 1d) + .put(Math.sqrt(4), 2d) + .clamp(clamp); + assertTransform(interpolator, Math.sqrt(2), 1); + assertTransform(interpolator, Math.sqrt(3), 1.5); + assertTransform(interpolator, Math.sqrt(4), 2); + // clamp + assertClose(clamp ? 1 : 0.5, interpolator.applyAsDouble(1)); + assertClose(clamp ? Math.sqrt(2) : 1, interpolator.invert().applyAsDouble(0.5)); + assertClose(clamp ? 2 : 2.5, interpolator.applyAsDouble(Math.sqrt(5))); + assertClose(clamp ? Math.sqrt(4) : Math.sqrt(5), interpolator.invert().applyAsDouble(2.5)); + } + + @Test + void testNegativePower() { + var interpolator = Scales.power(-1) + .put(1, 0d) + .put(2, 1d); + assertTransform(interpolator, 1, 0); + assertTransform(interpolator, 1.5, 2d / 3); + assertTransform(interpolator, 2, 1); + } + + @ParameterizedTest + @ValueSource(booleans = {false, true}) + void testSqrt(boolean clamp) { + var interpolator = Scales.sqrt() + .put(4, 1d) + .put(16, 2d) + .clamp(clamp); + assertTransform(interpolator, 4, 1); + assertTransform(interpolator, 9, 1.5); + assertTransform(interpolator, 16, 2); + // clamp + assertClose(clamp ? 1 : 0.5, interpolator.applyAsDouble(1)); + assertClose(clamp ? 4 : 1, interpolator.invert().applyAsDouble(0.5)); + assertClose(clamp ? 2 : 2.5, interpolator.applyAsDouble(25)); + assertClose(clamp ? 16 : 25, interpolator.invert().applyAsDouble(2.5)); + } + + @ParameterizedTest + @ValueSource(booleans = {false, true}) + void testMultipartDescending(boolean clamp) { + var interpolator = Scales.linear() + .put(4, 1d) + .put(2, 2d) + .put(1, 4d) + .clamp(clamp); + assertTransform(interpolator, 4, 1); + assertTransform(interpolator, 3, 1.5); + assertTransform(interpolator, 2, 2); + assertTransform(interpolator, 1.5, 3); + assertTransform(interpolator, 1, 4); + // clamp + assertClose(clamp ? 1 : 0.5, interpolator.applyAsDouble(5)); + assertClose(clamp ? 4 : 5, interpolator.invert().applyAsDouble(0.5)); + assertClose(clamp ? 4 : 6, interpolator.applyAsDouble(0)); + assertClose(clamp ? 1 : 0, interpolator.invert().applyAsDouble(6)); + } + + @ParameterizedTest + @ValueSource(booleans = {false, true}) + void testMultipartAscending(boolean clamp) { + var interpolator = Scales.linear() + .put(1, 4d) + .put(2, 2d) + .put(4, 1d) + .clamp(clamp); + assertTransform(interpolator, 1, 4); + assertTransform(interpolator, 1.5, 3); + assertTransform(interpolator, 2, 2); + assertTransform(interpolator, 3, 1.5); + assertTransform(interpolator, 4, 1); + // clamp + assertClose(clamp ? 1 : 0.5, interpolator.invert().applyAsDouble(5)); + assertClose(clamp ? 4 : 5, interpolator.applyAsDouble(0.5)); + assertClose(clamp ? 4 : 6, interpolator.invert().applyAsDouble(0)); + assertClose(clamp ? 1 : 0, interpolator.applyAsDouble(6)); + } + + @Test + void testUnknown() { + var interpolator = Scales.linear() + .put(0, 10d) + .put(1, 20d) + .defaultValue(100d); + assertEquals(100, interpolator.applyAsDouble(Double.NaN)); + } + + @Test + void testThreshold() { + var scale = Scales.threshold("a") + .putAbove(1.5, "b") + .putAbove(3, "c"); + assertEquals("a", scale.apply(1.4)); + assertEquals("b", scale.apply(1.5)); + assertEquals("b", scale.apply(2.99)); + assertEquals("c", scale.apply(3)); + assertEquals("c", scale.apply(3.1)); + + testInvert(scale, "a", Double.NEGATIVE_INFINITY, 1.5); + testInvert(scale, "b", 1.5, 3); + testInvert(scale, "c", 3, Double.POSITIVE_INFINITY); + testInvert(scale, "d", Double.NaN, Double.NaN); + } + + @ParameterizedTest + @ValueSource(booleans = {false, true}) + void testQuantize(boolean list) { + var scale = list ? + Scales.quantize(0, 3, List.of("a", "b", "c")) : + Scales.quantize(0, 3, "a", "b", "c"); + assertEquals("a", scale.apply(0)); + assertEquals("a", scale.apply(0.99)); + assertEquals("b", scale.apply(1)); + assertEquals("b", scale.apply(1.01)); + assertEquals("b", scale.apply(1.99)); + assertEquals("c", scale.apply(2)); + assertEquals("c", scale.apply(2.01)); + assertEquals("c", scale.apply(99)); + } + + @Test + void testExponential() { + var scale = Scales.exponential(2) + .put(1, 2) + .put(3, 6) + .clamp(true); + + assertClose(2, scale.apply(0)); + assertClose(2, scale.apply(1)); + assertClose(3.3333333333, scale.apply(2)); + assertClose(6, scale.apply(3)); + assertClose(6, scale.apply(5)); + } + + @Test + void testBezier() { + var scale = Scales.bezier(0.42, 0, 0.58, 1) + .put(0, 0d) + .put(100, 100d) + .clamp(true); + + assertEquals(0, scale.apply(0), 1e-4); + assertEquals(1.97224, scale.apply(10), 1e-4); + assertEquals(8.16597, scale.apply(20), 1e-4); + assertEquals(18.7395, scale.apply(30), 1e-4); + assertEquals(33.1883, scale.apply(40), 1e-4); + assertEquals(50, scale.apply(50), 1e-4); + assertEquals(66.8116, scale.apply(60), 1e-4); + assertEquals(81.2604, scale.apply(70), 1e-4); + assertEquals(91.834, scale.apply(80), 1e-4); + assertEquals(98.0277, scale.apply(90), 1e-4); + assertEquals(100, scale.apply(100), 1e-4); + } + + private static void testInvert(Scales.ThresholdScale scale, V val, double min, double max) { + assertEquals(min, scale.invertMin(val)); + assertEquals(max, scale.invertMax(val)); + assertEquals(max, scale.invertMax(val)); + } + + private static void assertClose(double expected, double actual) { + assertEquals(expected, actual, 1e-10); + } + + private static void assertTransform(Scales.DoubleContinuous interp, double x, double y) { + assertClose(y, interp.applyAsDouble(x)); + assertClose(x, interp.invert().applyAsDouble(y)); + } +}