Geometric Harnessing of Precipitation Records: Reexamining Four Storms from Iowa City Online Appendix

Huai-Hsien Huang, Carlos E. Puente*, Andrea Cortis
Stochastic Environmental Research and Risk Assessment (2012)
DOI 10.1007/s00477-012-0617-6

This page contains supplementary information deemed not important enough to include in the paper but still of relevance to the interested reader. Here one will find the values of parameters used in making the figures, example MATLAB source code, and interactive demos of the fractal-multifractal encoding method.

Parameters

These tables only list free parameters, rounded to two decimal places. One can mouse over a table cell to see a slightly more precise value (four decimal places).

All values are normalized between 0 and 1 for simplicity when running the General PSO search and thus shown here as such. A, p, r(1), r(2), x live in [0, 1]; y, z in [−5, 5]; d in [−1, 1]; θ(1), θ(2), and ω in [−π/2, π/2].

The first and last of the x-, y-, and p-values are fixed at 0 and 1, respectively. Additionally the first z-value is fixed at 0 for the higher-dimensonal marginals. Note that the p-values delineate the bounds of the probability of using a map. For example, the probabilities associated with the three maps in the first wire example are 26.44%, 30.70%, and 42.86%.

Tables marked in gray denote results that were not shown in the paper.

Iowa A

FM-Wire (maps: 3)
x y d p
0.3569 0.4562 0.4186 0.2644
0.8120 0.3669 0.6954 0.5714
0.7420
FM-Wire (maps: 4)
x y d p
0.4986 0.0000 0.5203 0.2250
0.6935 0.7123 0.2775 0.4844
0.8469 0.7898 0.4189 0.7506
0.1157
FM-Wire (maps: 4)
x y d p
0.6965 0.7877 0.5251 0.2764
0.8864 0.0041 0.7075 0.6026
0.9509 0.4341 0.4624 0.7140
0.1175
FM-Wire (maps: 4)
x y d p
0.4388 0.1077 0.4010 0.2032
0.6957 0.4323 0.5243 0.3931
0.8925 0.4120 0.0473 0.6941
0.3588
Extension 1 (maps: 4)
x y d p
0.5134 0.0047 0.2993 0.4388
0.2350 0.8299 0.1086 0.6736
0.5541 0.5485 0.4299 0.9121
0.4277 0.7382 0.4497
0.8277 0.8291
0.7765 0.7752
Extension 2 (maps: 3)
x y d ω A p
0.3886 0.5850 0.3098 1.0000 0.4051 0.3509
0.3517 0.7122 0.2797 0.9647 0.8092 0.5647
0.7310 0.6591 0.1513 0.8953 0.6706
0.5692 0.1678
Extension 3 (maps: 3)
x y z r(1) r(2) θ1 θ2 p
0.4094 0.0065 0.9795 0.6462 0.5001 0.6030 0.3840 0.3736
0.5767 0.3894 0.0316 0.7908 0.6071 0.0000 0.8713 0.6527
1.0000 0.6715 0.6789 0.5417 0.8932 0.6263 0.0089
0.1515 0.5428 0.5709
0.5130

Iowa B

FM-Wire (maps: 4)
x y d p
0.0488 0.5670 0.9804 0.3652
0.5487 0.0692 0.1346 0.5255
0.6768 0.9337 0.4494 0.7139
0.2687
Extension 1 (maps: 3)
x y d p
0.6813 0.1379 0.3486 0.5818
0.6338 0.2176 0.2289 0.8564
0.8481 0.4266 0.7987
0.7021 0.5818
Extension 1 (maps: 3)
x y d p
0.0001 0.7476 0.6791 0.3615
0.5940 0.0932 0.7372 0.7156
0.7509 0.8297 0.1429
0.7961 0.1167
Extension 1 (maps: 4)
x y d p
0.3282 0.5422 0.3205 0.3588
0.4858 0.7192 0.8440 0.5259
0.9669 0.4414 0.3642 0.7991
0.7317 0.9149 0.2924
0.9514 0.0000
0.3663 0.1624
Extension 1 (maps: 4)
x y d p
0.3006 0.8274 0.7851 0.4763
0.6647 0.8115 0.3732 0.8040
0.9102 0.1029 0.3750 0.9496
0.6626 0.3256 0.1430
0.8874 0.9305
0.1671 0.4528
Extension 2 (maps: 3)
x y d ω A p
0.3394 0.6668 0.6447 0.1075 0.0054 0.1754
0.0000 0.8390 0.3442 0.2192 0.5009 0.5987
0.7382 0.3361 0.2228 0.4745 0.7035
0.7594 0.3396
Extension 3 (maps: 3)
x y z r(1) r(2) θ1 θ2 p
0.5532 0.9756 0.0000 0.7227 0.7036 0.4695 0.0715 0.3653
0.0000 0.4301 0.0000 0.8181 0.6019 0.0090 0.3781 0.6035
0.7758 0.9549 0.6095 0.8125 0.9144 0.0900 0.0632
0.5912 0.0014 0.5173
0.6902

Iowa C

FM-Wire (maps: 4)
x y d p
0.6606 0.2674 0.6476 0.3899
0.8850 0.4225 0.3292 0.6065
0.9330 0.6413 0.0911 0.7241
0.7079
Extension 1 (maps: 3)
x y d p
0.4241 0.5149 0.6235 0.3347
0.4071 0.7321 0.5128 0.5913
0.9478 0.5679 0.5761 0.6814
0.7434
Extension 2 (maps: 2)
x y d ω A p
0.5226 0.8105 0.5786 0.5759 0.1242 0.5944
0.4871 0.5164 0.1663 0.5299 0.1890
Extension 2 (maps: 3)
x y d ω A p
0.4364 1.0000 0.7984 0.8822 0.5941 0.2219
0.2226 0.3974 0.5925 0.2639 0.0190 0.5348
0.5585 0.2710 0.1503 1.0000 0.3383
0.6633 0.3692
Extension 2 (maps: 3)
x y d ω A p
0.0283 0.8144 0.2486 0.5081 0.3864 0.1004
0.0147 0.4412 0.6202 0.0922 0.1642 0.5972
0.7944 0.2242 0.7329 0.8240 0.4581
0.6054 0.6405
Extension 2 (maps: 3)
x y d ω A p
0.6723 0.5581 0.5076 0.4175 0.4342 0.2837
0.2041 0.7414 0.2873 0.1399 0.4530 0.5785
0.5560 0.6145 0.0574 0.1897 0.7604
0.2401 0.8743
Extension 3 (maps: 2)
x y z r(1) r(2) θ1 θ2 p
0.0146 0.8228 0.5411 0.8443 0.8441 0.3302 0.4247 0.7410
0.8758 0.8144 0.3237 0.9825 0.7542 0.4336 0.2307
0.4069

Iowa D

FM-Wire (maps: 4)
x y d p
0.0324 0.7167 0.2417 0.4307
0.6012 0.4157 0.2608 0.7860
0.9490 0.5543 0.3748 0.9269
0.4320
Extension 1 (maps: 4)
x y d p
0.1016 0.7296 0.1891 0.2576
0.4569 0.8795 0.3869 0.3994
0.6600 0.3016 0.3816 0.6128
0.0428 0.3692 0.3960
0.2839 0.9016
0.0671 0.0035
Extension 2 (maps: 3)
x y d ω A p
0.6784 0.4228 0.3857 0.2828 0.1740 0.2413
0.8648 0.5145 0.1162 0.2158 0.0876 0.6383
1.0000 0.8131 0.3137 0.0000 0.2479
0.4686 0.7033
Extension 3 (maps: 2)
x y z r(1) r(2) θ1 θ2 p
0.7159 0.3931 0.1837 0.7902 0.9781 0.8896 0.7006 0.5846
0.5400 0.2011 0.0000 0.5430 0.7530 0.0000 0.1992
0.9503
Extension 3 (maps: 3)
x y z r(1) r(2) θ1 θ2 p
0.9993 0.0000 1.0000 0.8468 0.6512 0.0349 0.1765 0.5541
0.1402 0.7465 0.7566 0.7049 0.8335 0.8116 0.4028 0.8043
0.4111 1.0000 0.0712 0.8388 0.5692 0.9130 0.2270
0.8701 0.5646 0.3806
0.8562
Extension 3 (maps: 3)
x y z r(1) r(2) θ1 θ2 p
0.2116 0.3123 0.8713 0.8214 0.8446 0.9937 1.0000 0.6550
0.6598 0.9999 0.9712 0.9730 0.7659 0.9916 1.0000 0.9332
0.9825 0.9460 0.3455 0.5093 0.5697 0.0000 0.8860
0.9999 0.9999 0.6566
0.5818
Extension 3 (maps: 3)
x y z r(1) r(2) θ1 θ2 p
0.3574 0.2793 0.8624 0.8898 0.5418 0.5016 0.0545 0.4673
0.3887 0.9355 0.7481 0.9816 0.8238 0.0719 0.3034 0.8218
0.8967 0.7042 0.5661 0.5605 0.7010 0.5070 0.5337
0.6378 0.3455 0.5709
0.8781

MATLAB Code

Original FM (wire.m)

  1. function dy = wire(p)
  2. % This subroutine accepts a vector that contains all the
  3. % parameters used in the wire. The number of interpolating
  4. % points is automatically calculated from the vector length.
  5. ​
  6. % Examples:
  7. % pattern_3p = wire([x y d1 d2 p]);
  8. % pattern_4p = wire([x1 x2 y1 y2 d1 d2 d3 p1 p2]);
  9. % pattern_5p = wire([x1 x2 x3 y1 y2 y3 d1 d2 d3 d4 p1 p2 p3]);
  10. ​
  11. % resolution of the pattern (# bins in histogram)
  12. nbins = 2^8;
  13. ​
  14. % duration of chaos game
  15. ndots = 2^14;
  16. ​
  17. % # of interpolating points
  18. intp = (7+length(p))/4;
  19. ​
  20. % y in [-5,5]
  21. M = 5;
  22. ​
  23. % d in [-1,1]
  24. D = 1;
  25. ​
  26. % Put the parameters into neat little vectors.
  27. xs = [0 p(1:(intp-2)) 1];
  28. ys = [0 M*(2*p((intp-1):(2*intp-4))-1) 1];
  29. d = D*(2*p((2*intp-3):(3*intp-5))-1);
  30. ps = [0 p((3*intp-4):(4*intp-7)) 1];
  31. ​
  32. % Validation required when >3 interpolating pts
  33. if intp>3
  34. % x-values must be in increasing order
  35. for vx=2:(length(xs)-2)
  36. if xs(vx) >= xs(vx+1)
  37. dy = ones(1,nbins) / nbins;
  38. return
  39. end
  40. end
  41. % p-values must be in increasing order
  42. for vp=2:(length(ps)-2)
  43. if ps(vp) >= ps(vp+1)
  44. dy = ones(1,nbins) / nbins;
  45. return
  46. end
  47. end
  48. end
  49. ​
  50. % Calculate the other parameters of the map.
  51. h = xs(end)-xs(1);
  52. a =(xs(2:end)-xs(1:end-1))/h;
  53. e = xs(1:end-1)-xs(1)*a;
  54. c =(ys(2:end)-ys(1:end-1)-d*(ys(end) - ys(1)))/h;
  55. f = ys(1:end-1)-d*ys(1)-c*xs(1);
  56. ​
  57. % Determine when to use which map during the chaos game.
  58. W = zeros(ndots,1);
  59. prob = rand(ndots,1);
  60. for k=1:(length(ps)-1)
  61. W = W + and(prob<ps(k+1),prob>ps(k))*k;
  62. end
  63. kfun = sum(W,2);
  64. y = zeros(ndots,1);
  65. xold = xs(2);
  66. yold = ys(2);
  67. y(1) = yold;
  68. for i=1:ndots-1
  69. k = max(min(kfun(i),numel(a)),1);
  70. xnew = a(k)*xold + e(k);
  71. ynew = c(k)*xold + d(k)*yold + f(k);
  72. xold = xnew;
  73. yold = ynew;
  74. y(i+1) = ynew;
  75. end
  76. ​
  77. % Create histogram.
  78. dy = hist(y, nbins);
  79. ​
  80. % Normalize histogram so that the area underneath = 1
  81. dy = dy/sum(dy);
  82. ​
  83. % End.

Extension 1 (overlap.m)

  1. function dy = overlap(p)
  2. % This subroutine accepts a vector that contains all the
  3. % parameters used in the wire. The number of mappings is
  4. % automatically calculated from the vector length.
  5. ​
  6. % resolution of the pattern (# bins in histogram)
  7. nbins = 2^8;
  8. ​
  9. % duration of chaos game
  10. ndots = 2^14;
  11. ​
  12. % # of maps
  13. nobj = (5+length(p))/6;
  14. ​
  15. % y in [-5,5]
  16. M = 5;
  17. ​
  18. % d in [-1,1]
  19. D = 1;
  20. ​
  21. % Put the parameters into neat little vectors.
  22. xs = [0 p(1:(2*nobj-2)) 1];
  23. ys = [0 M*(2*p((2*nobj-1):(4*nobj-4))-1) 1];
  24. d = D*(2*p((4*nobj-3):(5*nobj-4))-1);
  25. ps = [0 p((5*nobj-3):(6*nobj-5)) 1];
  26. ​
  27. % Validation required when >2 maps
  28. if nobj>2
  29. % interval spans must be positive
  30. for vx=3:2:(length(xs)-3)
  31. if xs(vx) >= xs(vx+1)
  32. dy = ones(1,nbins) / nbins;
  33. return
  34. end
  35. end
  36. % p-values must be in increasing order
  37. for vp=2:(length(ps)-2)
  38. if ps(vp) >= ps(vp+1)
  39. dy = ones(1,nbins) / nbins;
  40. return
  41. end
  42. end
  43. end
  44. ​
  45. % Calculate the other parameters of the map.
  46. h = xs(end)-xs(1);
  47. a =(xs(2:2:end)-xs(1:2:end-1))/h;
  48. e = xs(1:2:end-1)-xs(1)*a;
  49. c =(ys(2:2:end)-ys(1:2:end-1)-d*(ys(end)-ys(1)))/h;
  50. f = ys(1:2:end-1)-d*ys(1)-c*xs(1);
  51. ​
  52. % Determine when to use which map during the chaos game.
  53. W = zeros(ndots,1);
  54. prob = rand(ndots,1);
  55. for k=1:(length(ps)-1)
  56. W = W + and(prob<ps(k+1),prob>ps(k))*k;
  57. end
  58. kfun = sum(W,2);
  59. y = zeros(ndots,1);
  60. xold = xs(2);
  61. yold = ys(2);
  62. y(1) = yold;
  63. for i=1:ndots-1
  64. k = max(min(kfun(i),numel(a)),1);
  65. xnew = a(k)*xold + e(k);
  66. ynew = c(k)*xold + d(k)*yold + f(k);
  67. xold = xnew;
  68. yold = ynew;
  69. y(i+1) = ynew;
  70. end
  71. ​
  72. % Create histogram.
  73. dy = hist(y, nbins);
  74. ​
  75. % Normalize histogram so that the area underneath = 1
  76. dy = dy/sum(dy);
  77. ​
  78. % End.

Extension 2 (nonlinear.m)

  1. function dy = nonlinear(p)
  2. % This subroutine accepts a vector that contains all the
  3. % parameters used in the wire. The number of mappings is
  4. % automatically calculated from the vector length.
  5. ​
  6. % resolution of the pattern (# bins in histogram)
  7. nbins = 2^8;
  8. ​
  9. % duration of chaos game
  10. ndots = 2^14;
  11. ​
  12. % # of maps
  13. nobj = (5+length(p))/8;
  14. ​
  15. % y in [-5,5]
  16. M = 5;
  17. ​
  18. % d in [-1,1]
  19. D = 1;
  20. ​
  21. % period in [-pi/2,pi/2]
  22. T = pi/2;
  23. ​
  24. xs = [0 p(1:(2*nobj-2)) 1];
  25. ys = [0 M*(2*p((2*nobj-1):(4*nobj-4))-1) 1];
  26. d = D*(2*p((4*nobj-3):(5*nobj-4))-1);
  27. amp = p((5*nobj-3):(6*nobj-4));
  28. per = T*(2*p((6*nobj-3):(7*nobj-4))-1);
  29. ps = [0 p((7*nobj-3):(8*nobj-5)) 1];
  30. ​
  31. % Validation required when >2 maps
  32. if nobj>2
  33. % interval spans must be positive
  34. for vx=3:2:(length(xs)-3)
  35. if xs(vx) >= xs(vx+1)
  36. dy = ones(1,nbins) / nbins;
  37. return
  38. end
  39. end
  40. % p-values must be in increasing order
  41. for vp=2:(length(ps)-2)
  42. if ps(vp) >= ps(vp+1)
  43. dy = ones(1,nbins) / nbins;
  44. return
  45. end
  46. end
  47. end
  48. ​
  49. % Calculate the other parameters of the map.
  50. hx = xs(end)-xs(1);
  51. hy = ys(end)-ys(1);
  52. a = (xs(2:2:end)-xs(1:2:end-1))/hx;
  53. e = xs(1:2:end-1) - xs(1)*a;
  54. c = (ys(2:2:end)-ys(1:2:end-1) - d*(hy) + d.*amp.*(cos(per*ys(end))-cos(per*ys(1))))/hx;
  55. f = ys(1:2:end-1) - d*ys(1) - d.*amp.*cos(per*ys(1)) - c*xs(1);
  56. ​
  57. % Determine when to use which map during the chaos game.
  58. W = zeros(ndots,1);
  59. prob = rand(ndots,1);
  60. for k=1:(length(ps)-1)
  61. W = W + and(prob<ps(k+1),prob>ps(k))*k;
  62. end
  63. kfun = sum(W,2);
  64. y = zeros(ndots,1);
  65. xold = xs(2);
  66. yold = ys(2);
  67. y(1) = yold;
  68. for i=1:ndots-1
  69. k = max(min(kfun(i),numel(a)),1);
  70. xnew = a(k)*xold + e(k);
  71. ynew = c(k)*xold + d(k)*(yold+amp(k)*cos(per(k)*yold)) + f(k);
  72. xold = xnew;
  73. yold = ynew;
  74. y(i+1) = ynew;
  75. end
  76. ​
  77. % Create histogram.
  78. dy = hist(y, nbins);
  79. ​
  80. % Normalize histogram so that the area underneath = 1
  81. dy = dy/sum(dy);
  82. ​
  83. % End.

Extension 3 (marginal.m)

  1. function res = marginal(p,map_dim)
  2. % This subroutine accepts a vector of length 16 (2 maps)
  3. % or 27 (three maps) as the parameters, plus the dimension/s
  4. % of interest (y, z, or yz; default: y).
  5. ​
  6. % resolution of the pattern (# bins in histogram)
  7. nbins = 2^8;
  8. ​
  9. % resolution of the 2D pattern (if map_dim set to 'yz')
  10. nbins2 = 2^5 -1;
  11. ​
  12. % duration of chaos game
  13. ndots = 2^14;
  14. ​
  15. % y in [-5,5]
  16. M = 5;
  17. ​
  18. % theta in [-pi/2,pi/2]
  19. T = pi/2;
  20. ​
  21. % Put the parameters into neat little vectors.
  22. switch length(p)
  23. case 16 % 2 maps
  24. xs = [0 p(1:2) 1];
  25. ys = [0 M*(2*p(3:4)-1) 1];
  26. zs = [0 M*(2*p(5:7)-1)];
  27. r1 = 2*p(8:9)-1;
  28. r2 = 2*p(10:11)-1;
  29. t1 = T*(2*p(12:13)-1);
  30. t2 = T*(2*p(14:15)-1);
  31. ps = [0 p(16) 1];
  32. case 27 % 3 maps
  33. if ~all([p(2)<p(3) p(26)<p(27)])
  34. res = ones(1,nbins) / nbins;
  35. return
  36. end
  37. xs = [0 p(1:4) 1];
  38. ys = [0 M*(2*p(5:8)-1) 1];
  39. zs = [0 M*(2*p(9:13)-1)];
  40. r1 = 2*p(14:16)-1;
  41. r2 = 2*p(17:19)-1;
  42. t1 = T*(2*p(20:22)-1);
  43. t2 = T*(2*p(23:25)-1);
  44. ps = [0 p(26:27) 1];
  45. otherwise
  46. res = ones(1,nbins) / nbins;
  47. return
  48. end
  49. ​
  50. % Calculate the other parameters of the map.
  51. e = r1.*cos(t1);
  52. f = -r2.*sin(t2);
  53. h = r1.*sin(t1);
  54. k = r2.*cos(t2);
  55. hx = xs(end)-xs(1);
  56. hy = ys(end)-ys(1);
  57. hz = zs(end)-zs(1);
  58. a = (xs(2:2:end)-xs(1:2:end-1)) / hx;
  59. l = xs(1:2:end-1) - xs(1)*a;
  60. d = (ys(2:2:end)-ys(1:2:end-1) - e*hy - f*hz) / hx;
  61. g = (zs(2:2:end)-zs(1:2:end-1) - h*hy - k*hz) / hx;
  62. m = (ys(1:2:end-1) - d*xs(1) - e*ys(1) - f*zs(1));
  63. n = (zs(1:2:end-1) - g*xs(1) - h*ys(1) - k*zs(1));
  64. ​
  65. % Determine when to use which map during the chaos game.
  66. kfun = zeros(ndots,1);
  67. prob = rand(ndots,1);
  68. for kind=1:(length(ps)-1)
  69. kfun = kfun + and(prob<ps(kind+1),prob>ps(kind)) * kind;
  70. end
  71. y = zeros(ndots,1);
  72. z = zeros(ndots,1);
  73. yz = zeros(nbins2+1,nbins2+1);
  74. xold = xs(2);
  75. yold = ys(2);
  76. zold = zs(2);
  77. y(1) = yold;
  78. z(1) = zold;
  79. for i=1:ndots-1
  80. kk = max(min(kfun(i),numel(a)),1);
  81. xnew = a(kk)*xold + l(kk);
  82. ynew = d(kk)*xold + e(kk)*yold + f(kk)*zold + m(kk);
  83. znew = g(kk)*xold + h(kk)*yold + k(kk)*zold + n(kk);
  84. xold = xnew;
  85. yold = ynew;
  86. zold = znew;
  87. y(i+1) = ynew;
  88. z(i+1) = znew;
  89. end
  90. ​
  91. % Decide on what to send back.
  92. switch map_dim
  93. case 'y' % Marginal over y.
  94. dy = hist(y, nbins);
  95. dy = dy/sum(dy);
  96. res = dy;
  97. case 'z' % Marginal over z.
  98. dz = hist(z, nbins);
  99. dz = dz/sum(dz);
  100. res = dz;
  101. case 'yz' % 2D histogram over yz.
  102. y1 = min(y)-eps; y2 = max(y)+eps;
  103. z1 = min(z)-eps; z2 = max(z)+eps;
  104. for i=1:ndots
  105. iy = floor((y(i) - y1)/(y2-y1)*nbins2)+1;
  106. iz = floor((z(i) - z1)/(z2-z1)*nbins2)+1;
  107. yz(iy,iz) = yz(iy,iz) + 1;
  108. end
  109. dyz = yz/sum(yz(:));
  110. res = dyz;
  111. otherwise % Marginal over y by default.
  112. dy = hist(y, nbins);
  113. dy = dy/sum(dy);
  114. res = dy;
  115. end
  116. ​
  117. % End.

Online Demo

Below are live demonstrations of the wire subroutine, using normalized parameters. The other procedures have been omitted as matrix manipulation is not very efficient in JavaScript.

Wire with 3 maps

X-values
Y-values
Scalings
Probability

Wire with 4 maps

X-values
Y-values
Scalings
Probability