Cesium-Examples/examples/cesiumEx/kriging.js

/*!
 * author: [object Object] 
 * @sakitam-gis/kriging v0.1.0
 * build-time: 2019-7-6 20:41
 * LICENSE: MIT
 * (c) 2019-2019 https://github.com/sakitam-gis/kriging.js
 */
(function (global, factory) {
  typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports) :
  typeof define === 'function' && define.amd ? define(['exports'], factory) :
  (global = global || self, factory(global.kriging = {}));
}(this, function (exports) { 'use strict';

  function max(source) {
      return Math.max.apply(null, source);
  }
  function min(source) {
      return Math.min.apply(null, source);
  }
  function rep(source, n) {
      var array = [];
      for (var i = 0; i < n; i++) {
          array.push(source);
      }
      return array;
  }
  function pip(source, x, y) {
      var i = 0;
      var j = source.length - 1;
      var c = false;
      var length = source.length;
      for (; i < length; j = i++) {
          if (((source[i][1] > y) !== (source[j][1] > y))
              && (x < (source[j][0] - source[i][0]) * (y - source[i][1]) / (source[j][1] - source[i][1]) + source[i][0])) {
              c = !c;
          }
      }
      return c;
  }
  function matrixDiag(c, n) {
      var i = 0;
      var Z = rep(0, n * n);
      for (; i < n; i++) {
          Z[i * n + i] = c;
      }
      return Z;
  }
  function matrixTranspose(X, n, m) {
      var i = 0;
      var j;
      var Z = Array(m * n);
      for (; i < n; i++) {
          j = 0;
          for (; j < m; j++) {
              Z[j * n + i] = X[i * m + j];
          }
      }
      return Z;
  }
  function matrixAdd(X, Y, n, m) {
      var i = 0;
      var j;
      var Z = Array(n * m);
      for (; i < n; i++) {
          j = 0;
          for (; j < m; j++) {
              Z[i * m + j] = X[i * m + j] + Y[i * m + j];
          }
      }
      return Z;
  }
  function matrixMultiply(X, Y, n, m, p) {
      var i = 0;
      var j;
      var k;
      var Z = Array(n * p);
      for (; i < n; i++) {
          j = 0;
          for (; j < p; j++) {
              Z[i * p + j] = 0;
              k = 0;
              for (; k < m; k++) {
                  Z[i * p + j] += X[i * m + k] * Y[k * p + j];
              }
          }
      }
      return Z;
  }
  function matrixChol(X, n) {
      var i;
      var j;
      var k;
      var p = Array(n);
      for (i = 0; i < n; i++)
          p[i] = X[i * n + i];
      for (i = 0; i < n; i++) {
          for (j = 0; j < i; j++)
              p[i] -= X[i * n + j] * X[i * n + j];
          if (p[i] <= 0)
              return false;
          p[i] = Math.sqrt(p[i]);
          for (j = i + 1; j < n; j++) {
              for (k = 0; k < i; k++)
                  X[j * n + i] -= X[j * n + k] * X[i * n + k];
              X[j * n + i] /= p[i];
          }
      }
      for (i = 0; i < n; i++)
          X[i * n + i] = p[i];
      return true;
  }
  function matrixChol2inv(X, n) {
      var i;
      var j;
      var k;
      var sum;
      for (i = 0; i < n; i++) {
          X[i * n + i] = 1 / X[i * n + i];
          for (j = i + 1; j < n; j++) {
              sum = 0;
              for (k = i; k < j; k++)
                  sum -= X[j * n + k] * X[k * n + i];
              X[j * n + i] = sum / X[j * n + j];
          }
      }
      for (i = 0; i < n; i++)
          for (j = i + 1; j < n; j++)
              X[i * n + j] = 0;
      for (i = 0; i < n; i++) {
          X[i * n + i] *= X[i * n + i];
          for (k = i + 1; k < n; k++)
              X[i * n + i] += X[k * n + i] * X[k * n + i];
          for (j = i + 1; j < n; j++)
              for (k = j; k < n; k++)
                  X[i * n + j] += X[k * n + i] * X[k * n + j];
      }
      for (i = 0; i < n; i++)
          for (j = 0; j < i; j++)
              X[i * n + j] = X[j * n + i];
  }
  function matrixSolve(X, n) {
      var m = n;
      var b = Array(n * n);
      var indxc = Array(n);
      var indxr = Array(n);
      var ipiv = Array(n);
      var i;
      var icol = 0;
      var irow = 0;
      var j;
      var k;
      var l;
      var ll;
      var big;
      var dum;
      var pivinv;
      var temp;
      for (i = 0; i < n; i++) {
          for (j = 0; j < n; j++) {
              if (i === j)
                  b[i * n + j] = 1;
              else
                  b[i * n + j] = 0;
          }
      }
      for (j = 0; j < n; j++)
          ipiv[j] = 0;
      for (i = 0; i < n; i++) {
          big = 0;
          for (j = 0; j < n; j++) {
              if (ipiv[j] !== 1) {
                  for (k = 0; k < n; k++) {
                      if (ipiv[k] === 0) {
                          if (Math.abs(X[j * n + k]) >= big) {
                              big = Math.abs(X[j * n + k]);
                              irow = j;
                              icol = k;
                          }
                      }
                  }
              }
          }
          ++(ipiv[icol]);
          if (irow !== icol) {
              for (l = 0; l < n; l++) {
                  temp = X[irow * n + l];
                  X[irow * n + l] = X[icol * n + l];
                  X[icol * n + l] = temp;
              }
              for (l = 0; l < m; l++) {
                  temp = b[irow * n + l];
                  b[irow * n + l] = b[icol * n + l];
                  b[icol * n + l] = temp;
              }
          }
          indxr[i] = irow;
          indxc[i] = icol;
          if (X[icol * n + icol] === 0)
              return false;
          pivinv = 1 / X[icol * n + icol];
          X[icol * n + icol] = 1;
          for (l = 0; l < n; l++)
              X[icol * n + l] *= pivinv;
          for (l = 0; l < m; l++)
              b[icol * n + l] *= pivinv;
          for (ll = 0; ll < n; ll++) {
              if (ll !== icol) {
                  dum = X[ll * n + icol];
                  X[ll * n + icol] = 0;
                  for (l = 0; l < n; l++)
                      X[ll * n + l] -= X[icol * n + l] * dum;
                  for (l = 0; l < m; l++)
                      b[ll * n + l] -= b[icol * n + l] * dum;
              }
          }
      }
      for (l = (n - 1); l >= 0; l--) {
          if (indxr[l] !== indxc[l]) {
              for (k = 0; k < n; k++) {
                  temp = X[k * n + indxr[l]];
                  X[k * n + indxr[l]] = X[k * n + indxc[l]];
                  X[k * n + indxc[l]] = temp;
              }
          }
      }
      return true;
  }
  function variogramGaussian(h, nugget, range, sill, A) {
      return nugget + ((sill - nugget) / range)
          * (1.0 - Math.exp(-(1.0 / A) * Math.pow(h / range, 2)));
  }
  function variogramExponential(h, nugget, range, sill, A) {
      return nugget + ((sill - nugget) / range)
          * (1.0 - Math.exp(-(1.0 / A) * (h / range)));
  }
  function variogramSpherical(h, nugget, range, sill) {
      if (h > range)
          return nugget + (sill - nugget) / range;
      return nugget + ((sill - nugget) / range)
          * (1.5 * (h / range) - 0.5 * Math.pow(h / range, 3));
  }

  function train(t, x, y, model, sigma2, alpha) {
      var variogram = {
          t: t,
          x: x,
          y: y,
          nugget: 0.0,
          range: 0.0,
          sill: 0.0,
          A: 1 / 3,
          n: 0,
          model: variogramExponential,
          K: [],
          M: [],
      };
      switch (model) {
          case 'gaussian':
              variogram.model = variogramGaussian;
              break;
          case 'exponential':
              variogram.model = variogramExponential;
              break;
          case 'spherical':
              variogram.model = variogramSpherical;
              break;
          default:
              variogram.model = variogramExponential;
      }
      var i;
      var j;
      var k;
      var l;
      var n = t.length;
      var distance = Array((n * n - n) / 2);
      for (i = 0, k = 0; i < n; i++) {
          for (j = 0; j < i; j++, k++) {
              distance[k] = Array(2);
              distance[k][0] = Math.pow(Math.pow(x[i] - x[j], 2)
                  + Math.pow(y[i] - y[j], 2), 0.5);
              distance[k][1] = Math.abs(t[i] - t[j]);
          }
      }
      distance.sort(function (a, b) { return a[0] - b[0]; });
      variogram.range = distance[(n * n - n) / 2 - 1][0];
      var lags = ((n * n - n) / 2) > 30 ? 30 : (n * n - n) / 2;
      var tolerance = variogram.range / lags;
      var lag = rep(0, lags);
      var semi = rep(0, lags);
      if (lags < 30) {
          for (l = 0; l < lags; l++) {
              lag[l] = distance[l][0];
              semi[l] = distance[l][1];
          }
      }
      else {
          for (i = 0, j = 0, k = 0, l = 0; i < lags && j < ((n * n - n) / 2); i++, k = 0) {
              while (distance[j][0] <= ((i + 1) * tolerance)) {
                  lag[l] += distance[j][0];
                  semi[l] += distance[j][1];
                  j++;
                  k++;
                  if (j >= ((n * n - n) / 2))
                      break;
              }
              if (k > 0) {
                  lag[l] /= k;
                  semi[l] /= k;
                  l++;
              }
          }
          if (l < 2)
              return variogram;
      }
      n = l;
      variogram.range = lag[n - 1] - lag[0];
      var X = rep(1, 2 * n);
      var Y = Array(n);
      var A = variogram.A;
      for (i = 0; i < n; i++) {
          switch (model) {
              case 'gaussian':
                  X[i * 2 + 1] = 1.0 - Math.exp(-(1.0 / A) * Math.pow(lag[i] / variogram.range, 2));
                  break;
              case 'exponential':
                  X[i * 2 + 1] = 1.0 - Math.exp(-(1.0 / A) * lag[i] / variogram.range);
                  break;
              case 'spherical':
                  X[i * 2 + 1] = 1.5 * (lag[i] / variogram.range)
                      - 0.5 * Math.pow(lag[i] / variogram.range, 3);
                  break;
          }
          Y[i] = semi[i];
      }
      var Xt = matrixTranspose(X, n, 2);
      var Z = matrixMultiply(Xt, X, 2, n, 2);
      Z = matrixAdd(Z, matrixDiag(1 / alpha, 2), 2, 2);
      var cloneZ = Z.slice(0);
      if (matrixChol(Z, 2)) {
          matrixChol2inv(Z, 2);
      }
      else {
          matrixSolve(cloneZ, 2);
          Z = cloneZ;
      }
      var W = matrixMultiply(matrixMultiply(Z, Xt, 2, 2, n), Y, 2, n, 1);
      variogram.nugget = W[0];
      variogram.sill = W[1] * variogram.range + variogram.nugget;
      variogram.n = x.length;
      n = x.length;
      var K = Array(n * n);
      for (i = 0; i < n; i++) {
          for (j = 0; j < i; j++) {
              K[i * n + j] = variogram.model(Math.pow(Math.pow(x[i] - x[j], 2)
                  + Math.pow(y[i] - y[j], 2), 0.5), variogram.nugget, variogram.range, variogram.sill, variogram.A);
              K[j * n + i] = K[i * n + j];
          }
          K[i * n + i] = variogram.model(0, variogram.nugget, variogram.range, variogram.sill, variogram.A);
      }
      var C = matrixAdd(K, matrixDiag(sigma2, n), n, n);
      var cloneC = C.slice(0);
      if (matrixChol(C, n)) {
          matrixChol2inv(C, n);
      }
      else {
          matrixSolve(cloneC, n);
          C = cloneC;
      }
      var K1 = C.slice(0);
      var M = matrixMultiply(C, t, n, n, 1);
      variogram.K = K1;
      variogram.M = M;
      return variogram;
  }
  function predict(x, y, variogram) {
      var i;
      var k = Array(variogram.n);
      for (i = 0; i < variogram.n; i++) {
          k[i] = variogram.model(Math.pow(Math.pow(x - variogram.x[i], 2)
              + Math.pow(y - variogram.y[i], 2), 0.5), variogram.nugget, variogram.range, variogram.sill, variogram.A);
      }
      return matrixMultiply(k, variogram.M, 1, variogram.n, 1)[0];
  }
  function variance(x, y, variogram) {
      var i;
      var k = Array(variogram.n);
      for (i = 0; i < variogram.n; i++) {
          k[i] = variogram.model(Math.pow(Math.pow(x - variogram.x[i], 2) + Math.pow(y - variogram.y[i], 2), 0.5), variogram.nugget, variogram.range, variogram.sill, variogram.A);
      }
      var val = matrixMultiply(matrixMultiply(k, variogram.K, 1, variogram.n, variogram.n), k, 1, variogram.n, 1)[0];
      return variogram.model(0, variogram.nugget, variogram.range, variogram.sill, variogram.A) + val;
  }
  function grid(polygons, variogram, width) {
      var i;
      var j;
      var k;
      var n = polygons.length;
      if (n === 0)
          return;
      var xlim = [polygons[0][0][0], polygons[0][0][0]];
      var ylim = [polygons[0][0][1], polygons[0][0][1]];
      for (i = 0; i < n; i++) {
          for (j = 0; j < polygons[i].length; j++) {
              if (polygons[i][j][0] < xlim[0])
                  xlim[0] = polygons[i][j][0];
              if (polygons[i][j][0] > xlim[1])
                  xlim[1] = polygons[i][j][0];
              if (polygons[i][j][1] < ylim[0])
                  ylim[0] = polygons[i][j][1];
              if (polygons[i][j][1] > ylim[1])
                  ylim[1] = polygons[i][j][1];
          }
      }
      var xtarget;
      var ytarget;
      var a = Array(2);
      var b = Array(2);
      var lxlim = Array(2);
      var lylim = Array(2);
      var x = Math.ceil((xlim[1] - xlim[0]) / width);
      var y = Math.ceil((ylim[1] - ylim[0]) / width);
      var A = Array(x + 1);
      for (i = 0; i <= x; i++)
          A[i] = Array(y + 1);
      for (i = 0; i < n; i++) {
          lxlim[0] = polygons[i][0][0];
          lxlim[1] = lxlim[0];
          lylim[0] = polygons[i][0][1];
          lylim[1] = lylim[0];
          for (j = 1; j < polygons[i].length; j++) {
              if (polygons[i][j][0] < lxlim[0])
                  lxlim[0] = polygons[i][j][0];
              if (polygons[i][j][0] > lxlim[1])
                  lxlim[1] = polygons[i][j][0];
              if (polygons[i][j][1] < lylim[0])
                  lylim[0] = polygons[i][j][1];
              if (polygons[i][j][1] > lylim[1])
                  lylim[1] = polygons[i][j][1];
          }
          a[0] = Math.floor(((lxlim[0] - ((lxlim[0] - xlim[0]) % width)) - xlim[0]) / width);
          a[1] = Math.ceil(((lxlim[1] - ((lxlim[1] - xlim[1]) % width)) - xlim[0]) / width);
          b[0] = Math.floor(((lylim[0] - ((lylim[0] - ylim[0]) % width)) - ylim[0]) / width);
          b[1] = Math.ceil(((lylim[1] - ((lylim[1] - ylim[1]) % width)) - ylim[0]) / width);
          for (j = a[0]; j <= a[1]; j++) {
              for (k = b[0]; k <= b[1]; k++) {
                  xtarget = xlim[0] + j * width;
                  ytarget = ylim[0] + k * width;
                  if (pip(polygons[i], xtarget, ytarget)) {
                      A[j][k] = predict(xtarget, ytarget, variogram);
                  }
              }
          }
      }
      return {
          xlim: xlim,
          ylim: ylim,
          width: width,
          data: A,
          zlim: [min(variogram.t), max(variogram.t)],
      };
  }
  function plot(canvas, grid, xlim, ylim, colors) {
      var ctx = canvas.getContext('2d');
      var data = grid.data, zlim = grid.zlim, width = grid.width;
      if (ctx) {
          ctx.clearRect(0, 0, canvas.width, canvas.height);
          var range = [xlim[1] - xlim[0], ylim[1] - ylim[0], zlim[1] - zlim[0]];
          var i = void 0;
          var j = void 0;
          var x = void 0;
          var y = void 0;
          var z = void 0;
          var n = data.length;
          var m = data[0].length;
          var wx = Math.ceil(width * canvas.width / (xlim[1] - xlim[0]));
          var wy = Math.ceil(width * canvas.height / (ylim[1] - ylim[0]));
          for (i = 0; i < n; i++) {
              for (j = 0; j < m; j++) {
                  if (data[i][j] === undefined)
                      continue;
                  x = canvas.width * (i * width + grid.xlim[0] - xlim[0]) / range[0];
                  y = canvas.height * (1 - (j * width + grid.ylim[0] - ylim[0]) / range[1]);
                  z = (data[i][j] - zlim[0]) / range[2];
                  if (z < 0.0)
                      z = 0.0;
                  if (z > 1.0)
                      z = 1.0;
                  ctx.fillStyle = colors[Math.floor((colors.length - 1) * z)];
                  ctx.fillRect(Math.round(x - wx / 2), Math.round(y - wy / 2), wx, wy);
              }
          }
      }
  }
  var index = {
      train: train,
      predict: predict,
      variance: variance,
      grid: grid,
      plot: plot,
      max: max,
      min: min,
      pip: pip,
      rep: rep,
      matrixDiag: matrixDiag,
      matrixTranspose: matrixTranspose,
      matrixAdd: matrixAdd,
      matrixMultiply: matrixMultiply,
      matrixChol: matrixChol,
      matrixChol2inv: matrixChol2inv,
      matrixSolve: matrixSolve,
      variogramGaussian: variogramGaussian,
      variogramExponential: variogramExponential,
      variogramSpherical: variogramSpherical,
  };

  exports.default = index;
  exports.grid = grid;
  exports.matrixAdd = matrixAdd;
  exports.matrixChol = matrixChol;
  exports.matrixChol2inv = matrixChol2inv;
  exports.matrixDiag = matrixDiag;
  exports.matrixMultiply = matrixMultiply;
  exports.matrixSolve = matrixSolve;
  exports.matrixTranspose = matrixTranspose;
  exports.max = max;
  exports.min = min;
  exports.pip = pip;
  exports.plot = plot;
  exports.predict = predict;
  exports.rep = rep;
  exports.train = train;
  exports.variance = variance;
  exports.variogramExponential = variogramExponential;
  exports.variogramGaussian = variogramGaussian;
  exports.variogramSpherical = variogramSpherical;

  Object.defineProperty(exports, '__esModule', { value: true });

}));
//# sourceMappingURL=kriging.js.map