Deterministic simulation of mildly intermittent hydrologic recordsOnline Appendix

Mahesh L. Maskey, Carlos E. Puente*, Bellie Sivakumar

This page contains the following information: (a) MATLAB source codes for the original Fractal-Muiltifractal method and an overlap extension, and (b) interactive demos that allow changing FM patterns to explore other FM sets using as a basis Figures 1 and 2 in the paper. This information fully empowers the interested reader in computing FM sets based on the notion of vertical thresholds.

MATLAB Code

Original FM

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

Extension to overlap

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

Online Demo

Wire with 2 maps

Parameters below may be modified to draw FM outcomes dys. When the parameter Smoothness equals 1 the raw Fractal-Multifractal shadow dy is found.

X-values
Y-values
Scalings
Proportion
Smoothness

Overlap with 2 maps

Parameters below may be modified to draw FM outcomes dys. When the parameter Smoothness equals 1 the raw Fractal-Multifractal shadow dy is found.

X-values
Y-values
Scalings
Proportion
Smoothness