sharpetronics/Map-Projections
1
0
Fork 0
mirror of https://github.com/csharpee/Map-Projections.git synced 2025-08-31 00:00:32 -04:00
Code Issues Projects Releases Wiki Activity

Good documentation is its own reward

Browse Source 
I tidied up the documentation on the BFGS stuff and updated my maps and
defaults to be informed by the new optimizations.
This commit is contained in:
Justin Kunimune 2018-01-31 13:22:22 -10:00
parent b51f632cf5
commit 5e01e339b9
9 changed files with 65 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
output/Unbiased Political.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 207 KiB

BIN
output/UnbiasedPolitical.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 268 KiB

BIN
output/graph - optimizer.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 17 KiB

After

Width:  |  Height:  |  Size: 27 KiB

BIN
output/graph - plotter.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 46 KiB

After

Width:  |  Height:  |  Size: 46 KiB

32
output/parameters.txt
View File

@ -1,12 +1,44 @@
We got the best Tobler Hyperelliptical projections using:
t0=30.506558786136164; t1=0.3975383525206906; t2=3.997720064276123; (2.5245792765063137E-4, 0.49042998348577777)
t0=30.53431423862674; t1=0.4062974245671901; t2=3.951968307864351; (2.3634131999763095E-4, 0.490413331660892)
t0=30.534996054570122; t1=0.40982217940437315; t2=3.9388293231446183; (2.3084286442176676E-4, 0.49040541967506746)
t0=30.526861993190156; t1=0.41289331825953984; t2=3.948490584936325; (2.2823612925236138E-4, 0.4903943344342276)
t0=30.171506742208667; t1=0.5447325922534694; t2=3.3045726206106174; (7.472092115181078E-5, 0.4905648364049878)
t0=30.472407744588843; t1=0.4624100943196009; t2=3.7058554158579327; (1.6026863026488797E-4, 0.4903554923556594)
t0=30.772855981645655; t1=0.3505145046575816; t2=4.142699303717857; (2.0497084855872877E-4, 0.4905873681746431)
t0=30.597300348860273; t1=0.34295825057460816; t2=4.058216809727585; (2.1112594263567147E-4, 0.49068831665768275)
t0=0.0; t1=1.0; t2=1.0; (2.3111608934529733E-8, 0.6137565646410882)
We got the best Winkel Tripel projections using:
t0=17.957933140628256; (0.25932019508548493, 0.36444804618226867)
t0=19.578229322577233; (0.2586497441950308, 0.3644980527902062)
t0=21.663672947605303; (0.2576915304376235, 0.36472868044150314)
t0=24.552760247498448; (0.2561800301642924, 0.36543738102508005)
t0=28.935732586336282; (0.2534534735464852, 0.367639315334747)
t0=36.48524369490269; (0.2473866667743297, 0.3758665329854081)
t0=89.64383463941283; (0.07078405279816037, 0.8459966547063894)
t0=90.20684487459667; (0.06743617153928232, 0.8587166159021128)
t0=90.31210519110017; (0.06728093742152608, 0.8611571986133976)
We got the best TetraPower projections using:
t0=0.7251510971705017; t1=0.8799415807400688; t2=0.8572850133622626; (0.422750395775632, 0.13826331216471252)
t0=0.7442516235907949; t1=0.9031123519136701; t2=0.8650260337289639; (0.41666407130907396, 0.13643981232990932)
t0=0.7995356266967016; t1=0.877198483919957; t2=0.8775689219350017; (0.4061846891910591, 0.13545117616364158)
t0=0.843813906814002; t1=0.9055863648561882; t2=0.9068219650428548; (0.39168587855607145, 0.13590934103565014)
t0=1.0887347058906243; t1=1.1906232589287455; t2=1.0068156358843527; (0.3136201572259122, 0.20723377252523384)
t0=1.8578831572159344; t1=1.3230427148676347; t2=1.3541168911984933; (0.09659313794029525, 0.5053885132535845)
t0=1.8609065719598639; t1=1.2417845624360808; t2=1.4437856176727957; (0.07501424763978838, 0.5451384929588722)
t0=1.8892521212625646; t1=1.1793474367429049; t2=1.5043789437785882; (0.06455989310517463, 0.5775973717174138)
t0=1.9077945281918272; t1=1.1571046889875514; t2=1.5308152214270465; (0.062421585098500404, 0.5926588448712382)
We got the best AuthaPower projections using:
t0=0.5500000000000003; (0.41312749977668417, 0.08828966383834298)
t0=0.5500000000000003; (0.41312749977668417, 0.08828966383834298)
t0=0.5500000000000003; (0.41312749977668417, 0.08828966383834298)
t0=0.5750000000000003; (0.3876223544028035, 0.0968718822112939)
t0=0.6000000000000003; (0.3621275092593102, 0.1161981251407432)
t0=0.9250000000000006; (0.04143197931427, 0.5374514903103988)
t0=0.9500000000000006; (0.030342030927343812, 0.5653681878864117)
t0=0.9500000000000006; (0.030342030927343812, 0.5653681878864117)
t0=0.9500000000000006; (0.030342030927343812, 0.5653681878864117)

66
src/apps/MapOptimizer.java
View File

@ -57,14 +57,13 @@ import utils.linalg.Vector;
public class MapOptimizer extends Application {
private static final Projection[] EXISTING_PROJECTIONS = { Cylindrical.BEHRMANN,
Arbitrary.ROBINSON, Cylindrical.GALL_STEREOGRAPHIC, Misc.PEIRCE_QUINCUNCIAL };
Arbitrary.ROBINSON, Cylindrical.PLATE_CARREE, Cylindrical.GALL_STEREOGRAPHIC,
Misc.PEIRCE_QUINCUNCIAL };
private static final Projection[] PROJECTIONS_TO_OPTIMIZE = { Tobler.TOBLER,
WinkelTripel.WINKEL_TRIPEL, Polyhedral.TETRAPOWER, Polyhedral.AUTHAPOWER };
// private static final double[] WEIGHTS = { 0., .125, .25, .375, .5, .625, .75, .875, 1. };
private static final double[] WEIGHTS = {0};
private static final int NUM_SAMPLE = 1;
private static final double[] WEIGHTS = { 0., .125, .25, .375, .5, .625, .75, .875, 1. };
private static final int NUM_BRUTE_FORCE = 30;
private static final int NUM_BFGS_ITERATE = 6;
private static final int NUM_BFGS_ITERATE = 7;
private static final double GOLDSTEIN_C = 0.5;
private static final double BACKTRACK_TAU = 0.5;
private static final double BACKTRACK_ALF0 = 4;
@ -126,17 +125,8 @@ public class MapOptimizer extends Application {
}
private static final double weighDistortion(Projection proj, double[] params, double weight) {
double areaDist = 0, anglDist = 0;
for (int i = 0; i < NUM_SAMPLE; i ++) {
double[] mcParams = new double[params.length];
for (int j = 0; j < mcParams.length; j ++)
mcParams[j] = params[j];// + (Math.random()-.5)*1e-3; //is it weird that I'm introducing stochasticity into this?
double[] distortions = proj.avgDistortion(GLOBE, mcParams);
areaDist += distortions[0];
anglDist += distortions[1];
}
return areaDist/NUM_SAMPLE*weight + anglDist/NUM_SAMPLE*(1-weight);
private static final double weighDistortion(double[] distortions, double weight) {
return distortions[0]*weight + distortions[1]*(1-weight);
}
@ -149,10 +139,10 @@ public class MapOptimizer extends Application {
for (int k = 0; k < WEIGHTS.length; k ++) {
final double weighFactor = WEIGHTS[k];
double[] currentBest = bruteForceMinimise(
(params) -> weighDistortion(proj, params, weighFactor),
(params) -> weighDistortion(proj.avgDistortion(GLOBE, params), weighFactor),
proj.getParameterValues());
best[k] = bfgsMinimise(
(params) -> weighDistortion(proj, params, weighFactor),
(params) -> weighDistortion(proj.avgDistortion(GLOBE, params), weighFactor),
currentBest);
}
@ -181,7 +171,6 @@ public class MapOptimizer extends Application {
* parameter sweep.
* @param func The function to minimise
* @param bounds Parameter limits for each argument
* @param numTries The approximate number of samples to take
* @return An array containing the best input to func that it found
*/
private static double[] bruteForceMinimise(Function<double[], Double> func, double[][] bounds) {
@ -222,6 +211,13 @@ public class MapOptimizer extends Application {
}
/**
* Calculates the set of parameters that minimises the function using BFGS optimisation with
* a backtracking line search
* @param arrFunction The function that takes a parameter array and returns a double value
* @param x0 The initial guess
* @return The array of parameters that mimimise arrFunction
*/
private static double[] bfgsMinimise(Function<double[], Double> arrFunction, double[] x0) { //The Broyden-Fletcher-Goldfarb-Shanno algorithm
System.out.println("BFGS = [");
final int n = x0.length;
@ -230,23 +226,9 @@ public class MapOptimizer extends Application {
Vector xk = new Vector(x0); //initial variable values
double fxk = func.apply(xk);
// System.out.println(hessian(func, xk, fxk));
// System.out.println(grad(func, xk, fxk));
// Matrix H = hessian(func, xk, fxk);
Matrix Binv = hessian(func, xk, fxk).inverse();
// Matrix Binv = Matrix.identity(n); //initial approximate Hessian inverse
Vector gradFxk = grad(func, xk, fxk); //function at current location
// System.out.println("anetauson!");
// System.out.print("A"+NUM_SAMPLE+" = [");
// for (double d = 0; d <= 1e-1; d += 5e-4) {
// xk = new Vector(17.8, -17.8);
// Vector dx = Vector.unit(0, n).times(d);
// double f = func.apply(xk.minus(new Vector(-1,1)).plus(dx));
// System.out.println(xk.minus(new Vector(-1,1)).plus(dx).getElement(0)+", "+f+";");
// }
// System.out.println("];");
for (int k = 0; k < NUM_BFGS_ITERATE; k ++) { //(I'm not sure how to test for convergence here, so I'm just running a set number of iterations)
Vector pk = Vector.fromMatrix(Binv.times(gradFxk)); //apply Newton's method for initial step direction
pk = pk.times(-Math.signum(pk.dot(gradFxk))); //but make sure it points downhill
@ -271,7 +253,6 @@ public class MapOptimizer extends Application {
Matrix a = I.minus(sk.times(yk.T()).times(1/yk.dot(sk)));
Matrix b = sk.times(sk.T()).times(1/yk.dot(sk));
// Binv = hessian(func, xkp1, fxkp1).inverse();
Binv = a.times(Binv).times(a.T()).plus(b); //update Binv
xk = xkp1;
@ -283,6 +264,13 @@ public class MapOptimizer extends Application {
}
/**
* Calculates the gradient of f at x
* @param f The function to differentiate
* @param x The point at which to differentiate
* @param fx The value of f(x), to speed computations
* @return
*/
private static Vector grad(Function<Vector, Double> f, Vector x, double fx) {
final int n = x.getLength();
Vector gradF = new Vector(n); //compute the gradient
@ -294,12 +282,18 @@ public class MapOptimizer extends Application {
}
for (double d: x.asArray())
System.out.print(d+", ");
// System.out.println(gradF.getElement(0)+", "+gradF.getElement(1)+";");
System.out.println(fx+";");
return gradF;
}
/**
* Computes the Hessian matrix of f at x
* @param f The function to differentiate
* @param x The point at which to differentiate
* @param fx The value of f(x), to aid in computation
* @return The Jacobian of the gradient, a symmetric Matrix of second derivatives
*/
private static Matrix hessian(Function<Vector, Double> f, Vector x, double fx) {
final int n = x.getLength();
@ -316,7 +310,6 @@ public class MapOptimizer extends Application {
values[k] = f.apply(x.plus(dx));
}
// System.out.println(Arrays.toString(values));
Matrix h = new Matrix(n, n);
for (int i = 0; i < n; i ++) { //compute the derivatives and fill the matrix
for (int j = i; j < n; j ++) {
@ -324,7 +317,6 @@ public class MapOptimizer extends Application {
int dxj = (int)Math.pow(3, j);
double dfdx0 = (values[dxi] - values[0])/DEL_X;
double dfdx1 = (values[dxi+dxj] - values[dxj])/DEL_X;
// System.out.println(dfdx0+", "+dfdx1);
double d2fdx2 = (dfdx1 - dfdx0)/DEL_X;
h.setElement(i, j, d2fdx2);
h.setElement(j, i, d2fdx2);

2
src/maps/Polyhedral.java
View File

@ -236,7 +236,7 @@ public class Polyhedral {
"TetraPower", "A parameterised tetrahedral projection that I invented.", 0b1111,
Configuration.TETRAHEDRON_WIDE_FACE, Property.COMPROMISE, 2,
new String[] {"k1","k2","k3"},
new double[][] {{.01,2.,.98},{.01,2.,1.2},{.01,2.,.98}}) {
new double[][] {{.01,2.,1.01},{.01,2.,1.19},{.01,2.,1.01}}) {
private double k1, k2, k3;

2
src/maps/Tobler.java
View File

@ -42,7 +42,7 @@ public class Tobler {
"Tobler Hyperelliptical", "An equal-area projection shaped like a hyperellipse.",
2*Math.PI, 0, 0b1001, Type.PSEUDOCYLINDRICAL, Property.EQUAL_AREA, 4,
new String[]{"Std. Parallel","alpha","K"},
new double[][] {{0,89,37}, {0,1,0}, {1,5,3.6}}) { //optimal parameters are 30.6,.50,3.63, but these defaults are more recognizably Tobler
new double[][] {{0,89,30}, {0,1,.46}, {1,5,3.7}}) { //optimal parameters are 30.6,.50,3.63, but these defaults are more recognizably Tobler
private static final int N = 20000;
private double alpha, kappa, epsilon; //epsilon is related to gamma, but defined somewhat differently

4
src/maps/WinkelTripel.java
View File

@ -48,7 +48,7 @@ public final class WinkelTripel {
public static final Projection WINKEL_TRIPEL =
new Projection("Winkel Tripel", "National Geographic's compromise projection of choice.",
0, 2*Math.PI, 0b1011, Type.OTHER, Property.COMPROMISE, 3,
new String[] {"Std. Parallel"}, new double[][] {{0,90,50.4598}}) {
new String[] {"Std. Parallel"}, new double[][] {{0,90,40}}) {
private double stdParallel;
@ -66,7 +66,7 @@ public final class WinkelTripel {
x, y,
y/2, x*(1 + Math.cos(y*Math.PI/2))/(2 + 2*Math.cos(stdParallel)), //inital guess is Eckert V
this::f1pX, this::f2pY,
this::df1dphi, this::df1dlam, this::df2dphi, this::df2dlam, .002);
this::df1dphi, this::df1dlam, this::df2dphi, this::df2dlam, .001);
}
private double f1pX(double phi, double lam) {