diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml new file mode 100644 index 0000000..3397381 --- /dev/null +++ b/.github/workflows/main.yml @@ -0,0 +1,132 @@ +name: Continuous Integration + +on: + push: + branches: + - "master" + - "develop" + tags: + - "*" + pull_request: + branches: + - "develop" + # Allows you to run this workflow manually from the Actions tab + workflow_dispatch: + +jobs: + format_check: + name: format check + runs-on: ubuntu-latest + strategy: + fail-fast: false + + steps: + - uses: actions/checkout@v2 + + - name: Set up Python 3.8 + uses: actions\setup-python@v2 + with: + python-version: 3.8 + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install black + + - name: black check + run: | + python -m black --check pykrige/ examples/ tests/ + + build_sdist: + name: sdist and coveralls + runs-on: ubuntu-latest + strategy: + fail-fast: false + + steps: + - uses: actions/checkout@v2 + with: + fetch-depth: '0' + + - name: Set up Python 3.8 + uses: actions\setup-python@v2 + with: + python-version: 3.8 + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -r requirements_setup.txt + pip install -r requirements.txt + pip install -r requirements_test.txt + pip install coveralls>=3.0.0 + + - name: Build sdist + run: | + python setup.py sdist -d dist + python setup.py build_ext --inplace + + - name: Run tests + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: | + python -m pytest --cov pykrige --cov-report term-missing -v tests/ + python -m coveralls --service=github + + - uses: actions/upload-artifact@v2 + with: + path: dist/*.tar.gz + + build_wheels: + name: wheels on ${{matrix.os}} + runs-on: ${{matrix.os}} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, windows-latest, macos-latest] + + steps: + - uses: actions/checkout@v2 + with: + fetch-depth: '0' + + - name: Build wheels + uses: joerick/cibuildwheel@v1.10.0 + env: + CIBW_BUILD: cp36-* cp37-* cp38-* cp39-* + CIBW_BEFORE_BUILD: pip install numpy==1.19.* scipy==1.5.* cython==0.29.* setuptools + CIBW_TEST_REQUIRES: pytest scikit-learn + CIBW_TEST_COMMAND: pytest -v {project}/tests + with: + output-dir: dist + + - uses: actions/upload-artifact@v2 + with: + path: ./dist/*.whl + + upload_to_pypi: + needs: [build_wheels, build_sdist] + runs-on: ubuntu-latest + + steps: + - uses: actions/download-artifact@v2 + with: + name: artifact + path: dist + + - name: Publish to Test PyPI + if: github.ref == 'refs/heads/master' || github.ref == 'refs/heads/develop' + uses: pypa/gh-action-pypi-publish@release/v1 + with: + user: __token__ + password: ${{ secrets.test_pypi_password }} + repository_url: https://test.pypi.org/legacy/ + skip_existing: true + + - name: Publish to PyPI + # only if tagged + if: startsWith(github.ref, 'refs/tags') + uses: pypa/gh-action-pypi-publish@release/v1 + with: + user: __token__ + password: ${{ secrets.pypi_password }} diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 4e41a65..0000000 --- a/.travis.yml +++ /dev/null @@ -1,111 +0,0 @@ -# do not build pull request from base repo twice (only build PRs of forks) -if: type != pull_request OR fork = true - -language: python -python: 3.8 - -# setuptools-scm needs all tags in order to obtain a proper version -git: - depth: false - -env: - global: - - TWINE_USERNAME=geostatframework - - CIBW_BEFORE_BUILD="pip install numpy==1.17.3 scipy==1.3.2 cython==0.29.14 setuptools" - - CIBW_TEST_REQUIRES="pytest scikit-learn gstools" - - CIBW_TEST_COMMAND="pytest -v {project}/tests" - - OMP_NUM_THREADS=2 - -before_install: - - | - if [[ "$TRAVIS_OS_NAME" = windows ]]; then - choco install python --version 3.8.0 - export PATH="/c/Python38:/c/Python38/Scripts:$PATH" - # make sure it's on PATH as 'python3' - ln -s /c/Python38/python.exe /c/Python38/python3.exe - fi - -install: python3 -m pip install cibuildwheel==1.3.0 - -script: python3 -m cibuildwheel --output-dir dist - -after_success: - - | - if [[ $TRAVIS_PULL_REQUEST == 'false' ]]; then - python3 -m pip install twine - python3 -m twine upload --verbose --skip-existing --repository-url https://test.pypi.org/legacy/ dist/* - if [[ $TRAVIS_TAG ]]; then python3 -m twine upload --verbose --skip-existing dist/*; fi - fi - -notifications: - email: - recipients: - - info@geostat-framework.org - -jobs: - include: - - name: "black check" - services: docker - install: python3 -m pip install black - script: python3 -m black --check pykrige/ examples/ tests/ - after_success: echo "all files formatted correctly" - after_failure: echo "some files not formatted correctly" - - - name: "sdist and coverage" - services: docker - install: - - python3 -m pip install -r requirements_setup.txt - - python3 -m pip install -r requirements_test.txt - - python3 -m pip install -r requirements.txt - script: - - python3 setup.py sdist -d dist - - python3 setup.py build_ext --inplace - - python3 -m pytest --cov pykrige --cov-report term-missing -v tests/ - - python3 -m coveralls - - - name: "Linux py35" - services: docker - env: CIBW_BUILD="cp35-*" - - name: "Linux py36" - services: docker - env: CIBW_BUILD="cp36-*" - - name: "Linux py37" - services: docker - env: CIBW_BUILD="cp37-*" - - name: "Linux py38" - services: docker - env: CIBW_BUILD="cp38-*" - - - name: "MacOS py35" - os: osx - language: shell - env: CIBW_BUILD="cp35-*" - - name: "MacOS py36" - os: osx - language: shell - env: CIBW_BUILD="cp36-*" - - name: "MacOS py37" - os: osx - language: shell - env: CIBW_BUILD="cp37-*" - - name: "MacOS py38" - os: osx - language: shell - env: CIBW_BUILD="cp38-*" - - - name: "Win py35" - os: windows - language: shell - env: CIBW_BUILD="cp35-*" - - name: "Win py36" - os: windows - language: shell - env: CIBW_BUILD="cp36-*" - - name: "Win py37" - os: windows - language: shell - env: CIBW_BUILD="cp37-*" - - name: "Win py38" - os: windows - language: shell - env: CIBW_BUILD="cp38-*" diff --git a/.zenodo.json b/.zenodo.json index 4e384e2..9600bfc 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -28,6 +28,10 @@ { "type": "Other", "name": "Daniel Mej\u00eda Raigosa" + }, + { + "type": "Other", + "name": "Marcelo Albuquerque" } ], "language": "eng", @@ -37,6 +41,7 @@ "universal kriging", "external drift kriging", "regression kriging", + "classification kriging", "variogram", "geostatistics", "Python", diff --git a/CHANGELOG.md b/CHANGELOG.md index 76f7ba1..6c36d69 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,23 @@ Changelog ========= +Version 1.6.0 +------------- +*April 04, 2021* + +**New features** + +* added Classification Kriging ([#165](https://github.com/GeoStat-Framework/PyKrige/pull/165), [#184](https://github.com/GeoStat-Framework/PyKrige/pull/184)) +* added wheels for Python 3.9 ([#175](https://github.com/GeoStat-Framework/PyKrige/pull/175)) + +**Changes** + +* moved scikit-learn compat-class `Krige` to `pykrige.compat` ([#165](https://github.com/GeoStat-Framework/PyKrige/pull/165)) +* dropped Python 3.5 support ([#183](https://github.com/GeoStat-Framework/PyKrige/pull/183)) +* moved CI to GitHub-Actions ([#175](https://github.com/GeoStat-Framework/PyKrige/pull/175), [#183](https://github.com/GeoStat-Framework/PyKrige/pull/183)) +* Fixed Typo in `02_kriging3D.py` example ([#182](https://github.com/GeoStat-Framework/PyKrige/pull/182)) + + Version 1.5.1 ------------- *August 20, 2020* diff --git a/README.md b/README.md new file mode 100644 index 0000000..766617a --- /dev/null +++ b/README.md @@ -0,0 +1,118 @@ +# PyKrige + +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3738604.svg)](https://doi.org/10.5281/zenodo.3738604) +[![PyPI version](https://badge.fury.io/py/PyKrige.svg)](https://badge.fury.io/py/PyKrige) +[![Conda Version](https://img.shields.io/conda/vn/conda-forge/pykrige.svg)](https://anaconda.org/conda-forge/pykrige) +[![Build Status](https://github.com/GeoStat-Framework/PyKrige/workflows/Continuous%20Integration/badge.svg?branch=develop)](https://github.com/GeoStat-Framework/PyKrige/actions) +[![Coverage Status](https://coveralls.io/repos/github/GeoStat-Framework/PyKrige/badge.svg?branch=develop)](https://coveralls.io/github/GeoStat-Framework/PyKrige?branch=develop) +[![Documentation Status](https://readthedocs.org/projects/pykrige/badge/?version=stable)](http://pykrige.readthedocs.io/en/stable/?badge=stable) +[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black) + + +

+PyKrige-LOGO +

+ +Kriging Toolkit for Python. + +## Purpose + +The code supports 2D and 3D ordinary and universal kriging. Standard +variogram models (linear, power, spherical, gaussian, exponential) are +built in, but custom variogram models can also be used. The 2D universal +kriging code currently supports regional-linear, point-logarithmic, and +external drift terms, while the 3D universal kriging code supports a +regional-linear drift term in all three spatial dimensions. Both +universal kriging classes also support generic 'specified' and +'functional' drift capabilities. With the 'specified' drift capability, +the user may manually specify the values of the drift(s) at each data +point and all grid points. With the 'functional' drift capability, the +user may provide callable function(s) of the spatial coordinates that +define the drift(s). The package includes a module that contains +functions that should be useful in working with ASCII grid files (`\*.asc`). + +See the documentation at for more +details and examples. + +## Installation + +PyKrige requires Python 3.5+ as well as numpy, scipy. It can be +installed from PyPi with, + +``` bash +pip install pykrige +``` + +scikit-learn is an optional dependency needed for parameter tuning and +regression kriging. matplotlib is an optional dependency needed for +plotting. + +If you use conda, PyKrige can be installed from the conda-forge channel with, + +``` bash +conda install -c conda-forge pykrige +``` + +## Features + +### Kriging algorithms + +- `OrdinaryKriging`: 2D ordinary kriging with estimated mean +- `UniversalKriging`: 2D universal kriging providing drift terms +- `OrdinaryKriging3D`: 3D ordinary kriging +- `UniversalKriging3D`: 3D universal kriging +- `RegressionKriging`: An implementation of Regression-Kriging +- `ClassificationKriging`: An implementation of Simplicial Indicator + Kriging + +### Wrappers + +- `rk.Krige`: A scikit-learn wrapper class for Ordinary and Universal + Kriging + +### Tools + +- `kriging_tools.write_asc_grid`: Writes gridded data to ASCII grid + file (`\*.asc`) +- `kriging_tools.read_asc_grid`: Reads ASCII grid file (`\*.asc`) + +### Kriging Parameters Tuning + +A scikit-learn compatible API for parameter tuning by cross-validation +is exposed in +[sklearn.model\_selection.GridSearchCV](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html). +See the [Krige +CV](http://pykrige.readthedocs.io/en/latest/examples/08_krige_cv.html#sphx-glr-examples-08-krige-cv-py) +example for a more practical illustration. + +### Regression Kriging + +[Regression kriging](https://en.wikipedia.org/wiki/Regression-Kriging) +can be performed with +[pykrige.rk.RegressionKriging](http://pykrige.readthedocs.io/en/latest/examples/07_regression_kriging2d.html). +This class takes as parameters a scikit-learn regression model, and +details of either either the `OrdinaryKriging` or the `UniversalKriging` +class, and performs a correction steps on the ML regression prediction. + +A demonstration of the regression kriging is provided in the +[corresponding +example](http://pykrige.readthedocs.io/en/latest/examples/07_regression_kriging2d.html#sphx-glr-examples-07-regression-kriging2d-py). + +### Classification Kriging + +[Simplifical Indicator +kriging](https://www.sciencedirect.com/science/article/abs/pii/S1002070508600254) +can be performed with +[pykrige.rk.ClassificationKriging](http://pykrige.readthedocs.io/en/latest/examples/10_classification_kriging2d.html). +This class takes as parameters a scikit-learn classification model, and +details of either the `OrdinaryKriging` or the `UniversalKriging` class, +and performs a correction steps on the ML classification prediction. + +A demonstration of the classification kriging is provided in the +[corresponding +example](http://pykrige.readthedocs.io/en/latest/examples/10_classification_kriging2d.html#sphx-glr-examples-10-classification-kriging2d-py). + +## License + +PyKrige uses the BSD 3-Clause License. diff --git a/README.rst b/README.rst deleted file mode 100644 index 822de88..0000000 --- a/README.rst +++ /dev/null @@ -1,116 +0,0 @@ -PyKrige -======= - -.. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.3738604.svg - :target: https://doi.org/10.5281/zenodo.3738604 -.. image:: https://badge.fury.io/py/PyKrige.svg - :target: https://badge.fury.io/py/PyKrige -.. image:: https://img.shields.io/conda/vn/conda-forge/pykrige.svg - :target: https://anaconda.org/conda-forge/pykrige -.. image:: https://travis-ci.com/GeoStat-Framework/PyKrige.svg?branch=master - :target: https://travis-ci.com/GeoStat-Framework/PyKrige -.. image:: https://coveralls.io/repos/github/GeoStat-Framework/PyKrige/badge.svg?branch=master - :target: https://coveralls.io/github/GeoStat-Framework/PyKrige?branch=master -.. image:: https://readthedocs.org/projects/pykrige/badge/?version=stable - :target: http://pykrige.readthedocs.io/en/stable/?badge=stable - :alt: Documentation Status -.. image:: https://img.shields.io/badge/code%20style-black-000000.svg - :target: https://github.com/psf/black - - -.. figure:: https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/docs/source/pics/PyKrige_250.png - :align: center - :alt: PyKrige - :figclass: align-center - - -Kriging Toolkit for Python. - - -Purpose -^^^^^^^ - -The code supports 2D and 3D ordinary and universal kriging. Standard variogram models -(linear, power, spherical, gaussian, exponential) are built in, but custom variogram models can also be used. -The 2D universal kriging code currently supports regional-linear, point-logarithmic, and external drift terms, -while the 3D universal kriging code supports a regional-linear drift term in all three spatial dimensions. -Both universal kriging classes also support generic 'specified' and 'functional' drift capabilities. -With the 'specified' drift capability, the user may manually specify the values of the drift(s) at each data -point and all grid points. With the 'functional' drift capability, the user may provide callable function(s) -of the spatial coordinates that define the drift(s). The package includes a module that contains functions -that should be useful in working with ASCII grid files (`*.asc`). - -See the documentation at `http://pykrige.readthedocs.io/ `_ -for more details and examples. - - -Installation -^^^^^^^^^^^^ - -PyKrige requires Python 3.5+ as well as numpy, scipy. It can be installed from PyPi with, - -.. code:: bash - - pip install pykrige - -scikit-learn is an optional dependency needed for parameter tuning and regression kriging. -matplotlib is an optional dependency needed for plotting. - -If you use conda, PyKrige can be installed from the `conda-forge` channel with, - -.. code:: bash - - conda install -c conda-forge pykrige - - -Features -^^^^^^^^ - -Kriging algorithms ------------------- - -* ``OrdinaryKriging``: 2D ordinary kriging with estimated mean -* ``UniversalKriging``: 2D universal kriging providing drift terms -* ``OrdinaryKriging3D``: 3D ordinary kriging -* ``UniversalKriging3D``: 3D universal kriging -* ``RegressionKriging``: An implementation of Regression-Kriging - - -Wrappers --------- - -* ``rk.Krige``: A scikit-learn wrapper class for Ordinary and Universal Kriging - - -Tools ------ - -* ``kriging_tools.write_asc_grid``: Writes gridded data to ASCII grid file (\*.asc) -* ``kriging_tools.read_asc_grid``: Reads ASCII grid file (\*.asc) - - -Kriging Parameters Tuning -------------------------- - -A scikit-learn compatible API for parameter tuning by cross-validation is exposed in -`sklearn.model_selection.GridSearchCV `_. -See the `Krige CV `_ -example for a more practical illustration. - - -Regression Kriging ------------------- - -`Regression kriging `_ can be performed -with `pykrige.rk.RegressionKriging `_. -This class takes as parameters a scikit-learn regression model, and details of either either -the ``OrdinaryKriging`` or the ``UniversalKriging`` class, and performs a correction steps on the ML regression prediction. - -A demonstration of the regression kriging is provided in the -`corresponding example `_. - - -License -^^^^^^^ - -PyKrige uses the BSD 3-Clause License. diff --git a/benchmarks/kriging_benchmarks.py b/benchmarks/kriging_benchmarks.py index 9f03688..927c333 100644 --- a/benchmarks/kriging_benchmarks.py +++ b/benchmarks/kriging_benchmarks.py @@ -3,11 +3,11 @@ from time import time import numpy as np from pykrige.ok import OrdinaryKriging + np.random.seed(19999) -VARIOGRAM_MODELS = ['power', 'gaussian', 'spherical', - 'exponential', 'linear'] -BACKENDS = ['vectorized', 'loop', 'C'] +VARIOGRAM_MODELS = ["power", "gaussian", "spherical", "exponential", "linear"] +BACKENDS = ["vectorized", "loop", "C"] N_MOVING_WINDOW = [None, 10, 50, 100] @@ -36,23 +36,32 @@ def make_benchark(n_train, n_test, n_dim=2): for variogram_model in VARIOGRAM_MODELS: tic = time() - OK = OrdinaryKriging(X_train[:, 0], X_train[:, 1], y_train, - variogram_model='linear', - verbose=False, enable_plotting=False) - res['t_train_{}'.format(variogram_model)] = time() - tic + OK = OrdinaryKriging( + X_train[:, 0], + X_train[:, 1], + y_train, + variogram_model="linear", + verbose=False, + enable_plotting=False, + ) + res["t_train_{}".format(variogram_model)] = time() - tic # All the following tests are performed with the linear variogram model for backend in BACKENDS: for n_closest_points in N_MOVING_WINDOW: - if backend == 'vectorized' and n_closest_points is not None: + if backend == "vectorized" and n_closest_points is not None: continue # this is not supported tic = time() - OK.execute('points', X_test[:, 0], X_test[:, 1], - backend=backend, - n_closest_points=n_closest_points) - res['t_test_{}_{}'.format(backend, n_closest_points)] = time() - tic + OK.execute( + "points", + X_test[:, 0], + X_test[:, 1], + backend=backend, + n_closest_points=n_closest_points, + ) + res["t_test_{}_{}".format(backend, n_closest_points)] = time() - tic return res @@ -71,34 +80,40 @@ def print_benchmark(n_train, n_test, n_dim, res): res : dict a dictionary with the timing results """ - print('='*80) - print(' '*10, 'N_dim={}, N_train={}, N_test={}'.format(n_dim, - n_train, n_test)) - print('='*80) - print('\n', '# Training the model', '\n') - print('|'.join(['{:>11} '.format(el) for el in ['t_train (s)'] + - VARIOGRAM_MODELS])) - print('-' * (11 + 2) * (len(VARIOGRAM_MODELS) + 1)) - print('|'.join(['{:>11} '.format('Training')] + - ['{:>11.2} '.format(el) for el in - [res['t_train_{}'.format(mod)] - for mod in VARIOGRAM_MODELS]])) - - print('\n', '# Predicting kriging points', '\n') - print('|'.join(['{:>11} '.format(el) for el in ['t_test (s)'] + BACKENDS])) - print('-' * (11 + 2) * (len(BACKENDS) + 1)) + print("=" * 80) + print(" " * 10, "N_dim={}, N_train={}, N_test={}".format(n_dim, n_train, n_test)) + print("=" * 80) + print("\n", "# Training the model", "\n") + print("|".join(["{:>11} ".format(el) for el in ["t_train (s)"] + VARIOGRAM_MODELS])) + print("-" * (11 + 2) * (len(VARIOGRAM_MODELS) + 1)) + print( + "|".join( + ["{:>11} ".format("Training")] + + [ + "{:>11.2} ".format(el) + for el in [res["t_train_{}".format(mod)] for mod in VARIOGRAM_MODELS] + ] + ) + ) + + print("\n", "# Predicting kriging points", "\n") + print("|".join(["{:>11} ".format(el) for el in ["t_test (s)"] + BACKENDS])) + print("-" * (11 + 2) * (len(BACKENDS) + 1)) for n_closest_points in N_MOVING_WINDOW: - timing_results = [res.get( - 't_test_{}_{}'.format(mod, n_closest_points), '') - for mod in BACKENDS] - print('|'.join(['{:>11} '.format('N_nn=' + str(n_closest_points))] + - ['{:>11.2} '.format(el) for el in timing_results])) - - -if __name__ == '__main__': - for no_train, no_test in [(400, 1000), - (400, 2000), - (800, 2000)]: + timing_results = [ + res.get("t_test_{}_{}".format(mod, n_closest_points), "") + for mod in BACKENDS + ] + print( + "|".join( + ["{:>11} ".format("N_nn=" + str(n_closest_points))] + + ["{:>11.2} ".format(el) for el in timing_results] + ) + ) + + +if __name__ == "__main__": + for no_train, no_test in [(400, 1000), (400, 2000), (800, 2000)]: results = make_benchark(no_train, no_test) print_benchmark(no_train, no_test, 2, results) diff --git a/docs/requirements_doc.txt b/docs/requirements_doc.txt index 34ee608..2c4ec45 100644 --- a/docs/requirements_doc.txt +++ b/docs/requirements_doc.txt @@ -2,7 +2,7 @@ scikit-learn>=0.19 gstools>=1.1.1 sphinx pillow -recommonmark sphinx_rtd_theme sphinxcontrib-napoleon sphinx-gallery +m2r2 \ No newline at end of file diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 4bb9e62..558bef5 100755 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -1,2 +1,2 @@ -.. include:: ../../CHANGELOG.md +.. mdinclude:: ../../CHANGELOG.md diff --git a/docs/source/conf.py b/docs/source/conf.py index b36d737..9ddd224 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -48,7 +48,7 @@ "sphinxcontrib.napoleon", "sphinx_gallery.gen_gallery", "sphinx.ext.linkcode", - "recommonmark", + "m2r2", ] @@ -244,6 +244,8 @@ # -- Options for LaTeX output --------------------------------------------- +latex_logo = "pics/PyKrige_150.png" + latex_elements = { # The paper size ('letterpaper' or 'a4paper'). #'papersize': 'letterpaper', diff --git a/docs/source/index.rst b/docs/source/index.rst index 4d2c74b..8a2fb51 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -1,2 +1,2 @@ -.. include:: ../../README.rst +.. mdinclude:: ../../README.md diff --git a/examples/02_kriging3D.py b/examples/02_kriging3D.py index 7cf14ea..09438db 100755 --- a/examples/02_kriging3D.py +++ b/examples/02_kriging3D.py @@ -98,7 +98,7 @@ ax2.set_title("regional lin. drift") ax3.imshow(k3d3[:, :, 0], origin="lower") ax3.set_title("specified drift") -ax4.imshow(k3d3[:, :, 0], origin="lower") +ax4.imshow(k3d4[:, :, 0], origin="lower") ax4.set_title("functional drift") plt.tight_layout() plt.show() diff --git a/examples/07_regression_kriging2d.py b/examples/07_regression_kriging2d.py index 19b927d..3716412 100644 --- a/examples/07_regression_kriging2d.py +++ b/examples/07_regression_kriging2d.py @@ -11,9 +11,9 @@ from sklearn.ensemble import RandomForestRegressor from sklearn.linear_model import LinearRegression from sklearn.datasets import fetch_california_housing +from sklearn.model_selection import train_test_split from pykrige.rk import RegressionKriging -from pykrige.compat import train_test_split svr_model = SVR(C=0.1, gamma="auto") rf_model = RandomForestRegressor(n_estimators=100) diff --git a/examples/10_classification_kriging2d.py b/examples/10_classification_kriging2d.py new file mode 100644 index 0000000..c625754 --- /dev/null +++ b/examples/10_classification_kriging2d.py @@ -0,0 +1,50 @@ +""" +Classification kriging +---------------------- + +An example of classification kriging +""" + +import sys + +from sklearn.svm import SVC +from sklearn.ensemble import RandomForestClassifier +from sklearn.linear_model import LogisticRegression +from sklearn.datasets import fetch_california_housing +from sklearn.preprocessing import KBinsDiscretizer +from sklearn.model_selection import train_test_split + +from pykrige.ck import ClassificationKriging + +svc_model = SVC(C=0.1, gamma="auto", probability=True) +rf_model = RandomForestClassifier(n_estimators=100) +lr_model = LogisticRegression(max_iter=10000) + +models = [svc_model, rf_model, lr_model] + +try: + housing = fetch_california_housing() +except PermissionError: + # this dataset can occasionally fail to download on Windows + sys.exit(0) + +# take the first 5000 as Kriging is memory intensive +p = housing["data"][:5000, :-2] +x = housing["data"][:5000, -2:] +target = housing["target"][:5000] +discretizer = KBinsDiscretizer(encode="ordinal") +target = discretizer.fit_transform(target.reshape(-1, 1)) + +p_train, p_test, x_train, x_test, target_train, target_test = train_test_split( + p, x, target, test_size=0.3, random_state=42 +) + +for m in models: + print("=" * 40) + print("classification model:", m.__class__.__name__) + m_ck = ClassificationKriging(classification_model=m, n_closest_points=10) + m_ck.fit(p_train, x_train, target_train) + print( + "Classification Score: ", m_ck.classification_model.score(p_test, target_test) + ) + print("CK score: ", m_ck.score(p_test, x_test, target_test)) diff --git a/pykrige/ck.py b/pykrige/ck.py new file mode 100644 index 0000000..ea9ea50 --- /dev/null +++ b/pykrige/ck.py @@ -0,0 +1,290 @@ +# coding: utf-8 +"""Classification Kriging.""" +import numpy as np +from pykrige.compat import Krige, validate_sklearn, check_sklearn_model + +validate_sklearn() + +from sklearn.metrics import accuracy_score +from sklearn.svm import SVC +from sklearn.preprocessing import OneHotEncoder +from scipy.linalg import helmert + + +class ClassificationKriging: + """ + An implementation of Simplicial Indicator Kriging applied to classification ilr transformed residuals. + + Parameters + ---------- + classification_model: machine learning model instance from sklearn + method: str, optional + type of kriging to be performed + variogram_model: str, optional + variogram model to be used during Kriging + n_closest_points: int + number of closest points to be used during Ordinary Kriging + nlags: int + see OK/UK class description + weight: bool + see OK/UK class description + verbose: bool + see OK/UK class description + exact_values : bool + see OK/UK class description + variogram_parameters : list or dict + see OK/UK class description + variogram_function : callable + see OK/UK class description + anisotropy_scaling : tuple + single value for 2D (UK/OK) and two values in 3D (UK3D/OK3D) + anisotropy_angle : tuple + single value for 2D (UK/OK) and three values in 3D (UK3D/OK3D) + enable_statistics : bool + see OK class description + coordinates_type : str + see OK/UK class description + drift_terms : list of strings + see UK/UK3D class description + point_drift : array_like + see UK class description + ext_drift_grid : tuple + Holding the three values external_drift, external_drift_x and + external_drift_z for the UK class + functional_drift : list of callable + see UK/UK3D class description + """ + + def __init__( + self, + classification_model=SVC(), + method="ordinary", + variogram_model="linear", + n_closest_points=10, + nlags=6, + weight=False, + verbose=False, + exact_values=True, + pseudo_inv=False, + pseudo_inv_type="pinv", + variogram_parameters=None, + variogram_function=None, + anisotropy_scaling=(1.0, 1.0), + anisotropy_angle=(0.0, 0.0, 0.0), + enable_statistics=False, + coordinates_type="euclidean", + drift_terms=None, + point_drift=None, + ext_drift_grid=(None, None, None), + functional_drift=None, + ): + check_sklearn_model(classification_model, task="classification") + self.classification_model = classification_model + self.n_closest_points = n_closest_points + self._kriging_kwargs = dict( + method=method, + variogram_model=variogram_model, + nlags=nlags, + weight=weight, + n_closest_points=n_closest_points, + verbose=verbose, + exact_values=exact_values, + pseudo_inv=pseudo_inv, + pseudo_inv_type=pseudo_inv_type, + variogram_parameters=variogram_parameters, + variogram_function=variogram_function, + anisotropy_scaling=anisotropy_scaling, + anisotropy_angle=anisotropy_angle, + enable_statistics=enable_statistics, + coordinates_type=coordinates_type, + drift_terms=drift_terms, + point_drift=point_drift, + ext_drift_grid=ext_drift_grid, + functional_drift=functional_drift, + ) + + def fit(self, p, x, y): + """ + Fit the classification method and also krige the residual. + + Parameters + ---------- + p: ndarray + (Ns, d) array of predictor variables (Ns samples, d dimensions) + for classification + x: ndarray + ndarray of (x, y) points. Needs to be a (Ns, 2) array + corresponding to the lon/lat, for example 2d classification kriging. + array of Points, (x, y, z) pairs of shape (N, 3) for 3d kriging + y: ndarray + array of targets (Ns, ) + """ + self.classification_model.fit(p, y.ravel()) + print("Finished learning classification model") + self.classes_ = self.classification_model.classes_ + + self.krige = [] + for i in range(len(self.classes_) - 1): + self.krige.append(Krige(**self._kriging_kwargs)) + + ml_pred = self.classification_model.predict_proba(p) + ml_pred_ilr = ilr_transformation(ml_pred) + + self.onehotencode = OneHotEncoder(categories=[self.classes_]) + y_ohe = np.array(self.onehotencode.fit_transform(y).todense()) + y_ohe_ilr = ilr_transformation(y_ohe) + + for i in range(len(self.classes_) - 1): + self.krige[i].fit(x=x, y=y_ohe_ilr[:, i] - ml_pred_ilr[:, i]) + + print("Finished kriging residuals") + + def predict(self, p, x, **kwargs): + """ + Predict. + + Parameters + ---------- + p: ndarray + (Ns, d) array of predictor variables (Ns samples, d dimensions) + for classification + x: ndarray + ndarray of (x, y) points. Needs to be a (Ns, 2) array + corresponding to the lon/lat, for example. + array of Points, (x, y, z) pairs of shape (N, 3) for 3d kriging + + Returns + ------- + pred: ndarray + The expected value of ys for the query inputs, of shape (Ns,). + + """ + + ml_pred = self.classification_model.predict_proba(p) + ml_pred_ilr = ilr_transformation(ml_pred) + + pred_proba_ilr = self.krige_residual(x, **kwargs) + ml_pred_ilr + pred_proba = inverse_ilr_transformation(pred_proba_ilr) + + return np.argmax(pred_proba, axis=1) + + def krige_residual(self, x, **kwargs): + """ + Calculate the residuals. + + Parameters + ---------- + x: ndarray + ndarray of (x, y) points. Needs to be a (Ns, 2) array + corresponding to the lon/lat, for example. + + Returns + ------- + residual: ndarray + kriged residual values + """ + + krig_pred = [ + self.krige[i].predict(x=x, **kwargs) for i in range(len(self.classes_) - 1) + ] + + return np.vstack(krig_pred).T + + def score(self, p, x, y, sample_weight=None, **kwargs): + """ + Overloading default classification score method. + + Parameters + ---------- + p: ndarray + (Ns, d) array of predictor variables (Ns samples, d dimensions) + for classification + x: ndarray + ndarray of (x, y) points. Needs to be a (Ns, 2) array + corresponding to the lon/lat, for example. + array of Points, (x, y, z) pairs of shape (N, 3) for 3d kriging + y: ndarray + array of targets (Ns, ) + """ + return accuracy_score( + y_pred=self.predict(p, x, **kwargs), y_true=y, sample_weight=sample_weight + ) + + +def closure(data, k=1.0): + """Apply closure to data, sample-wise. + Adapted from https://github.com/ofgulban/compoda. + + Parameters + ---------- + data : 2d numpy array, shape [n_samples, n_measurements] + Data to be closed to a certain constant. Do not forget to deal with + zeros in the data before this operation. + k : float, positive + Sum of the measurements will be equal to this number. + + Returns + ------- + data : 2d numpy array, shape [n_samples, n_measurements] + Closed data. + + Reference + --------- + [1] Pawlowsky-Glahn, V., Egozcue, J. J., & Tolosana-Delgado, R. + (2015). Modelling and Analysis of Compositional Data, pg. 9. + Chichester, UK: John Wiley & Sons, Ltd. + DOI: 10.1002/9781119003144 + """ + + return k * data / np.sum(data, axis=1)[:, np.newaxis] + + +def ilr_transformation(data): + """Isometric logratio transformation (not vectorized). + Adapted from https://github.com/ofgulban/compoda. + + Parameters + ---------- + data : 2d numpy array, shape [n_samples, n_coordinates] + Barycentric coordinates (closed) in simplex space. + + Returns + ------- + out : 2d numpy array, shape [n_samples, n_coordinates-1] + Coordinates in real space. + + Reference + --------- + [1] Pawlowsky-Glahn, V., Egozcue, J. J., & Tolosana-Delgado, R. + (2015). Modelling and Analysis of Compositional Data, pg. 37. + Chichester, UK: John Wiley & Sons, Ltd. + DOI: 10.1002/9781119003144 + """ + data = np.maximum(data, np.finfo(float).eps) + + return np.einsum("ij,jk->ik", np.log(data), -helmert(data.shape[1]).T) + + +def inverse_ilr_transformation(data): + """Inverse isometric logratio transformation (not vectorized). + Adapted from https://github.com/ofgulban/compoda. + + Parameters + ---------- + data : 2d numpy array, shape [n_samples, n_coordinates] + Isometric log-ratio transformed coordinates in real space. + + Returns + ------- + out : 2d numpy array, shape [n_samples, n_coordinates+1] + Barycentric coordinates (closed) in simplex space. + + Reference + --------- + [1] Pawlowsky-Glahn, V., Egozcue, J. J., & Tolosana-Delgado, R. + (2015). Modelling and Analysis of Compositional Data, pg. 37. + Chichester, UK: John Wiley & Sons, Ltd. + DOI: 10.1002/9781119003144 + """ + + return closure(np.exp(np.einsum("ij,jk->ik", data, -helmert(data.shape[1] + 1)))) diff --git a/pykrige/compat.py b/pykrige/compat.py index 73434aa..d092606 100644 --- a/pykrige/compat.py +++ b/pykrige/compat.py @@ -1,26 +1,307 @@ # coding: utf-8 # pylint: disable= invalid-name, unused-import -"""For compatibility""" -from functools import partial - +"""For compatibility.""" +from pykrige.uk3d import UniversalKriging3D +from pykrige.ok3d import OrdinaryKriging3D +from pykrige.uk import UniversalKriging +from pykrige.ok import OrdinaryKriging # sklearn try: - from sklearn.model_selection import GridSearchCV + # keep train_test_split here for backward compatibility from sklearn.model_selection import train_test_split + from sklearn.base import RegressorMixin, ClassifierMixin, BaseEstimator SKLEARN_INSTALLED = True except ImportError: SKLEARN_INSTALLED = False + train_test_split = None + + class RegressorMixin: + """Mock RegressorMixin.""" + + class ClassifierMixin: + """Mock ClassifierMixin.""" + + class BaseEstimator: + """Mock BaseEstimator.""" + + +krige_methods = { + "ordinary": OrdinaryKriging, + "universal": UniversalKriging, + "ordinary3d": OrdinaryKriging3D, + "universal3d": UniversalKriging3D, +} + +threed_krige = ("ordinary3d", "universal3d") + +krige_methods_kws = { + "ordinary": [ + "anisotropy_scaling", + "anisotropy_angle", + "enable_statistics", + "coordinates_type", + ], + "universal": [ + "anisotropy_scaling", + "anisotropy_angle", + "drift_terms", + "point_drift", + "external_drift", + "external_drift_x", + "external_drift_y", + "functional_drift", + ], + "ordinary3d": [ + "anisotropy_scaling_y", + "anisotropy_scaling_z", + "anisotropy_angle_x", + "anisotropy_angle_y", + "anisotropy_angle_z", + ], + "universal3d": [ + "anisotropy_scaling_y", + "anisotropy_scaling_z", + "anisotropy_angle_x", + "anisotropy_angle_y", + "anisotropy_angle_z", + "drift_terms", + "functional_drift", + ], +} + class SklearnException(Exception): - pass + """Exception for missing scikit-learn.""" + + +def validate_method(method): + """Validate the kriging method in use.""" + if method not in krige_methods.keys(): + raise ValueError( + "Kriging method must be one of {}".format(krige_methods.keys()) + ) def validate_sklearn(): + """Validate presence of scikit-learn.""" if not SKLEARN_INSTALLED: raise SklearnException( "sklearn needs to be installed in order to use this module" ) + + +class Krige(RegressorMixin, BaseEstimator): + """ + A scikit-learn wrapper class for Ordinary and Universal Kriging. + + This works with both Grid/RandomSearchCv for finding the best + Krige parameters combination for a problem. + + Parameters + ---------- + method: str, optional + type of kriging to be performed + variogram_model: str, optional + variogram model to be used during Kriging + nlags: int + see OK/UK class description + weight: bool + see OK/UK class description + n_closest_points: int + number of closest points to be used during Ordinary Kriging + verbose: bool + see OK/UK class description + exact_values : bool + see OK/UK class description + variogram_parameters : list or dict + see OK/UK class description + variogram_function : callable + see OK/UK class description + anisotropy_scaling : tuple + single value for 2D (UK/OK) and two values in 3D (UK3D/OK3D) + anisotropy_angle : tuple + single value for 2D (UK/OK) and three values in 3D (UK3D/OK3D) + enable_statistics : bool + see OK class description + coordinates_type : str + see OK/UK class description + drift_terms : list of strings + see UK/UK3D class description + point_drift : array_like + see UK class description + ext_drift_grid : tuple + Holding the three values external_drift, external_drift_x and + external_drift_z for the UK class + functional_drift : list of callable + see UK/UK3D class description + """ + + def __init__( + self, + method="ordinary", + variogram_model="linear", + nlags=6, + weight=False, + n_closest_points=10, + verbose=False, + exact_values=True, + pseudo_inv=False, + pseudo_inv_type="pinv", + variogram_parameters=None, + variogram_function=None, + anisotropy_scaling=(1.0, 1.0), + anisotropy_angle=(0.0, 0.0, 0.0), + enable_statistics=False, + coordinates_type="euclidean", + drift_terms=None, + point_drift=None, + ext_drift_grid=(None, None, None), + functional_drift=None, + ): + validate_method(method) + self.variogram_model = variogram_model + self.variogram_parameters = variogram_parameters + self.variogram_function = variogram_function + self.nlags = nlags + self.weight = weight + self.verbose = verbose + self.exact_values = exact_values + self.pseudo_inv = pseudo_inv + self.pseudo_inv_type = pseudo_inv_type + self.anisotropy_scaling = anisotropy_scaling + self.anisotropy_angle = anisotropy_angle + self.enable_statistics = enable_statistics + self.coordinates_type = coordinates_type + self.drift_terms = drift_terms + self.point_drift = point_drift + self.ext_drift_grid = ext_drift_grid + self.functional_drift = functional_drift + self.model = None # not trained + self.n_closest_points = n_closest_points + self.method = method + + def fit(self, x, y, *args, **kwargs): + """ + Fit the current model. + + Parameters + ---------- + x: ndarray + array of Points, (x, y) pairs of shape (N, 2) for 2d kriging + array of Points, (x, y, z) pairs of shape (N, 3) for 3d kriging + y: ndarray + array of targets (N, ) + """ + val_kw = "val" if self.method in threed_krige else "z" + setup = dict( + variogram_model=self.variogram_model, + variogram_parameters=self.variogram_parameters, + variogram_function=self.variogram_function, + nlags=self.nlags, + weight=self.weight, + verbose=self.verbose, + exact_values=self.exact_values, + pseudo_inv=self.pseudo_inv, + pseudo_inv_type=self.pseudo_inv_type, + ) + add_setup = dict( + anisotropy_scaling=self.anisotropy_scaling[0], + anisotropy_angle=self.anisotropy_angle[0], + enable_statistics=self.enable_statistics, + coordinates_type=self.coordinates_type, + anisotropy_scaling_y=self.anisotropy_scaling[0], + anisotropy_scaling_z=self.anisotropy_scaling[1], + anisotropy_angle_x=self.anisotropy_angle[0], + anisotropy_angle_y=self.anisotropy_angle[1], + anisotropy_angle_z=self.anisotropy_angle[2], + drift_terms=self.drift_terms, + point_drift=self.point_drift, + external_drift=self.ext_drift_grid[0], + external_drift_x=self.ext_drift_grid[1], + external_drift_y=self.ext_drift_grid[2], + functional_drift=self.functional_drift, + ) + for kw in krige_methods_kws[self.method]: + setup[kw] = add_setup[kw] + input_kw = self._dimensionality_check(x) + input_kw.update(setup) + input_kw[val_kw] = y + self.model = krige_methods[self.method](**input_kw) + + def _dimensionality_check(self, x, ext=""): + if self.method in ("ordinary", "universal"): + if x.shape[1] != 2: + raise ValueError("2d krige can use only 2d points") + else: + return {"x" + ext: x[:, 0], "y" + ext: x[:, 1]} + if self.method in ("ordinary3d", "universal3d"): + if x.shape[1] != 3: + raise ValueError("3d krige can use only 3d points") + else: + return { + "x" + ext: x[:, 0], + "y" + ext: x[:, 1], + "z" + ext: x[:, 2], + } + + def predict(self, x, *args, **kwargs): + """ + Predict. + + Parameters + ---------- + x: ndarray + array of Points, (x, y) pairs of shape (N, 2) for 2d kriging + array of Points, (x, y, z) pairs of shape (N, 3) for 3d kriging + Returns + ------- + Prediction array + """ + if not self.model: + raise Exception("Not trained. Train first") + points = self._dimensionality_check(x, ext="points") + return self.execute(points, *args, **kwargs)[0] + + def execute(self, points, *args, **kwargs): + # TODO array of Points, (x, y) pairs of shape (N, 2) + """ + Execute. + + Parameters + ---------- + points: dict + + Returns + ------- + Prediction array + Variance array + """ + default_kw = dict(style="points", backend="loop") + default_kw.update(kwargs) + points.update(default_kw) + if isinstance(self.model, (OrdinaryKriging, OrdinaryKriging3D)): + points.update(dict(n_closest_points=self.n_closest_points)) + else: + print("n_closest_points will be ignored for UniversalKriging") + prediction, variance = self.model.execute(**points) + return prediction, variance + + +def check_sklearn_model(model, task="regression"): + """Check the sklearn method in use.""" + if task == "regression": + if not (isinstance(model, BaseEstimator) and isinstance(model, RegressorMixin)): + raise RuntimeError( + "Needs to supply an instance of a scikit-learn regression class." + ) + elif task == "classification": + if not ( + isinstance(model, BaseEstimator) and isinstance(model, ClassifierMixin) + ): + raise RuntimeError( + "Needs to supply an instance of a scikit-learn classification class." + ) diff --git a/pykrige/ok.py b/pykrige/ok.py index 358520d..92d42a7 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -958,7 +958,12 @@ def execute( zvalues, sigmasq = self._exec_loop_moving_window(a, bd, mask, bd_idx) elif backend == "C": zvalues, sigmasq = _c_exec_loop_moving_window( - a, bd, mask.astype("int8"), bd_idx, self.X_ADJUSTED.shape[0], c_pars + a, + bd, + mask.astype("int8"), + bd_idx.astype(int), + self.X_ADJUSTED.shape[0], + c_pars, ) else: raise ValueError( diff --git a/pykrige/rk.py b/pykrige/rk.py index e7b4760..70fcc14 100644 --- a/pykrige/rk.py +++ b/pykrige/rk.py @@ -1,271 +1,11 @@ # coding: utf-8 -from pykrige.compat import validate_sklearn +"""Regression Kriging.""" +from pykrige.compat import Krige, validate_sklearn, check_sklearn_model validate_sklearn() -from pykrige.ok import OrdinaryKriging -from pykrige.uk import UniversalKriging -from pykrige.ok3d import OrdinaryKriging3D -from pykrige.uk3d import UniversalKriging3D -from sklearn.base import RegressorMixin, BaseEstimator -from sklearn.svm import SVR -from sklearn.metrics import r2_score - -krige_methods = { - "ordinary": OrdinaryKriging, - "universal": UniversalKriging, - "ordinary3d": OrdinaryKriging3D, - "universal3d": UniversalKriging3D, -} - -threed_krige = ("ordinary3d", "universal3d") - -krige_methods_kws = { - "ordinary": [ - "anisotropy_scaling", - "anisotropy_angle", - "enable_statistics", - "coordinates_type", - ], - "universal": [ - "anisotropy_scaling", - "anisotropy_angle", - "drift_terms", - "point_drift", - "external_drift", - "external_drift_x", - "external_drift_y", - "functional_drift", - ], - "ordinary3d": [ - "anisotropy_scaling_y", - "anisotropy_scaling_z", - "anisotropy_angle_x", - "anisotropy_angle_y", - "anisotropy_angle_z", - ], - "universal3d": [ - "anisotropy_scaling_y", - "anisotropy_scaling_z", - "anisotropy_angle_x", - "anisotropy_angle_y", - "anisotropy_angle_z", - "drift_terms", - "functional_drift", - ], -} - - -def validate_method(method): - """Validate the kriging method in use.""" - if method not in krige_methods.keys(): - raise ValueError( - "Kriging method must be one of {}".format(krige_methods.keys()) - ) - - -class Krige(RegressorMixin, BaseEstimator): - """ - A scikit-learn wrapper class for Ordinary and Universal Kriging. - - This works with both Grid/RandomSearchCv for finding the best - Krige parameters combination for a problem. - - Parameters - ---------- - method: str, optional - type of kriging to be performed - variogram_model: str, optional - variogram model to be used during Kriging - nlags: int - see OK/UK class description - weight: bool - see OK/UK class description - n_closest_points: int - number of closest points to be used during Ordinary Kriging - verbose: bool - see OK/UK class description - exact_values : bool - see OK/UK class description - variogram_parameters : list or dict - see OK/UK class description - variogram_function : callable - see OK/UK class description - anisotropy_scaling : tuple - single value for 2D (UK/OK) and two values in 3D (UK3D/OK3D) - anisotropy_angle : tuple - single value for 2D (UK/OK) and three values in 3D (UK3D/OK3D) - enable_statistics : bool - see OK class description - coordinates_type : str - see OK/UK class description - drift_terms : list of strings - see UK/UK3D class description - point_drift : array_like - see UK class description - ext_drift_grid : tuple - Holding the three values external_drift, external_drift_x and - external_drift_z for the UK class - functional_drift : list of callable - see UK/UK3D class description - """ - def __init__( - self, - method="ordinary", - variogram_model="linear", - nlags=6, - weight=False, - n_closest_points=10, - verbose=False, - exact_values=True, - pseudo_inv=False, - pseudo_inv_type="pinv", - variogram_parameters=None, - variogram_function=None, - anisotropy_scaling=(1.0, 1.0), - anisotropy_angle=(0.0, 0.0, 0.0), - enable_statistics=False, - coordinates_type="euclidean", - drift_terms=None, - point_drift=None, - ext_drift_grid=(None, None, None), - functional_drift=None, - ): - validate_method(method) - self.variogram_model = variogram_model - self.variogram_parameters = variogram_parameters - self.variogram_function = variogram_function - self.nlags = nlags - self.weight = weight - self.verbose = verbose - self.exact_values = exact_values - self.pseudo_inv = pseudo_inv - self.pseudo_inv_type = pseudo_inv_type - self.anisotropy_scaling = anisotropy_scaling - self.anisotropy_angle = anisotropy_angle - self.enable_statistics = enable_statistics - self.coordinates_type = coordinates_type - self.drift_terms = drift_terms - self.point_drift = point_drift - self.ext_drift_grid = ext_drift_grid - self.functional_drift = functional_drift - self.model = None # not trained - self.n_closest_points = n_closest_points - self.method = method - - def fit(self, x, y, *args, **kwargs): - """ - Fit the current model. - - Parameters - ---------- - x: ndarray - array of Points, (x, y) pairs of shape (N, 2) for 2d kriging - array of Points, (x, y, z) pairs of shape (N, 3) for 3d kriging - y: ndarray - array of targets (N, ) - """ - val_kw = "val" if self.method in threed_krige else "z" - setup = dict( - variogram_model=self.variogram_model, - variogram_parameters=self.variogram_parameters, - variogram_function=self.variogram_function, - nlags=self.nlags, - weight=self.weight, - verbose=self.verbose, - exact_values=self.exact_values, - pseudo_inv=self.pseudo_inv, - pseudo_inv_type=self.pseudo_inv_type, - ) - add_setup = dict( - anisotropy_scaling=self.anisotropy_scaling[0], - anisotropy_angle=self.anisotropy_angle[0], - enable_statistics=self.enable_statistics, - coordinates_type=self.coordinates_type, - anisotropy_scaling_y=self.anisotropy_scaling[0], - anisotropy_scaling_z=self.anisotropy_scaling[1], - anisotropy_angle_x=self.anisotropy_angle[0], - anisotropy_angle_y=self.anisotropy_angle[1], - anisotropy_angle_z=self.anisotropy_angle[2], - drift_terms=self.drift_terms, - point_drift=self.point_drift, - external_drift=self.ext_drift_grid[0], - external_drift_x=self.ext_drift_grid[1], - external_drift_y=self.ext_drift_grid[2], - functional_drift=self.functional_drift, - ) - for kw in krige_methods_kws[self.method]: - setup[kw] = add_setup[kw] - input_kw = self._dimensionality_check(x) - input_kw.update(setup) - input_kw[val_kw] = y - self.model = krige_methods[self.method](**input_kw) - - def _dimensionality_check(self, x, ext=""): - if self.method in ("ordinary", "universal"): - if x.shape[1] != 2: - raise ValueError("2d krige can use only 2d points") - else: - return {"x" + ext: x[:, 0], "y" + ext: x[:, 1]} - if self.method in ("ordinary3d", "universal3d"): - if x.shape[1] != 3: - raise ValueError("3d krige can use only 3d points") - else: - return { - "x" + ext: x[:, 0], - "y" + ext: x[:, 1], - "z" + ext: x[:, 2], - } - - def predict(self, x, *args, **kwargs): - """ - Predict. - - Parameters - ---------- - x: ndarray - array of Points, (x, y) pairs of shape (N, 2) for 2d kriging - array of Points, (x, y, z) pairs of shape (N, 3) for 3d kriging - Returns - ------- - Prediction array - """ - if not self.model: - raise Exception("Not trained. Train first") - points = self._dimensionality_check(x, ext="points") - return self.execute(points, *args, **kwargs)[0] - - def execute(self, points, *args, **kwargs): - # TODO array of Points, (x, y) pairs of shape (N, 2) - """ - Execute. - - Parameters - ---------- - points: dict - - Returns - ------- - Prediction array - Variance array - """ - default_kw = dict(style="points", backend="loop") - default_kw.update(kwargs) - points.update(default_kw) - if isinstance(self.model, (OrdinaryKriging, OrdinaryKriging3D)): - points.update(dict(n_closest_points=self.n_closest_points)) - else: - print("n_closest_points will be ignored for UniversalKriging") - prediction, variance = self.model.execute(**points) - return prediction, variance - - -def check_sklearn_model(model): - """Check the sklearn method in use.""" - if not (isinstance(model, BaseEstimator) and isinstance(model, RegressorMixin)): - raise RuntimeError( - "Needs to supply an instance of a scikit-learn regression class." - ) +from sklearn.metrics import r2_score +from sklearn.svm import SVR class RegressionKriging: diff --git a/pykrige/uk.py b/pykrige/uk.py index 707cf34..6b6863d 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -969,7 +969,8 @@ def _exec_vector(self, a, bd, xy, xy_orig, mask, n_withdrifts, spec_drift_grids) i += 1 if i != n_withdrifts: warnings.warn( - "Error in setting up kriging system. Kriging may fail.", RuntimeWarning, + "Error in setting up kriging system. Kriging may fail.", + RuntimeWarning, ) if self.UNBIAS: b[:, n_withdrifts, 0] = 1.0 diff --git a/pykrige/uk3d.py b/pykrige/uk3d.py index 3fe8813..115893c 100644 --- a/pykrige/uk3d.py +++ b/pykrige/uk3d.py @@ -778,7 +778,8 @@ def _exec_vector(self, a, bd, xyz, mask, n_withdrifts, spec_drift_grids): i += 1 if i != n_withdrifts: warnings.warn( - "Error in setting up kriging system. Kriging may fail.", RuntimeWarning, + "Error in setting up kriging system. Kriging may fail.", + RuntimeWarning, ) if self.UNBIAS: b[:, n_withdrifts, 0] = 1.0 diff --git a/setup.py b/setup.py index 0962fdf..e9a2364 100755 --- a/setup.py +++ b/setup.py @@ -35,7 +35,7 @@ # setup ####################################################################### -with open(os.path.join(HERE, "README.rst"), encoding="utf-8") as f: +with open(os.path.join(HERE, "README.md"), encoding="utf-8") as f: README = f.read() with open(os.path.join(HERE, "requirements.txt"), encoding="utf-8") as f: REQ = f.read().splitlines() @@ -64,10 +64,10 @@ "Operating System :: Unix", "Programming Language :: Python", "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.5", "Programming Language :: Python :: 3.6", "Programming Language :: Python :: 3.7", "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3 :: Only", "Topic :: Scientific/Engineering", "Topic :: Scientific/Engineering :: GIS", @@ -78,7 +78,7 @@ name="PyKrige", description=DOCLINE, long_description=README, - long_description_content_type="text/x-rst", + long_description_content_type="text/markdown", author="Benjamin S. Murphy", author_email="bscott.murphy@gmail.com", maintainer="Sebastian Mueller, Roman Yurchak", @@ -88,7 +88,7 @@ classifiers=CLASSIFIERS, platforms=["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"], include_package_data=True, - python_requires=">=3.5", + python_requires=">=3.6", use_scm_version={ "relative_to": __file__, "write_to": "pykrige/_version.py", diff --git a/tests/test_api.py b/tests/test_api.py index 98d2de0..6c033ca 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -3,8 +3,8 @@ import numpy as np import pytest -from pykrige.rk import Krige -from pykrige.rk import threed_krige +from pykrige.compat import Krige +from pykrige.compat import threed_krige def _method_and_vergiogram(): diff --git a/tests/test_classification_krige.py b/tests/test_classification_krige.py new file mode 100644 index 0000000..3e5e93a --- /dev/null +++ b/tests/test_classification_krige.py @@ -0,0 +1,101 @@ +from itertools import product +import pytest + +import numpy as np + +from pykrige.ck import ClassificationKriging + +try: + from sklearn.svm import SVC + from sklearn.datasets import fetch_california_housing + from sklearn.ensemble import RandomForestClassifier + from sklearn.preprocessing import KBinsDiscretizer + from sklearn.model_selection import train_test_split + + SKLEARN_INSTALLED = True +except ImportError: + SKLEARN_INSTALLED = False + + +def _methods(): + krige_methods = ["ordinary", "universal"] + ml_methods = [ + SVC(C=0.01, gamma="auto", probability=True), + RandomForestClassifier(n_estimators=50), + ] + return product(ml_methods, krige_methods) + + +@pytest.mark.skipif(not SKLEARN_INSTALLED, reason="requires scikit-learn") +def test_classification_krige(): + np.random.seed(1) + x = np.linspace(-1.0, 1.0, 100) + # create a feature matrix with 5 features + X = np.tile(x, reps=(5, 1)).T + y = ( + 1 + + 5 * X[:, 0] + - 2 * X[:, 1] + - 2 * X[:, 2] + + 3 * X[:, 3] + + 4 * X[:, 4] + + 2 * (np.random.rand(100) - 0.5) + ) + + # create lat/lon array + lon = np.linspace(-180.0, 180.0, 10) + lat = np.linspace(-90.0, 90.0, 10) + lon_lat = np.array(list(product(lon, lat))) + + discretizer = KBinsDiscretizer(encode="ordinal") + y = discretizer.fit_transform(y.reshape(-1, 1)) + + X_train, X_test, y_train, y_test, lon_lat_train, lon_lat_test = train_test_split( + X, y, lon_lat, train_size=0.7, random_state=10 + ) + + for ml_model, krige_method in _methods(): + class_model = ClassificationKriging( + classification_model=ml_model, method=krige_method, n_closest_points=2 + ) + class_model.fit(X_train, lon_lat_train, y_train) + assert class_model.score(X_test, lon_lat_test, y_test) > 0.25 + + +@pytest.mark.skipif(not SKLEARN_INSTALLED, reason="requires scikit-learn") +def test_krige_classification_housing(): + import ssl + import urllib + + try: + housing = fetch_california_housing() + except (ssl.SSLError, urllib.error.URLError): + ssl._create_default_https_context = ssl._create_unverified_context + try: + housing = fetch_california_housing() + except PermissionError: + # This can raise permission error on Appveyor + pytest.skip("Failed to load california housing dataset") + ssl._create_default_https_context = ssl.create_default_context + + # take only first 1000 + p = housing["data"][:1000, :-2] + x = housing["data"][:1000, -2:] + target = housing["target"][:1000] + discretizer = KBinsDiscretizer(encode="ordinal") + target = discretizer.fit_transform(target.reshape(-1, 1)) + + p_train, p_test, y_train, y_test, x_train, x_test = train_test_split( + p, target, x, train_size=0.7, random_state=10 + ) + + for ml_model, krige_method in _methods(): + + class_model = ClassificationKriging( + classification_model=ml_model, method=krige_method, n_closest_points=2 + ) + class_model.fit(p_train, x_train, y_train) + if krige_method == "ordinary": + assert class_model.score(p_test, x_test, y_test) > 0.5 + else: + assert class_model.score(p_test, x_test, y_test) > 0.0 diff --git a/tests/test_core.py b/tests/test_core.py index cba4e00..6b07f7d 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -439,7 +439,12 @@ def test_core_krige_3d(): def test_non_exact(): # custom data for this test data = np.array( - [[0.0, 0.0, 0.47], [1.5, 1.5, 0.56], [3, 3, 0.74], [4.5, 4.5, 1.47],] + [ + [0.0, 0.0, 0.47], + [1.5, 1.5, 0.56], + [3, 3, 0.74], + [4.5, 4.5, 1.47], + ] ) # construct grid points so diagonal diff --git a/tests/test_regression_krige.py b/tests/test_regression_krige.py index 1eeeb21..631af56 100644 --- a/tests/test_regression_krige.py +++ b/tests/test_regression_krige.py @@ -11,7 +11,7 @@ from sklearn.linear_model import ElasticNet, Lasso from sklearn.ensemble import RandomForestRegressor from sklearn.linear_model import LinearRegression - from pykrige.compat import train_test_split + from sklearn.model_selection import train_test_split SKLEARN_INSTALLED = True except ImportError: