diff --git a/.github/workflows/build-docs.yaml b/.github/workflows/build-docs.yaml index f5f778f29..33b95200f 100644 --- a/.github/workflows/build-docs.yaml +++ b/.github/workflows/build-docs.yaml @@ -12,6 +12,8 @@ jobs: steps: - name: Clone repo uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7 + with: + submodules: "true" - name: Clone docs repo uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7 with: diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 4d9a1aa3e..a386130e4 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -14,15 +14,85 @@ env: NIXTLA_NUMBA_CACHE: "1" jobs: + build_wheels: + name: Build wheels for cp${{ matrix.python-version }}-${{ matrix.os-platform[1] }} + runs-on: ${{ matrix.os-platform[0] }} + strategy: + fail-fast: false + matrix: + python-version: [38, 39, 310, 311, 312] + os-platform: + [ + [ubuntu-latest, manylinux_x86_64], + [ubuntu-latest, manylinux_aarch64], + [windows-latest, win_amd64], + [macos-12, macosx_x86_64], + [macos-14, macosx_arm64], + ] + steps: + - uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7 + with: + submodules: "true" + + - name: Set up QEMU + if: matrix.os-platform[1] == 'manylinux_aarch64' + uses: docker/setup-qemu-action@49b3bc8e6bdd4a60e6116a5414239cba5943d3cf # v3.2.0 + with: + platforms: arm64 + + - name: Build wheels + uses: pypa/cibuildwheel@bd033a44476646b606efccdd5eed92d5ea1d77ad # v2.20.0 + env: + CIBW_BUILD: cp${{ matrix.python-version }}-${{ matrix.os-platform[1] }} + + - uses: actions/upload-artifact@834a144ee995460fba8ed112a2fc961b36a5ec5a # v4.3.6 + with: + name: artifact-cp${{ matrix.python-version }}-${{ matrix.os-platform[1] }} + path: wheelhouse/*.whl + retention-days: 1 + + publish-nightly-wheels: + if: github.event_name == 'push' + needs: [build_wheels] + env: + AWS_ACCESS_KEY_ID: ${{ secrets.AWS_ACCESS_KEY_ID_NIXTLA_PACKAGES }} + AWS_SECRET_ACCESS_KEY: ${{ secrets.AWS_SECRET_ACCESS_KEY_NIXTLA_PACKAGES }} + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7 + + - name: Download wheels + uses: actions/download-artifact@fa0a91b85d4f404e444e00e005971372dc801d16 # v4.1.8 + with: + path: dist + pattern: "artifact-*" + merge-multiple: true + + - uses: actions/setup-python@39cd14951b08e74b54015e9e001cdefcf80e669f # v5.1.1 + with: + python-version: "3.10" + + - name: Publish wheels + run: | + pip install awscli beautifulsoup4 + aws s3 cp s3://nixtla-packages/statsforecast/index.html . + python action_files/update_index.py + aws s3 sync dist s3://nixtla-packages/statsforecast/ + run-local-tests: - runs-on: windows-latest + needs: [build_wheels] + runs-on: ${{ matrix.os }} + timeout-minutes: 60 strategy: fail-fast: false matrix: + os: [macos-13, macos-14, windows-latest] python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] steps: - name: Clone repo uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7 + with: + submodules: "true" - name: Set up environment uses: actions/setup-python@39cd14951b08e74b54015e9e001cdefcf80e669f # v5.1.1 @@ -30,22 +100,34 @@ jobs: python-version: ${{ matrix.python-version }} - name: Install pip requirements - run: pip install uv && uv pip install --system ".[dev,plotly]" + run: pip install uv && uv pip install --system -r setup.py --extra dev --extra plotly + + - name: Download wheels + uses: actions/download-artifact@fa0a91b85d4f404e444e00e005971372dc801d16 # v4.1.8 + with: + path: dist + pattern: "artifact-*" + merge-multiple: true + + - name: Install the library + run: pip install --no-index --find-links dist statsforecast - name: Run local tests - run: nbdev_test --skip_file_re '(models|distributed|ets|prophet).*.ipynb' --pause 1.0 --do_print --timing + run: nbdev_test --skip_file_re '(distributed|prophet).*.ipynb' --pause 1.0 --do_print --timing run-tests: - runs-on: ${{ matrix.os }} + needs: [build_wheels] + runs-on: ubuntu-latest timeout-minutes: 60 strategy: fail-fast: false matrix: - os: [ubuntu-latest, macos-13] python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] steps: - name: Clone repo uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7 + with: + submodules: "true" - name: Set up environment uses: actions/setup-python@39cd14951b08e74b54015e9e001cdefcf80e669f # v5.1.1 @@ -53,7 +135,17 @@ jobs: python-version: ${{ matrix.python-version }} - name: Install pip requirements - run: pip install uv && uv pip install --system ".[all]" + run: pip install uv && uv pip install --system -r setup.py --extra all + + - name: Download wheels + uses: actions/download-artifact@fa0a91b85d4f404e444e00e005971372dc801d16 # v4.1.8 + with: + path: dist + pattern: "artifact-*" + merge-multiple: true + + - name: Install the library + run: pip install --no-index --find-links dist statsforecast - name: Run tests run: nbdev_test --skip_file_re '(distributed|prophet).*.ipynb' --pause 1.0 --do_print --timing @@ -64,6 +156,7 @@ jobs: pytest --durations=0 action_files test-m3-performance: + needs: [build_wheels] runs-on: ubuntu-latest steps: - name: Clone repo @@ -74,8 +167,18 @@ jobs: with: python-version: "3.10" - - name: Install library and extra deps - run: pip install uv && uv pip install --system ".[dev]" + - name: Install pip requirements + run: pip install uv && uv pip install --system -r setup.py --extra dev + + - name: Download wheels + uses: actions/download-artifact@fa0a91b85d4f404e444e00e005971372dc801d16 # v4.1.8 + with: + path: dist + pattern: "artifact-*" + merge-multiple: true + + - name: Install the library + run: pip install --no-index --find-links dist statsforecast - name: Run M3 experiment run: | diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml index c8f175034..15a3096d4 100644 --- a/.github/workflows/lint.yaml +++ b/.github/workflows/lint.yaml @@ -16,10 +16,10 @@ jobs: - name: Set up python uses: actions/setup-python@39cd14951b08e74b54015e9e001cdefcf80e669f # v5.1.1 with: - python-version: '3.10' + python-version: "3.10" - name: Install linters - run: pip install black nbdev pre-commit + run: pip install black 'nbdev<2.3.26' pre-commit - name: Run pre-commit - run: pre-commit run --files statsforecast/* + run: pre-commit run --files python/statsforecast/* diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml index 481a49894..a4b096376 100644 --- a/.github/workflows/python-publish.yml +++ b/.github/workflows/python-publish.yml @@ -1,25 +1,77 @@ name: Upload Python Package on: - release: - types: [published] + push: + tags: ["v*"] jobs: - deploy: + build_wheels: + name: Build wheels for cp${{ matrix.python-version }}-${{ matrix.os-platform[1] }} + runs-on: ${{ matrix.os-platform[0] }} + strategy: + fail-fast: false + matrix: + python-version: [38, 39, 310, 311, 312] + os-platform: + [ + [ubuntu-latest, manylinux_x86_64], + [ubuntu-latest, manylinux_aarch64], + [windows-latest, win_amd64], + [macos-12, macosx_x86_64], + [macos-14, macosx_arm64], + ] + steps: + - uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7 + with: + submodules: "true" + + - name: Set up QEMU + if: matrix.os-platform[1] == 'manylinux_aarch64' + uses: docker/setup-qemu-action@49b3bc8e6bdd4a60e6116a5414239cba5943d3cf # v3.2.0 + with: + platforms: arm64 + + - name: Build wheels + uses: pypa/cibuildwheel@bd033a44476646b606efccdd5eed92d5ea1d77ad # v2.20.0 + env: + CIBW_BUILD: cp${{ matrix.python-version }}-${{ matrix.os-platform[1] }} + + - uses: actions/upload-artifact@834a144ee995460fba8ed112a2fc961b36a5ec5a # v4.3.6 + with: + name: artifact-cp${{ matrix.python-version }}-${{ matrix.os-platform[1] }} + path: wheelhouse/*.whl + retention-days: 1 + + build_sdist: + name: Build sdist + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7 + + - uses: actions/setup-python@39cd14951b08e74b54015e9e001cdefcf80e669f # v5.1.1 + with: + python-version: "3.10" + + - name: Build sdist + run: | + python -m pip install build + python -m build --sdist --outdir dist + - uses: actions/upload-artifact@834a144ee995460fba8ed112a2fc961b36a5ec5a # v4.3.6 + with: + path: dist/*.tar.gz + + upload_to_pypi: + name: Upload to PyPI + needs: [build_wheels, build_sdist] runs-on: ubuntu-latest permissions: id-token: write - steps: - - uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7 - - name: Set up Python - uses: actions/setup-python@39cd14951b08e74b54015e9e001cdefcf80e669f # v5.1.1 + - uses: actions/download-artifact@fa0a91b85d4f404e444e00e005971372dc801d16 # v4.1.8 with: - python-version: '3.10' - - name: Install dependencies - run: python -m pip install build wheel - - name: Build package - run: python -m build -sw - - name: Publish package + name: artifact + path: dist + + - name: Publish package to PyPI uses: pypa/gh-action-pypi-publish@ec4db0b4ddc65acdf4bff5fa45ac92d78b56bdf0 # v1.9.0 diff --git a/.gitignore b/.gitignore index ec552379a..4032fd610 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,4 @@ nbs/.last_checked .idea mlruns/ .luarc.json +*.so diff --git a/.gitmodules b/.gitmodules index e69de29bb..bac27244e 100644 --- a/.gitmodules +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "external_libs/eigen"] + path = external_libs/eigen + url = https://gitlab.com/libeigen/eigen diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 000000000..7ada6bd54 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,6 @@ +{ + "editor.formatOnSave": true, + "C_Cpp.formatting": "clangFormat", + "C_Cpp.clang_format_sortIncludes": true, + "C_Cpp.default.cppStandard": "c++17" +} diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index b645c8ee0..65eab7a9f 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -2,34 +2,33 @@ > This document contains instructions for collaborating on the different libraries of Nixtla. -Sometimes, diving into a new technology can be challenging and overwhelming. We've been there too, and we're more than ready to assist you with any issues you may encounter while following these steps. Don't hesitate to reach out to us on [Slack](https://join.slack.com/t/nixtlacommunity/shared_invite/zt-1pmhan9j5-F54XR20edHk0UtYAPcW4KQ). Just give fede a ping, and she'll be glad to assist you. +Sometimes, diving into a new technology can be challenging and overwhelming. We've been there too, and we're more than ready to assist you with any issues you may encounter while following these steps. Don't hesitate to reach out to us on [Slack](https://join.slack.com/t/nixtlacommunity/shared_invite/zt-1pmhan9j5-F54XR20edHk0UtYAPcW4KQ). Just give Azul a ping, and she'll be glad to assist you. ## Table of Contents 📚 1. [Prerequisites](#prerequisites) 2. [Git `fork-and-pull` worklow](#git-fork-and-pull-worklow) -3. [Set Up a Conda Environment](#set-up-a-conda-environment) +3. [Set Up a Virtual Environment](#set-up-a-virtual-environment) 4. [Install required libraries for development](#install-required-libraries-for-development) 5. [Start editable mode](#start-editable-mode) 6. [Set Up your Notebook based development environment](#set-up-your-notebook-based-development-environment) 7. [Start Coding](#start-coding) 8. [Example with Screen-shots](#example-with-screen-shots) -## Prerequisites +## Prerequisites -- *GitHub*: You should already have a GitHub account and a basic understanding of its functionalities. Alternatively check [this guide](https://docs.github.com/en/get-started). -- *Python*: Python should be installed on your system. Alternatively check [this guide](https://www.python.org/downloads/). -- *conda*: You need to have conda installed, along with a good grasp of fundamental operations such as creating environments, and activating and deactivating them. Alternatively check [this guide](https://conda.io/projects/conda/en/latest/user-guide/install/index.html). +- **GitHub**: You should already have a GitHub account and a basic understanding of its functionalities. Alternatively check [this guide](https://docs.github.com/en/get-started). +- **uv**: You need to have `uv` installed. You can refer to the [docs](https://docs.astral.sh/uv/getting-started/installation/) in order to install it. ## Git `fork-and-pull` worklow -**1. Fork the Project:** +**1. Fork the Project:** Start by forking the Nixtla repository to your own GitHub account. This creates a personal copy of the project where you can make changes without affecting the main repository. **2. Clone the Forked Repository** -Clone the forked repository to your local machine using `git clone https://github.com//nixtla.git`. This allows you to work with the code directly on your system. +Clone the forked repository to your local machine using `git clone --recursive https://github.com//statsforecast.git`. This allows you to work with the code directly on your system. -**3. Create a Branch:** +**3. Create a Branch:** Branching in GitHub is a key strategy for effectively managing and isolating changes to your project. It allows you to segregate work on different features, fixes, and issues without interfering with the main, production-ready codebase. @@ -43,33 +42,36 @@ Branching in GitHub is a key strategy for effectively managing and isolating cha After testing, branches are merged back into the main branch via a pull request, and then typically deleted to maintain a clean repository. You can read more about github and branching [here](https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/creating-and-deleting-branches-within-your-repository). -## Set Up a Conda Environment +## Set Up a Virtual Environment > If you want to use Docker or Codespaces, let us know opening an issue and we will set you up. -Next, you'll need to set up a [Conda](https://docs.conda.io/en/latest/) environment. Conda is an open-source package management and environment management system that runs on Windows, macOS, and Linux. It allows you to create separate environments containing files, packages, and dependencies that will not interact with each other. +Next, you'll need to create a virtual environment. uv is an open-source package management and environment management system that runs on Windows, macOS, and Linux. It allows you to create isolated environments with the dependencies required for a specific project. -First, ensure you have Anaconda or Miniconda installed on your system. Alternatively checkout these guides: [Anaconda](https://www.anaconda.com/), [Miniconda](https://docs.conda.io/en/latest/miniconda.html), and [Mamba](https://mamba.readthedocs.io/en/latest/). +First, ensure you have uv installed on your system. Alternatively, refer to the [installation docs](https://docs.astral.sh/uv/getting-started/installation/). -Then, you can create a new environment using `conda create -n nixtla-env python=3.10`. +Then, you can create a new environment using `uv venv`. This will use your default python interpreter, you can also define a specific python version (which will be downloaded if you don't have it already) by providing the `--python` flag. For example, `uv venv --python 3.12`. -You can also use mamba for creating the environment (mamba is faster than Conda) using `mamba create -n nixtla-env python=3.10`. - -Activate your new environment with `conda activate nixtla-env`. +Activate your new environment with `source .venv/bin/activate` for MacOS and Linux or `.\.venv\Scripts\activate` for Windows. ## Install required libraries for development -The `environment.yml` file contains all the dependencies required for the project. To install these dependencies, use the `mamba` package manager, which offers faster package installation and environment resolution than Conda. If you haven't installed `mamba` yet, you can do so using `conda install mamba -c conda-forge`. Run the following command to install the dependencies: +The `setup.py` file contains all the dependencies required for the project. To install these dependencies you can use `uv pip install -r setup.py --extra dev` + +### Setup pre-commit hooks + +We use [pre-commit](https://pre-commit.com/) to ease the development process, which run some checks before you make a commit to have a faster feedback loop. -``` -mamba env update -f environment.yml -``` +To setup the pre-commit hooks run: `pre-commit install` ## Start editable mode -Install the library in editable mode using `pip install -e ".[dev]"`. +Install the library in editable mode using `uv pip install --no-build-isolation -e .` (this requires a C++ compiler). -This means the package is linked directly to the source code, allowing any changes made to the source code to be immediately reflected in your Python environment without the need to reinstall the package. This is useful for testing changes during package development. +By using the `-e` flag the package is linked directly to the source code, allowing any changes made to the source code to be immediately reflected in your Python environment without the need to reinstall the package. This is useful for testing changes during package development. + +### Re-compiling the shared library +If you're working on the C++ code, you'll need to re-compile the shared library, which can be done with: `python setup.py build_ext --inplace` (this will compile it into the `build` directory and copy it to the python package location). ## Set Up your Notebook based development environment @@ -98,7 +100,7 @@ Open a jupyter notebook using `jupyter lab` (or VS Code). 2. **Commit Your Changes:** Add the changed files using `git add [your_modified_file_0.ipynb] [your_modified_file_1.ipynb]`, then commit these changes using `git commit -m ": "`. Please use [Conventional Commits](https://www.conventionalcommits.org/en/v1.0.0/) -3. **Push Your Changes:** +3. **Push Your Changes:** Push your changes to the remote repository on GitHub with `git push origin feature/your-feature-name`. 4. **Open a Pull Request:** Open a pull request from your new branch on the Nixtla repository on GitHub. Provide a thorough description of your changes when creating the pull request. @@ -107,7 +109,7 @@ Push your changes to the remote repository on GitHub with `git push origin featu Remember, contributing to open-source projects is a collaborative effort. Respect the work of others, welcome feedback, and always strive to improve. Happy coding! -> Nixtla offers the possibility of assisting with stipends for computing infrastructure for our contributors. If you are interested, please join our [slack](https://nixtlacommunity.slack.com/join/shared_invite/zt-1pmhan9j5-F54XR20edHk0UtYAPcW4KQ#/shared-invite/email) and write to fede or Max. +> Nixtla offers the possibility of assisting with stipends for computing infrastructure for our contributors. If you are interested, please join our [slack](https://nixtlacommunity.slack.com/join/shared_invite/zt-1pmhan9j5-F54XR20edHk0UtYAPcW4KQ#/shared-invite/email) and write to Azul or Max. You can find a detailed step by step buide with screen-shots below. @@ -125,7 +127,7 @@ In that repository, you can make your changes and then request to have them adde ### 2. Clone the repository -In this tutorial, we are using Mac (also compatible with other Linux distributions). If you are a collaborator of Nixtla, you can request an AWS instance to collaborate from there. If this is the case, please reach out to Max or Fede on [Slack](https://join.slack.com/t/nixtlacommunity/shared_invite/zt-1pmhan9j5-F54XR20edHk0UtYAPcW4KQ) to receive the appropriate access. We also use Visual Studio Code, which you can download from [here](https://code.visualstudio.com/download). +In this tutorial, we are using Mac (also compatible with other Linux distributions). If you are a collaborator of Nixtla, you can request an AWS instance to collaborate from there. If this is the case, please reach out to Max or Azul on [Slack](https://join.slack.com/t/nixtlacommunity/shared_invite/zt-1pmhan9j5-F54XR20edHk0UtYAPcW4KQ) to receive the appropriate access. We also use Visual Studio Code, which you can download from [here](https://code.visualstudio.com/download). Once the repository is created, you need to clone it to your own computer. Simply copy the repository URL from GitHub as shown below: @@ -147,42 +149,13 @@ You will end up with something like this: ![image](https://github.com/Nixtla/how-to-contribute-nixtlaverse/assets/10517170/ea4aed6f-2000-4ec8-a242-36b9dfd68d26) -### 3. Create the Conda environment - -Open a terminal within Visual Studio Code, as shown in the image: - -image - -You can use conda but we highly recommend using Mamba to speed up the creation of the Conda environment. To install it, simply use `conda install mamba -c conda-forge` in the terminal you just opened: - -![image](https://github.com/Nixtla/how-to-contribute-nixtlaverse/assets/10517170/08482b00-9434-46f0-9452-c3f4920eca6d) - -Create an empty environment named `mlforecast` with the following command: `mamba create -n mlforecast python=3.10`: -![image](https://github.com/Nixtla/how-to-contribute-nixtlaverse/assets/10517170/5e9032e8-3f5b-4a1c-93e7-3d390d5f73f1) - -Activate the newly created environment using `conda activate mlforecast`: - -![image](https://github.com/Nixtla/how-to-contribute-nixtlaverse/assets/10517170/803ae2b7-8369-4a24-9b7c-9326d52c13ef) - -Install the libraries within the environment file `environment.yml` using `mamba env update -f environment.yml`: - -![image](https://github.com/Nixtla/how-to-contribute-nixtlaverse/assets/10517170/e9672d58-b477-4963-9751-277c944a4d8a) - -Now install the library to make interactive changes and other additional dependencies using `pip install -e ".[dev]"`: - -![image](https://github.com/Nixtla/how-to-contribute-nixtlaverse/assets/10517170/501c8223-862d-40a9-8f2d-ecdaceaeaedb) - -Finally, setup the pre-commit hooks, which will run some checks when you commit your changes: `pre-commit install` - -### 4. Make the changes you want. +### 3. Make the changes you want. In this section, we assume that we want to increase the default number of windows used to create prediction intervals from 2 to 3. The first thing we need to do is create a specific branch for that change using `git checkout -b [new_branch]` like this: ![image](https://github.com/Nixtla/how-to-contribute-nixtlaverse/assets/10517170/7884f89a-ecc6-4200-8176-6a9b9f7c0aa2) -Once created, open the notebook you want to modify. In this case, it's `nbs/utils.ipynb`, which contains the metadata for the prediction intervals. After opening it, click on the environment you want to use (top right) and select the `mlforecast` environment: - -![image](https://github.com/Nixtla/how-to-contribute-nixtlaverse/assets/10517170/0a0a8285-9344-471e-b699-8bc13159e3a8) +Once created, open the notebook you want to modify. In this case, it's `nbs/utils.ipynb`, which contains the metadata for the prediction intervals. After opening it, click on the environment you want to use (top right) and select the `.venv` environment. Next, execute the notebook and make the necessary changes. In this case, we want to modify the `PredictionIntervals` class: @@ -209,7 +182,7 @@ Finally, push your changes using `git push`: ![image](https://github.com/Nixtla/how-to-contribute-nixtlaverse/assets/10517170/49c6851c-949b-4ca7-ac38-6b17ec103437) -### 5. Create a pull request. +### 4. Create a pull request. In GitHub, open your repository that contains your fork of the original repo. Once inside, you will see the changes you just pushed. Click on "Compare and pull request": diff --git a/MANIFEST.in b/MANIFEST.in index fc209e5a2..92bbd62a8 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -4,8 +4,10 @@ include README.md include LICENSE include settings.ini -recursive-include statsforecast * +recursive-include python/statsforecast * +recursive-include include/statsforecast * +recursive-include external_libs * ## Exclude List @@ -22,3 +24,4 @@ recursive-exclude docs * recursive-exclude examples * recursive-exclude experiments * recursive-exclude nbs * +recursive-exclude tests * diff --git a/action_files/test_dask.py b/action_files/test_dask.py index 5e6d94a8e..15749e07a 100644 --- a/action_files/test_dask.py +++ b/action_files/test_dask.py @@ -7,22 +7,34 @@ def to_distributed(df): return dd.from_pandas(df, npartitions=2) + @pytest.fixture() def sample_data(local_data): series, X_df = local_data return to_distributed(series), to_distributed(X_df) + def test_dask_flow(horizon, sample_data, n_series): - renamer = {'unique_id': 'id', 'ds': 'time', 'y': 'value'} + renamer = {"unique_id": "id", "ds": "time", "y": "value"} series, X_df = sample_data series = series.rename(columns=renamer) X_df = X_df.rename(columns=renamer) - pipeline(series, X_df, n_series, horizon, id_col='id', time_col='time', target_col='value') + pipeline( + series, + X_df, + n_series, + horizon, + id_col="id", + time_col="time", + target_col="value", + ) + def test_dask_flow_with_level(horizon, sample_data, n_series): pipeline_with_level(*sample_data, n_series, horizon) -@pytest.mark.parametrize('use_x', [True, False]) + +@pytest.mark.parametrize("use_x", [True, False]) def test_dask_flow_with_fitted(horizon, use_x, sample_data): series, X_df = sample_data if not use_x: diff --git a/action_files/test_spark.py b/action_files/test_spark.py index a0d2338e0..e067351e4 100644 --- a/action_files/test_spark.py +++ b/action_files/test_spark.py @@ -10,23 +10,28 @@ def spark(): return SparkSession.builder.getOrCreate() + def to_distributed(spark, df): - return spark.createDataFrame(df).repartition(2, 'unique_id') + return spark.createDataFrame(df).repartition(2, "unique_id") + @pytest.fixture() def sample_data(spark, local_data): series, X_df = local_data return to_distributed(spark, series), to_distributed(spark, X_df) + @pytest.mark.skipif(sys.version_info < (3, 8), reason="requires python >= 3.8") def test_spark_flow(horizon, sample_data, n_series): pipeline(*sample_data, n_series, horizon) + @pytest.mark.skipif(sys.version_info < (3, 8), reason="requires python >= 3.8") def test_spark_flow_with_level(horizon, sample_data, n_series): pipeline_with_level(*sample_data, n_series, horizon) -@pytest.mark.parametrize('use_x', [True, False]) + +@pytest.mark.parametrize("use_x", [True, False]) def test_spark_flow_with_fitted(horizon, use_x, sample_data): series, X_df = sample_data if not use_x: diff --git a/action_files/update_index.py b/action_files/update_index.py new file mode 100644 index 000000000..3afb872bc --- /dev/null +++ b/action_files/update_index.py @@ -0,0 +1,14 @@ +from pathlib import Path +from bs4 import BeautifulSoup + +index = Path("index.html").read_text() +soup = BeautifulSoup(index, "html.parser") +wheels = [f.name for f in Path("dist").glob("*.whl")] +existing_wheels = [a["href"] for a in soup.find_all("a")] +to_add = set(wheels) - set(existing_wheels) +for wheel in to_add: + new_link = soup.new_tag("a", href=wheel) + new_link.string = wheel + soup.body.append(new_link) + soup.body.append(soup.new_tag("br")) +Path("dist/index.html").write_text(soup.prettify()) diff --git a/external_libs/eigen b/external_libs/eigen new file mode 160000 index 000000000..9df21dc8b --- /dev/null +++ b/external_libs/eigen @@ -0,0 +1 @@ +Subproject commit 9df21dc8b4b576a7aa5c0094daa8d7e8b8be60f0 diff --git a/include/statsforecast/nelder_mead.h b/include/statsforecast/nelder_mead.h new file mode 100644 index 000000000..e79b41314 --- /dev/null +++ b/include/statsforecast/nelder_mead.h @@ -0,0 +1,135 @@ +#pragma once + +#include +#include + +namespace nm { +using Eigen::VectorXd; +using RowMajorMatrixXd = + Eigen::Matrix; + +auto Clamp(const VectorXd &x, const VectorXd &lower, const VectorXd &upper) { + return x.cwiseMax(lower).cwiseMin(upper); +} + +double StandardDeviation(const VectorXd &x) { + return std::sqrt((x.array() - x.mean()).square().mean()); +} + +Eigen::VectorX ArgSort(const VectorXd &v) { + Eigen::VectorX indices(v.size()); + std::iota(indices.begin(), indices.end(), 0); + std::sort(indices.begin(), indices.end(), + [&v](Eigen::Index i1, Eigen::Index i2) { return v(i1) < v(i2); }); + return indices; +} + +template +std::tuple NelderMead(Func F, const VectorXd &x, const VectorXd &lower, + const VectorXd upper, double init_step, + double zero_pert, double alpha, double gamma, + double rho, double sigma, int max_iter, double tol_std, + bool adaptive, Args &&...args) { + auto x0 = Clamp(x, lower, upper); + auto n = x0.size(); + if (adaptive) { + gamma = 1.0 + 2.0 / n; + rho = 0.75 - 1.0 / (2 * n); + sigma = 1.0 - 1.0 / n; + } + RowMajorMatrixXd simplex = x0.transpose().replicate(n + 1, 1); + // perturb simplex using init_step + for (Eigen::Index i = 0; i < n; ++i) { + double val = simplex(i, i); + if (val == 0) { + val = zero_pert; + } else { + val *= 1.0 + init_step; + } + simplex(i, i) = std::clamp(val, lower(i), upper(i)); + } + // array of the value of f + auto f_simplex = VectorXd(simplex.rows()); + for (Eigen::Index i = 0; i < simplex.rows(); ++i) { + f_simplex(i) = F(simplex.row(i), std::forward(args)...); + } + int i; + Eigen::Index best_idx = 0; + for (i = 0; i < max_iter; ++i) { + // check whether method should stop + if (StandardDeviation(f_simplex) < tol_std) { + break; + } + + // Step1: order of f_simplex + Eigen::VectorX order_f = ArgSort(f_simplex); + best_idx = order_f(0); + auto worst_idx = order_f(n); + auto second_worst_idx = order_f(n - 1); + + // calculate centroid as the col means removing the row with the max fval + VectorXd x_o = (simplex.colwise().sum() - simplex.row(worst_idx)) / n; + + // Step2: Reflection, Compute reflected point + VectorXd x_r = x_o + alpha * (x_o - simplex.row(worst_idx).transpose()); + x_r = Clamp(x_r, lower, upper); + double f_r = F(x_r, std::forward(args)...); + if (f_simplex(best_idx) <= f_r && f_r < f_simplex(second_worst_idx)) { + // accept reflection point + simplex(worst_idx, Eigen::all) = x_r; + f_simplex(worst_idx) = f_r; + continue; + } + + // Step3: Expansion, reflected point is the best point so far + if (f_r < f_simplex(best_idx)) { + VectorXd x_e = x_o + gamma * (x_r - x_o); + x_e = Clamp(x_e, lower, upper); + double f_e = F(x_e, std::forward(args)...); + if (f_e < f_r) { + // accept expansion point + simplex(worst_idx, Eigen::all) = x_e; + f_simplex(worst_idx) = f_e; + } else { + // accept reflection point + simplex(worst_idx, Eigen::all) = x_r; + f_simplex(worst_idx) = f_r; + } + continue; + } + + // Step4: Contraction + if (f_simplex(second_worst_idx) <= f_r && f_r < f_simplex(worst_idx)) { + VectorXd x_oc = x_o + rho * (x_r - x_o); + x_oc = Clamp(x_oc, lower, upper); + double f_oc = F(x_oc, std::forward(args)...); + if (f_oc <= f_r) { + // accept contraction point + simplex(worst_idx, Eigen::all) = x_oc; + f_simplex(worst_idx) = f_oc; + continue; + } + } else { + // Step5: Inside contraction + VectorXd x_ic = x_o - rho * (x_r - x_o); + x_ic = Clamp(x_ic, lower, upper); + double f_ic = F(x_ic, std::forward(args)...); + if (f_ic < f_simplex(worst_idx)) { + // accept inside contraction point + simplex(worst_idx, Eigen::all) = x_ic; + f_simplex(worst_idx) = f_ic; + continue; + } + } + + // Step6: Shrink + for (Eigen::Index j = 1; j < simplex.rows(); ++j) { + simplex.row(j) = simplex.row(best_idx) + + sigma * (simplex.row(j) - simplex.row(best_idx)); + simplex.row(j) = Clamp(simplex.row(j), lower, upper); + f_simplex(j) = F(simplex.row(j), std::forward(args)...); + } + } + return {simplex.row(best_idx), f_simplex(best_idx), i + 1}; +} +} // namespace nm diff --git a/nbs/docs/getting-started/0_Installation.ipynb b/nbs/docs/getting-started/0_Installation.ipynb index 62172ac14..12744cf1c 100644 --- a/nbs/docs/getting-started/0_Installation.ipynb +++ b/nbs/docs/getting-started/0_Installation.ipynb @@ -12,20 +12,19 @@ ] }, { - "attachments": {}, "cell_type": "markdown", - "id": "0f1d1483-6da7-4372-8390-84c9c280109e", + "id": "6b6943b9-4738-4fef-b3d1-1f18664b4a6b", "metadata": {}, "source": [ "You can install the *released version* of `StatsForecast` from the [Python package index](https://pypi.org) with:\n", "\n", - "```python\n", + "```shell\n", "pip install statsforecast\n", "```\n", "\n", "or \n", "\n", - "```python\n", + "```shell\n", "conda install -c conda-forge statsforecast\n", "``` \n", "\n", @@ -35,44 +34,17 @@ "\n", ":::{.callout-tip}\n", "We recommend installing your libraries inside a python virtual or [conda environment](https://docs.conda.io/projects/conda/en/latest/user-guide/install/macos.html).\n", - ":::\n", - "\n", - "#### User our env (optional)\n", - "\n", - "If you don't have a Conda environment and need tools like Numba, Pandas, NumPy, Jupyter, StatsModels, and Nbdev you can use ours by following these steps:\n", - "\n", - "1. Clone the StatsForecast repo: \n", - "\n", - "```bash \n", - "$ git clone https://github.com/Nixtla/statsforecast.git && cd statsforecast\n", - "```\n", - "\n", - "2. Create the environment using the `environment.yml` file: \n", - "\n", - "```bash \n", - "$ conda env create -f environment.yml\n", - "```\n", - "\n", - "3. Activate the environment:\n", - "```bash\n", - "$ conda activate statsforecast\n", - "```" - ] - }, - { - "cell_type": "markdown", - "id": "58909f00-5456-4398-80d6-72f8afbdd0e1", - "metadata": {}, - "source": [ - "## Extras" + ":::" ] }, { "attachments": {}, "cell_type": "markdown", - "id": "ecb9b7ac-36a3-48f9-a847-0fcba777d921", + "id": "02b9c8ff-5416-4197-b3ce-e416f27378b5", "metadata": {}, "source": [ + "#### Extras\n", + "\n", "The following features can also be installed by specifying the extra inside the install command, e.g. `pip install 'statsforecast[extra1,extra2]'`\n", "\n", "* **polars**: provide polars dataframes to StatsForecast.\n", @@ -81,6 +53,19 @@ "* **spark**: perform distributed forecasting with spark.\n", "* **ray**: perform distributed forecasting with ray." ] + }, + { + "cell_type": "markdown", + "id": "b5b913e0-8013-445e-bd4c-1a0e404b7d19", + "metadata": {}, + "source": [ + "#### Development version\n", + "\n", + "If you want to try out a new feature that hasn’t made it into a release yet you have the following options:\n", + "\n", + "* Install from our nightly wheels: `pip install --extra-index-url=http://nixtla-packages.s3-website.us-east-2.amazonaws.com --trusted-host nixtla-packages.s3-website.us-east-2.amazonaws.com statsforecast`\n", + "* Install from github: `pip install git+https://github.com/Nixtla/statsforecast`. This requires that you have a C++ compiler installed, so we encourage you to try the previous option first." + ] } ], "metadata": { diff --git a/nbs/src/ces.ipynb b/nbs/src/ces.ipynb index 5fa55071d..e8b1ba0ba 100644 --- a/nbs/src/ces.ipynb +++ b/nbs/src/ces.ipynb @@ -25,8 +25,7 @@ "from numba import njit\n", "from statsmodels.tsa.seasonal import seasonal_decompose\n", "\n", - "from statsforecast.ets import restrict_to_bounds, results\n", - "from statsforecast.utils import CACHE, NOGIL" + "from statsforecast.utils import CACHE, NOGIL, restrict_to_bounds, results" ] }, { diff --git a/nbs/src/ets.ipynb b/nbs/src/ets.ipynb index a82cf5d49..3199cc958 100644 --- a/nbs/src/ets.ipynb +++ b/nbs/src/ets.ipynb @@ -10,6 +10,18 @@ "#| default_exp ets" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "a0154a55-0a2a-48bb-b96d-8de1c02da878", + "metadata": {}, + "outputs": [], + "source": [ + "#| hide\n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, { "cell_type": "code", "execution_count": null, @@ -19,15 +31,12 @@ "source": [ "#| export\n", "import math\n", - "from collections import namedtuple\n", - "from typing import Tuple\n", "\n", "import numpy as np\n", - "from numba import njit\n", - "from numba.typed import List\n", "from statsmodels.tsa.seasonal import seasonal_decompose\n", "\n", - "from statsforecast.utils import _calculate_intervals, CACHE, NOGIL" + "from statsforecast._lib import ets as _ets\n", + "from statsforecast.utils import _calculate_intervals, results" ] }, { @@ -50,14 +59,6 @@ "# ETS Model" ] }, - { - "cell_type": "markdown", - "id": "03e43d2e-6341-4efd-81d6-b4ee8cf268e4", - "metadata": {}, - "source": [ - "## etscalc" - ] - }, { "cell_type": "code", "execution_count": null, @@ -65,211 +66,13 @@ "metadata": {}, "outputs": [], "source": [ - "#| export\n", - "# Global variables \n", - "NONE = 0\n", - "ADD = 1\n", - "MULT = 2\n", - "DAMPED = 1\n", - "TOL = 1.0e-10\n", - "HUGEN = 1.0e10\n", - "NA = -99999.0\n", - "smalno = np.finfo(float).eps\n", + "#| exporti\n", + "# Global variables\n", + "_smalno = np.finfo(float).eps\n", "_PHI_LOWER = 0.8\n", "_PHI_UPPER = 0.98" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "467c7ef6-267d-4d6a-9268-55a0c1ca24f6", - "metadata": {}, - "outputs": [], - "source": [ - "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def etscalc(y, n, x, m, \n", - " error, trend, season, \n", - " alpha, beta, \n", - " gamma, phi, e, \n", - " amse, nmse):\n", - " oldb = 0.\n", - " olds = np.zeros(max(24, m))\n", - " s = np.zeros(max(24, m))\n", - " f = np.zeros(30)\n", - " denom = np.zeros(30)\n", - " if m < 1:\n", - " m = 1 \n", - " if nmse > 30:\n", - " nmse = 30 \n", - " nstates = m * (season > NONE) + 1 + (trend > NONE) \n", - " #Copy initial state components \n", - " l = x[0]\n", - " if trend > NONE:\n", - " b = x[1]\n", - " else:\n", - " b = 0.\n", - " if season > NONE:\n", - " for j in range(m):\n", - " s[j] = x[(trend > NONE) + j + 1]\n", - " lik = 0.\n", - " lik2 = 0.\n", - " for j in range(nmse):\n", - " amse[j] = 0.\n", - " denom[j] = 0.\n", - " for i in range(n):\n", - " # Copy previous state\n", - " oldl = l \n", - " if trend > NONE:\n", - " oldb = b\n", - " if season > NONE:\n", - " for j in range(m):\n", - " olds[j] = s[j]\n", - " # one step forecast \n", - " forecast(oldl, oldb, olds, m, trend, season, phi, f, nmse)\n", - " if math.fabs(f[0] - NA) < TOL:\n", - " lik = NA\n", - " return lik\n", - " if error == ADD:\n", - " e[i] = y[i] - f[0]\n", - " else:\n", - " if math.fabs(f[0]) < TOL:\n", - " f_0 = f[0] + TOL\n", - " else:\n", - " f_0 = f[0]\n", - " e[i] = (y[i] - f[0]) / f_0\n", - " for j in range(nmse):\n", - " if (i + j) < n:\n", - " denom[j] += 1.\n", - " tmp = y[i + j] - f[j]\n", - " amse[j] = (amse[j] * (denom[j] - 1.0) + (tmp * tmp)) / denom[j]\n", - " # update state\n", - " l, b, s = update(oldl, l, oldb, b, olds, s, m, trend, season, alpha, beta, gamma, phi, y[i])\n", - " # store new state\n", - " x[nstates * (i + 1)] = l \n", - " if trend > NONE:\n", - " x[nstates * (i + 1) + 1] = b \n", - " if season > NONE:\n", - " for j in range(m):\n", - " x[nstates * (i + 1) + (trend > NONE) + j + 1] = s[j]\n", - " lik = lik + e[i] * e[i]\n", - " val = math.fabs(f[0])\n", - " if val > 0.:\n", - " lik2 += math.log(val)\n", - " else:\n", - " lik2 += math.log(val + 1e-8)\n", - " if lik > 0.:\n", - " lik = n * math.log(lik)\n", - " else:\n", - " lik = n * math.log(lik + 1e-8)\n", - " if error == MULT:\n", - " lik += 2 * lik2 \n", - " return lik" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e2d34a68-2641-4eeb-b025-17f07472bbdc", - "metadata": {}, - "outputs": [], - "source": [ - "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def forecast(l, b, s, m, \n", - " trend, season, phi, f, h):\n", - " #f is modified and it is mutable\n", - " phistar = phi \n", - " #forecasts\n", - " for i in range(h):\n", - " #print(phistar)\n", - " if trend == NONE:\n", - " f[i] = l\n", - " elif trend == ADD:\n", - " f[i] = l + phistar * b \n", - " elif b < 0:\n", - " f[i] = NA\n", - " else:\n", - " f[i] = l * math.pow(b, phistar)\n", - " j = m - 1 - i \n", - " while j < 0:\n", - " j += m\n", - " if season == ADD:\n", - " f[i] = f[i] + s[j]\n", - " elif season == MULT:\n", - " f[i] = f[i] * s[j]\n", - " if i < h - 1:\n", - " if math.fabs(phi - 1.0) < TOL:\n", - " phistar = phistar + 1.0 \n", - " else:\n", - " phistar = phistar + math.pow(phi, i + 1)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5a88b4c7-fad0-489f-b7c3-eb0e6db34f12", - "metadata": {}, - "outputs": [], - "source": [ - "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def update(oldl, l, oldb, b, \n", - " olds, s, \n", - " m, trend, season, \n", - " alpha, beta, gamma, \n", - " phi, y):\n", - " # New Level \n", - " if trend == NONE:\n", - " q = oldl # l(t - 1)\n", - " phib = 0 \n", - " elif trend == ADD:\n", - " phib = phi * oldb\n", - " q = oldl + phib #l(t - 1) + phi * b(t - 1)\n", - " elif math.fabs(phi - 1.0) < TOL:\n", - " phib = oldb \n", - " q = oldl * oldb #l(t - 1) * b(t - 1)\n", - " else:\n", - " phib = math.pow(oldb, phi)\n", - " q = oldl * phib #l(t - 1) * b(t - 1)^phi\n", - " # season\n", - " if season == NONE:\n", - " p = y \n", - " elif season == ADD:\n", - " p = y - olds[m - 1] #y[t] - s[t - m]\n", - " else:\n", - " if math.fabs(olds[m - 1]) < TOL:\n", - " p = HUGEN \n", - " else:\n", - " p = y / olds[m - 1] #y[t] / s[t - m]\n", - " l = q + alpha * (p - q)\n", - " # New Growth \n", - " if trend > NONE:\n", - " if trend == ADD:\n", - " r = l - oldl #l[t] - l[t-1]\n", - " else: #if(trend == MULT)\n", - " if math.fabs(oldl) < TOL:\n", - " r = HUGEN\n", - " else:\n", - " r = l / oldl #l[t] / l[t-1]\n", - " b = phib + (beta / alpha) * (r - phib) \n", - " # b[t] = phi*b[t-1] + beta*(r - phi*b[t-1])\n", - " # b[t] = b[t-1]^phi + beta*(r - b[t-1]^phi)\n", - " # New Seasonal\n", - " if season > NONE:\n", - " if season == ADD:\n", - " t = y - q \n", - " else: #if(season == MULT)\n", - " if math.fabs(q) < TOL:\n", - " t = HUGEN \n", - " else:\n", - " t = y / q \n", - " s[0] = olds[m - 1] + gamma * (t - olds[m - 1]) # s[t] = s[t - m] + gamma * (t - s[t - m])\n", - " for j in range(1, m):\n", - " s[j] = olds[j - 1] # s[t] = s[t]\n", - " return l, b, s" - ] - }, { "cell_type": "code", "execution_count": null, @@ -278,46 +81,79 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def etssimulate(x, m, error, trend, \n", - " season, alpha, beta, \n", - " gamma, phi, h, \n", - " y, e):\n", + "def etssimulate(\n", + " x: np.ndarray,\n", + " m: int,\n", + " error: _ets.Component,\n", + " trend: _ets.Component, \n", + " season: _ets.Component,\n", + " alpha: float,\n", + " beta: float, \n", + " gamma: float,\n", + " phi: float,\n", + " h: int, \n", + " y: np.ndarray,\n", + " e: np.ndarray\n", + ") -> None:\n", " oldb = 0.\n", " olds = np.zeros(24)\n", " s = np.zeros(24)\n", " f = np.zeros(10)\n", - " if m > 24 and season > NONE:\n", + " if m > 24 and season != _ets.Component.Nothing:\n", " return \n", " elif m < 1:\n", " m = 1 \n", - " #nstates = m * (season > NONE) + 1 + (trend > NONE)\n", " # Copy initial state components \n", " l = x[0]\n", - " if trend > NONE:\n", + " if trend != _ets.Component.Nothing:\n", " b = x[1]\n", - " if season > NONE:\n", + " if season != _ets.Component.Nothing:\n", " for j in range(m):\n", - " s[j] = x[(trend > NONE) + j + 1]\n", + " s[j] = x[(trend != _ets.Component.Nothing) + j + 1]\n", " for i in range(h):\n", " # Copy previous state\n", " oldl = l \n", - " if trend > NONE:\n", + " if trend != _ets.Component.Nothing:\n", " oldb = b \n", - " if season > NONE:\n", + " if season != _ets.Component.Nothing:\n", " for j in range(m):\n", " olds[j] = s[j]\n", " # one step forecast\n", - " forecast(oldl, oldb, olds, m, trend, season, phi, f, 1)\n", - " if math.fabs(f[0] - NA) < TOL:\n", - " y[0] = NA\n", + " _ets.forecast(\n", + " f,\n", + " oldl,\n", + " oldb,\n", + " olds,\n", + " m,\n", + " trend,\n", + " season,\n", + " phi,\n", + " 1,\n", + " )\n", + " if math.fabs(f[0] - _ets.NA) < _ets.TOL:\n", + " y[0] = _ets.NA\n", " return \n", - " if error == ADD:\n", + " if error == _ets.Component.Additive:\n", " y[i] = f[0] + e[i]\n", " else:\n", " y[i] = f[0] * (1.0 + e[i])\n", - " # Update state \n", - " l, b, s = update(oldl, l, oldb, b, olds, s, m, trend, season, alpha, beta, gamma, phi, y[i])" + " # Update state\n", + " l, b = _ets.update(\n", + " s,\n", + " l,\n", + " b,\n", + " oldl,\n", + " oldb,\n", + " olds,\n", + " m,\n", + " trend,\n", + " season,\n", + " alpha,\n", + " beta,\n", + " gamma,\n", + " phi,\n", + " y[i],\n", + " )" ] }, { @@ -328,29 +164,44 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def etsforecast(x, m, trend, season, \n", - " phi, h, f):\n", - " s = np.zeros(m)\n", + "def etsforecast(\n", + " x: np.ndarray,\n", + " m: int,\n", + " trend: _ets.Component,\n", + " season: _ets.Component, \n", + " phi: float,\n", + " h: int,\n", + " f: np.ndarray\n", + ") -> None:\n", + " s = np.zeros(m, dtype=np.float64)\n", " if m < 1:\n", " m = 1 \n", " # Copy initial state components\n", " l = x[0]\n", " b = 0.0\n", - " if trend > NONE:\n", + " has_trend = trend != _ets.Component.Nothing\n", + " if has_trend:\n", " b = x[1]\n", - " if season > NONE:\n", - " for j in range(m):\n", - " s[j] = x[(trend > NONE) + j + 1]\n", - "\n", + " if season != _ets.Component.Nothing:\n", + " s[:m] = x[has_trend + 1 : has_trend + 1 + m]\n", " # compute forecasts\n", - " forecast(l, b, s, m, trend, season, phi, f, h) " + " _ets.forecast(\n", + " f,\n", + " l,\n", + " b,\n", + " s,\n", + " m,\n", + " trend,\n", + " season,\n", + " phi,\n", + " h,\n", + " ) " ] }, { "cell_type": "code", "execution_count": null, - "id": "2585519c-6c14-45d3-9b28-ac4ef984d55e", + "id": "2bf5e405-d1bb-4473-bcd3-415f1f275fcb", "metadata": {}, "outputs": [], "source": [ @@ -363,27 +214,49 @@ "beta = 0.001766333 \n", "gamma = 0.\n", "phi = 0.\n", - "init_states = np.zeros(len(ap) * (2 + 1))\n", + "init_states = np.zeros(len(ap) * (2 + 1), dtype=np.float64)\n", "init_states[0] = 118.466667\n", "init_states[1] = 2.060606\n", - "etscalc(ap, len(ap), \n", - " init_states, 12, 1, 1, 0, \n", - " alpha, beta, gamma, phi, \n", - " e_, amse_, 3)" + "_ets.calc(\n", + " init_states,\n", + " e_,\n", + " amse_,\n", + " nmse_,\n", + " ap,\n", + " _ets.Component.Additive,\n", + " _ets.Component.Additive,\n", + " _ets.Component.Nothing,\n", + " alpha,\n", + " beta,\n", + " gamma,\n", + " phi,\n", + " 12,\n", + ")" ] }, { "cell_type": "code", "execution_count": null, - "id": "cf836cad-8071-48f8-bfe1-c0f716f2ab71", + "id": "b7dfb3f6-948d-46bd-b8c8-8633cdbcc0f5", "metadata": {}, "outputs": [], "source": [ "#| hide\n", - "etscalc(ap, len(ap), \n", - " init_states, 12, 0, 0, 0, \n", - " alpha, beta, gamma, phi, \n", - " e_, amse_, nmse_)" + "_ets.calc(\n", + " init_states,\n", + " e_,\n", + " amse_,\n", + " nmse_,\n", + " ap,\n", + " _ets.Component.Nothing,\n", + " _ets.Component.Nothing,\n", + " _ets.Component.Nothing, \n", + " alpha,\n", + " beta,\n", + " gamma,\n", + " phi, \n", + " 12,\n", + ")" ] }, { @@ -394,7 +267,6 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def initparam(alpha: float, beta: float, gamma: float, phi: float, \n", " trendtype: str, seasontype: str, \n", " damped: bool, \n", @@ -523,23 +395,6 @@ " return True" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "9395a3ef-430f-41bf-bb40-c75107ab53af", - "metadata": {}, - "outputs": [], - "source": [ - "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def sinpi(x):\n", - " return np.sin(np.pi * x)\n", - "\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def cospi(x):\n", - " return np.cos(np.pi * x)" - ] - }, { "cell_type": "code", "execution_count": null, @@ -548,7 +403,6 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def fourier(x, period, K, h=None):\n", " if h is None:\n", " times = np.arange(1, len(x) + 1)\n", @@ -556,10 +410,7 @@ " times = np.arange(len(x) + 1, len(x) + h + 1)\n", " # compute periods of all fourier terms\n", " # numba doesnt support list comprehension\n", - " len_p = 0\n", - " for k in K:\n", - " if k > 0:\n", - " len_p += k\n", + " len_p = sum(K)\n", " p = np.full(len_p, fill_value=np.nan)\n", " idx = 0\n", " for j, p_ in enumerate(period):\n", @@ -568,13 +419,13 @@ " idx += K[j]\n", " p = np.unique(p)\n", " # Remove columns where sinpi=0\n", - " k = np.abs(2 * p - np.round(2 * p, 0, np.empty_like(p))) > smalno\n", + " k = np.abs(2 * p - np.round(2 * p, 0, np.empty_like(p))) > _smalno\n", " # Compute matrix of fourier terms\n", " X = np.full((len(times), 2 * len(p)), fill_value=np.nan)\n", " for j in range(len(p)):\n", " if k[j]:\n", - " X[:, 2 * j - 1] = sinpi(2 * p[j] * times)\n", - " X[:, 2 * j] = cospi(2 * p[j] * times)\n", + " X[:, 2 * j - 1] = np.sin(2 * np.pi * p[j] * times)\n", + " X[:, 2 * j] = np.cos(2 * np.pi * p[j] * times)\n", " X = X[:, ~np.isnan(X.sum(axis=0))]\n", " return X" ] @@ -588,7 +439,7 @@ "source": [ "#| hide\n", "period = 12\n", - "fourier_terms = fourier(ap, List(x for x in [period]), List(x for x in [1]))" + "fourier_terms = fourier(ap, [period], [1])" ] }, { @@ -605,7 +456,7 @@ " if n < 4:\n", " raise ValueError(\"You've got to be joking (not enough data).\")\n", " elif n < 3 * m: #fit simple Fourier model\n", - " fouriery = fourier(y, List(x for x in [m]), List(x for x in [1]))\n", + " fouriery = fourier(y, [m], [1])\n", " X_fourier = np.full((n, 4), fill_value=np.nan)\n", " X_fourier[:, 0] = np.ones(n)\n", " X_fourier[:, 1] = np.arange(1, n + 1)\n", @@ -718,9 +569,14 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def switch(x: str):\n", - " return {'N': 0, 'A': 1, 'M': 2}[x]" + "def switch(x: str) -> _ets.Component:\n", + " if x == 'N':\n", + " return _ets.Component.Nothing\n", + " if x == 'A':\n", + " return _ets.Component.Additive\n", + " if x == 'M':\n", + " return _ets.Component.Multiplicative\n", + " raise ValueError(f'Unknown component {x}')" ] }, { @@ -734,6 +590,39 @@ "switch('A')" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "211d683b-32cd-4837-bfbf-19f539e899ad", + "metadata": {}, + "outputs": [], + "source": [ + "#| exporti\n", + "def switch_criterion(x: str) -> _ets.Criterion:\n", + " if x == 'lik':\n", + " return _ets.Criterion.Likelihood\n", + " if x == 'mse':\n", + " return _ets.Criterion.MSE\n", + " if x == 'amse':\n", + " return _ets.Criterion.AMSE\n", + " if x == 'sigma':\n", + " return _ets.Criterion.Sigma\n", + " if x == 'mae':\n", + " return _ets.Criterion.MAE\n", + " raise ValueError(f'Unknown crtierion {x}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57b3d637-f195-40f8-b65c-0b620d184fc7", + "metadata": {}, + "outputs": [], + "source": [ + "#| hide\n", + "switch_criterion('sigma')" + ] + }, { "cell_type": "code", "execution_count": null, @@ -742,7 +631,6 @@ "outputs": [], "source": [ "#| exporti\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", "def pegelsresid_C(y: np.ndarray, \n", " m: int, \n", " init_state: np.ndarray, \n", @@ -764,12 +652,21 @@ " if seasontype == 'N':\n", " gamma = 0.\n", " amse = np.full(nmse, fill_value=np.nan)\n", - " lik = etscalc(y=y, n=n, x=x, m=m, \n", - " error=switch(errortype),\n", - " trend=switch(trendtype), \n", - " season=switch(seasontype),\n", - " alpha=alpha, beta=beta, gamma=gamma, phi=phi,\n", - " e=e, amse=amse, nmse=nmse)\n", + " lik = _ets.calc(\n", + " x,\n", + " e,\n", + " amse,\n", + " nmse,\n", + " y,\n", + " switch(errortype),\n", + " switch(trendtype),\n", + " switch(seasontype),\n", + " alpha,\n", + " beta,\n", + " gamma,\n", + " phi,\n", + " m,\n", + " )\n", " x = x.reshape((n + 1, p))\n", " if not np.isnan(lik):\n", " if np.abs(lik + 99999) < 1e-7:\n", @@ -777,230 +674,6 @@ " return amse, e, x, lik" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "61ad54e9-0510-42f4-b34a-f947d6823c8c", - "metadata": {}, - "outputs": [], - "source": [ - "#| exporti\n", - "results = namedtuple('results', 'x fn nit simplex')\n", - "\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def restrict_to_bounds(x, lower, upper):\n", - " new_x = np.full_like(x, fill_value=np.nan, dtype=x.dtype)\n", - " for i in range(x.size):\n", - " lo = lower[i]\n", - " up = upper[i]\n", - " if x[i] < lo:\n", - " new_x[i] = lo\n", - " elif x[i] > up:\n", - " new_x[i] = up\n", - " else:\n", - " new_x[i] = x[i]\n", - " return new_x\n", - " \n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def nelder_mead_ets(\n", - " x0: np.ndarray, \n", - " args: Tuple = (), \n", - " lower: np.ndarray = np.empty(0), \n", - " upper: np.ndarray = np.empty(0), \n", - " init_step: float = 0.05,\n", - " zero_pert: float = 0.0001,\n", - " alpha: float = 1.,\n", - " gamma: float = 2.,\n", - " rho: float = 0.5,\n", - " sigma: float = 0.5,\n", - " max_iter: int = 2_000,\n", - " tol_std: float = 1e-10,\n", - " adaptive: bool = False,\n", - " ):\n", - " #We are trying to minimize the function fn(x, args)\n", - " #with initial point x0.\n", - " #Step 0:\n", - " # get x1, ..., x_{n+1}\n", - " # the original article suggested a simplex where an initial point is given as x0\n", - " # with the others generated with a fixed step along each dimension in turn.\n", - " bounds = len(lower) and len(upper)\n", - " if bounds:\n", - " x0 = restrict_to_bounds(x0, lower, upper)\n", - " \n", - " n = x0.size\n", - " if adaptive:\n", - " gamma = 1. + 2. / n\n", - " rho = 0.75 - 1. / (2. * n)\n", - " sigma = 1. - 1. / n\n", - " simplex = np.full((n + 1, n), fill_value=np.nan, dtype=np.float64) #each row is x_j\n", - " simplex[:] = x0\n", - " # perturb simplex using `init_step`\n", - " diag = np.copy(np.diag(simplex))\n", - " diag[diag == 0.] = zero_pert\n", - " diag[diag != 0.] *= (1 + init_step)\n", - " np.fill_diagonal(simplex, diag)\n", - " # restrict simplex to bounds if passed\n", - " if bounds:\n", - " for j in range(n + 1):\n", - " simplex[j] = restrict_to_bounds(simplex[j], lower, upper)\n", - " # array of the value of f\n", - " f_simplex = np.full(n + 1, fill_value=np.nan)\n", - " for j in range(n + 1):\n", - " f_simplex[j] = ets_target_fn(simplex[j], *args)\n", - " for it in range(max_iter):\n", - " #Step1: order of f_simplex\n", - " #print(simplex)\n", - " #print(f_simplex)\n", - " order_f = f_simplex.argsort()\n", - " best_idx = order_f[0]\n", - " worst_idx = order_f[-1]\n", - " second_worst_idx = order_f[-2]\n", - " #Check whether method should stop.\n", - " if np.std(f_simplex) < tol_std:\n", - " break\n", - " #calculate centroid except argmax f_simplex\n", - " x_o = simplex[np.delete(order_f, -1)].sum(axis=0) / n\n", - " #Step2: Reflection, Compute reflected point\n", - " x_r = x_o + alpha * (x_o - simplex[worst_idx])\n", - " # restrict x_r to bounds if passed\n", - " if bounds:\n", - " x_r = restrict_to_bounds(x_r, lower, upper)\n", - " f_r = ets_target_fn(x_r, *args)\n", - " if f_simplex[best_idx] <= f_r < f_simplex[second_worst_idx]:\n", - " simplex[worst_idx] = x_r\n", - " f_simplex[worst_idx] = f_r\n", - " continue\n", - " #Step3: Expansion, reflected point is the best point so far\n", - " if f_r < f_simplex[best_idx]:\n", - " x_e = x_o + gamma * (x_r - x_o)\n", - " # restrict x_e to bounds if passed\n", - " if bounds:\n", - " x_e = restrict_to_bounds(x_e, lower, upper)\n", - " f_e = ets_target_fn(x_e, *args)\n", - " if f_e < f_r:\n", - " simplex[worst_idx] = x_e\n", - " f_simplex[worst_idx] = f_e\n", - " else:\n", - " simplex[worst_idx] = x_r\n", - " f_simplex[worst_idx] = f_r\n", - " continue\n", - " #Step4: outside Contraction\n", - " if f_simplex[second_worst_idx] <= f_r < f_simplex[worst_idx]:\n", - " x_oc = x_o + rho * (x_r - x_o)\n", - " if bounds:\n", - " x_oc = restrict_to_bounds(x_oc, lower, upper)\n", - " f_oc = ets_target_fn(x_oc, *args)\n", - " if f_oc <= f_r:\n", - " simplex[worst_idx] = x_oc\n", - " f_simplex[worst_idx] = f_oc\n", - " continue\n", - " #step 5 inside contraction\n", - " else:\n", - " x_ic = x_o - rho * (x_r - x_o)\n", - " # restrict x_c to bounds if passed\n", - " if bounds:\n", - " x_ic = restrict_to_bounds(x_ic, lower, upper)\n", - " f_ic = ets_target_fn(x_ic, *args)\n", - " if f_ic < f_simplex[worst_idx]:\n", - " simplex[worst_idx] = x_ic\n", - " f_simplex[worst_idx] = f_ic\n", - " continue\n", - " #step 6: shrink\n", - " simplex[np.delete(order_f, 0)] = simplex[best_idx] + sigma * (simplex[np.delete(order_f, 0)] - simplex[best_idx])\n", - " for i in np.delete(order_f, 0):\n", - " simplex[i] = restrict_to_bounds(simplex[i], lower, upper)\n", - " f_simplex[i] = ets_target_fn(simplex[i], *args)\n", - " return results(simplex[best_idx], f_simplex[best_idx], it + 1, simplex)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6ffc4052-c1d0-4350-a55d-1482c164907c", - "metadata": {}, - "outputs": [], - "source": [ - "#| export\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", - "def ets_target_fn(\n", - " par,\n", - " p_y, p_nstate, \n", - " p_errortype, p_trendtype, p_seasontype, p_damped, \n", - " p_lower, p_upper, p_opt_crit, p_nmse, p_bounds, p_m, \n", - " p_optAlpha, p_optBeta, p_optGamma, p_optPhi, \n", - " p_givenAlpha, p_givenBeta, p_givenGamma, p_givenPhi, \n", - " alpha, beta, gamma, phi\n", - " ):\n", - " #p_par_length = len(par)\n", - " j = 0 \n", - " if p_optAlpha: \n", - " alpha = par[j]\n", - " j += 1\n", - " if p_optBeta:\n", - " beta = par[j]\n", - " j += 1\n", - " if p_optGamma:\n", - " gamma = par[j]\n", - " j += 1\n", - " if p_optPhi:\n", - " phi = par[j]\n", - " #size of state\n", - " p = p_nstate + (1 if p_seasontype != 0 else 0)\n", - " n = len(p_y)\n", - " state = np.zeros(p * (n + 1))\n", - " state[:p_nstate] = par[-p_nstate:]\n", - " # Add extra state\n", - " if p_seasontype != 0:\n", - " start = 1 + (1 if p_trendtype != 0 else 0)\n", - " sum_ = state[start:p_nstate].sum()\n", - " new_state = p_m * (1 if p_seasontype == 2 else 0) - sum_\n", - " state[p_nstate] = new_state\n", - " # check states \n", - " if p_seasontype == 2 and state[start:].min() < 0:\n", - " return np.inf\n", - " amse = np.zeros(30)\n", - " e = np.zeros(n)\n", - " lik = etscalc(\n", - " #we dont need to change these param\n", - " p_y, n, state, p_m, \n", - " p_errortype, p_trendtype, p_seasontype, \n", - " #inmutable, wee need to change them?\n", - " alpha, beta, gamma, phi, \n", - " #mutable\n", - " e, amse, \n", - " #inmutable\n", - " p_nmse\n", - " )\n", - " if lik < -1e10: \n", - " lik = -1e10 \n", - " if math.isnan(lik): \n", - " lik = -np.inf\n", - " if math.fabs(lik + 99999) < 1e-7: \n", - " lik = -np.inf\n", - " if p_opt_crit == \"lik\":\n", - " objval = lik \n", - " elif p_opt_crit == \"mse\": \n", - " objval = amse[0]\n", - " elif p_opt_crit == \"amse\":\n", - " mean = 0.\n", - " for i in range(p_nmse):\n", - " mean += amse[i] / p_nmse \n", - " objval = mean \n", - " elif p_opt_crit == \"sigma\":\n", - " mean = 0.\n", - " ne = len(e)\n", - " for i in range(ne):\n", - " mean += e[i] * e[i] / ne \n", - " objval = mean \n", - " elif p_opt_crit == \"mae\":\n", - " mean = 0\n", - " ne = len(e)\n", - " for i in range(ne):\n", - " mean += math.fabs(e[i]) / ne \n", - " objval = mean\n", - " return objval" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1043,46 +716,46 @@ " optGamma = not np.isnan(gamma)\n", " optPhi = not np.isnan(phi)\n", " \n", - " givenAlpha = False\n", - " givenBeta = False\n", - " givenGamma = False\n", - " givenPhi = False\n", - " \n", " if not np.isnan(par_noopt['alpha']):\n", " optAlpha = False\n", - " givenAlpha = True\n", " if not np.isnan(par_noopt['beta']):\n", " optBeta = False\n", - " givenBeta = True\n", " if not np.isnan(par_noopt['gamma']):\n", " optGamma = False\n", - " givenGamma = True\n", " if not np.isnan(par_noopt['phi']):\n", " optPhi = False\n", - " givenPhi = True\n", - " \n", + "\n", " if not damped:\n", " phi = 1.\n", " if trendtype == 'N':\n", " beta = 0.\n", " if seasontype == 'N':\n", " gamma = 0.\n", - " res = nelder_mead_ets(\n", - " x0, \n", - " args=(\n", - " y, nstate, switch(errortype), switch(trendtype), switch(seasontype),\n", - " damped, lowerb, upperb, opt_crit, nmse, bounds, m, \n", - " optAlpha, optBeta, optGamma, optPhi, \n", - " givenAlpha, givenBeta, givenGamma, givenPhi,\n", - " alpha, beta, gamma, phi\n", - " ),\n", - " lower=lowerb,\n", - " upper=upperb,\n", - " tol_std=1e-4, \n", - " max_iter=1_000,\n", - " adaptive=True,\n", + " opt_res = _ets.optimize(\n", + " x0,\n", + " y,\n", + " nstate,\n", + " switch(errortype),\n", + " switch(trendtype),\n", + " switch(seasontype),\n", + " switch_criterion(opt_crit),\n", + " nmse,\n", + " m,\n", + " optAlpha,\n", + " optBeta,\n", + " optGamma,\n", + " optPhi,\n", + " alpha,\n", + " beta,\n", + " gamma,\n", + " phi,\n", + " lowerb,\n", + " upperb,\n", + " 1e-4,\n", + " 1_000,\n", + " True,\n", " )\n", - " return res" + " return results(*opt_res, None)" ] }, { @@ -1144,7 +817,6 @@ " if np_ >= len(y) - 1:\n", " return dict(aic=np.inf, bic=np.inf, aicc=np.inf, mse=np.inf,\n", " amse=np.inf, fit=None, par=par, states=init_state)\n", - " \n", " fred = optimize_ets_target_fn(\n", " x0=par, par=par_, y=y, nstate=nstate, \n", " errortype=errortype, trendtype=trendtype,\n", @@ -1240,8 +912,7 @@ "metadata": {}, "outputs": [], "source": [ - "#| export\n", - "@njit(nogil=NOGIL, cache=CACHE)\n", + "#| exporti\n", "def is_constant(x):\n", " return np.all(x[0] == x)" ] @@ -1264,7 +935,7 @@ "metadata": {}, "outputs": [], "source": [ - "#| exporti\n", + "#| export\n", "def ets_f(y, m, model='ZZZ', \n", " damped=None, alpha=None, beta=None, gamma=None, phi=None,\n", " additive_only=None, blambda=None, biasadj=None, \n", @@ -1273,8 +944,7 @@ " ic='aicc', restrict=True, allow_multiplicative_trend=False,\n", " use_initial_values=False, \n", " maxit=2_000):\n", - " if y.dtype not in (np.float32, np.float64):\n", - " y = y.astype(np.float32)\n", + " y = y.astype(np.float64, copy=False)\n", " # converting params to floats \n", " # to improve numba compilation\n", " if alpha is None:\n", @@ -1471,24 +1141,15 @@ "outputs": [], "source": [ "#| exporti\n", - "#@njit(nogil=NOGIL, cache=CACHE)\n", "def _compute_sigmah(pf, h, sigma, cvals):\n", - " \n", " theta = np.full(h, np.nan)\n", " theta[0] = pf[0]**2\n", " \n", " for k in range(1,h): \n", - " sum_val = 0 \n", - " for j in range(1,k+1): \n", - " val = cvals[j-1]**2*theta[k-j]\n", - " sum_val = sum_val+val \n", - " theta[k] = pf[k]**2+sigma*sum_val\n", - " \n", - " sigmah = np.full(h, np.nan) \n", - " for k in range(0,h): \n", - " sigmah[k] = (1+sigma)*theta[k]-pf[k]**2\n", - " \n", - " return sigmah " + " sum_val = cvals[:k]**2 @ theta[:k][::-1] \n", + " theta[k] = pf[k]**2 + sigma * sum_val\n", + "\n", + " return (1 + sigma) * theta - pf**2" ] }, { @@ -1541,9 +1202,9 @@ " mu = np.zeros(h)\n", " var = np.zeros(h) \n", "\n", - " for i in range(0,h): \n", - " mu[i] = np.matmul(H1, np.matmul(Mh, np.transpose(H2)))\n", - " var[i] = (1+sigma)*np.matmul(H21, np.matmul(Vh, np.transpose(H21)))+sigma*mu[i]**2\n", + " for i in range(0,h):\n", + " mu[i] = (H1 @ (Mh @ H2.T)).item()\n", + " var[i] = (1 + sigma) * (H21 @ (Vh @ H21.T)).item() + sigma * mu[i] ** 2\n", " vecMh = Mh.flatten()\n", " exp1 = np.matmul(F21, np.matmul(Vh, np.transpose(F21)))\n", " exp2 = np.matmul(F21, np.matmul(Vh, np.transpose(G21)))\n", @@ -1746,7 +1407,7 @@ "outputs": [], "source": [ "#| hide\n", - "import matplotlib.pyplot as plt\n" + "import matplotlib.pyplot as plt" ] }, { @@ -1959,7 +1620,7 @@ "outputs": [], "source": [ "#| hide\n", - "test_series = np.array([ \n", + "test_series = np.array([\n", " 7., 17., 7., 38., 89., 269., 33., 13., 17., 10., 21.,\n", " 61., 18., 35., 19., 23., 75., 188., 62., 31., 39., 10.,\n", " 26., 16., 35., 40., 26., 43., 72., 301., 69., 24.\n", @@ -2038,7 +1699,26 @@ "test_data = sales[train_size:]\n", "\n", "mod = ets_f(train_data, m=52, model='ZAZ')\n", - "assert not np.isnan(forecast_ets(mod, len(test_data))['mean']).any()" + "assert not np.isnan(forecast_ets(mod, len(test_data), level=[80])['mean']).any()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5114a5f9-bd89-4ca5-96ea-160f4f388b6f", + "metadata": {}, + "outputs": [], + "source": [ + "#| hide\n", + "# test all model combinations\n", + "for error in ['M', 'A', 'Z']:\n", + " for trend in ['N', 'M', 'A', 'Z']:\n", + " for seasonal in ['N', 'M', 'A', 'Z']:\n", + " model = f'{error}{trend}{seasonal}'\n", + " mod = ets_f(train_data, m=52, model=model, restrict=False)\n", + " forecasts = forecast_ets(mod, len(test_data), level=[80])\n", + " mape = np.abs(forecasts['mean'] / test_data - 1).mean()\n", + " assert mape < 0.3" ] } ], diff --git a/nbs/src/theta.ipynb b/nbs/src/theta.ipynb index d26ea83f9..f9bc07259 100644 --- a/nbs/src/theta.ipynb +++ b/nbs/src/theta.ipynb @@ -27,8 +27,14 @@ "from statsmodels.tsa.seasonal import seasonal_decompose\n", "from statsmodels.tsa.stattools import acf\n", "\n", - "from statsforecast.ets import restrict_to_bounds, results\n", - "from statsforecast.utils import _seasonal_naive, _repeat_val_seas, CACHE, NOGIL" + "from statsforecast.utils import (\n", + " CACHE,\n", + " NOGIL,\n", + " _repeat_val_seas,\n", + " _seasonal_naive,\n", + " restrict_to_bounds,\n", + " results,\n", + ")" ] }, { diff --git a/nbs/src/utils.ipynb b/nbs/src/utils.ipynb index c933ded36..417156136 100644 --- a/nbs/src/utils.ipynb +++ b/nbs/src/utils.ipynb @@ -41,11 +41,13 @@ "import math\n", "import os\n", "import warnings\n", + "from collections import namedtuple\n", "from functools import wraps\n", "from typing import Dict\n", "\n", "import numpy as np\n", "import pandas as pd\n", + "from numba import njit\n", "from scipy.stats import norm\n", "\n", "from utilsforecast.compat import DataFrame\n", @@ -68,7 +70,30 @@ " 'Please set that one instead.',\n", " FutureWarning,\n", " )\n", - "CACHE = bool(os.getenv('NIXTLA_NUMBA_CACHE', '')) or LEGACY_CACHE" + "CACHE = bool(os.getenv('NIXTLA_NUMBA_CACHE', '')) or LEGACY_CACHE\n", + "results = namedtuple(\"results\", \"x fn nit simplex\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#| exporti\n", + "@njit(nogil=NOGIL, cache=CACHE)\n", + "def restrict_to_bounds(x, lower, upper):\n", + " new_x = np.full_like(x, fill_value=np.nan, dtype=x.dtype)\n", + " for i in range(x.size):\n", + " lo = lower[i]\n", + " up = upper[i]\n", + " if x[i] < lo:\n", + " new_x[i] = lo\n", + " elif x[i] > up:\n", + " new_x[i] = up\n", + " else:\n", + " new_x[i] = x[i]\n", + " return new_x" ] }, { diff --git a/pyproject.toml b/pyproject.toml index af5c7b40e..375c37033 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,7 +1,16 @@ +[build-system] +requires = ["setuptools>=42", "pybind11>=2.10"] +build-backend = "setuptools.build_meta" + +[tool.cibuildwheel] +archs = "all" +build-verbosity = 1 +test-skip = "*" + [tool.ruff] target-version = "py38" line-length = 88 lint.select = ["F"] [tool.mypy] -ignore_missing_imports = true \ No newline at end of file +ignore_missing_imports = true diff --git a/statsforecast/__init__.py b/python/statsforecast/__init__.py similarity index 79% rename from statsforecast/__init__.py rename to python/statsforecast/__init__.py index 7c4ca7076..3ec14b09c 100644 --- a/statsforecast/__init__.py +++ b/python/statsforecast/__init__.py @@ -1,4 +1,4 @@ -__version__ = "1.7.6" +__version__ = "1.7.6.99" __all__ = ["StatsForecast"] from .core import StatsForecast from .distributed import fugue # noqa diff --git a/statsforecast/_modidx.py b/python/statsforecast/_modidx.py similarity index 99% rename from statsforecast/_modidx.py rename to python/statsforecast/_modidx.py index fbefefec3..246699083 100644 --- a/statsforecast/_modidx.py +++ b/python/statsforecast/_modidx.py @@ -4,7 +4,7 @@ 'doc_baseurl': '/statsforecast/', 'doc_host': 'https://Nixtla.github.io', 'git_url': 'https://github.com/Nixtla/statsforecast/', - 'lib_path': 'statsforecast'}, + 'lib_path': 'python/statsforecast'}, 'syms': { 'statsforecast.adapters.prophet': { 'statsforecast.adapters.prophet.AutoARIMAProphet': ( 'src/adapters.prophet.html#autoarimaprophet', 'statsforecast/adapters/prophet.py'), 'statsforecast.adapters.prophet.AutoARIMAProphet.__init__': ( 'src/adapters.prophet.html#autoarimaprophet.__init__', @@ -250,29 +250,22 @@ 'statsforecast.ets._compute_sigmah': ('src/ets.html#_compute_sigmah', 'statsforecast/ets.py'), 'statsforecast.ets.admissible': ('src/ets.html#admissible', 'statsforecast/ets.py'), 'statsforecast.ets.check_param': ('src/ets.html#check_param', 'statsforecast/ets.py'), - 'statsforecast.ets.cospi': ('src/ets.html#cospi', 'statsforecast/ets.py'), 'statsforecast.ets.ets_f': ('src/ets.html#ets_f', 'statsforecast/ets.py'), - 'statsforecast.ets.ets_target_fn': ('src/ets.html#ets_target_fn', 'statsforecast/ets.py'), - 'statsforecast.ets.etscalc': ('src/ets.html#etscalc', 'statsforecast/ets.py'), 'statsforecast.ets.etsforecast': ('src/ets.html#etsforecast', 'statsforecast/ets.py'), 'statsforecast.ets.etsmodel': ('src/ets.html#etsmodel', 'statsforecast/ets.py'), 'statsforecast.ets.etssimulate': ('src/ets.html#etssimulate', 'statsforecast/ets.py'), - 'statsforecast.ets.forecast': ('src/ets.html#forecast', 'statsforecast/ets.py'), 'statsforecast.ets.forecast_ets': ('src/ets.html#forecast_ets', 'statsforecast/ets.py'), 'statsforecast.ets.forward_ets': ('src/ets.html#forward_ets', 'statsforecast/ets.py'), 'statsforecast.ets.fourier': ('src/ets.html#fourier', 'statsforecast/ets.py'), 'statsforecast.ets.initparam': ('src/ets.html#initparam', 'statsforecast/ets.py'), 'statsforecast.ets.initstate': ('src/ets.html#initstate', 'statsforecast/ets.py'), 'statsforecast.ets.is_constant': ('src/ets.html#is_constant', 'statsforecast/ets.py'), - 'statsforecast.ets.nelder_mead_ets': ('src/ets.html#nelder_mead_ets', 'statsforecast/ets.py'), 'statsforecast.ets.optimize_ets_target_fn': ( 'src/ets.html#optimize_ets_target_fn', 'statsforecast/ets.py'), 'statsforecast.ets.pegelsfcast_C': ('src/ets.html#pegelsfcast_c', 'statsforecast/ets.py'), 'statsforecast.ets.pegelsresid_C': ('src/ets.html#pegelsresid_c', 'statsforecast/ets.py'), - 'statsforecast.ets.restrict_to_bounds': ('src/ets.html#restrict_to_bounds', 'statsforecast/ets.py'), - 'statsforecast.ets.sinpi': ('src/ets.html#sinpi', 'statsforecast/ets.py'), 'statsforecast.ets.switch': ('src/ets.html#switch', 'statsforecast/ets.py'), - 'statsforecast.ets.update': ('src/ets.html#update', 'statsforecast/ets.py')}, + 'statsforecast.ets.switch_criterion': ('src/ets.html#switch_criterion', 'statsforecast/ets.py')}, 'statsforecast.feature_engineering': { 'statsforecast.feature_engineering.mstl_decomposition': ( 'src/feature_engineering.html#mstl_decomposition', 'statsforecast/feature_engineering.py')}, 'statsforecast.garch': { 'statsforecast.garch.garch_cons': ('src/garch.html#garch_cons', 'statsforecast/garch.py'), @@ -814,4 +807,6 @@ 'statsforecast.utils._repeat_val': ('src/utils.html#_repeat_val', 'statsforecast/utils.py'), 'statsforecast.utils._repeat_val_seas': ('src/utils.html#_repeat_val_seas', 'statsforecast/utils.py'), 'statsforecast.utils._seasonal_naive': ('src/utils.html#_seasonal_naive', 'statsforecast/utils.py'), - 'statsforecast.utils.generate_series': ('src/utils.html#generate_series', 'statsforecast/utils.py')}}} + 'statsforecast.utils.generate_series': ('src/utils.html#generate_series', 'statsforecast/utils.py'), + 'statsforecast.utils.restrict_to_bounds': ( 'src/utils.html#restrict_to_bounds', + 'statsforecast/utils.py')}}} diff --git a/statsforecast/adapters/__init__.py b/python/statsforecast/adapters/__init__.py similarity index 100% rename from statsforecast/adapters/__init__.py rename to python/statsforecast/adapters/__init__.py diff --git a/statsforecast/adapters/prophet.py b/python/statsforecast/adapters/prophet.py similarity index 98% rename from statsforecast/adapters/prophet.py rename to python/statsforecast/adapters/prophet.py index dbac58460..ce3d3a667 100644 --- a/statsforecast/adapters/prophet.py +++ b/python/statsforecast/adapters/prophet.py @@ -1,9 +1,9 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/adapters.prophet.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../../nbs/src/adapters.prophet.ipynb. # %% auto 0 __all__ = ['AutoARIMAProphet'] -# %% ../../nbs/src/adapters.prophet.ipynb 4 +# %% ../../../nbs/src/adapters.prophet.ipynb 4 import sys import pandas as pd @@ -32,7 +32,7 @@ ) raise ModuleNotFoundError(msg) from e -# %% ../../nbs/src/adapters.prophet.ipynb 7 +# %% ../../../nbs/src/adapters.prophet.ipynb 7 class AutoARIMAProphet(Prophet): """AutoARIMAProphet adapter. diff --git a/statsforecast/arima.py b/python/statsforecast/arima.py similarity index 98% rename from statsforecast/arima.py rename to python/statsforecast/arima.py index 7341f5f09..42c9b5255 100644 --- a/statsforecast/arima.py +++ b/python/statsforecast/arima.py @@ -1,10 +1,10 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/arima.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/arima.ipynb. # %% auto 0 __all__ = ['predict_arima', 'arima_string', 'forecast_arima', 'fitted_arima', 'auto_arima_f', 'print_statsforecast_ARIMA', 'ARIMASummary', 'AutoARIMA'] -# %% ../nbs/src/arima.ipynb 4 +# %% ../../nbs/src/arima.ipynb 4 import math import warnings from collections import namedtuple @@ -21,10 +21,10 @@ from .mstl import mstl from .utils import CACHE, NOGIL -# %% ../nbs/src/arima.ipynb 6 +# %% ../../nbs/src/arima.ipynb 6 OptimResult = namedtuple("OptimResult", "success status x fun hess_inv") -# %% ../nbs/src/arima.ipynb 7 +# %% ../../nbs/src/arima.ipynb 7 @njit(nogil=NOGIL, cache=CACHE) def partrans(p, raw, new): if p > 100: @@ -39,7 +39,7 @@ def partrans(p, raw, new): work[k] -= a * new[j - k - 1] new[:j] = work[:j] -# %% ../nbs/src/arima.ipynb 8 +# %% ../../nbs/src/arima.ipynb 8 @njit(nogil=NOGIL, cache=CACHE) def arima_gradtrans(x, arma): eps = 1e-3 @@ -71,7 +71,7 @@ def arima_gradtrans(x, arma): w1[i] -= eps return y -# %% ../nbs/src/arima.ipynb 10 +# %% ../../nbs/src/arima.ipynb 10 @njit(nogil=NOGIL, cache=CACHE) def arima_undopars(x, arma): mp, mq, msp = arma[:3] @@ -83,7 +83,7 @@ def arima_undopars(x, arma): partrans(msp, x[v:], res[v:]) return res -# %% ../nbs/src/arima.ipynb 12 +# %% ../../nbs/src/arima.ipynb 12 @njit(nogil=NOGIL, cache=CACHE) def tsconv(a, b): na = len(a) @@ -98,7 +98,7 @@ def tsconv(a, b): return ab -# %% ../nbs/src/arima.ipynb 14 +# %% ../../nbs/src/arima.ipynb 14 @njit(nogil=NOGIL, cache=CACHE) def inclu2(np_, xnext, xrow, ynext, d, rbar, thetab): for i in range(np_): @@ -127,7 +127,7 @@ def inclu2(np_, xnext, xrow, ynext, d, rbar, thetab): else: ithisr = ithisr + np_ - i - 1 -# %% ../nbs/src/arima.ipynb 15 +# %% ../../nbs/src/arima.ipynb 15 @njit(nogil=NOGIL, cache=CACHE) def invpartrans(p, phi, new): if p > 100: @@ -145,7 +145,7 @@ def invpartrans(p, phi, new): for j in range(p): new[j] = math.atanh(new[j]) -# %% ../nbs/src/arima.ipynb 16 +# %% ../../nbs/src/arima.ipynb 16 @njit(nogil=NOGIL, cache=CACHE) def ARIMA_invtrans(x, arma): mp, mq, msp = arma[:3] @@ -157,7 +157,7 @@ def ARIMA_invtrans(x, arma): invpartrans(msp, x[v:], y[v:]) return y -# %% ../nbs/src/arima.ipynb 18 +# %% ../../nbs/src/arima.ipynb 18 @njit(nogil=NOGIL, cache=CACHE) def getQ0(phi, theta): p = len(phi) @@ -287,7 +287,7 @@ def getQ0(phi, theta): res = res.reshape((r, r)) return res -# %% ../nbs/src/arima.ipynb 20 +# %% ../../nbs/src/arima.ipynb 20 @njit(nogil=NOGIL, cache=CACHE) def arima_transpar(params_in, arma, trans): # TODO check trans=True results @@ -326,7 +326,7 @@ def arima_transpar(params_in, arma, trans): return phi, theta -# %% ../nbs/src/arima.ipynb 23 +# %% ../../nbs/src/arima.ipynb 23 @njit(nogil=NOGIL, cache=CACHE) def arima_css(y, arma, phi, theta, ncond): n = len(y) @@ -370,7 +370,7 @@ def arima_css(y, arma, phi, theta, ncond): return res, resid -# %% ../nbs/src/arima.ipynb 25 +# %% ../../nbs/src/arima.ipynb 25 @njit(nogil=NOGIL, cache=CACHE) def _make_arima(phi, theta, delta, kappa=1e6, tol=np.finfo(float).eps): # check nas phi @@ -423,7 +423,7 @@ def make_arima(phi, theta, delta, kappa=1e6, tol=np.finfo(np.float64).eps): res = _make_arima(phi, theta, delta, kappa, tol) return dict(zip(keys, res)) -# %% ../nbs/src/arima.ipynb 27 +# %% ../../nbs/src/arima.ipynb 27 @njit(nogil=NOGIL, cache=CACHE) def arima_like(y, phi, theta, delta, a, P, Pn, up, use_resid): n = len(y) @@ -560,7 +560,7 @@ def arima_like(y, phi, theta, delta, a, P, Pn, up, use_resid): rsResid = None return ssq, sumlog, nu, rsResid -# %% ../nbs/src/arima.ipynb 29 +# %% ../../nbs/src/arima.ipynb 29 @njit(nogil=NOGIL, cache=CACHE) def diff1d(x, lag, differences): y = x.copy() @@ -592,7 +592,7 @@ def diff(x, lag, differences): raise ValueError(x.ndim) return y[~nan_mask] -# %% ../nbs/src/arima.ipynb 30 +# %% ../../nbs/src/arima.ipynb 30 def fixed_params_from_dict( fixed_dict: dict, order: tuple, seasonal: dict, intercept: bool, n_ex: int ): @@ -616,7 +616,7 @@ def fixed_params_from_dict( ) # prevent adding non-existing keys return list(full_dict.values()) -# %% ../nbs/src/arima.ipynb 32 +# %% ../../nbs/src/arima.ipynb 32 def arima( x: np.ndarray, order=(0, 0, 0), @@ -1058,7 +1058,7 @@ def arma_css_op(p, x): } return ans -# %% ../nbs/src/arima.ipynb 40 +# %% ../../nbs/src/arima.ipynb 40 @njit(nogil=NOGIL, cache=CACHE) def kalman_forecast(n, Z, a, P, T, V, h): p = len(a) @@ -1100,13 +1100,13 @@ def kalman_forecast(n, Z, a, P, T, V, h): return forecasts, se -# %% ../nbs/src/arima.ipynb 43 +# %% ../../nbs/src/arima.ipynb 43 def checkarima(obj): if obj["var_coef"] is None: return False return any(np.isnan(np.sqrt(np.diag(obj["var_coef"])))) -# %% ../nbs/src/arima.ipynb 44 +# %% ../../nbs/src/arima.ipynb 44 def predict_arima(model, n_ahead, newxreg=None, se_fit=True): myNCOL = lambda x: x.shape[1] if x is not None else 0 @@ -1165,7 +1165,7 @@ def predict_arima(model, n_ahead, newxreg=None, se_fit=True): return pred -# %% ../nbs/src/arima.ipynb 48 +# %% ../../nbs/src/arima.ipynb 48 def convert_coef_name(name, inverse=False): if not inverse: if "ex" in name: @@ -1187,13 +1187,13 @@ def convert_coef_name(name, inverse=False): else: return name -# %% ../nbs/src/arima.ipynb 49 +# %% ../../nbs/src/arima.ipynb 49 def change_drift_name(model_coef, inverse=False): return { convert_coef_name(name, inverse): value for name, value in model_coef.items() } -# %% ../nbs/src/arima.ipynb 50 +# %% ../../nbs/src/arima.ipynb 50 def myarima( x, order=(0, 0, 0), @@ -1292,7 +1292,7 @@ def myarima( raise e return {"ic": math.inf} -# %% ../nbs/src/arima.ipynb 53 +# %% ../../nbs/src/arima.ipynb 53 def search_arima( x, d=0, @@ -1386,7 +1386,7 @@ def search_arima( ) return best_fit -# %% ../nbs/src/arima.ipynb 55 +# %% ../../nbs/src/arima.ipynb 55 def arima2(x, model, xreg, method): m = model["arma"][4] # 5 use_drift = "drift" in model["coef"].keys() @@ -1453,7 +1453,7 @@ def arima2(x, model, xreg, method): refit["coef"] = change_drift_name(refit["coef"]) return refit -# %% ../nbs/src/arima.ipynb 56 +# %% ../../nbs/src/arima.ipynb 56 def Arima( x, order=(0, 0, 0), @@ -1547,7 +1547,7 @@ def Arima( tmp["sigma2"] = np.nansum(tmp["residuals"] ** 2) / (nstar - npar + 1) return tmp -# %% ../nbs/src/arima.ipynb 64 +# %% ../../nbs/src/arima.ipynb 64 def arima_string(model, padding=False): order = tuple(model["arma"][i] for i in [0, 5, 1, 2, 6, 3, 4]) m = order[6] @@ -1577,11 +1577,11 @@ def arima_string(model, padding=False): return result -# %% ../nbs/src/arima.ipynb 67 +# %% ../../nbs/src/arima.ipynb 67 def is_constant(x): return np.all(x[0] == x) -# %% ../nbs/src/arima.ipynb 68 +# %% ../../nbs/src/arima.ipynb 68 def forecast_arima( model, h=None, @@ -1676,7 +1676,7 @@ def forecast_arima( return ans -# %% ../nbs/src/arima.ipynb 75 +# %% ../../nbs/src/arima.ipynb 75 def fitted_arima(model, h=1): """Returns h-step forecasts for the data used in fitting the model.""" if h == 1: @@ -1692,7 +1692,7 @@ def fitted_arima(model, h=1): else: raise NotImplementedError("h > 1") -# %% ../nbs/src/arima.ipynb 80 +# %% ../../nbs/src/arima.ipynb 80 def seas_heuristic(x, period): # nperiods = period > 1 season = math.nan @@ -1704,7 +1704,7 @@ def seas_heuristic(x, period): season = max(0, min(1, 1 - vare / np.var(remainder + seasonal, ddof=1))) return season -# %% ../nbs/src/arima.ipynb 82 +# %% ../../nbs/src/arima.ipynb 82 def nsdiffs(x, test="seas", alpha=0.05, period=1, max_D=1, **kwargs): D = 0 if alpha < 0.01: @@ -1767,7 +1767,7 @@ def run_tests(x, test, alpha): dodiff = False return D -# %% ../nbs/src/arima.ipynb 84 +# %% ../../nbs/src/arima.ipynb 84 def ndiffs(x, alpha=0.05, test="kpss", kind="level", max_d=2): x = x[~np.isnan(x)] d = 0 @@ -1812,13 +1812,13 @@ def run_tests(x, test, alpha): return d - 1 return d -# %% ../nbs/src/arima.ipynb 86 +# %% ../../nbs/src/arima.ipynb 86 def newmodel(p, d, q, P, D, Q, constant, results): curr = np.array([p, d, q, P, D, Q, constant]) in_results = (curr == results[:, :7]).all(1).any() return not in_results -# %% ../nbs/src/arima.ipynb 88 +# %% ../../nbs/src/arima.ipynb 88 def auto_arima_f( x, d=None, @@ -2356,11 +2356,11 @@ def try_params(p, d, q, P, D, Q, constant, k, bestfit): return bestfit -# %% ../nbs/src/arima.ipynb 90 +# %% ../../nbs/src/arima.ipynb 90 def forward_arima(fitted_model, y, xreg=None, method="CSS-ML"): return Arima(x=y, model=fitted_model, xreg=xreg, method=method) -# %% ../nbs/src/arima.ipynb 99 +# %% ../../nbs/src/arima.ipynb 99 def print_statsforecast_ARIMA(model, digits=3, se=True): print(arima_string(model, padding=False)) if model["lambda"] is not None: @@ -2390,7 +2390,7 @@ def print_statsforecast_ARIMA(model, digits=3, se=True): if not np.isnan(model["aic"]): print(f'AIC={round(model["aic"], 2)}') -# %% ../nbs/src/arima.ipynb 101 +# %% ../../nbs/src/arima.ipynb 101 class ARIMASummary: """ARIMA Summary.""" @@ -2403,7 +2403,7 @@ def __repr__(self): def summary(self): return print_statsforecast_ARIMA(self.model) -# %% ../nbs/src/arima.ipynb 102 +# %% ../../nbs/src/arima.ipynb 102 class AutoARIMA: """An AutoARIMA estimator. diff --git a/statsforecast/ces.py b/python/statsforecast/ces.py similarity index 96% rename from statsforecast/ces.py rename to python/statsforecast/ces.py index 40bae2743..39a0edba7 100644 --- a/statsforecast/ces.py +++ b/python/statsforecast/ces.py @@ -1,9 +1,9 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/ces.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/ces.ipynb. # %% auto 0 __all__ = ['ces_target_fn'] -# %% ../nbs/src/ces.ipynb 1 +# %% ../../nbs/src/ces.ipynb 1 import math from typing import Tuple @@ -11,10 +11,9 @@ from numba import njit from statsmodels.tsa.seasonal import seasonal_decompose -from .ets import restrict_to_bounds, results -from .utils import CACHE, NOGIL +from .utils import CACHE, NOGIL, restrict_to_bounds, results -# %% ../nbs/src/ces.ipynb 4 +# %% ../../nbs/src/ces.ipynb 4 # Global variables NONE = 0 SIMPLE = 1 @@ -25,7 +24,7 @@ NA = -99999.0 smalno = np.finfo(float).eps -# %% ../nbs/src/ces.ipynb 6 +# %% ../../nbs/src/ces.ipynb 6 def initstate(y, m, seasontype): n = len(y) components = 2 + (seasontype == "P") + 2 * (seasontype == "F") @@ -53,7 +52,7 @@ def initstate(y, m, seasontype): return states -# %% ../nbs/src/ces.ipynb 8 +# %% ../../nbs/src/ces.ipynb 8 @njit(nogil=NOGIL, cache=CACHE) def cescalc( y: np.ndarray, @@ -153,7 +152,7 @@ def cescalc( lik = n * math.log(lik) return lik -# %% ../nbs/src/ces.ipynb 9 +# %% ../../nbs/src/ces.ipynb 9 @njit(nogil=NOGIL, cache=CACHE) def cesfcst(states, i, m, season, f, h, alpha_0, alpha_1, beta_0, beta_1): # obs: @@ -174,7 +173,7 @@ def cesfcst(states, i, m, season, f, h, alpha_0, alpha_1, beta_0, beta_1): ) return new_states -# %% ../nbs/src/ces.ipynb 10 +# %% ../../nbs/src/ces.ipynb 10 @njit(nogil=NOGIL, cache=CACHE) def cesupdate( states, i, m, season, alpha_0, alpha_1, beta_0, beta_1, y # kind of season @@ -220,7 +219,7 @@ def cesupdate( states[i - m, 2] + (1 - beta_0) * states[i - m, 3] + (beta_0 + beta_1) * e ) -# %% ../nbs/src/ces.ipynb 11 +# %% ../../nbs/src/ces.ipynb 11 @njit(nogil=NOGIL, cache=CACHE) def cesforecast(states, n, m, season, f, h, alpha_0, alpha_1, beta_0, beta_1): # compute forecasts @@ -239,7 +238,7 @@ def cesforecast(states, n, m, season, f, h, alpha_0, alpha_1, beta_0, beta_1): ) return new_states -# %% ../nbs/src/ces.ipynb 20 +# %% ../../nbs/src/ces.ipynb 20 @njit(nogil=NOGIL, cache=CACHE) def initparamces( alpha_0: float, alpha_1: float, beta_0: float, beta_1: float, seasontype: str @@ -290,12 +289,12 @@ def initparamces( "optimize_beta_1": optimize_beta_1, } -# %% ../nbs/src/ces.ipynb 22 +# %% ../../nbs/src/ces.ipynb 22 @njit(nogil=NOGIL, cache=CACHE) def switch_ces(x: str): return {"N": 0, "S": 1, "P": 2, "F": 3}[x] -# %% ../nbs/src/ces.ipynb 24 +# %% ../../nbs/src/ces.ipynb 24 @njit(nogil=NOGIL, cache=CACHE) def pegelsresid_ces( y: np.ndarray, @@ -332,7 +331,7 @@ def pegelsresid_ces( lik = np.nan return amse, e, states, lik -# %% ../nbs/src/ces.ipynb 25 +# %% ../../nbs/src/ces.ipynb 25 @njit(nogil=NOGIL, cache=CACHE) def ces_target_fn( optimal_param, @@ -402,7 +401,7 @@ def ces_target_fn( lik = -np.inf return lik -# %% ../nbs/src/ces.ipynb 26 +# %% ../../nbs/src/ces.ipynb 26 @njit(nogil=NOGIL, cache=CACHE) def nelder_mead_ces( x0: np.ndarray, @@ -518,7 +517,7 @@ def nelder_mead_ces( f_simplex[i] = ces_target_fn(simplex[i], *args) return results(simplex[best_idx], f_simplex[best_idx], it + 1, simplex) -# %% ../nbs/src/ces.ipynb 27 +# %% ../../nbs/src/ces.ipynb 27 def optimize_ces_target_fn( init_par, optimize_params, y, m, init_states, n_components, seasontype, nmse ): @@ -563,7 +562,7 @@ def optimize_ces_target_fn( ) return res -# %% ../nbs/src/ces.ipynb 28 +# %% ../../nbs/src/ces.ipynb 28 def cesmodel( y: np.ndarray, m: int, @@ -660,7 +659,7 @@ def cesmodel( sigma2=sigma2, ) -# %% ../nbs/src/ces.ipynb 30 +# %% ../../nbs/src/ces.ipynb 30 def pegelsfcast_C(h, obj, npaths=None, level=None, bootstrap=None): forecast = np.full(h, fill_value=np.nan) m = obj["m"] @@ -677,7 +676,7 @@ def pegelsfcast_C(h, obj, npaths=None, level=None, bootstrap=None): ) return forecast -# %% ../nbs/src/ces.ipynb 31 +# %% ../../nbs/src/ces.ipynb 31 def _simulate_pred_intervals(model, h, level): np.random.seed(1) @@ -708,7 +707,7 @@ def _simulate_pred_intervals(model, h, level): return pi -# %% ../nbs/src/ces.ipynb 32 +# %% ../../nbs/src/ces.ipynb 32 def forecast_ces(obj, h, level=None): fcst = pegelsfcast_C(h, obj) out = {"mean": fcst} @@ -718,7 +717,7 @@ def forecast_ces(obj, h, level=None): out = {**out, **pi} return out -# %% ../nbs/src/ces.ipynb 34 +# %% ../../nbs/src/ces.ipynb 34 def auto_ces( y, m, @@ -783,7 +782,7 @@ def auto_ces( raise Exception("no model able to be fitted") return model -# %% ../nbs/src/ces.ipynb 37 +# %% ../../nbs/src/ces.ipynb 37 def forward_ces(fitted_model, y): m = fitted_model["m"] model = fitted_model["seasontype"] diff --git a/statsforecast/core.py b/python/statsforecast/core.py similarity index 99% rename from statsforecast/core.py rename to python/statsforecast/core.py index 8bd0a8189..a02355e2a 100644 --- a/statsforecast/core.py +++ b/python/statsforecast/core.py @@ -1,9 +1,9 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/core/core.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/core/core.ipynb. # %% auto 0 __all__ = ['StatsForecast'] -# %% ../nbs/src/core/core.ipynb 6 +# %% ../../nbs/src/core/core.ipynb 6 import datetime as dt import errno import inspect @@ -35,7 +35,7 @@ from .utils import ConformalIntervals -# %% ../nbs/src/core/core.ipynb 7 +# %% ../../nbs/src/core/core.ipynb 7 if __name__ == "__main__": logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s: %(message)s", @@ -43,7 +43,7 @@ ) logger = logging.getLogger(__name__) -# %% ../nbs/src/core/core.ipynb 10 +# %% ../../nbs/src/core/core.ipynb 10 class GroupedArray(BaseGroupedArray): def __eq__(self, other): @@ -461,7 +461,7 @@ def _single_threaded_cross_validation( target_col=target_col, ) -# %% ../nbs/src/core/core.ipynb 24 +# %% ../../nbs/src/core/core.ipynb 24 def _get_n_jobs(n_groups, n_jobs): if n_jobs == -1 or (n_jobs is None): actual_n_jobs = os.cpu_count() @@ -469,7 +469,7 @@ def _get_n_jobs(n_groups, n_jobs): actual_n_jobs = n_jobs return min(n_groups, actual_n_jobs) -# %% ../nbs/src/core/core.ipynb 27 +# %% ../../nbs/src/core/core.ipynb 27 def _warn_df_constructor(): warnings.warn( "The `df` argument of the StatsForecast constructor as well as reusing stored " @@ -500,7 +500,7 @@ def _warn_id_as_idx(): def _id_as_idx() -> bool: return not bool(os.getenv("NIXTLA_ID_AS_COL", "")) -# %% ../nbs/src/core/core.ipynb 28 +# %% ../../nbs/src/core/core.ipynb 28 _param_descriptions = { "freq": """freq : str or int Frequency of the data. Must be a valid pandas or polars offset alias, or an integer.""", @@ -547,7 +547,7 @@ def _id_as_idx() -> bool: If int, train the models every `refit` windows.""", } -# %% ../nbs/src/core/core.ipynb 29 +# %% ../../nbs/src/core/core.ipynb 29 class _StatsForecast: """The `StatsForecast` class allows you to efficiently fit multiple `StatsForecast` models for large sets of time series. It operates on a DataFrame `df` with at least three columns @@ -1579,7 +1579,7 @@ def __repr__(self): _StatsForecast.plot.__doc__ = _StatsForecast.plot.__doc__.format(**_param_descriptions) # type: ignore[union-attr] -# %% ../nbs/src/core/core.ipynb 30 +# %% ../../nbs/src/core/core.ipynb 30 class ParallelBackend: def forecast( self, @@ -1660,7 +1660,7 @@ def cross_validation( def make_backend(obj: Any, *args: Any, **kwargs: Any) -> ParallelBackend: return ParallelBackend() -# %% ../nbs/src/core/core.ipynb 31 +# %% ../../nbs/src/core/core.ipynb 31 class StatsForecast(_StatsForecast): def forecast( self, diff --git a/statsforecast/distributed/__init__.py b/python/statsforecast/distributed/__init__.py similarity index 100% rename from statsforecast/distributed/__init__.py rename to python/statsforecast/distributed/__init__.py diff --git a/statsforecast/distributed/fugue.py b/python/statsforecast/distributed/fugue.py similarity index 98% rename from statsforecast/distributed/fugue.py rename to python/statsforecast/distributed/fugue.py index 470a92f7f..dbb9c0315 100644 --- a/statsforecast/distributed/fugue.py +++ b/python/statsforecast/distributed/fugue.py @@ -1,9 +1,9 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/core/distributed.fugue.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../../nbs/src/core/distributed.fugue.ipynb. # %% auto 0 __all__ = ['FugueBackend'] -# %% ../../nbs/src/core/distributed.fugue.ipynb 4 +# %% ../../../nbs/src/core/distributed.fugue.ipynb 4 import inspect from typing import Any, Dict, Iterable, List, Optional, Tuple, Union @@ -25,7 +25,7 @@ ) from ..utils import ConformalIntervals -# %% ../../nbs/src/core/distributed.fugue.ipynb 5 +# %% ../../../nbs/src/core/distributed.fugue.ipynb 5 def _cotransform( df1: Any, df2: Any, @@ -54,7 +54,7 @@ def _cotransform( return result return result.as_pandas() if result.is_local else result.native # type:ignore -# %% ../../nbs/src/core/distributed.fugue.ipynb 6 +# %% ../../../nbs/src/core/distributed.fugue.ipynb 6 class FugueBackend(ParallelBackend): """FugueBackend for Distributed Computation. [Source code](https://github.com/Nixtla/statsforecast/blob/main/statsforecast/distributed/fugue.py). diff --git a/statsforecast/distributed/multiprocess.py b/python/statsforecast/distributed/multiprocess.py similarity index 87% rename from statsforecast/distributed/multiprocess.py rename to python/statsforecast/distributed/multiprocess.py index be62f3d1c..56f23f948 100644 --- a/statsforecast/distributed/multiprocess.py +++ b/python/statsforecast/distributed/multiprocess.py @@ -1,14 +1,14 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/distributed.multiprocess.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../../nbs/src/distributed.multiprocess.ipynb. # %% auto 0 __all__ = ['MultiprocessBackend'] -# %% ../../nbs/src/distributed.multiprocess.ipynb 4 +# %% ../../../nbs/src/distributed.multiprocess.ipynb 4 from typing import Any from ..core import _StatsForecast, ParallelBackend -# %% ../../nbs/src/distributed.multiprocess.ipynb 5 +# %% ../../../nbs/src/distributed.multiprocess.ipynb 5 # This parent class holds common `forecast` and `cross_validation` methods # from `core.StatsForecast` to enable the `FugueBackend` and the `RayBackend`. diff --git a/statsforecast/ets.py b/python/statsforecast/ets.py similarity index 66% rename from statsforecast/ets.py rename to python/statsforecast/ets.py index 2dcde3a9b..2dfa99fbb 100644 --- a/statsforecast/ets.py +++ b/python/statsforecast/ets.py @@ -1,256 +1,133 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/ets.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/ets.ipynb. # %% auto 0 -__all__ = ['NONE', 'ADD', 'MULT', 'DAMPED', 'TOL', 'HUGEN', 'NA', 'smalno', 'ets_target_fn', 'is_constant'] +__all__ = ['ets_f'] -# %% ../nbs/src/ets.ipynb 1 +# %% ../../nbs/src/ets.ipynb 2 import math -from collections import namedtuple -from typing import Tuple import numpy as np -from numba import njit -from numba.typed import List from statsmodels.tsa.seasonal import seasonal_decompose -from .utils import _calculate_intervals, CACHE, NOGIL +from ._lib import ets as _ets +from .utils import _calculate_intervals, results -# %% ../nbs/src/ets.ipynb 5 +# %% ../../nbs/src/ets.ipynb 5 # Global variables -NONE = 0 -ADD = 1 -MULT = 2 -DAMPED = 1 -TOL = 1.0e-10 -HUGEN = 1.0e10 -NA = -99999.0 -smalno = np.finfo(float).eps +_smalno = np.finfo(float).eps _PHI_LOWER = 0.8 _PHI_UPPER = 0.98 -# %% ../nbs/src/ets.ipynb 6 -@njit(nogil=NOGIL, cache=CACHE) -def etscalc(y, n, x, m, error, trend, season, alpha, beta, gamma, phi, e, amse, nmse): - oldb = 0.0 - olds = np.zeros(max(24, m)) - s = np.zeros(max(24, m)) - f = np.zeros(30) - denom = np.zeros(30) - if m < 1: - m = 1 - if nmse > 30: - nmse = 30 - nstates = m * (season > NONE) + 1 + (trend > NONE) - # Copy initial state components - l = x[0] - if trend > NONE: - b = x[1] - else: - b = 0.0 - if season > NONE: - for j in range(m): - s[j] = x[(trend > NONE) + j + 1] - lik = 0.0 - lik2 = 0.0 - for j in range(nmse): - amse[j] = 0.0 - denom[j] = 0.0 - for i in range(n): - # Copy previous state - oldl = l - if trend > NONE: - oldb = b - if season > NONE: - for j in range(m): - olds[j] = s[j] - # one step forecast - forecast(oldl, oldb, olds, m, trend, season, phi, f, nmse) - if math.fabs(f[0] - NA) < TOL: - lik = NA - return lik - if error == ADD: - e[i] = y[i] - f[0] - else: - if math.fabs(f[0]) < TOL: - f_0 = f[0] + TOL - else: - f_0 = f[0] - e[i] = (y[i] - f[0]) / f_0 - for j in range(nmse): - if (i + j) < n: - denom[j] += 1.0 - tmp = y[i + j] - f[j] - amse[j] = (amse[j] * (denom[j] - 1.0) + (tmp * tmp)) / denom[j] - # update state - l, b, s = update( - oldl, l, oldb, b, olds, s, m, trend, season, alpha, beta, gamma, phi, y[i] - ) - # store new state - x[nstates * (i + 1)] = l - if trend > NONE: - x[nstates * (i + 1) + 1] = b - if season > NONE: - for j in range(m): - x[nstates * (i + 1) + (trend > NONE) + j + 1] = s[j] - lik = lik + e[i] * e[i] - val = math.fabs(f[0]) - if val > 0.0: - lik2 += math.log(val) - else: - lik2 += math.log(val + 1e-8) - if lik > 0.0: - lik = n * math.log(lik) - else: - lik = n * math.log(lik + 1e-8) - if error == MULT: - lik += 2 * lik2 - return lik - -# %% ../nbs/src/ets.ipynb 7 -@njit(nogil=NOGIL, cache=CACHE) -def forecast(l, b, s, m, trend, season, phi, f, h): - # f is modified and it is mutable - phistar = phi - # forecasts - for i in range(h): - # print(phistar) - if trend == NONE: - f[i] = l - elif trend == ADD: - f[i] = l + phistar * b - elif b < 0: - f[i] = NA - else: - f[i] = l * math.pow(b, phistar) - j = m - 1 - i - while j < 0: - j += m - if season == ADD: - f[i] = f[i] + s[j] - elif season == MULT: - f[i] = f[i] * s[j] - if i < h - 1: - if math.fabs(phi - 1.0) < TOL: - phistar = phistar + 1.0 - else: - phistar = phistar + math.pow(phi, i + 1) - -# %% ../nbs/src/ets.ipynb 8 -@njit(nogil=NOGIL, cache=CACHE) -def update(oldl, l, oldb, b, olds, s, m, trend, season, alpha, beta, gamma, phi, y): - # New Level - if trend == NONE: - q = oldl # l(t - 1) - phib = 0 - elif trend == ADD: - phib = phi * oldb - q = oldl + phib # l(t - 1) + phi * b(t - 1) - elif math.fabs(phi - 1.0) < TOL: - phib = oldb - q = oldl * oldb # l(t - 1) * b(t - 1) - else: - phib = math.pow(oldb, phi) - q = oldl * phib # l(t - 1) * b(t - 1)^phi - # season - if season == NONE: - p = y - elif season == ADD: - p = y - olds[m - 1] # y[t] - s[t - m] - else: - if math.fabs(olds[m - 1]) < TOL: - p = HUGEN - else: - p = y / olds[m - 1] # y[t] / s[t - m] - l = q + alpha * (p - q) - # New Growth - if trend > NONE: - if trend == ADD: - r = l - oldl # l[t] - l[t-1] - else: # if(trend == MULT) - if math.fabs(oldl) < TOL: - r = HUGEN - else: - r = l / oldl # l[t] / l[t-1] - b = phib + (beta / alpha) * (r - phib) - # b[t] = phi*b[t-1] + beta*(r - phi*b[t-1]) - # b[t] = b[t-1]^phi + beta*(r - b[t-1]^phi) - # New Seasonal - if season > NONE: - if season == ADD: - t = y - q - else: # if(season == MULT) - if math.fabs(q) < TOL: - t = HUGEN - else: - t = y / q - s[0] = olds[m - 1] + gamma * ( - t - olds[m - 1] - ) # s[t] = s[t - m] + gamma * (t - s[t - m]) - for j in range(1, m): - s[j] = olds[j - 1] # s[t] = s[t] - return l, b, s - -# %% ../nbs/src/ets.ipynb 9 -@njit(nogil=NOGIL, cache=CACHE) -def etssimulate(x, m, error, trend, season, alpha, beta, gamma, phi, h, y, e): +# %% ../../nbs/src/ets.ipynb 6 +def etssimulate( + x: np.ndarray, + m: int, + error: _ets.Component, + trend: _ets.Component, + season: _ets.Component, + alpha: float, + beta: float, + gamma: float, + phi: float, + h: int, + y: np.ndarray, + e: np.ndarray, +) -> None: oldb = 0.0 olds = np.zeros(24) s = np.zeros(24) f = np.zeros(10) - if m > 24 and season > NONE: + if m > 24 and season != _ets.Component.Nothing: return elif m < 1: m = 1 - # nstates = m * (season > NONE) + 1 + (trend > NONE) # Copy initial state components l = x[0] - if trend > NONE: + if trend != _ets.Component.Nothing: b = x[1] - if season > NONE: + if season != _ets.Component.Nothing: for j in range(m): - s[j] = x[(trend > NONE) + j + 1] + s[j] = x[(trend != _ets.Component.Nothing) + j + 1] for i in range(h): # Copy previous state oldl = l - if trend > NONE: + if trend != _ets.Component.Nothing: oldb = b - if season > NONE: + if season != _ets.Component.Nothing: for j in range(m): olds[j] = s[j] # one step forecast - forecast(oldl, oldb, olds, m, trend, season, phi, f, 1) - if math.fabs(f[0] - NA) < TOL: - y[0] = NA + _ets.forecast( + f, + oldl, + oldb, + olds, + m, + trend, + season, + phi, + 1, + ) + if math.fabs(f[0] - _ets.NA) < _ets.TOL: + y[0] = _ets.NA return - if error == ADD: + if error == _ets.Component.Additive: y[i] = f[0] + e[i] else: y[i] = f[0] * (1.0 + e[i]) # Update state - l, b, s = update( - oldl, l, oldb, b, olds, s, m, trend, season, alpha, beta, gamma, phi, y[i] + l, b = _ets.update( + s, + l, + b, + oldl, + oldb, + olds, + m, + trend, + season, + alpha, + beta, + gamma, + phi, + y[i], ) -# %% ../nbs/src/ets.ipynb 10 -@njit(nogil=NOGIL, cache=CACHE) -def etsforecast(x, m, trend, season, phi, h, f): - s = np.zeros(m) +# %% ../../nbs/src/ets.ipynb 7 +def etsforecast( + x: np.ndarray, + m: int, + trend: _ets.Component, + season: _ets.Component, + phi: float, + h: int, + f: np.ndarray, +) -> None: + s = np.zeros(m, dtype=np.float64) if m < 1: m = 1 # Copy initial state components l = x[0] b = 0.0 - if trend > NONE: + has_trend = trend != _ets.Component.Nothing + if has_trend: b = x[1] - if season > NONE: - for j in range(m): - s[j] = x[(trend > NONE) + j + 1] - + if season != _ets.Component.Nothing: + s[:m] = x[has_trend + 1 : has_trend + 1 + m] # compute forecasts - forecast(l, b, s, m, trend, season, phi, f, h) + _ets.forecast( + f, + l, + b, + s, + m, + trend, + season, + phi, + h, + ) -# %% ../nbs/src/ets.ipynb 13 -@njit(nogil=NOGIL, cache=CACHE) +# %% ../../nbs/src/ets.ipynb 10 def initparam( alpha: float, beta: float, @@ -294,7 +171,7 @@ def initparam( phi = upper[3] - 1e-3 return {"alpha": alpha, "beta": beta, "gamma": gamma, "phi": phi} -# %% ../nbs/src/ets.ipynb 15 +# %% ../../nbs/src/ets.ipynb 12 def admissible(alpha: float, beta: float, gamma: float, phi: float, m: int): if np.isnan(phi): phi = 1 @@ -331,7 +208,7 @@ def admissible(alpha: float, beta: float, gamma: float, phi: float, m: int): # passed all tests return True -# %% ../nbs/src/ets.ipynb 16 +# %% ../../nbs/src/ets.ipynb 13 def check_param( alpha: float, beta: float, @@ -360,18 +237,7 @@ def check_param( return False return True -# %% ../nbs/src/ets.ipynb 17 -@njit(nogil=NOGIL, cache=CACHE) -def sinpi(x): - return np.sin(np.pi * x) - - -@njit(nogil=NOGIL, cache=CACHE) -def cospi(x): - return np.cos(np.pi * x) - -# %% ../nbs/src/ets.ipynb 18 -@njit(nogil=NOGIL, cache=CACHE) +# %% ../../nbs/src/ets.ipynb 14 def fourier(x, period, K, h=None): if h is None: times = np.arange(1, len(x) + 1) @@ -379,10 +245,7 @@ def fourier(x, period, K, h=None): times = np.arange(len(x) + 1, len(x) + h + 1) # compute periods of all fourier terms # numba doesnt support list comprehension - len_p = 0 - for k in K: - if k > 0: - len_p += k + len_p = sum(K) p = np.full(len_p, fill_value=np.nan) idx = 0 for j, p_ in enumerate(period): @@ -391,24 +254,24 @@ def fourier(x, period, K, h=None): idx += K[j] p = np.unique(p) # Remove columns where sinpi=0 - k = np.abs(2 * p - np.round(2 * p, 0, np.empty_like(p))) > smalno + k = np.abs(2 * p - np.round(2 * p, 0, np.empty_like(p))) > _smalno # Compute matrix of fourier terms X = np.full((len(times), 2 * len(p)), fill_value=np.nan) for j in range(len(p)): if k[j]: - X[:, 2 * j - 1] = sinpi(2 * p[j] * times) - X[:, 2 * j] = cospi(2 * p[j] * times) + X[:, 2 * j - 1] = np.sin(2 * np.pi * p[j] * times) + X[:, 2 * j] = np.cos(2 * np.pi * p[j] * times) X = X[:, ~np.isnan(X.sum(axis=0))] return X -# %% ../nbs/src/ets.ipynb 20 +# %% ../../nbs/src/ets.ipynb 16 def initstate(y, m, trendtype, seasontype): n = len(y) if seasontype != "N": if n < 4: raise ValueError("You've got to be joking (not enough data).") elif n < 3 * m: # fit simple Fourier model - fouriery = fourier(y, List(x for x in [m]), List(x for x in [1])) + fouriery = fourier(y, [m], [1]) X_fourier = np.full((n, 4), fill_value=np.nan) X_fourier[:, 0] = np.ones(n) X_fourier[:, 1] = np.arange(1, n + 1) @@ -479,13 +342,31 @@ def initstate(y, m, trendtype, seasontype): b0 = max(y_sa[1] / div, 1e-3) return np.concatenate([[l0, b0], init_seas]) -# %% ../nbs/src/ets.ipynb 24 -@njit(nogil=NOGIL, cache=CACHE) -def switch(x: str): - return {"N": 0, "A": 1, "M": 2}[x] - -# %% ../nbs/src/ets.ipynb 26 -@njit(nogil=NOGIL, cache=CACHE) +# %% ../../nbs/src/ets.ipynb 20 +def switch(x: str) -> _ets.Component: + if x == "N": + return _ets.Component.Nothing + if x == "A": + return _ets.Component.Additive + if x == "M": + return _ets.Component.Multiplicative + raise ValueError(f"Unknown component {x}") + +# %% ../../nbs/src/ets.ipynb 22 +def switch_criterion(x: str) -> _ets.Criterion: + if x == "lik": + return _ets.Criterion.Likelihood + if x == "mse": + return _ets.Criterion.MSE + if x == "amse": + return _ets.Criterion.AMSE + if x == "sigma": + return _ets.Criterion.Sigma + if x == "mae": + return _ets.Criterion.MAE + raise ValueError(f"Unknown crtierion {x}") + +# %% ../../nbs/src/ets.ipynb 24 def pegelsresid_C( y: np.ndarray, m: int, @@ -512,21 +393,20 @@ def pegelsresid_C( if seasontype == "N": gamma = 0.0 amse = np.full(nmse, fill_value=np.nan) - lik = etscalc( - y=y, - n=n, - x=x, - m=m, - error=switch(errortype), - trend=switch(trendtype), - season=switch(seasontype), - alpha=alpha, - beta=beta, - gamma=gamma, - phi=phi, - e=e, - amse=amse, - nmse=nmse, + lik = _ets.calc( + x, + e, + amse, + nmse, + y, + switch(errortype), + switch(trendtype), + switch(seasontype), + alpha, + beta, + gamma, + phi, + m, ) x = x.reshape((n + 1, p)) if not np.isnan(lik): @@ -534,248 +414,7 @@ def pegelsresid_C( lik = np.nan return amse, e, x, lik -# %% ../nbs/src/ets.ipynb 27 -results = namedtuple("results", "x fn nit simplex") - - -@njit(nogil=NOGIL, cache=CACHE) -def restrict_to_bounds(x, lower, upper): - new_x = np.full_like(x, fill_value=np.nan, dtype=x.dtype) - for i in range(x.size): - lo = lower[i] - up = upper[i] - if x[i] < lo: - new_x[i] = lo - elif x[i] > up: - new_x[i] = up - else: - new_x[i] = x[i] - return new_x - - -@njit(nogil=NOGIL, cache=CACHE) -def nelder_mead_ets( - x0: np.ndarray, - args: Tuple = (), - lower: np.ndarray = np.empty(0), - upper: np.ndarray = np.empty(0), - init_step: float = 0.05, - zero_pert: float = 0.0001, - alpha: float = 1.0, - gamma: float = 2.0, - rho: float = 0.5, - sigma: float = 0.5, - max_iter: int = 2_000, - tol_std: float = 1e-10, - adaptive: bool = False, -): - # We are trying to minimize the function fn(x, args) - # with initial point x0. - # Step 0: - # get x1, ..., x_{n+1} - # the original article suggested a simplex where an initial point is given as x0 - # with the others generated with a fixed step along each dimension in turn. - bounds = len(lower) and len(upper) - if bounds: - x0 = restrict_to_bounds(x0, lower, upper) - - n = x0.size - if adaptive: - gamma = 1.0 + 2.0 / n - rho = 0.75 - 1.0 / (2.0 * n) - sigma = 1.0 - 1.0 / n - simplex = np.full( - (n + 1, n), fill_value=np.nan, dtype=np.float64 - ) # each row is x_j - simplex[:] = x0 - # perturb simplex using `init_step` - diag = np.copy(np.diag(simplex)) - diag[diag == 0.0] = zero_pert - diag[diag != 0.0] *= 1 + init_step - np.fill_diagonal(simplex, diag) - # restrict simplex to bounds if passed - if bounds: - for j in range(n + 1): - simplex[j] = restrict_to_bounds(simplex[j], lower, upper) - # array of the value of f - f_simplex = np.full(n + 1, fill_value=np.nan) - for j in range(n + 1): - f_simplex[j] = ets_target_fn(simplex[j], *args) - for it in range(max_iter): - # Step1: order of f_simplex - # print(simplex) - # print(f_simplex) - order_f = f_simplex.argsort() - best_idx = order_f[0] - worst_idx = order_f[-1] - second_worst_idx = order_f[-2] - # Check whether method should stop. - if np.std(f_simplex) < tol_std: - break - # calculate centroid except argmax f_simplex - x_o = simplex[np.delete(order_f, -1)].sum(axis=0) / n - # Step2: Reflection, Compute reflected point - x_r = x_o + alpha * (x_o - simplex[worst_idx]) - # restrict x_r to bounds if passed - if bounds: - x_r = restrict_to_bounds(x_r, lower, upper) - f_r = ets_target_fn(x_r, *args) - if f_simplex[best_idx] <= f_r < f_simplex[second_worst_idx]: - simplex[worst_idx] = x_r - f_simplex[worst_idx] = f_r - continue - # Step3: Expansion, reflected point is the best point so far - if f_r < f_simplex[best_idx]: - x_e = x_o + gamma * (x_r - x_o) - # restrict x_e to bounds if passed - if bounds: - x_e = restrict_to_bounds(x_e, lower, upper) - f_e = ets_target_fn(x_e, *args) - if f_e < f_r: - simplex[worst_idx] = x_e - f_simplex[worst_idx] = f_e - else: - simplex[worst_idx] = x_r - f_simplex[worst_idx] = f_r - continue - # Step4: outside Contraction - if f_simplex[second_worst_idx] <= f_r < f_simplex[worst_idx]: - x_oc = x_o + rho * (x_r - x_o) - if bounds: - x_oc = restrict_to_bounds(x_oc, lower, upper) - f_oc = ets_target_fn(x_oc, *args) - if f_oc <= f_r: - simplex[worst_idx] = x_oc - f_simplex[worst_idx] = f_oc - continue - # step 5 inside contraction - else: - x_ic = x_o - rho * (x_r - x_o) - # restrict x_c to bounds if passed - if bounds: - x_ic = restrict_to_bounds(x_ic, lower, upper) - f_ic = ets_target_fn(x_ic, *args) - if f_ic < f_simplex[worst_idx]: - simplex[worst_idx] = x_ic - f_simplex[worst_idx] = f_ic - continue - # step 6: shrink - simplex[np.delete(order_f, 0)] = simplex[best_idx] + sigma * ( - simplex[np.delete(order_f, 0)] - simplex[best_idx] - ) - for i in np.delete(order_f, 0): - simplex[i] = restrict_to_bounds(simplex[i], lower, upper) - f_simplex[i] = ets_target_fn(simplex[i], *args) - return results(simplex[best_idx], f_simplex[best_idx], it + 1, simplex) - -# %% ../nbs/src/ets.ipynb 28 -@njit(nogil=NOGIL, cache=CACHE) -def ets_target_fn( - par, - p_y, - p_nstate, - p_errortype, - p_trendtype, - p_seasontype, - p_damped, - p_lower, - p_upper, - p_opt_crit, - p_nmse, - p_bounds, - p_m, - p_optAlpha, - p_optBeta, - p_optGamma, - p_optPhi, - p_givenAlpha, - p_givenBeta, - p_givenGamma, - p_givenPhi, - alpha, - beta, - gamma, - phi, -): - # p_par_length = len(par) - j = 0 - if p_optAlpha: - alpha = par[j] - j += 1 - if p_optBeta: - beta = par[j] - j += 1 - if p_optGamma: - gamma = par[j] - j += 1 - if p_optPhi: - phi = par[j] - # size of state - p = p_nstate + (1 if p_seasontype != 0 else 0) - n = len(p_y) - state = np.zeros(p * (n + 1)) - state[:p_nstate] = par[-p_nstate:] - # Add extra state - if p_seasontype != 0: - start = 1 + (1 if p_trendtype != 0 else 0) - sum_ = state[start:p_nstate].sum() - new_state = p_m * (1 if p_seasontype == 2 else 0) - sum_ - state[p_nstate] = new_state - # check states - if p_seasontype == 2 and state[start:].min() < 0: - return np.inf - amse = np.zeros(30) - e = np.zeros(n) - lik = etscalc( - # we dont need to change these param - p_y, - n, - state, - p_m, - p_errortype, - p_trendtype, - p_seasontype, - # inmutable, wee need to change them? - alpha, - beta, - gamma, - phi, - # mutable - e, - amse, - # inmutable - p_nmse, - ) - if lik < -1e10: - lik = -1e10 - if math.isnan(lik): - lik = -np.inf - if math.fabs(lik + 99999) < 1e-7: - lik = -np.inf - if p_opt_crit == "lik": - objval = lik - elif p_opt_crit == "mse": - objval = amse[0] - elif p_opt_crit == "amse": - mean = 0.0 - for i in range(p_nmse): - mean += amse[i] / p_nmse - objval = mean - elif p_opt_crit == "sigma": - mean = 0.0 - ne = len(e) - for i in range(ne): - mean += e[i] * e[i] / ne - objval = mean - elif p_opt_crit == "mae": - mean = 0 - ne = len(e) - for i in range(ne): - mean += math.fabs(e[i]) / ne - objval = mean - return objval - -# %% ../nbs/src/ets.ipynb 29 +# %% ../../nbs/src/ets.ipynb 25 def optimize_ets_target_fn( x0, par, @@ -823,23 +462,14 @@ def optimize_ets_target_fn( optGamma = not np.isnan(gamma) optPhi = not np.isnan(phi) - givenAlpha = False - givenBeta = False - givenGamma = False - givenPhi = False - if not np.isnan(par_noopt["alpha"]): optAlpha = False - givenAlpha = True if not np.isnan(par_noopt["beta"]): optBeta = False - givenBeta = True if not np.isnan(par_noopt["gamma"]): optGamma = False - givenGamma = True if not np.isnan(par_noopt["phi"]): optPhi = False - givenPhi = True if not damped: phi = 1.0 @@ -847,43 +477,33 @@ def optimize_ets_target_fn( beta = 0.0 if seasontype == "N": gamma = 0.0 - res = nelder_mead_ets( + opt_res = _ets.optimize( x0, - args=( - y, - nstate, - switch(errortype), - switch(trendtype), - switch(seasontype), - damped, - lowerb, - upperb, - opt_crit, - nmse, - bounds, - m, - optAlpha, - optBeta, - optGamma, - optPhi, - givenAlpha, - givenBeta, - givenGamma, - givenPhi, - alpha, - beta, - gamma, - phi, - ), - lower=lowerb, - upper=upperb, - tol_std=1e-4, - max_iter=1_000, - adaptive=True, + y, + nstate, + switch(errortype), + switch(trendtype), + switch(seasontype), + switch_criterion(opt_crit), + nmse, + m, + optAlpha, + optBeta, + optGamma, + optPhi, + alpha, + beta, + gamma, + phi, + lowerb, + upperb, + 1e-4, + 1_000, + True, ) - return res + return results(*opt_res, None) -# %% ../nbs/src/ets.ipynb 30 +# %% ../../nbs/src/ets.ipynb 26 def etsmodel( y: np.ndarray, m: int, @@ -957,7 +577,6 @@ def etsmodel( par=par, states=init_state, ) - fred = optimize_ets_target_fn( x0=par, par=par_, @@ -1056,12 +675,11 @@ def etsmodel( n_params=np_, ) -# %% ../nbs/src/ets.ipynb 32 -@njit(nogil=NOGIL, cache=CACHE) +# %% ../../nbs/src/ets.ipynb 28 def is_constant(x): return np.all(x[0] == x) -# %% ../nbs/src/ets.ipynb 34 +# %% ../../nbs/src/ets.ipynb 30 def ets_f( y, m, @@ -1085,8 +703,7 @@ def ets_f( use_initial_values=False, maxit=2_000, ): - if y.dtype not in (np.float32, np.float64): - y = y.astype(np.float32) + y = y.astype(np.float64, copy=False) # converting params to floats # to improve numba compilation if alpha is None: @@ -1306,7 +923,7 @@ def ets_f( model["method"] = f"ETS({best_e},{best_t}{'d' if best_d else ''},{best_s})" return model -# %% ../nbs/src/ets.ipynb 35 +# %% ../../nbs/src/ets.ipynb 31 def pegelsfcast_C(h, obj, npaths=None, level=None, bootstrap=None): forecast = np.full(h, fill_value=np.nan) states = obj["states"][-1, :] @@ -1316,27 +933,18 @@ def pegelsfcast_C(h, obj, npaths=None, level=None, bootstrap=None): etsforecast(x=states, m=m, trend=ttype, season=stype, phi=phi, h=h, f=forecast) return forecast -# %% ../nbs/src/ets.ipynb 36 -# @njit(nogil=NOGIL, cache=CACHE) +# %% ../../nbs/src/ets.ipynb 32 def _compute_sigmah(pf, h, sigma, cvals): - theta = np.full(h, np.nan) theta[0] = pf[0] ** 2 for k in range(1, h): - sum_val = 0 - for j in range(1, k + 1): - val = cvals[j - 1] ** 2 * theta[k - j] - sum_val = sum_val + val + sum_val = cvals[:k] ** 2 @ theta[:k][::-1] theta[k] = pf[k] ** 2 + sigma * sum_val - sigmah = np.full(h, np.nan) - for k in range(0, h): - sigmah[k] = (1 + sigma) * theta[k] - pf[k] ** 2 - - return sigmah + return (1 + sigma) * theta - pf**2 -# %% ../nbs/src/ets.ipynb 37 +# %% ../../nbs/src/ets.ipynb 33 def _class3models( h, sigma, @@ -1406,10 +1014,8 @@ def _class3models( var = np.zeros(h) for i in range(0, h): - mu[i] = np.matmul(H1, np.matmul(Mh, np.transpose(H2))) - var[i] = (1 + sigma) * np.matmul( - H21, np.matmul(Vh, np.transpose(H21)) - ) + sigma * mu[i] ** 2 + mu[i] = (H1 @ (Mh @ H2.T)).item() + var[i] = (1 + sigma) * (H21 @ (Vh @ H21.T)).item() + sigma * mu[i] ** 2 vecMh = Mh.flatten() exp1 = np.matmul(F21, np.matmul(Vh, np.transpose(F21))) exp2 = np.matmul(F21, np.matmul(Vh, np.transpose(G21))) @@ -1439,7 +1045,7 @@ def _class3models( return var -# %% ../nbs/src/ets.ipynb 38 +# %% ../../nbs/src/ets.ipynb 34 def _compute_pred_intervals(model, forecasts, h, level): sigma = model["sigma2"] season_length = model["m"] @@ -1630,7 +1236,7 @@ def _compute_pred_intervals(model, forecasts, h, level): return pi -# %% ../nbs/src/ets.ipynb 39 +# %% ../../nbs/src/ets.ipynb 35 def forecast_ets(obj, h, level=None): fcst = pegelsfcast_C(h, obj) out = {"mean": fcst} @@ -1641,6 +1247,6 @@ def forecast_ets(obj, h, level=None): out = {**out, **pi} return out -# %% ../nbs/src/ets.ipynb 46 +# %% ../../nbs/src/ets.ipynb 42 def forward_ets(fitted_model, y): return ets_f(y=y, m=fitted_model["m"], model=fitted_model) diff --git a/statsforecast/feature_engineering.py b/python/statsforecast/feature_engineering.py similarity index 93% rename from statsforecast/feature_engineering.py rename to python/statsforecast/feature_engineering.py index af753d1db..a3f04c14e 100644 --- a/statsforecast/feature_engineering.py +++ b/python/statsforecast/feature_engineering.py @@ -1,9 +1,9 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/feature_engineering.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/feature_engineering.ipynb. # %% auto 0 __all__ = ['mstl_decomposition'] -# %% ../nbs/src/feature_engineering.ipynb 3 +# %% ../../nbs/src/feature_engineering.ipynb 3 from typing import Tuple import pandas as pd @@ -20,7 +20,7 @@ from .core import _id_as_idx from .models import MSTL, _predict_mstl_components -# %% ../nbs/src/feature_engineering.ipynb 4 +# %% ../../nbs/src/feature_engineering.ipynb 4 def mstl_decomposition( df: DataFrame, model: MSTL, diff --git a/statsforecast/garch.py b/python/statsforecast/garch.py similarity index 92% rename from statsforecast/garch.py rename to python/statsforecast/garch.py index 56bdc4b0b..3cc706733 100644 --- a/statsforecast/garch.py +++ b/python/statsforecast/garch.py @@ -1,16 +1,16 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/garch.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/garch.ipynb. # %% auto 0 __all__ = ['garch_model', 'garch_forecast'] -# %% ../nbs/src/garch.ipynb 4 +# %% ../../nbs/src/garch.ipynb 4 import numpy as np from numba import njit from scipy.optimize import minimize from .utils import CACHE, NOGIL -# %% ../nbs/src/garch.ipynb 7 +# %% ../../nbs/src/garch.ipynb 7 @njit(nogil=NOGIL, cache=CACHE) def generate_garch_data(n, w, alpha, beta): @@ -54,7 +54,7 @@ def generate_garch_data(n, w, alpha, beta): return y -# %% ../nbs/src/garch.ipynb 12 +# %% ../../nbs/src/garch.ipynb 12 @njit(nogil=NOGIL, cache=CACHE) def garch_sigma2(x0, x, p, q): @@ -77,14 +77,14 @@ def garch_sigma2(x0, x, p, q): return sigma2 -# %% ../nbs/src/garch.ipynb 14 +# %% ../../nbs/src/garch.ipynb 14 @njit(nogil=NOGIL, cache=CACHE) def garch_cons(x0): # Constraints for GARCH model # alpha+beta < 1 return 1 - (x0[1:].sum()) -# %% ../nbs/src/garch.ipynb 16 +# %% ../../nbs/src/garch.ipynb 16 @njit(nogil=NOGIL, cache=CACHE) def garch_loglik(x0, x, p, q): @@ -101,7 +101,7 @@ def garch_loglik(x0, x, p, q): return -loglik -# %% ../nbs/src/garch.ipynb 18 +# %% ../../nbs/src/garch.ipynb 18 def garch_model(x, p, q): np.random.seed(1) @@ -132,7 +132,7 @@ def garch_model(x, p, q): return res -# %% ../nbs/src/garch.ipynb 22 +# %% ../../nbs/src/garch.ipynb 22 def garch_forecast(mod, h): np.random.seed(1) diff --git a/statsforecast/mfles.py b/python/statsforecast/mfles.py similarity index 99% rename from statsforecast/mfles.py rename to python/statsforecast/mfles.py index f8a885479..eb33be432 100644 --- a/statsforecast/mfles.py +++ b/python/statsforecast/mfles.py @@ -1,9 +1,9 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/mfles.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/mfles.ipynb. # %% auto 0 __all__ = ['MFLES'] -# %% ../nbs/src/mfles.ipynb 3 +# %% ../../nbs/src/mfles.ipynb 3 import itertools import warnings @@ -14,7 +14,7 @@ from .utils import _ensure_float -# %% ../nbs/src/mfles.ipynb 4 +# %% ../../nbs/src/mfles.ipynb 4 # utility functions def calc_mse(y_true, y_pred): sq_err = (y_true - y_pred) ** 2 @@ -345,7 +345,7 @@ class Zeros: def predict(self, X): return np.zeros(X.shape[0]) -# %% ../nbs/src/mfles.ipynb 5 +# %% ../../nbs/src/mfles.ipynb 5 class MFLES: def __init__(self, verbose=1, robust=None): self.penalty = None diff --git a/statsforecast/models.py b/python/statsforecast/models.py similarity index 98% rename from statsforecast/models.py rename to python/statsforecast/models.py index f1623f299..4476d48de 100644 --- a/statsforecast/models.py +++ b/python/statsforecast/models.py @@ -1,4 +1,4 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/core/models.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/core/models.ipynb. # %% auto 0 __all__ = ['AutoARIMA', 'AutoETS', 'ETS', 'AutoCES', 'AutoTheta', 'ARIMA', 'AutoRegressive', 'SimpleExponentialSmoothing', @@ -9,7 +9,7 @@ 'DynamicTheta', 'DynamicOptimizedTheta', 'GARCH', 'ARCH', 'SklearnModel', 'MFLES', 'AutoMFLES', 'ConstantModel', 'ZeroModel', 'NaNModel'] -# %% ../nbs/src/core/models.ipynb 5 +# %% ../../nbs/src/core/models.ipynb 5 import warnings from math import trunc from typing import Any, Dict, List, Optional, Sequence, Tuple, Union @@ -55,7 +55,7 @@ NOGIL, ) -# %% ../nbs/src/core/models.ipynb 9 +# %% ../../nbs/src/core/models.ipynb 9 def _add_fitted_pi(res, se, level): level = sorted(level) level = np.asarray(level) @@ -68,7 +68,7 @@ def _add_fitted_pi(res, se, level): res = {**res, **lo, **hi} return res -# %% ../nbs/src/core/models.ipynb 10 +# %% ../../nbs/src/core/models.ipynb 10 def _add_conformal_distribution_intervals( fcst: Dict, cs: np.ndarray, @@ -97,7 +97,7 @@ def _add_conformal_distribution_intervals( fcst[col] = quantiles[i] return fcst -# %% ../nbs/src/core/models.ipynb 11 +# %% ../../nbs/src/core/models.ipynb 11 def _get_conformal_method(method: str): available_methods = { "conformal_distribution": _add_conformal_distribution_intervals, @@ -110,7 +110,7 @@ def _get_conformal_method(method: str): ) return available_methods[method] -# %% ../nbs/src/core/models.ipynb 12 +# %% ../../nbs/src/core/models.ipynb 12 class _TS: uses_exog = False @@ -171,7 +171,7 @@ def _add_conformal_intervals(self, fcst, y, X, level): def _add_predict_conformal_intervals(self, fcst, level): return self._add_conformal_intervals(fcst=fcst, y=None, X=None, level=level) -# %% ../nbs/src/core/models.ipynb 17 +# %% ../../nbs/src/core/models.ipynb 17 class AutoARIMA(_TS): r"""AutoARIMA model. @@ -596,7 +596,7 @@ def forward( res = _add_fitted_pi(res=res, se=se, level=level) return res -# %% ../nbs/src/core/models.ipynb 33 +# %% ../../nbs/src/core/models.ipynb 33 class AutoETS(_TS): r"""Automatic Exponential Smoothing model. @@ -858,7 +858,7 @@ def forward( res = _add_fitted_pi(res=res, se=se, level=level) return res -# %% ../nbs/src/core/models.ipynb 48 +# %% ../../nbs/src/core/models.ipynb 48 class ETS(AutoETS): @classmethod def _warn(cls): @@ -887,7 +887,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../nbs/src/core/models.ipynb 53 +# %% ../../nbs/src/core/models.ipynb 53 class AutoCES(_TS): r"""Complex Exponential Smoothing model. @@ -1139,7 +1139,7 @@ def forward( res = _add_fitted_pi(res=res, se=se, level=level) return res -# %% ../nbs/src/core/models.ipynb 71 +# %% ../../nbs/src/core/models.ipynb 71 class AutoTheta(_TS): r"""AutoTheta model. @@ -1356,7 +1356,7 @@ def forward( res = _add_fitted_pi(res=res, se=se, level=level) return res -# %% ../nbs/src/core/models.ipynb 87 +# %% ../../nbs/src/core/models.ipynb 87 class ARIMA(_TS): r"""ARIMA model. @@ -1656,7 +1656,7 @@ def forward( res = _add_fitted_pi(res=res, se=se, level=level) return res -# %% ../nbs/src/core/models.ipynb 102 +# %% ../../nbs/src/core/models.ipynb 102 class AutoRegressive(ARIMA): r"""Simple Autoregressive model. @@ -1728,7 +1728,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../nbs/src/core/models.ipynb 117 +# %% ../../nbs/src/core/models.ipynb 117 @njit(nogil=NOGIL, cache=CACHE) def _ses_fcst_mse(x: np.ndarray, alpha: float) -> Tuple[float, float, np.ndarray]: r"""Perform simple exponential smoothing on a series. @@ -1799,7 +1799,7 @@ def _chunk_sums(array: np.ndarray, chunk_size: int) -> np.ndarray: n_elems = n_chunks * chunk_size return array[:n_elems].reshape(n_chunks, chunk_size).sum(axis=1) -# %% ../nbs/src/core/models.ipynb 118 +# %% ../../nbs/src/core/models.ipynb 118 def _ses( y: np.ndarray, # time series h: int, # forecasting horizon @@ -1812,7 +1812,7 @@ def _ses( fcst["fitted"] = fitted_vals return fcst -# %% ../nbs/src/core/models.ipynb 119 +# %% ../../nbs/src/core/models.ipynb 119 class SimpleExponentialSmoothing(_TS): r"""SimpleExponentialSmoothing model. @@ -1972,7 +1972,7 @@ def forecast( raise Exception("You must pass `prediction_intervals` to " "compute them.") return res -# %% ../nbs/src/core/models.ipynb 131 +# %% ../../nbs/src/core/models.ipynb 131 def _ses_optimized( y: np.ndarray, # time series h: int, # forecasting horizon @@ -1985,7 +1985,7 @@ def _ses_optimized( fcst["fitted"] = fitted_vals return fcst -# %% ../nbs/src/core/models.ipynb 132 +# %% ../../nbs/src/core/models.ipynb 132 class SimpleExponentialSmoothingOptimized(_TS): r"""SimpleExponentialSmoothing model. @@ -2139,7 +2139,7 @@ def forecast( raise Exception("You must pass `prediction_intervals` to compute them.") return res -# %% ../nbs/src/core/models.ipynb 144 +# %% ../../nbs/src/core/models.ipynb 144 def _seasonal_exponential_smoothing( y: np.ndarray, # time series h: int, # forecasting horizon @@ -2163,7 +2163,7 @@ def _seasonal_exponential_smoothing( fcst["fitted"] = fitted_vals return fcst -# %% ../nbs/src/core/models.ipynb 145 +# %% ../../nbs/src/core/models.ipynb 145 class SeasonalExponentialSmoothing(_TS): r"""SeasonalExponentialSmoothing model. @@ -2338,7 +2338,7 @@ def forecast( raise Exception("You must pass `prediction_intervals` to compute them.") return res -# %% ../nbs/src/core/models.ipynb 160 +# %% ../../nbs/src/core/models.ipynb 160 def _seasonal_ses_optimized( y: np.ndarray, # time series h: int, # forecasting horizon @@ -2361,7 +2361,7 @@ def _seasonal_ses_optimized( fcst["fitted"] = fitted_vals return fcst -# %% ../nbs/src/core/models.ipynb 161 +# %% ../../nbs/src/core/models.ipynb 161 class SeasonalExponentialSmoothingOptimized(_TS): def __init__( @@ -2533,7 +2533,7 @@ def forecast( raise Exception("You must pass `prediction_intervals` to compute them.") return res -# %% ../nbs/src/core/models.ipynb 174 +# %% ../../nbs/src/core/models.ipynb 174 class Holt(AutoETS): r"""Holt's method. @@ -2575,7 +2575,7 @@ def __init__( season_length, model, alias=alias, prediction_intervals=prediction_intervals ) -# %% ../nbs/src/core/models.ipynb 188 +# %% ../../nbs/src/core/models.ipynb 188 class HoltWinters(AutoETS): r"""Holt-Winters' method. @@ -2615,7 +2615,7 @@ def __init__( season_length, model, alias=alias, prediction_intervals=prediction_intervals ) -# %% ../nbs/src/core/models.ipynb 203 +# %% ../../nbs/src/core/models.ipynb 203 def _historic_average( y: np.ndarray, # time series h: int, # forecasting horizon @@ -2627,7 +2627,7 @@ def _historic_average( fcst["fitted"] = fitted_vals return fcst -# %% ../nbs/src/core/models.ipynb 204 +# %% ../../nbs/src/core/models.ipynb 204 class HistoricAverage(_TS): def __init__( @@ -2803,7 +2803,7 @@ def forecast( return res -# %% ../nbs/src/core/models.ipynb 217 +# %% ../../nbs/src/core/models.ipynb 217 class Naive(_TS): def __init__( @@ -3010,7 +3010,7 @@ def forward( ) return res -# %% ../nbs/src/core/models.ipynb 233 +# %% ../../nbs/src/core/models.ipynb 233 def _random_walk_with_drift( y: np.ndarray, # time series h: int, # forecasting horizon @@ -3029,7 +3029,7 @@ def _random_walk_with_drift( fcst["fitted"] = fitted_vals return fcst -# %% ../nbs/src/core/models.ipynb 234 +# %% ../../nbs/src/core/models.ipynb 234 class RandomWalkWithDrift(_TS): def __init__( @@ -3205,7 +3205,7 @@ def forecast( return res -# %% ../nbs/src/core/models.ipynb 249 +# %% ../../nbs/src/core/models.ipynb 249 class SeasonalNaive(_TS): def __init__( @@ -3391,7 +3391,7 @@ def forecast( return res -# %% ../nbs/src/core/models.ipynb 264 +# %% ../../nbs/src/core/models.ipynb 264 def _window_average( y: np.ndarray, # time series h: int, # forecasting horizon @@ -3406,7 +3406,7 @@ def _window_average( mean = _repeat_val(val=wavg, h=h) return {"mean": mean} -# %% ../nbs/src/core/models.ipynb 265 +# %% ../../nbs/src/core/models.ipynb 265 class WindowAverage(_TS): def __init__( @@ -3562,7 +3562,7 @@ def forecast( raise Exception("You must pass `prediction_intervals` to compute them.") return res -# %% ../nbs/src/core/models.ipynb 276 +# %% ../../nbs/src/core/models.ipynb 276 def _seasonal_window_average( y: np.ndarray, h: int, @@ -3579,7 +3579,7 @@ def _seasonal_window_average( out = _repeat_val_seas(season_vals=season_avgs, h=h) return {"mean": out} -# %% ../nbs/src/core/models.ipynb 277 +# %% ../../nbs/src/core/models.ipynb 277 class SeasonalWindowAverage(_TS): def __init__( @@ -3749,7 +3749,7 @@ def forecast( raise Exception("You must pass `prediction_intervals` to compute them.") return res -# %% ../nbs/src/core/models.ipynb 289 +# %% ../../nbs/src/core/models.ipynb 289 def _chunk_forecast(y, aggregation_level): lost_remainder_data = len(y) % aggregation_level y_cut = y[lost_remainder_data:] @@ -3832,7 +3832,7 @@ def _adida( res["fitted"] = np.append(np.nan, sums_fitted / fitted_aggregation_levels) return res -# %% ../nbs/src/core/models.ipynb 290 +# %% ../../nbs/src/core/models.ipynb 290 class ADIDA(_TS): def __init__( @@ -3999,7 +3999,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../nbs/src/core/models.ipynb 302 +# %% ../../nbs/src/core/models.ipynb 302 def _croston_classic( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4027,7 +4027,7 @@ def _croston_classic( out["fitted"] = ydf / yif return out -# %% ../nbs/src/core/models.ipynb 303 +# %% ../../nbs/src/core/models.ipynb 303 class CrostonClassic(_TS): def __init__( @@ -4189,7 +4189,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../nbs/src/core/models.ipynb 314 +# %% ../../nbs/src/core/models.ipynb 314 def _croston_optimized( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4231,7 +4231,7 @@ def _croston_optimized( out["fitted"] = ydf / yif return out -# %% ../nbs/src/core/models.ipynb 315 +# %% ../../nbs/src/core/models.ipynb 315 class CrostonOptimized(_TS): def __init__( @@ -4390,7 +4390,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../nbs/src/core/models.ipynb 326 +# %% ../../nbs/src/core/models.ipynb 326 def _croston_sba( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4402,7 +4402,7 @@ def _croston_sba( out["fitted"] *= 0.95 return out -# %% ../nbs/src/core/models.ipynb 327 +# %% ../../nbs/src/core/models.ipynb 327 class CrostonSBA(_TS): def __init__( @@ -4566,7 +4566,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../nbs/src/core/models.ipynb 338 +# %% ../../nbs/src/core/models.ipynb 338 def _imapa( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4600,7 +4600,7 @@ def _imapa( res["fitted"] = fitted_vals return res -# %% ../nbs/src/core/models.ipynb 339 +# %% ../../nbs/src/core/models.ipynb 339 class IMAPA(_TS): def __init__( @@ -4763,7 +4763,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../nbs/src/core/models.ipynb 350 +# %% ../../nbs/src/core/models.ipynb 350 def _tsb( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4788,7 +4788,7 @@ def _tsb( res["fitted"] = ypft * ydft return res -# %% ../nbs/src/core/models.ipynb 351 +# %% ../../nbs/src/core/models.ipynb 351 class TSB(_TS): def __init__( @@ -4961,7 +4961,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../nbs/src/core/models.ipynb 363 +# %% ../../nbs/src/core/models.ipynb 363 def _predict_mstl_components(mstl_ob, h, season_length): seasoncolumns = mstl_ob.filter(regex="seasonal*").columns nseasons = len(seasoncolumns) @@ -4982,7 +4982,7 @@ def _predict_mstl_seas(mstl_ob, h, season_length): seascomp = _predict_mstl_components(mstl_ob, h, season_length) return seascomp.sum(axis=1) -# %% ../nbs/src/core/models.ipynb 364 +# %% ../../nbs/src/core/models.ipynb 364 class MSTL(_TS): r"""MSTL model. @@ -5259,7 +5259,7 @@ def forward( } return res -# %% ../nbs/src/core/models.ipynb 380 +# %% ../../nbs/src/core/models.ipynb 380 class TBATS(_TS): r"""Trigonometric Box-Cox transform, ARMA errors, Trend and Seasonal components (TBATS) model. @@ -5471,7 +5471,7 @@ def forecast( res_trans = res return res_trans -# %% ../nbs/src/core/models.ipynb 388 +# %% ../../nbs/src/core/models.ipynb 388 class AutoTBATS(TBATS): r"""AutoTBATS model. @@ -5530,7 +5530,7 @@ def __init__( alias=alias, ) -# %% ../nbs/src/core/models.ipynb 398 +# %% ../../nbs/src/core/models.ipynb 398 class Theta(AutoTheta): r"""Standard Theta Method. @@ -5567,7 +5567,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../nbs/src/core/models.ipynb 412 +# %% ../../nbs/src/core/models.ipynb 412 class OptimizedTheta(AutoTheta): r"""Optimized Theta Method. @@ -5604,7 +5604,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../nbs/src/core/models.ipynb 426 +# %% ../../nbs/src/core/models.ipynb 426 class DynamicTheta(AutoTheta): r"""Dynamic Standard Theta Method. @@ -5641,7 +5641,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../nbs/src/core/models.ipynb 440 +# %% ../../nbs/src/core/models.ipynb 440 class DynamicOptimizedTheta(AutoTheta): r"""Dynamic Optimized Theta Method. @@ -5678,7 +5678,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../nbs/src/core/models.ipynb 455 +# %% ../../nbs/src/core/models.ipynb 455 class GARCH(_TS): r"""Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model. @@ -5870,7 +5870,7 @@ def forecast( res = _add_fitted_pi(res=res, se=se, level=level) return res -# %% ../nbs/src/core/models.ipynb 468 +# %% ../../nbs/src/core/models.ipynb 468 class ARCH(GARCH): r"""Autoregressive Conditional Heteroskedasticity (ARCH) model. @@ -5914,7 +5914,7 @@ def __init__( self.alias = alias super().__init__(p, q=0, alias=alias) -# %% ../nbs/src/core/models.ipynb 479 +# %% ../../nbs/src/core/models.ipynb 479 class SklearnModel(_TS): r"""scikit-learn model wrapper @@ -6120,7 +6120,7 @@ def forward( res = _add_fitted_pi(res=res, se=se, level=level) return res -# %% ../nbs/src/core/models.ipynb 489 +# %% ../../nbs/src/core/models.ipynb 489 class MFLES(_TS): r"""MFLES model. @@ -6398,7 +6398,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../nbs/src/core/models.ipynb 497 +# %% ../../nbs/src/core/models.ipynb 497 class AutoMFLES(_TS): r"""AutoMFLES @@ -6597,7 +6597,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../nbs/src/core/models.ipynb 501 +# %% ../../nbs/src/core/models.ipynb 501 class ConstantModel(_TS): def __init__(self, constant: float, alias: str = "ConstantModel"): @@ -6780,7 +6780,7 @@ def forward( ) return res -# %% ../nbs/src/core/models.ipynb 515 +# %% ../../nbs/src/core/models.ipynb 515 class ZeroModel(ConstantModel): def __init__(self, alias: str = "ZeroModel"): @@ -6795,7 +6795,7 @@ def __init__(self, alias: str = "ZeroModel"): """ super().__init__(constant=0, alias=alias) -# %% ../nbs/src/core/models.ipynb 529 +# %% ../../nbs/src/core/models.ipynb 529 class NaNModel(ConstantModel): def __init__(self, alias: str = "NaNModel"): diff --git a/statsforecast/mstl.py b/python/statsforecast/mstl.py similarity index 94% rename from statsforecast/mstl.py rename to python/statsforecast/mstl.py index 6814306d8..7d0cf4802 100644 --- a/statsforecast/mstl.py +++ b/python/statsforecast/mstl.py @@ -1,16 +1,16 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/mstl.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/mstl.ipynb. # %% auto 0 __all__ = ['mstl'] -# %% ../nbs/src/mstl.ipynb 3 +# %% ../../nbs/src/mstl.ipynb 3 from typing import Dict, List, Optional, Union import numpy as np import pandas as pd import statsmodels.api as sm -# %% ../nbs/src/mstl.ipynb 4 +# %% ../../nbs/src/mstl.ipynb 4 def mstl( x: np.ndarray, # time series period: Union[int, List[int]], # season length diff --git a/statsforecast/tbats.py b/python/statsforecast/tbats.py similarity index 97% rename from statsforecast/tbats.py rename to python/statsforecast/tbats.py index e1fb977e7..8b4eb436a 100644 --- a/statsforecast/tbats.py +++ b/python/statsforecast/tbats.py @@ -1,9 +1,9 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/tbats.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/tbats.ipynb. # %% auto 0 __all__ = ['tbats_model', 'tbats_selection', 'tbats_forecast'] -# %% ../nbs/src/tbats.ipynb 2 +# %% ../../nbs/src/tbats.ipynb 2 import warnings from functools import partial from itertools import product @@ -19,7 +19,7 @@ from .arima import auto_arima_f from .utils import NOGIL, CACHE -# %% ../nbs/src/tbats.ipynb 7 +# %% ../../nbs/src/tbats.ipynb 7 def find_harmonics(y, m): # Compute a 2 x m moving average to estimate the trend @@ -74,7 +74,7 @@ def find_harmonics(y, m): return num_harmonics, y - X[:, : 2 * num_harmonics] @ best_model -# %% ../nbs/src/tbats.ipynb 9 +# %% ../../nbs/src/tbats.ipynb 9 def initial_parameters(k_vector, use_trend, use_damped_trend, ar_coeffs, ma_coeffs): alpha = 0.09 @@ -119,7 +119,7 @@ def initial_parameters(k_vector, use_trend, use_damped_trend, ar_coeffs, ma_coef epsilon_vector, ) -# %% ../nbs/src/tbats.ipynb 11 +# %% ../../nbs/src/tbats.ipynb 11 def makeXMatrix(b, s_vector, d_vector, epsilon_vector): # x = (l_t, b_t, s_vector, d_vector, epsilon_vector) @@ -135,13 +135,13 @@ def makeXMatrix(b, s_vector, d_vector, epsilon_vector): return x -# %% ../nbs/src/tbats.ipynb 13 +# %% ../../nbs/src/tbats.ipynb 13 def findPQ(ar_coeffs, ma_coeffs): p = 0 if ar_coeffs is None else len(ar_coeffs) q = 0 if ma_coeffs is None else len(ma_coeffs) return p, q -# %% ../nbs/src/tbats.ipynb 15 +# %% ../../nbs/src/tbats.ipynb 15 def makeTBATSWMatrix(phi, k_vector, ar_coeffs, ma_coeffs, tau): # w_transpose = (1, phi, a, varphi, theta) p, q = findPQ(ar_coeffs, ma_coeffs) @@ -174,7 +174,7 @@ def makeTBATSWMatrix(phi, k_vector, ar_coeffs, ma_coeffs, tau): return w_transpose -# %% ../nbs/src/tbats.ipynb 17 +# %% ../../nbs/src/tbats.ipynb 17 def makeTBATSGMatrix( k_vector, alpha, adj_beta, beta, gamma_one_v, gamma_two_v, p, q, tau ): @@ -203,7 +203,7 @@ def makeTBATSGMatrix( return g, gamma_bold -# %% ../nbs/src/tbats.ipynb 19 +# %% ../../nbs/src/tbats.ipynb 19 @njit(nogil=NOGIL, cache=CACHE) def makeTBATSFMatrix( phi, tau, alpha, beta, ar_coeffs, ma_coeffs, gamma_bold, seasonal_periods, k_vector @@ -316,7 +316,7 @@ def makeTBATSFMatrix( return F -# %% ../nbs/src/tbats.ipynb 21 +# %% ../../nbs/src/tbats.ipynb 21 @njit(nogil=NOGIL, cache=CACHE) def calcTBATSFaster(y_trans, w_transpose, g, F, x_nought): @@ -337,7 +337,7 @@ def calcTBATSFaster(y_trans, w_transpose, g, F, x_nought): return yhat, e, x -# %% ../nbs/src/tbats.ipynb 23 +# %% ../../nbs/src/tbats.ipynb 23 def extract_params( params, use_boxcox, @@ -393,7 +393,7 @@ def extract_params( ma_coeffs, ) -# %% ../nbs/src/tbats.ipynb 25 +# %% ../../nbs/src/tbats.ipynb 25 def updateTBATSWMatrix(w_transpose, phi, tau, ar_coeffs, ma_coeffs, p, q): adjBeta = 0 @@ -413,7 +413,7 @@ def updateTBATSWMatrix(w_transpose, phi, tau, ar_coeffs, ma_coeffs, p, q): return w_transpose -# %% ../nbs/src/tbats.ipynb 27 +# %% ../../nbs/src/tbats.ipynb 27 def updateTBATSGMatrix(g, gamma_bold, alpha, beta, k_vector, gamma_one_v, gamma_two_v): # This function also updates gamma_bold adjBeta = 0 @@ -433,7 +433,7 @@ def updateTBATSGMatrix(g, gamma_bold, alpha, beta, k_vector, gamma_one_v, gamma_ return g -# %% ../nbs/src/tbats.ipynb 29 +# %% ../../nbs/src/tbats.ipynb 29 def updateTBATSFMatrix( F, phi, alpha, beta, gamma_bold, ar_coeffs, ma_coeffs, p, q, tau ): @@ -478,7 +478,7 @@ def updateTBATSFMatrix( return F -# %% ../nbs/src/tbats.ipynb 31 +# %% ../../nbs/src/tbats.ipynb 31 def checkAdmissibility( BoxCox_lambda, bc_lower_bound, bc_upper_bound, phi, ar_coeffs, ma_coeffs, D ): @@ -512,7 +512,7 @@ def checkAdmissibility( return np.all(abs(D_eigen_values) < 1 + 1e-2) -# %% ../nbs/src/tbats.ipynb 33 +# %% ../../nbs/src/tbats.ipynb 33 def calcLikelihoodTBATS( params, use_boxcox, @@ -579,7 +579,7 @@ def calcLikelihoodTBATS( else: return 10**20 -# %% ../nbs/src/tbats.ipynb 36 +# %% ../../nbs/src/tbats.ipynb 36 def tbats_model_generator( y, seasonal_periods, @@ -838,7 +838,7 @@ def tbats_model_generator( return res -# %% ../nbs/src/tbats.ipynb 38 +# %% ../../nbs/src/tbats.ipynb 38 def tbats_model( y, seasonal_periods, @@ -919,7 +919,7 @@ def tbats_model( return best_model -# %% ../nbs/src/tbats.ipynb 40 +# %% ../../nbs/src/tbats.ipynb 40 def tbats_selection( y, seasonal_periods, @@ -1004,7 +1004,7 @@ def tbats_selection( return mod -# %% ../nbs/src/tbats.ipynb 42 +# %% ../../nbs/src/tbats.ipynb 42 def tbats_forecast(mod, h): # this function is the same as bats_forecast fcst = np.zeros(h) xx = np.zeros((h, mod["x"].shape[1])) @@ -1022,7 +1022,7 @@ def tbats_forecast(mod, h): # this function is the same as bats_forecast return res -# %% ../nbs/src/tbats.ipynb 43 +# %% ../../nbs/src/tbats.ipynb 43 def _compute_sigmah(obj, h): """ Computes the sigmah requiered for prediction intervals diff --git a/statsforecast/theta.py b/python/statsforecast/theta.py similarity index 95% rename from statsforecast/theta.py rename to python/statsforecast/theta.py index c1668bd81..9999ae6ef 100644 --- a/statsforecast/theta.py +++ b/python/statsforecast/theta.py @@ -1,9 +1,9 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/theta.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/theta.ipynb. # %% auto 0 __all__ = ['theta_target_fn', 'is_constant'] -# %% ../nbs/src/theta.ipynb 1 +# %% ../../nbs/src/theta.ipynb 1 import math from typing import Tuple @@ -13,10 +13,16 @@ from statsmodels.tsa.seasonal import seasonal_decompose from statsmodels.tsa.stattools import acf -from .ets import restrict_to_bounds, results -from .utils import _seasonal_naive, _repeat_val_seas, CACHE, NOGIL +from statsforecast.utils import ( + CACHE, + NOGIL, + _repeat_val_seas, + _seasonal_naive, + restrict_to_bounds, + results, +) -# %% ../nbs/src/theta.ipynb 4 +# %% ../../nbs/src/theta.ipynb 4 # Global variables STM = 0 OTM = 1 @@ -27,7 +33,7 @@ NA = -99999.0 smalno = np.finfo(float).eps -# %% ../nbs/src/theta.ipynb 6 +# %% ../../nbs/src/theta.ipynb 6 @njit(nogil=NOGIL, cache=CACHE) def initstate(y, modeltype, initial_smoothed, alpha, theta): states = np.zeros((1, 5), dtype=np.float32) @@ -53,7 +59,7 @@ def initstate(y, modeltype, initial_smoothed, alpha, theta): return states -# %% ../nbs/src/theta.ipynb 8 +# %% ../../nbs/src/theta.ipynb 8 @njit(nogil=NOGIL, cache=CACHE) def thetacalc( y: np.ndarray, @@ -116,7 +122,7 @@ def thetacalc( mse = np.sum(e[3:] ** 2) / mean_y return mse -# %% ../nbs/src/theta.ipynb 9 +# %% ../../nbs/src/theta.ipynb 9 @njit(nogil=NOGIL, cache=CACHE) def thetafcst(states, i, modeltype, f, h, alpha, theta): # obs: @@ -137,7 +143,7 @@ def thetafcst(states, i, modeltype, f, h, alpha, theta): ) f[i_h] = new_states[i + i_h, 4] # mu is the forecast -# %% ../nbs/src/theta.ipynb 10 +# %% ../../nbs/src/theta.ipynb 10 @njit(nogil=NOGIL, cache=CACHE) def thetaupdate(states, i, modeltype, alpha, theta, y, usemu): # kind of model # states @@ -166,7 +172,7 @@ def thetaupdate(states, i, modeltype, alpha, theta, y, usemu): # kind of model states[i, 2] = An states[i, 3] = Bn -# %% ../nbs/src/theta.ipynb 11 +# %% ../../nbs/src/theta.ipynb 11 @njit(nogil=NOGIL, cache=CACHE) def thetaforecast(states, n, modeltype, f, h, alpha, theta): # compute forecasts @@ -175,7 +181,7 @@ def thetaforecast(states, n, modeltype, f, h, alpha, theta): ) return new_states -# %% ../nbs/src/theta.ipynb 16 +# %% ../../nbs/src/theta.ipynb 16 @njit(nogil=NOGIL, cache=CACHE) def initparamtheta( initial_smoothed: float, alpha: float, theta: float, y: np.ndarray, modeltype: str @@ -218,12 +224,12 @@ def initparamtheta( "optimize_theta": optimize_theta, } -# %% ../nbs/src/theta.ipynb 18 +# %% ../../nbs/src/theta.ipynb 18 @njit(nogil=NOGIL, cache=CACHE) def switch_theta(x: str): return {"STM": 0, "OTM": 1, "DSTM": 2, "DOTM": 3}[x] -# %% ../nbs/src/theta.ipynb 20 +# %% ../../nbs/src/theta.ipynb 20 @njit(nogil=NOGIL, cache=CACHE) def pegelsresid_theta( y: np.ndarray, @@ -252,7 +258,7 @@ def pegelsresid_theta( mse = np.nan return amse, e, states, mse -# %% ../nbs/src/theta.ipynb 21 +# %% ../../nbs/src/theta.ipynb 21 @njit(nogil=NOGIL, cache=CACHE) def theta_target_fn( optimal_param, @@ -307,7 +313,7 @@ def theta_target_fn( mse = -np.inf return mse -# %% ../nbs/src/theta.ipynb 22 +# %% ../../nbs/src/theta.ipynb 22 @njit(nogil=NOGIL, cache=CACHE) def nelder_mead_theta( x0: np.ndarray, @@ -423,7 +429,7 @@ def nelder_mead_theta( f_simplex[i] = theta_target_fn(simplex[i], *args) return results(simplex[best_idx], f_simplex[best_idx], it + 1, simplex) -# %% ../nbs/src/theta.ipynb 23 +# %% ../../nbs/src/theta.ipynb 23 def optimize_theta_target_fn(init_par, optimize_params, y, modeltype, nmse): x0 = [init_par[key] for key, val in optimize_params.items() if val] x0 = np.array(x0, dtype=np.float32) @@ -459,12 +465,12 @@ def optimize_theta_target_fn(init_par, optimize_params, y, modeltype, nmse): ) return res -# %% ../nbs/src/theta.ipynb 24 +# %% ../../nbs/src/theta.ipynb 24 @njit(nogil=NOGIL, cache=CACHE) def is_constant(x): return np.all(x[0] == x) -# %% ../nbs/src/theta.ipynb 26 +# %% ../../nbs/src/theta.ipynb 26 def thetamodel( y: np.ndarray, m: int, @@ -521,7 +527,7 @@ def thetamodel( mean_y=np.mean(y), ) -# %% ../nbs/src/theta.ipynb 28 +# %% ../../nbs/src/theta.ipynb 28 def compute_pi_samples( n, h, states, sigma, alpha, theta, mean_y, seed=0, n_samples=200 ): @@ -540,7 +546,7 @@ def compute_pi_samples( A = mean_y - B * (i + 2) / 2 return samples -# %% ../nbs/src/theta.ipynb 29 +# %% ../../nbs/src/theta.ipynb 29 def forecast_theta(obj, h, level=None): forecast = np.full(h, fill_value=np.nan) n = obj["n"] @@ -585,7 +591,7 @@ def forecast_theta(obj, h, level=None): res[key] = res[key] + seas_forecast return res -# %% ../nbs/src/theta.ipynb 31 +# %% ../../nbs/src/theta.ipynb 31 def auto_theta( y, m, @@ -685,7 +691,7 @@ def auto_theta( model["seas_forecast"] = dict(seas_forecast) return model -# %% ../nbs/src/theta.ipynb 41 +# %% ../../nbs/src/theta.ipynb 41 def forward_theta(fitted_model, y): m = fitted_model["m"] model = fitted_model["modeltype"] diff --git a/statsforecast/utils.py b/python/statsforecast/utils.py similarity index 90% rename from statsforecast/utils.py rename to python/statsforecast/utils.py index e20bf6f7b..22aad491e 100644 --- a/statsforecast/utils.py +++ b/python/statsforecast/utils.py @@ -1,24 +1,26 @@ -# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/utils.ipynb. +# AUTOGENERATED! DO NOT EDIT! File to edit: ../../nbs/src/utils.ipynb. # %% auto 0 __all__ = ['AirPassengers', 'AirPassengersDF', 'generate_series'] -# %% ../nbs/src/utils.ipynb 3 +# %% ../../nbs/src/utils.ipynb 3 import inspect import math import os import warnings +from collections import namedtuple from functools import wraps from typing import Dict import numpy as np import pandas as pd +from numba import njit from scipy.stats import norm from utilsforecast.compat import DataFrame from utilsforecast.data import generate_series as utils_generate_series -# %% ../nbs/src/utils.ipynb 4 +# %% ../../nbs/src/utils.ipynb 4 # Global variables NOGIL = bool(os.getenv("NIXTLA_NUMBA_RELEASE_GIL", "")) LEGACY_CACHE = bool(os.getenv("NUMBA_CACHE", "")) @@ -29,8 +31,24 @@ FutureWarning, ) CACHE = bool(os.getenv("NIXTLA_NUMBA_CACHE", "")) or LEGACY_CACHE +results = namedtuple("results", "x fn nit simplex") -# %% ../nbs/src/utils.ipynb 7 +# %% ../../nbs/src/utils.ipynb 5 +@njit(nogil=NOGIL, cache=CACHE) +def restrict_to_bounds(x, lower, upper): + new_x = np.full_like(x, fill_value=np.nan, dtype=x.dtype) + for i in range(x.size): + lo = lower[i] + up = upper[i] + if x[i] < lo: + new_x[i] = lo + elif x[i] > up: + new_x[i] = up + else: + new_x[i] = x[i] + return new_x + +# %% ../../nbs/src/utils.ipynb 8 def generate_series( n_series: int, freq: str = "D", @@ -82,7 +100,7 @@ def generate_series( seed=seed, ) -# %% ../nbs/src/utils.ipynb 11 +# %% ../../nbs/src/utils.ipynb 12 AirPassengers = np.array( [ 112.0, @@ -232,7 +250,7 @@ def generate_series( ] ) -# %% ../nbs/src/utils.ipynb 12 +# %% ../../nbs/src/utils.ipynb 13 AirPassengersDF = pd.DataFrame( { "unique_id": np.ones(len(AirPassengers)), @@ -243,7 +261,7 @@ def generate_series( } ) -# %% ../nbs/src/utils.ipynb 17 +# %% ../../nbs/src/utils.ipynb 18 def _repeat_val_seas(season_vals: np.ndarray, h: int) -> np.ndarray: repeats = math.ceil(h / season_vals.size) return np.tile(season_vals, repeats)[:h] @@ -292,7 +310,7 @@ def _ensure_float(x: np.ndarray) -> np.ndarray: x = x.astype(np.float32) return x -# %% ../nbs/src/utils.ipynb 19 +# %% ../../nbs/src/utils.ipynb 20 # Functions used for calculating prediction intervals def _quantiles(level): level = np.asarray(level) @@ -322,7 +340,7 @@ def _calculate_sigma(residuals, n): sigma = 0 return sigma -# %% ../nbs/src/utils.ipynb 20 +# %% ../../nbs/src/utils.ipynb 21 class ConformalIntervals: """Class for storing conformal intervals metadata information.""" @@ -343,7 +361,7 @@ def __init__( self.h = h self.method = method -# %% ../nbs/src/utils.ipynb 21 +# %% ../../nbs/src/utils.ipynb 22 def _old_kw_to_pos(old_names, new_positions): def decorator(f): @wraps(f) diff --git a/settings.ini b/settings.ini index b6ad4f962..83ac737ac 100644 --- a/settings.ini +++ b/settings.ini @@ -8,7 +8,7 @@ author = Nixtla author_email = business@nixtla.io copyright = Nixtla Inc. branch = main -version = 1.7.6 +version = 1.7.6.99 min_python = 3.8 audience = Developers language = English @@ -21,14 +21,14 @@ ray_requirements = fugue[ray]>=0.8.1 protobuf>=3.15.3,<4.0.0 numpy<2 pandas<2.2 dask_requirements = fugue[dask]>=0.8.1 spark_requirements = fugue[spark]>=0.8.1 plotly_requirements = plotly plotly-resampler -dev_requirements = black datasetsforecast fire nbdev nbformat nbdev_plotly pandas[plot] pmdarima polars[numpy] pre-commit prophet pyarrow scikit-learn setuptools<70 supersmoother +dev_requirements = black datasetsforecast fire nbdev nbformat nbdev_plotly pandas[plot] pmdarima polars[numpy] pre-commit prophet pyarrow pybind11 scikit-learn setuptools<70 supersmoother nbs_path = nbs doc_path = _docs recursive = True doc_host = https://%(user)s.github.io doc_baseurl = /%(lib_name)s/ git_url = https://github.com/%(user)s/statsforecast/ -lib_path = %(lib_name)s +lib_path = python/statsforecast title = %(lib_name)s black_formatting = True jupyter_hooks = True diff --git a/setup.py b/setup.py index 98298c50b..34987d13b 100644 --- a/setup.py +++ b/setup.py @@ -1,8 +1,9 @@ +import glob import sys -from pkg_resources import parse_version -from configparser import ConfigParser + import setuptools -assert parse_version(setuptools.__version__)>=parse_version('36.2') +from configparser import ConfigParser +from pybind11.setup_helpers import Pybind11Extension # note: all settings are in settings.ini; edit there, not here config = ConfigParser(delimiters=['=']) @@ -52,6 +53,15 @@ if sys.version_info < (3, 12): all_requirements.extend(ray_requirements) +ext_modules = [ + Pybind11Extension( + name="statsforecast._lib", + sources=glob.glob("src/*.cpp"), + include_dirs=["include/statsforecast", "external_libs/eigen"], + cxx_std=17, + ) +] + setuptools.setup( name = 'statsforecast', license = lic[0], @@ -61,7 +71,8 @@ 'Natural Language :: ' + cfg['language'].title(), ] + ['Programming Language :: Python :: '+o for o in py_versions[py_versions.index(min_python):]] + (['License :: ' + lic[1] ] if lic[1] else []), url = cfg['git_url'], - packages = setuptools.find_packages(), + package_dir={"": "python"}, + packages = setuptools.find_packages(where="python"), include_package_data = True, install_requires = requirements, extras_require={ @@ -80,6 +91,8 @@ zip_safe = False, entry_points = { 'console_scripts': cfg.get('console_scripts','').split(), - 'nbdev': [f'{cfg.get("lib_path")}={cfg.get("lib_path")}._modidx:d'] + 'nbdev': ['statsforecast=statsforecast._modidx:d'] }, - **setup_cfg) + ext_modules=ext_modules, + **setup_cfg +) diff --git a/src/ets.cpp b/src/ets.cpp new file mode 100644 index 000000000..bc577a1eb --- /dev/null +++ b/src/ets.cpp @@ -0,0 +1,327 @@ +#include + +#include "nelder_mead.h" + +namespace ets { +namespace py = pybind11; +using Eigen::VectorXd; + +enum class Component { + Nothing = 0, + Additive = 1, + Multiplicative = 2, +}; +enum class Criterion { + Likelihood = 0, + MSE = 1, + AMSE = 2, + Sigma = 3, + MAE = 4, +}; +constexpr double HUGE_N = 1e10; +constexpr double NA = -99999.0; +constexpr double TOL = 1e-10; + +template +std::tuple +Update(Ref s, double l, double b, double old_l, double old_b, CRef old_s, int m, + Component trend, Component season, double alpha, double beta, + double gamma, double phi, double y) { + double q, phi_b; + // new level + if (trend == Component::Nothing) { + q = old_l; + phi_b = 0.0; + } else if (trend == Component::Additive) { + phi_b = phi * old_b; + q = old_l + phi_b; + } else if (std::abs(phi - 1.0) < TOL) { + phi_b = old_b; + q = old_l * old_b; + } else { + phi_b = std::pow(old_b, phi); + q = old_l * phi_b; + } + // season + double p; + if (season == Component::Nothing) { + p = y; + } else if (season == Component::Additive) { + p = y - old_s[m - 1]; + } else { + if (std::abs(old_s[m - 1]) < TOL) { + p = HUGE_N; + } else { + p = y / old_s[m - 1]; + } + } + l = q + alpha * (p - q); + // new growth + double r; + if (trend != Component::Nothing) { + if (trend == Component::Additive) { + r = l - old_l; + } else { + if (std::abs(old_l) < TOL) { + r = HUGE_N; + } else { + r = l / old_l; + } + } + b = phi_b + (beta / alpha) * (r - phi_b); + } + // new seasonal + double t; + if (season != Component::Nothing) { + if (season == Component::Additive) { + t = y - q; + } else { + if (std::abs(q) < TOL) { + t = HUGE_N; + } else { + t = y / q; + } + } + s[0] = old_s[m - 1] + gamma * (t - old_s[m - 1]); + std::copy(old_s.data(), old_s.data() + m - 1, s.data() + 1); + } + return {l, b}; +} + +template +void Forecast(Ref f, double l, double b, CRef s, int m, Component trend, + Component season, double phi, int h) { + double phistar = phi; + for (int i = 0; i < h; ++i) { + if (trend == Component::Nothing) { + f[i] = l; + } else if (trend == Component::Additive) { + f[i] = l + phistar * b; + } else if (b < 0) { + f[i] = std::numeric_limits::quiet_NaN(); + } else { + f[i] = l * std::pow(b, phistar); + } + int j = m - 1 - i; + while (j < 0) { + j += m; + } + if (season == Component::Additive) { + f[i] += s[j]; + } else if (season == Component::Multiplicative) { + f[i] *= s[j]; + } + if (i < h - 1) { + if (std::abs(phi - 1.0) < TOL) { + phistar += 1.0; + } else { + phistar += std::pow(phi, i + 1); + } + } + } +} + +template +double Calc(Ref x, Ref e, Ref a_mse, int n_mse, CRef y, Component error, + Component trend, Component season, double alpha, double beta, + double gamma, double phi, int m) { + auto n = y.size(); + int n_s = std::max(m, 24); + m = std::max(m, 1); + n_mse = std::min(n_mse, 30); + int n_states = + m * (season != Component::Nothing) + (trend != Component::Nothing) + 1; + + // copy initial state components + double l = x[0]; + double b = (trend != Component::Nothing) ? x[1] : 0.0; + auto s = std::vector(n_s); + if (season != Component::Nothing) { + std::copy(x.data() + 1 + (trend != Component::Nothing), + x.data() + 1 + (trend != Component::Nothing) + m, s.data()); + } + + std::fill(a_mse.data(), a_mse.data() + n_mse, 0.0); + auto old_s = std::vector(n_s); + auto denom = std::vector(30); + auto f = std::vector(30); + double old_b = 0.0; + double lik = 0.0; + double lik2 = 0.0; + double f_0, old_l, val; + for (Eigen::Index i = 0; i < n; ++i) { + // copy previous state + old_l = l; + if (trend != Component::Nothing) { + old_b = b; + } + if (season != Component::Nothing) { + std::copy(s.begin(), s.begin() + m, old_s.begin()); + } + // one step forecast + Forecast &, const std::vector &>( + f, old_l, old_b, old_s, m, trend, season, phi, n_mse); + if (std::abs(f[0] - NA) < TOL) { + return NA; + } + if (error == Component::Additive) { + e[i] = y[i] - f[0]; + } else { + if (std::abs(f[0]) < TOL) { + f_0 = f[0] + TOL; + } else { + f_0 = f[0]; + } + e[i] = (y[i] - f[0]) / f_0; + } + for (Eigen::Index j = 0; j < n_mse; ++j) { + if (i + j < y.size()) { + denom[j] += 1.0; + double tmp = y[i + j] - f[j]; + a_mse[j] = (a_mse[j] * (denom[j] - 1.0) + tmp * tmp) / denom[j]; + } + } + // update state + std::tie(l, b) = Update &, const std::vector &>( + s, l, b, old_l, old_b, old_s, m, trend, season, alpha, beta, gamma, phi, + y[i]); + + // store new state + x[n_states * (i + 1)] = l; + if (trend != Component::Nothing) { + x[n_states * (i + 1) + 1] = b; + } + if (season != Component::Nothing) { + std::copy(s.data(), s.data() + m, + x.data() + n_states * (i + 1) + 1 + + (trend != Component::Nothing)); + } + lik += e[i] * e[i]; + val = std::abs(f[0]); + if (val > 0.0) { + lik2 += std::log(val); + } else { + lik2 += std::log(val + 1e-8); + } + } + if (lik > 0.0) { + lik = static_cast(n) * std::log(lik); + } else { + lik = static_cast(n) * std::log(lik + 1e-8); + } + if (error == Component::Multiplicative) { + lik += 2 * lik2; + } + return lik; +} + +double ObjectiveFunction(const VectorXd ¶ms, const VectorXd &y, int n_state, + Component error, Component trend, Component season, + Criterion opt_crit, int n_mse, int m, bool opt_alpha, + bool opt_beta, bool opt_gamma, bool opt_phi, + double alpha, double beta, double gamma, double phi) { + int j = 0; + if (opt_alpha) { + alpha = params(j++); + } + if (opt_beta) { + beta = params(j++); + } + if (opt_gamma) { + gamma = params(j++); + } + if (opt_phi) { + phi = params(j++); + } + auto n_params = params.size(); + auto n = y.size(); + int p = n_state + (season != Component::Nothing); + VectorXd state = VectorXd::Zero(p * (n + 1)); + std::copy(params.data() + n_params - n_state, params.data() + n_params, + state.data()); + if (season != Component::Nothing) { + // add extra state + int start = 1 + (trend != Component::Nothing); + double sum = state(Eigen::seq(start, n_state - 1)).sum(); + state(n_state) = + static_cast(m * (season == Component::Multiplicative)) - sum; + if (season == Component::Multiplicative && + state.tail(state.size() - start).minCoeff() < 0.0) { + return std::numeric_limits::infinity(); + } + } + VectorXd a_mse = VectorXd::Zero(30); + VectorXd e = VectorXd::Zero(n); + double lik = Calc(state, e, a_mse, n_mse, y, + error, trend, season, alpha, + beta, gamma, phi, m); + lik = std::max(lik, -1e10); + if (std::isnan(lik) || std::abs(lik + 99999.0) < 1e-7) { + lik = -std::numeric_limits::infinity(); + } + double obj_val = 0.0; + switch (opt_crit) { + case Criterion::Likelihood: + obj_val = lik; + break; + case Criterion::MSE: + obj_val = a_mse(0); + break; + case Criterion::AMSE: + obj_val = a_mse.head(n_mse).mean(); + break; + case Criterion::Sigma: + obj_val = e.array().square().mean(); + break; + case Criterion::MAE: + obj_val = e.array().abs().mean(); + break; + } + return obj_val; +} +std::tuple +Optimize(const Eigen::Ref &x0, + const Eigen::Ref &y, int n_state, Component error, + Component trend, Component season, Criterion opt_crit, int n_mse, + int m, bool opt_alpha, bool opt_beta, bool opt_gamma, bool opt_phi, + double alpha, double beta, double gamma, double phi, + const Eigen::Ref &lower, + const Eigen::Ref &upper, double tol_std, int max_iter, + bool adaptive) { + double init_step = 0.05; + double nm_alpha = 1.0; + double nm_gamma = 2.0; + double nm_rho = 0.5; + double nm_sigma = 0.5; + double zero_pert = 1.0e-4; + return nm::NelderMead(ObjectiveFunction, x0, lower, upper, init_step, + zero_pert, nm_alpha, nm_gamma, nm_rho, nm_sigma, + max_iter, tol_std, adaptive, y, n_state, error, trend, + season, opt_crit, n_mse, m, opt_alpha, opt_beta, + opt_gamma, opt_phi, alpha, beta, gamma, phi); +} + +void init(py::module_ &m) { + py::module_ ets = m.def_submodule("ets"); + ets.attr("HUGE_N") = HUGE_N; + ets.attr("NA") = NA; + ets.attr("TOL") = TOL; + py::enum_(ets, "Component") + .value("Nothing", Component::Nothing) + .value("Additive", Component::Additive) + .value("Multiplicative", Component::Multiplicative); + py::enum_(ets, "Criterion") + .value("Likelihood", Criterion::Likelihood) + .value("MSE", Criterion::MSE) + .value("AMSE", Criterion::AMSE) + .value("Sigma", Criterion::Sigma) + .value("MAE", Criterion::MAE); + ets.def("update", + &Update, const Eigen::Ref &>); + ets.def("forecast", + &Forecast, const Eigen::Ref &>); + ets.def("calc", + &Calc, const Eigen::Ref &>); + ets.def("optimize", &Optimize); +} +} // namespace ets diff --git a/src/statsforecast.cpp b/src/statsforecast.cpp new file mode 100644 index 000000000..f240b4a1f --- /dev/null +++ b/src/statsforecast.cpp @@ -0,0 +1,11 @@ +#include + +namespace py = pybind11; + +namespace ets { +void init(py::module_ &); +} + +PYBIND11_MODULE(_lib, m) { + ets::init(m); +}