diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml index a81379493..40fb6df9a 100644 --- a/.github/workflows/lint.yml +++ b/.github/workflows/lint.yml @@ -13,31 +13,31 @@ permissions: contents: read jobs: - super-linter: - name: Lint Code Base - runs-on: ubuntu-latest - permissions: - contents: read - pull-requests: write - actions: read - checks: read - statuses: write - steps: - - name: Checkout Code - uses: actions/checkout@v4 - with: - # Full git history is needed to get a proper - # list of changed files within `super-linter` - fetch-depth: 0 - - - name: Lint Code Base - uses: super-linter/super-linter/slim@v6 - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - VALIDATE_ALL_CODEBASE: false - DEFAULT_BRANCH: master - LINTER_RULES_PATH: . - MARKDOWN_CONFIG_FILE: .markdownlint.json +## super-linter: +## name: Lint Code Base +## runs-on: ubuntu-latest +## permissions: +## contents: read +## pull-requests: write +## actions: read +## checks: read +## statuses: write +## steps: +## - name: Checkout Code +## uses: actions/checkout@v4 +## with: +## # Full git history is needed to get a proper +## # list of changed files within `super-linter` +## fetch-depth: 0 +## +## - name: Lint Code Base +## uses: super-linter/super-linter/slim@v6.7.0 +## env: +## GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} +## VALIDATE_ALL_CODEBASE: false +## DEFAULT_BRANCH: main +## LINTER_RULES_PATH: . +## MARKDOWN_CONFIG_FILE: .markdownlint.json docs: runs-on: ubuntu-latest diff --git a/Dockerfile b/Dockerfile index 1dd45ed1a..790faeb28 100644 --- a/Dockerfile +++ b/Dockerfile @@ -77,3 +77,7 @@ EOF # install Docker tools (cli, buildx, compose) COPY --from=gloursdocker/docker / / + +USER ${USERNAME} +HEALTHCHECK CMD true + diff --git a/consts/README.md b/consts/README.md index 4f35eeb22..fa566cfdd 100644 --- a/consts/README.md +++ b/consts/README.md @@ -46,6 +46,12 @@ consts.mksa_plancks_constant_hbar Planck's constant divided by `2\pi`, `\hbar`. +```console +consts.mksa_planck_temperature +``` + +Planck temperature `T_p`. + ```console consts.num_avogadro ``` diff --git a/consts/cgs.v b/consts/cgs.v index 512694b8a..146d88977 100644 --- a/consts/cgs.v +++ b/consts/cgs.v @@ -8,6 +8,8 @@ pub const cgs_plancks_constant_h = 6.62606896e-27 // g cm^2 / s pub const cgs_plancks_constant_hbar = 1.05457162825e-27 // g cm^2 / s +pub const cgs_planck_temperature = 1.416785e+32 // Kelvin + pub const cgs_astronomical_unit = 1.49597870691e+13 // cm pub const cgs_light_year = 9.46053620707e+17 // cm diff --git a/consts/cgsm.v b/consts/cgsm.v index 8b2d4b0d0..7a2c18aff 100644 --- a/consts/cgsm.v +++ b/consts/cgsm.v @@ -8,6 +8,8 @@ pub const cgsm_plancks_constant_h = 6.62606896e-27 // g cm^2 / s pub const cgsm_plancks_constant_hbar = 1.05457162825e-27 // g cm^2 / s +pub const cgsm_planck_temperature = 1.416785e+32 // Kelvin + pub const cgsm_astronomical_unit = 1.49597870691e+13 // cm pub const cgsm_light_year = 9.46053620707e+17 // cm diff --git a/consts/mks.v b/consts/mks.v index 0de00d50d..8fedb081d 100644 --- a/consts/mks.v +++ b/consts/mks.v @@ -8,6 +8,8 @@ pub const mks_plancks_constant_h = 6.62606896e-34 // kg m^2 / s pub const mks_plancks_constant_hbar = 1.05457162825e-34 // kg m^2 / s +pub const mks_planck_temperature = 1.416785e+32 // Kelvin + pub const mks_astronomical_unit = 1.49597870691e+11 // m pub const mks_light_year = 9.46053620707e+15 // m diff --git a/consts/mksa.v b/consts/mksa.v index a8a48d7ef..3e2da5c8a 100644 --- a/consts/mksa.v +++ b/consts/mksa.v @@ -8,6 +8,8 @@ pub const mksa_plancks_constant_h = 6.62606896e-34 // kg m^2 / s pub const mksa_plancks_constant_hbar = 1.05457162825e-34 // kg m^2 / s +pub const mksa_planck_temperature = 1.416785e+32 // Kelvin + pub const mksa_astronomical_unit = 1.49597870691e+11 // m pub const mksa_light_year = 9.46053620707e+15 // m diff --git a/deriv/README.md b/deriv/README.md index 096613e8d..b66312a63 100644 --- a/deriv/README.md +++ b/deriv/README.md @@ -63,6 +63,18 @@ for values greater than `x`. πŸ“‰ This function is equivalent to calling `deriv.forward` with a negative step-size. +### `partial` + +```v ignore +pub fn partial(f fn ([]f64) f64, x []f64, variable int, h f64) (f64, f64) +``` + +This function computes the partial derivative of the function `f` with respect to +a specified variable at point `x` using step-size `h`. It returns the derivative +in `result` and an error estimate in `abserr`. The function `f` should take an array +of coordinates and return a single value. This method provides both the derivative +and its error estimate. + ## References and Further Reading This work is a spiritual descendent of the Differentiation module in [GSL](https://github.com/ampl/gsl). πŸ“– diff --git a/deriv/deriv.v b/deriv/deriv.v index 47a2b33cd..eb53c0758 100644 --- a/deriv/deriv.v +++ b/deriv/deriv.v @@ -1,6 +1,7 @@ module deriv import vsl.func +import vsl.errors import vsl.internal.prec import math @@ -113,3 +114,32 @@ pub fn forward(f func.Fn, x f64, h f64) (f64, f64) { pub fn backward(f func.Fn, x f64, h f64) (f64, f64) { return forward(f, x, -h) } + +/** +* partial is a function that computes the partial derivative of a multivariable function with respect to a specified variable. +* +* @param f The multivariable function for which the partial derivative is to be computed. +* @param x The point at which the partial derivative is to be computed, represented as an array of coordinates. +* @param variable The index of the variable with respect to which the partial derivative is to be computed. +* @param h The step size to be used in the central difference method. +* +* @return A tuple containing the value of the partial derivative and the estimated error. +*/ +pub fn partial(f fn ([]f64) f64, x []f64, variable int, h f64) (f64, f64) { + if variable < 0 || variable >= x.len { + errors.vsl_panic('Invalid variable index', .efailed) + } + + // Define a helper function that converts the multivariate function + // to a univariate function for the specified variable + partial_helper := func.Fn{ + f: fn [f, x, variable] (t f64, _ []f64) f64 { + mut x_new := x.clone() + x_new[variable] = t + return f(x_new) + } + } + + // Use the central difference method to compute the partial derivative + return central(partial_helper, x[variable], h) +} diff --git a/deriv/deriv_test.v b/deriv/deriv_test.v index f25a5573f..041c53d79 100644 --- a/deriv/deriv_test.v +++ b/deriv/deriv_test.v @@ -68,6 +68,18 @@ fn df6(x f64, _ []f64) f64 { return -1.0 / (x * x) } +fn f_multi(x []f64) f64 { + return x[0] * x[0] + x[1] * x[1] // f(x,y) = x^2 + y^2 +} + +fn df_multi_dx(x []f64) f64 { + return 2 * x[0] // βˆ‚f/βˆ‚x = 2x +} + +fn df_multi_dy(x []f64) f64 { + return 2 * x[1] // βˆ‚f/βˆ‚y = 2y +} + fn test_deriv() { f1_ := func.Fn.new(f: f1) df1_ := func.Fn.new(f: df1) @@ -100,6 +112,18 @@ fn test_deriv() { assert deriv_test('central', f6_, df6_, 10.0) assert deriv_test('forward', f6_, df6_, 10.0) assert deriv_test('backward', f6_, df6_, 10.0) + + // Partial derivative test + x := [2.0, 3.0] + h := 1e-5 + + // Partial derivative with respect to x + dx, _ := partial(f_multi, x, 0, h) + assert float64.tolerance(dx, df_multi_dx(x), 1e-5) + + // Partial derivative with respect to y + dy, _ := partial(f_multi, x, 1, h) + assert float64.tolerance(dy, df_multi_dy(x), 1e-5) } fn deriv_test(deriv_method string, f func.Fn, df func.Fn, x f64) bool { diff --git a/errors/errno.v b/errors/errno.v index accf1083c..d958ee0b6 100644 --- a/errors/errno.v +++ b/errors/errno.v @@ -2,73 +2,73 @@ module errors pub enum ErrorCode { // success - success = 0 + success = 0 // generic failure - failure = -1 + failure = -1 // iteration has not converged can_continue = -2 // input domain error, e.g sqrt(-1) - edom = 1 + edom = 1 // output range error, e.g. exp(1e+100) - erange = 2 + erange = 2 // invalid pointer - efault = 3 + efault = 3 // invalid argument supplied by user - einval = 4 + einval = 4 // generic failure - efailed = 5 + efailed = 5 // factorization failed - efactor = 6 + efactor = 6 // sanity check failed - shouldn't happen - esanity = 7 + esanity = 7 // malloc failed - enomem = 8 + enomem = 8 // problem with user-supplied function - ebadfunc = 9 + ebadfunc = 9 // iterative process is out of control - erunaway = 10 + erunaway = 10 // exceeded max number of iterations - emaxiter = 11 + emaxiter = 11 // tried to divide by zero - ezerodiv = 12 + ezerodiv = 12 // user specified an invalid tolerance - ebadtol = 13 + ebadtol = 13 // failed to reach the specified tolerance - etol = 14 + etol = 14 // underflow - eundrflw = 15 + eundrflw = 15 // overflow - eovrflw = 16 + eovrflw = 16 // loss of accuracy - eloss = 17 + eloss = 17 // failed because of roundoff error - eround = 18 + eround = 18 // matrix, vector lengths are not conformant - ebadlen = 19 + ebadlen = 19 // matrix not square - enotsqr = 20 + enotsqr = 20 // apparent singularity detected - esing = 21 + esing = 21 // integral or series is divergent - ediverge = 22 + ediverge = 22 // requested feature is not supported by the hardware - eunsup = 23 + eunsup = 23 // requested feature not (yet) implemented - eunimpl = 24 + eunimpl = 24 // cache limit exceeded - ecache = 25 + ecache = 25 // table limit exceeded - etable = 26 + etable = 26 // iteration is not making progress towards solution - enoprog = 27 + enoprog = 27 // jacobian evaluations are not improving the solution - enoprogj = 28 + enoprogj = 28 // cannot reach the specified tolerance in F - etolf = 29 + etolf = 29 // cannot reach the specified tolerance in X - etolx = 30 + etolx = 30 // cannot reach the specified tolerance in gradient - etolg = 31 + etolg = 31 // end of file - eof = 32 + eof = 32 } diff --git a/examples/README.md b/examples/README.md index deeb4f53f..0063b9e8a 100644 --- a/examples/README.md +++ b/examples/README.md @@ -1,6 +1,9 @@ # V Scientific Library - Examples Directory πŸ“š -This directory contains various examples demonstrating the capabilities and usage of the V Scientific Library. Each example showcases different functionalities, from machine learning algorithms to plotting and mathematical computations. Below is a summary of the examples included in this directory: +This directory contains various examples demonstrating the capabilities and usage +of the V Scientific Library. Each example showcases different functionalities, +from machine learning algorithms to plotting and mathematical computations. +Below is a summary of the examples included in this directory: ## Machine Learning Examples πŸ€– @@ -91,7 +94,11 @@ This directory contains various examples demonstrating the capabilities and usag ### Important Information ⚠️ -- Each example contains a `main.v` file with the V code demonstrating the specific functionality. -- Some examples include an additional `README.md` file. **You must read the `README.md` before running the example** to understand any prerequisites or specific instructions. +- Each example contains a `main.v` file with the V code demonstrating +the specific functionality. +- Some examples include an additional `README.md` file. +**You must read the `README.md` before running the example** to understand +any prerequisites or specific instructions. -To get started, navigate to the respective example directory and run the `main.v` file using the V compiler. +To get started, navigate to the respective example directory +and run the `main.v` file using the V compiler. diff --git a/examples/deriv_example/README.md b/examples/deriv_example/README.md index 3c28b51a0..1139a6e5a 100644 --- a/examples/deriv_example/README.md +++ b/examples/deriv_example/README.md @@ -1,6 +1,7 @@ # Example - deriv_example πŸ“˜ -This example demonstrates the usage of the V Scientific Library for demonstrating derivative calculation. +This example demonstrates the usage of the V Scientific Library +for demonstrating derivative calculation. ## Instructions diff --git a/examples/dist_histogram/README.md b/examples/dist_histogram/README.md index 60d2afc42..fd83f2e99 100644 --- a/examples/dist_histogram/README.md +++ b/examples/dist_histogram/README.md @@ -1,11 +1,13 @@ # Example - dist_histogram πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a distribution histogram. +This example demonstrates the usage of the V Scientific Library +for creating a distribution histogram. ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. You can do it +following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/fft_plot_example/README.md b/examples/fft_plot_example/README.md index a3a8e4842..6db29fc88 100644 --- a/examples/fft_plot_example/README.md +++ b/examples/fft_plot_example/README.md @@ -1,6 +1,7 @@ # Example - fft_plot_example πŸ“˜ -This example demonstrates the usage of the V Scientific Library for demonstrating Fast Fourier Transform (FFT) with plotting. +This example demonstrates the usage of the V Scientific Library +for demonstrating Fast Fourier Transform (FFT) with plotting. ## Instructions diff --git a/examples/io_h5_dataset/README.md b/examples/io_h5_dataset/README.md index e0be3c97f..bc1434206 100644 --- a/examples/io_h5_dataset/README.md +++ b/examples/io_h5_dataset/README.md @@ -1,11 +1,13 @@ # Example - io_h5_dataset πŸ“˜ -This example demonstrates the usage of the V Scientific Library for demonstrating HDF5 I/O for datasets. +This example demonstrates the usage of the V Scientific Library +for demonstrating HDF5 I/O for datasets. ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. You can do it +following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/io_h5_relax/README.md b/examples/io_h5_relax/README.md index fe9b5e8ae..b4618978c 100644 --- a/examples/io_h5_relax/README.md +++ b/examples/io_h5_relax/README.md @@ -1,6 +1,7 @@ # Example - io_h5_relax πŸ“˜ -This example demonstrates the usage of the V Scientific Library for demonstrating HDF5 I/O for relaxation data. +This example demonstrates the usage of the V Scientific Library +for demonstrating HDF5 I/O for relaxation data. ## Instructions diff --git a/examples/iter_lazy_generation/README.md b/examples/iter_lazy_generation/README.md index 3b0bc80ba..b8c1dda6c 100644 --- a/examples/iter_lazy_generation/README.md +++ b/examples/iter_lazy_generation/README.md @@ -1,6 +1,7 @@ # Example - iter_lazy_generation πŸ“˜ -This example demonstrates the usage of the V Scientific Library for demonstrating lazy generation using iterators. +This example demonstrates the usage of the V Scientific Library +for demonstrating lazy generation using iterators. ## Instructions diff --git a/examples/la_triplet01/README.md b/examples/la_triplet01/README.md index e987cd7a9..a45200f41 100644 --- a/examples/la_triplet01/README.md +++ b/examples/la_triplet01/README.md @@ -1,6 +1,7 @@ # Example - la_triplet01 πŸ“˜ -This example demonstrates the usage of the V Scientific Library for demonstrating linear algebra operations. +This example demonstrates the usage of the V Scientific Library +for demonstrating linear algebra operations. ## Instructions diff --git a/examples/ml_kmeans/README.md b/examples/ml_kmeans/README.md index 5f635d365..ff45c5d17 100644 --- a/examples/ml_kmeans/README.md +++ b/examples/ml_kmeans/README.md @@ -1,6 +1,7 @@ # Example - ml_kmeans πŸ“˜ -This example demonstrates the usage of the V Scientific Library for demonstrating the K-means clustering algorithm. +This example demonstrates the usage of the V Scientific Library +for demonstrating the K-means clustering algorithm. ## Instructions diff --git a/examples/ml_kmeans_plot/README.md b/examples/ml_kmeans_plot/README.md index 21855624f..d5cc605b9 100644 --- a/examples/ml_kmeans_plot/README.md +++ b/examples/ml_kmeans_plot/README.md @@ -1,6 +1,7 @@ # Example - ml_kmeans_plot πŸ“˜ -This example demonstrates the usage of the V Scientific Library for performing K-means clustering with plotting. +This example demonstrates the usage of the V Scientific Library +for performing K-means clustering with plotting. ## Instructions diff --git a/examples/ml_knn_plot/README.md b/examples/ml_knn_plot/README.md index 47472c8e6..4490f3dbb 100644 --- a/examples/ml_knn_plot/README.md +++ b/examples/ml_knn_plot/README.md @@ -1,6 +1,7 @@ # Example - ml_knn_plot πŸ“˜ -This example demonstrates the usage of the V Scientific Library for performing K-Nearest Neighbors algorithm with plotting. +This example demonstrates the usage of the V Scientific Library +for performing K-Nearest Neighbors algorithm with plotting. ## Instructions diff --git a/examples/ml_linreg01/README.md b/examples/ml_linreg01/README.md index 8e5955a4a..773b5b574 100644 --- a/examples/ml_linreg01/README.md +++ b/examples/ml_linreg01/README.md @@ -1,6 +1,7 @@ # Example - ml_linreg01 πŸ“˜ -This example demonstrates the usage of the V Scientific Library for performing a basic linear regression. +This example demonstrates the usage of the V Scientific Library +for performing a basic linear regression. ## Instructions diff --git a/examples/ml_linreg02/README.md b/examples/ml_linreg02/README.md index 29cc85fb2..7ccd49931 100644 --- a/examples/ml_linreg02/README.md +++ b/examples/ml_linreg02/README.md @@ -1,6 +1,7 @@ # Example - ml_linreg02 πŸ“˜ -This example demonstrates the usage of the V Scientific Library for performing an advanced linear regression. +This example demonstrates the usage of the V Scientific Library +for performing an advanced linear regression. ## Instructions diff --git a/examples/ml_linreg_plot/README.md b/examples/ml_linreg_plot/README.md index afee7f994..16c1b3da9 100644 --- a/examples/ml_linreg_plot/README.md +++ b/examples/ml_linreg_plot/README.md @@ -1,6 +1,7 @@ # Example - ml_linreg_plot πŸ“˜ -This example demonstrates the usage of the V Scientific Library for performing linear regression with plotting. +This example demonstrates the usage of the V Scientific Library +for performing linear regression with plotting. ## Instructions diff --git a/examples/ml_sentiment_analysis/README.md b/examples/ml_sentiment_analysis/README.md index aca5b2bf9..f73db9582 100644 --- a/examples/ml_sentiment_analysis/README.md +++ b/examples/ml_sentiment_analysis/README.md @@ -1,6 +1,7 @@ # Example - ml_sentiment_analysis πŸ“˜ -This example demonstrates the usage of the V Scientific Library for performing sentiment analysis using machine learning. +This example demonstrates the usage of the V Scientific Library +for performing sentiment analysis using machine learning. ## Instructions diff --git a/examples/mpi_basic_example/README.md b/examples/mpi_basic_example/README.md index 432504cd3..7a257bdbb 100644 --- a/examples/mpi_basic_example/README.md +++ b/examples/mpi_basic_example/README.md @@ -1,6 +1,7 @@ # Example - mpi_basic_example πŸ“˜ -This example demonstrates the usage of the V Scientific Library for demonstrating basic MPI functionality. +This example demonstrates the usage of the V Scientific Library +for demonstrating basic MPI functionality. ## Instructions diff --git a/examples/noise_fractal_2d/README.md b/examples/noise_fractal_2d/README.md new file mode 100644 index 000000000..ea94dd58f --- /dev/null +++ b/examples/noise_fractal_2d/README.md @@ -0,0 +1,16 @@ +# Example - noise_fractal_2d πŸ“˜ + +This example demonstrates the usage of the V Scientific Library for generating pink noise. + +## Instructions + +1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). +2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +3. Navigate to this directory. +4. Run the example using the following command: + +```sh +v run main.v +``` + +Enjoy exploring the capabilities of the V Scientific Library! 😊 diff --git a/examples/noise_fractal_2d/main.v b/examples/noise_fractal_2d/main.v new file mode 100644 index 000000000..a5ef30d25 --- /dev/null +++ b/examples/noise_fractal_2d/main.v @@ -0,0 +1,54 @@ +module main + +import rand +import vsl.noise +import vsl.plot + +fn main() { + // Creation of the noise generator. + // Other noise functions and dimensions count are available. + rand.seed([u32(3155200429), u32(3208395956)]) + mut generator := noise.Generator.new() + generator.randomize() + + mut x := []f64{} + mut y := []f64{} + mut z := []f64{} + + // 5 layers of simplex noise + octaves := 5 + persistence := 0.5 + + for xx in 0 .. 200 { + for yy in 0 .. 200 { + mut total := 0.0 + mut frequency := 1.0 + mut amplitude := 1.0 + mut max_value := 0.0 + + for _ in 0 .. octaves { + total += generator.simplex_2d(0.03 * xx * frequency, 0.03 * yy * frequency) * amplitude + max_value += amplitude + amplitude *= persistence + frequency *= 2.0 + } + + // Normalize to [-1, 1] + total /= max_value + x << xx + y << yy + z << total + } + } + + mut plt := plot.Plot.new() + plt.heatmap( + x: x + y: y + z: z + ) + plt.layout( + title: 'Pink `fractal` noise 2d example' + ) + plt.show()! +} diff --git a/examples/noise_simplex_2d/README.md b/examples/noise_simplex_2d/README.md new file mode 100644 index 000000000..f21dc9268 --- /dev/null +++ b/examples/noise_simplex_2d/README.md @@ -0,0 +1,16 @@ +# Example - noise_simple_2d πŸ“˜ + +This example demonstrates the usage of the V Scientific Library for generating simplex noise. + +## Instructions + +1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). +2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +3. Navigate to this directory. +4. Run the example using the following command: + +```sh +v run main.v +``` + +Enjoy exploring the capabilities of the V Scientific Library! 😊 diff --git a/examples/noise_simplex_2d/main.v b/examples/noise_simplex_2d/main.v new file mode 100644 index 000000000..031441cb3 --- /dev/null +++ b/examples/noise_simplex_2d/main.v @@ -0,0 +1,38 @@ +module main + +import rand +import vsl.noise +import vsl.plot + +fn main() { + // Creation of the noise generator. + // Other noise functions and dimensions count are available. + rand.seed([u32(3155200429), u32(3208395956)]) + mut generator := noise.Generator.new() + generator.randomize() + + mut x := []f64{} + mut y := []f64{} + mut z := []f64{} + + for xx in 0 .. 40 { + for yy in 0 .. 40 { + // We need to scale the coordinates to control the frequency of the noise. + val := generator.simplex_2d(0.03 * xx, 0.03 * yy) + x << xx + y << yy + z << val + } + } + + mut plt := plot.Plot.new() + plt.heatmap( + x: x + y: y + z: z + ) + plt.layout( + title: 'Simplex noise 2d example' + ) + plt.show()! +} diff --git a/examples/plot_heatmap_golden_ratio/README.md b/examples/plot_heatmap_golden_ratio/README.md index 39ec2c106..ecfe7129a 100644 --- a/examples/plot_heatmap_golden_ratio/README.md +++ b/examples/plot_heatmap_golden_ratio/README.md @@ -1,11 +1,13 @@ # Example - plot_heatmap_golden_ratio πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a heatmap with the golden ratio. +This example demonstrates the usage of the V Scientific Library +for creating a heatmap with the golden ratio. ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. You can do it +following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/plot_histogram/README.md b/examples/plot_histogram/README.md index 25e42a614..a519f8012 100644 --- a/examples/plot_histogram/README.md +++ b/examples/plot_histogram/README.md @@ -1,11 +1,13 @@ # Example - plot_histogram πŸ“˜ -This example demonstrates the usage of the V Scientific Library for showing how to create a histogram. +This example demonstrates the usage of the V Scientific Library +for showing how to create a histogram. ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. You can do it +following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/plot_line_axis_titles/README.md b/examples/plot_line_axis_titles/README.md index 8ec387e4b..bfcc17089 100644 --- a/examples/plot_line_axis_titles/README.md +++ b/examples/plot_line_axis_titles/README.md @@ -1,6 +1,7 @@ # Example - plot_line_axis_titles πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a line plot with axis titles. +This example demonstrates the usage of the V Scientific Library +for creating a line plot with axis titles. ## Instructions diff --git a/examples/plot_line_plot_with_areas/README.md b/examples/plot_line_plot_with_areas/README.md index 981174df5..71f8eaa66 100644 --- a/examples/plot_line_plot_with_areas/README.md +++ b/examples/plot_line_plot_with_areas/README.md @@ -1,11 +1,13 @@ # Example - plot_line_plot_with_areas πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a line plot with shaded areas. +This example demonstrates the usage of the V Scientific Library +for creating a line plot with shaded areas. ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. You can do it +following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/plot_scatter3d_1/README.md b/examples/plot_scatter3d_1/README.md index c1517cf6a..aed9f51d0 100644 --- a/examples/plot_scatter3d_1/README.md +++ b/examples/plot_scatter3d_1/README.md @@ -1,11 +1,13 @@ # Example - plot_scatter3d_1 πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a 3D scatter plot (example 1). +This example demonstrates the usage of the V Scientific Library for +creating a 3D scatter plot (example 1). ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. You can do it following +the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/plot_scatter3d_2/README.md b/examples/plot_scatter3d_2/README.md index 8b2c82c5e..d5853e956 100644 --- a/examples/plot_scatter3d_2/README.md +++ b/examples/plot_scatter3d_2/README.md @@ -1,6 +1,7 @@ # Example - plot_scatter3d_2 πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a 3D scatter plot (example 2). +This example demonstrates the usage of the V Scientific Library +for creating a 3D scatter plot (example 2). ## Instructions diff --git a/examples/plot_scatter3d_easing/README.md b/examples/plot_scatter3d_easing/README.md index 52e104d7f..7c8146dd2 100644 --- a/examples/plot_scatter3d_easing/README.md +++ b/examples/plot_scatter3d_easing/README.md @@ -1,6 +1,7 @@ # Example - plot_scatter3d_easing πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a 3D scatter plot with easing. +This example demonstrates the usage of the V Scientific Library +for creating a 3D scatter plot with easing. ## Instructions diff --git a/examples/plot_scatter_with_annotations/README.md b/examples/plot_scatter_with_annotations/README.md index bbe7e8164..9e750b95f 100644 --- a/examples/plot_scatter_with_annotations/README.md +++ b/examples/plot_scatter_with_annotations/README.md @@ -1,11 +1,13 @@ # Example - plot_scatter_with_annotations πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a scatter plot with annotations. +This example demonstrates the usage of the V Scientific Library for creating +a scatter plot with annotations. ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. +You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/plot_scatter_with_bars/README.md b/examples/plot_scatter_with_bars/README.md index 0e9f69a07..52dea8ee6 100644 --- a/examples/plot_scatter_with_bars/README.md +++ b/examples/plot_scatter_with_bars/README.md @@ -1,6 +1,7 @@ # Example - plot_scatter_with_bars πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a scatter plot with bars. +This example demonstrates the usage of the V Scientific Library +for creating a scatter plot with bars. ## Instructions diff --git a/examples/plot_scatter_with_histogram/README.md b/examples/plot_scatter_with_histogram/README.md index 648aa8d55..5082070c8 100644 --- a/examples/plot_scatter_with_histogram/README.md +++ b/examples/plot_scatter_with_histogram/README.md @@ -1,11 +1,13 @@ # Example - plot_scatter_with_histogram πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a scatter plot with histograms. +This example demonstrates the usage of the V Scientific Library +for creating a scatter plot with histograms. ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. You can do it +following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/plot_scatter_with_regression/README.md b/examples/plot_scatter_with_regression/README.md index da5588c28..8a830cf0a 100644 --- a/examples/plot_scatter_with_regression/README.md +++ b/examples/plot_scatter_with_regression/README.md @@ -1,6 +1,7 @@ # Example - plot_scatter_with_regression πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a scatter plot with regression line. +This example demonstrates the usage of the V Scientific Library +for creating a scatter plot with regression line. ## Instructions diff --git a/examples/plot_script_mode_ac_signal/README.md b/examples/plot_script_mode_ac_signal/README.md index d4873e78d..e4f5e61e3 100644 --- a/examples/plot_script_mode_ac_signal/README.md +++ b/examples/plot_script_mode_ac_signal/README.md @@ -1,6 +1,7 @@ # Example - plot_script_mode_ac_signal πŸ“˜ -This example demonstrates the usage of the V Scientific Library for plotting an AC signal in script mode. +This example demonstrates the usage of the V Scientific Library +for plotting an AC signal in script mode. ## Instructions diff --git a/examples/plot_script_mode_simple_plot/README.md b/examples/plot_script_mode_simple_plot/README.md index 9aaaea1d7..d56c092cc 100644 --- a/examples/plot_script_mode_simple_plot/README.md +++ b/examples/plot_script_mode_simple_plot/README.md @@ -1,6 +1,7 @@ # Example - plot_script_mode_simple_plot πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a simple plot in script mode. +This example demonstrates the usage of the V Scientific Library +for creating a simple plot in script mode. ## Instructions diff --git a/examples/plot_script_mode_three_phase_signal/README.md b/examples/plot_script_mode_three_phase_signal/README.md index fc6f8fad6..ddf2ea1c1 100644 --- a/examples/plot_script_mode_three_phase_signal/README.md +++ b/examples/plot_script_mode_three_phase_signal/README.md @@ -1,11 +1,13 @@ # Example - plot_script_mode_three_phase_signal πŸ“˜ -This example demonstrates the usage of the V Scientific Library for plotting a three-phase signal in script mode. +This example demonstrates the usage of the V Scientific Library +for plotting a three-phase signal in script mode. ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. You can do it +following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/plot_shaded_area_sin/README.md b/examples/plot_shaded_area_sin/README.md index 45ee03622..95ac176bc 100644 --- a/examples/plot_shaded_area_sin/README.md +++ b/examples/plot_shaded_area_sin/README.md @@ -1,11 +1,13 @@ # Example - plot_shaded_area_sin πŸ“˜ -This example demonstrates the usage of the V Scientific Library for creating a shaded area plot of the sine function. +This example demonstrates the usage of the V Scientific Library +for creating a shaded area plot of the sine function. ## Instructions 1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io). -2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! +2. Ensure you have the VSL installed. You can do it +following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)! 3. Navigate to this directory. 4. Run the example using the following command: diff --git a/examples/plot_sin_cos_surface/README.md b/examples/plot_sin_cos_surface/README.md index 9d6210252..77924ea8b 100644 --- a/examples/plot_sin_cos_surface/README.md +++ b/examples/plot_sin_cos_surface/README.md @@ -1,6 +1,7 @@ # Example - plot_sin_cos_surface πŸ“˜ -This example demonstrates the usage of the V Scientific Library for plotting the sine and cosine surface. +This example demonstrates the usage of the V Scientific Library +for plotting the sine and cosine surface. ## Instructions diff --git a/examples/plot_surface_easing/README.md b/examples/plot_surface_easing/README.md index a1820e673..f18e56d12 100644 --- a/examples/plot_surface_easing/README.md +++ b/examples/plot_surface_easing/README.md @@ -1,6 +1,7 @@ # Example - plot_surface_easing πŸ“˜ -This example demonstrates the usage of the V Scientific Library for plotting a surface with easing functions. +This example demonstrates the usage of the V Scientific Library +for plotting a surface with easing functions. ## Instructions diff --git a/examples/roots_bisection_solver/README.md b/examples/roots_bisection_solver/README.md index 7f8e212b5..3a0c7a709 100644 --- a/examples/roots_bisection_solver/README.md +++ b/examples/roots_bisection_solver/README.md @@ -1,6 +1,7 @@ # Example - roots_bisection_solver πŸ“˜ -This example demonstrates the usage of the V Scientific Library for finding roots using the bisection method. +This example demonstrates the usage of the V Scientific Library +for finding roots using the bisection method. ## Instructions diff --git a/examples/vcl_opencl_basic/README.md b/examples/vcl_opencl_basic/README.md index 74e1ed1d1..092cac13c 100644 --- a/examples/vcl_opencl_basic/README.md +++ b/examples/vcl_opencl_basic/README.md @@ -1,6 +1,7 @@ # Example - vcl_opencl_basic πŸ“˜ -This example demonstrates the usage of the V Scientific Library for demonstrating basic OpenCL functionality. +This example demonstrates the usage of the V Scientific Library +for demonstrating basic OpenCL functionality. ## Instructions diff --git a/examples/vcl_opencl_image_example/README.md b/examples/vcl_opencl_image_example/README.md index ac177d21b..4f3d8ad5b 100644 --- a/examples/vcl_opencl_image_example/README.md +++ b/examples/vcl_opencl_image_example/README.md @@ -8,4 +8,5 @@ This example shows how to do basic image processing with VCL using OpenCL as bac v run main.v ``` -after running it you should see the output image in the directory `output/` with the name `inverted.png`. +After running it, you should see the output image in the directory `output/` +with the name `inverted.png`. diff --git a/gm/bins.v b/gm/bins.v index c4f0a6026..5cffe7cf3 100644 --- a/gm/bins.v +++ b/gm/bins.v @@ -18,7 +18,7 @@ pub mut: @[heap] pub struct Bin { pub mut: - index int // index of bin + index int // index of bin entries []&BinEntry // entries } diff --git a/graph/graph.v b/graph/graph.v index 919785976..4d0ad390c 100644 --- a/graph/graph.v +++ b/graph/graph.v @@ -92,7 +92,8 @@ pub fn (g &Graph) get_edge(i int, j int) !int { */ pub fn (g &Graph) shortest_paths(method ShortestPaths) Graph { if method != .fw { - panic('shortest_paths works with FW (Floyd-Warshall) method only for now') + errors.vsl_panic('shortest_paths works with FW (Floyd-Warshall) method only for now', + .efailed) } g2 := g.calc_dist() nv := g2.dist.len @@ -176,7 +177,8 @@ pub fn (g &Graph) calc_dist() Graph { dist[i][j] = d next[i][j] = j if dist[i][j] < 0 { - panic('distance between vertices cannot be negative: ${dist}[i][j]') + errors.vsl_panic('distance between vertices cannot be negative: ${dist}[i][j]', + .efailed) } } return Graph{ diff --git a/iter/ranges.v b/iter/ranges.v index 23f28096c..27e48b6b6 100644 --- a/iter/ranges.v +++ b/iter/ranges.v @@ -129,8 +129,8 @@ pub: @[params] pub struct LinearIterParams { pub: - start f64 @[required] - stop f64 @[required] + start f64 @[required] + stop f64 @[required] len i64 = 50 endpoint bool = true } @@ -193,8 +193,8 @@ mut: @[params] pub struct LogIterParams { - start f64 @[required] - stop f64 @[required] + start f64 @[required] + stop f64 @[required] len i64 = 50 base f64 = 10.0 endpoint bool = true diff --git a/la/matrix_ops.v b/la/matrix_ops.v index 9e8cfa1f8..a541ed011 100644 --- a/la/matrix_ops.v +++ b/la/matrix_ops.v @@ -90,8 +90,7 @@ pub fn matrix_svd(mut s []f64, mut u Matrix[f64], mut vt Matrix[f64], mut a Matr if copy_a { acpy = a.clone() } - vlas.dgesvd(&char('A'.str), &char('A'.str), a.m, a.n, acpy.data, 1, s, u.data, a.m, - vt.data, a.n, superb) + vlas.dgesvd(`A`, `A`, a.m, a.n, acpy.data, 1, s, u.data, a.m, vt.data, a.n, superb) } // matrix_inv computes the inverse of a general matrix (square or not). It also computes the diff --git a/la/sparse_config.v b/la/sparse_config.v index 66ac022dc..9c2271d0a 100644 --- a/la/sparse_config.v +++ b/la/sparse_config.v @@ -7,10 +7,10 @@ import vsl.mpi pub struct SparseConfig { mut: communicator ?&mpi.Communicator // MPI communicator for parallel solvers [may be none] - mumps_ordering int // ICNTL(7) default = "" == "auto" - mumps_scaling int // Scaling type (check MUMPS solver) [may be empty] - symmetric bool // indicates symmetric system. NOTE: when using MUMPS, only the upper or lower part of the matrix must be provided - sym_pos_def bool // indicates symmetric-positive-defined system. NOTE: when using MUMPS, only the upper or lower part of the matrix must be provided + mumps_ordering int // ICNTL(7) default = "" == "auto" + mumps_scaling int // Scaling type (check MUMPS solver) [may be empty] + symmetric bool // indicates symmetric system. NOTE: when using MUMPS, only the upper or lower part of the matrix must be provided + sym_pos_def bool // indicates symmetric-positive-defined system. NOTE: when using MUMPS, only the upper or lower part of the matrix must be provided pub mut: verbose bool // run on verbose mode mumps_increase_of_working_space_pct int // MUMPS control parameters (check MUMPS solver manual) ICNTL(14) default = 100% diff --git a/ml/README.md b/ml/README.md new file mode 100644 index 000000000..e9c8d40a5 --- /dev/null +++ b/ml/README.md @@ -0,0 +1,80 @@ +# VSL Machine Learning (vsl.ml) + +VSL aims to provide a robust set of tools for scientific computing with an emphasis +on performance and ease of use. In the `vsl.ml` module, some machine learning +models are designed as observers of data, meaning they re-train automatically when +data changes, while others do not require this functionality. + +## Key Features + +- **Observers of Data**: Some machine learning models in VSL act as observers, + re-training automatically when data changes. +- **High Performance**: Leverages V’s performance optimizations and can integrate + with C and Fortran libraries like Open BLAS and LAPACK. +- **Versatile Algorithms**: Supports a variety of machine learning algorithms and + models. + +## Usage + +### Loading Data + +The `Data` struct in `vsl.ml` is designed to hold data in matrix format for machine +learning tasks. Here's a brief overview of how to use it: + +#### Creating a Data Object + +You can create a `Data` object using the following methods: + +- `Data.new`: Creates a new `Data` object with specified dimensions. +- `Data.from_raw_x`: Creates a `Data` object from raw x values (without y values). +- `Data.from_raw_xy`: Creates a `Data` object from raw x and y values combined in a single matrix. +- `Data.from_raw_xy_sep`: Creates a `Data` object from separate x and y raw values. + +### Data Methods + +The `Data` struct has several key methods to manage and manipulate data: + +- `set(x, y)`: Sets the x matrix and y vector and notifies observers. +- `set_y(y)`: Sets the y vector and notifies observers. +- `set_x(x)`: Sets the x matrix and notifies observers. +- `split(ratio)`: Splits the data into two parts based on the given ratio. +- `clone()`: Returns a deep copy of the Data object without observers. +- `clone_with_same_x()`: Returns a deep copy of the Data object but shares the same x reference. +- `add_observer(obs)`: Adds an observer to the data object. +- `notify_update()`: Notifies observers of data changes. + +### Stat Observer + +The `Stat` struct is an observer of `Data`, providing statistical analysis of the +data it observes. It automatically updates its statistics when the underlying data +changes. + +## Observer Models + +The following machine learning models in VSL are compatible with the `Observer` +pattern. This means they can observe data changes and automatically update +themselves. + +### K-Means Clustering + +K-Means Clustering is used for unsupervised learning to group data points into +clusters. As an observer model, it re-trains automatically when the data changes, +which is useful for dynamic datasets that require continuous updates. + +### K-Nearest Neighbors (KNN) + +K-Nearest Neighbors (KNN) is used for classification tasks where the target +variable is categorical. As an observer model, it re-trains automatically when the +data changes, which is beneficial for datasets that are frequently updated. + +## Non-Observer Models + +The following machine learning models in VSL do not require the observer pattern +and are trained once on a dataset without continuous updates. + +### Linear Regression + +Linear Regression is used for predicting a continuous target variable based on one +or more predictor variables. It is typically trained once on a dataset and used to +make predictions without requiring continuous updates. Hence, it is not implemented +as an observer model. diff --git a/ml/data.v b/ml/data.v index 273d3c47b..1a186bca2 100644 --- a/ml/data.v +++ b/ml/data.v @@ -26,7 +26,7 @@ pub mut: nb_samples int // number of data points (samples). number of rows in x and y nb_features int // number of features. number of columns in x x &la.Matrix[T] = unsafe { nil } // [nb_samples][nb_features] x values - y []T // [nb_samples] y values [optional] + y []T // [nb_samples] y values [optional] } // Data.new returns a new object to hold ML data diff --git a/ml/kmeans.v b/ml/kmeans.v index 9e869916a..87c39f583 100644 --- a/ml/kmeans.v +++ b/ml/kmeans.v @@ -9,12 +9,12 @@ import vsl.plot @[heap] pub struct Kmeans { mut: - name string // name of this "observer" + name string // name of this "observer" data &Data[f64] = unsafe { nil } // x data stat &Stat[f64] = unsafe { nil } // statistics about x (data) - nb_classes int // expected number of classes + nb_classes int // expected number of classes bins &gm.Bins = unsafe { nil } // "bins" to speed up searching for data points given their coordinates (2D or 3D only at the moment) - nb_iter int // number of iterations + nb_iter int // number of iterations pub mut: classes []int // [nb_samples] indices of classes of each sample centroids [][]f64 // [nb_classes][nb_features] coordinates of centroids diff --git a/ml/linreg.v b/ml/linreg.v index 65b0ae9e3..20b1cf98a 100644 --- a/ml/linreg.v +++ b/ml/linreg.v @@ -9,7 +9,7 @@ import vsl.util pub struct LinReg { mut: // main - name string // name of this "observer" + name string // name of this "observer" data &Data[f64] = unsafe { nil } // x-y data // workspace e []f64 // vector e = bβ‹…o + xβ‹…theta - y [nb_samples] diff --git a/ml/workspace.v b/ml/workspace.v index 52810d4bd..110d03640 100644 --- a/ml/workspace.v +++ b/ml/workspace.v @@ -12,19 +12,19 @@ import vsl.la pub struct Stat[T] { pub mut: data &Data[T] = unsafe { nil } // data - name string // name of this object - min_x []T // [n_features] min x values - max_x []T // [n_features] max x values - sum_x []T // [n_features] sum of x values - mean_x []T // [n_features] mean of x values - sig_x []T // [n_features] standard deviations of x - del_x []T // [n_features] difference: max(x) - min(x) - min_y T // min of y values - max_y T // max of y values - sum_y T // sum of y values - mean_y T // mean of y values - sig_y T // standard deviation of y - del_y T // difference: max(y) - min(y) + name string // name of this object + min_x []T // [n_features] min x values + max_x []T // [n_features] max x values + sum_x []T // [n_features] sum of x values + mean_x []T // [n_features] mean of x values + sig_x []T // [n_features] standard deviations of x + del_x []T // [n_features] difference: max(x) - min(x) + min_y T // min of y values + max_y T // max of y values + sum_y T // sum of y values + mean_y T // mean of y values + sig_y T // standard deviation of y + del_y T // difference: max(y) - min(y) } // stat returns a new Stat object diff --git a/noise/noise.v b/noise/noise.v new file mode 100644 index 000000000..f5c1aa994 --- /dev/null +++ b/noise/noise.v @@ -0,0 +1,20 @@ +module noise + +import rand + +// Generator is a struct holding the permutation table used in perlin and simplex noise +pub struct Generator { +mut: + perm []int = rand.shuffle_clone(permutations) or { panic(err) } +} + +// new is a function that return a new Generator struct +pub fn Generator.new() Generator { + return Generator{} +} + +// randomize is a function that shuffle the permutation set inside the Generator struct +// will not shuffle if rand.seed is not changed +pub fn (mut generator Generator) randomize() { + generator.perm = rand.shuffle_clone(permutations) or { panic(err) } +} diff --git a/noise/perlin2d.v b/noise/perlin2d.v index 4e5d2661e..0b60ffda0 100644 --- a/noise/perlin2d.v +++ b/noise/perlin2d.v @@ -1,26 +1,7 @@ module noise -import rand - -// Perlin is a struct that hold the permutation set for perlin noise -pub struct Perlin { -mut: - perm []int = rand.shuffle_clone(permutations) or { panic(err) } -} - -// new is a function that return a new Perlin struct -pub fn Perlin.new() Perlin { - return Perlin{} -} - -// randomize is a function that shuffle the permutation set inside the Perlin struct -// will not shuffle if rand.seed is not changed -pub fn (mut perlin Perlin) randomize() { - perlin.perm = rand.shuffle_clone(permutations) or { panic(err) } -} - // perlin2d is a function that return a single value of perlin noise for a given 2d position -pub fn (perlin Perlin) perlin2d(x f64, y f64) f64 { +pub fn (generator Generator) perlin2d(x f64, y f64) f64 { xi := int(x) & 0xFF yi := int(y) & 0xFF @@ -30,13 +11,13 @@ pub fn (perlin Perlin) perlin2d(x f64, y f64) f64 { u := fade(xf) v := fade(yf) - pxi := perlin.perm[xi] - pxi1 := perlin.perm[xi + 1] + pxi := generator.perm[xi] + pxi1 := generator.perm[xi + 1] - aa := perlin.perm[pxi + yi] - ab := perlin.perm[pxi + yi + 1] - ba := perlin.perm[pxi1 + yi] - bb := perlin.perm[pxi1 + yi + 1] + aa := generator.perm[pxi + yi] + ab := generator.perm[pxi + yi + 1] + ba := generator.perm[pxi1 + yi] + bb := generator.perm[pxi1 + yi + 1] x1 := lerp(grad2d(aa, xf, yf), grad2d(ba, xf - 1, yf), u) x2 := lerp(grad2d(ab, xf, yf - 1), grad2d(bb, xf - 1, yf - 1), u) diff --git a/noise/perlin2d_test.v b/noise/perlin2d_test.v index 6b4c81d83..f63a5d3e2 100644 --- a/noise/perlin2d_test.v +++ b/noise/perlin2d_test.v @@ -6,7 +6,7 @@ import vsl.float.float64 fn test_perlin2d() { rand.seed([u32(3155200429), u32(3208395956)]) - mut gen := Perlin.new() + mut gen := Generator.new() gen.randomize() result := gen.perlin2d(0.125, 0.125) diff --git a/noise/perlin3d.v b/noise/perlin3d.v index c32a6131e..9b77e820b 100644 --- a/noise/perlin3d.v +++ b/noise/perlin3d.v @@ -1,7 +1,7 @@ module noise // perlin3d is a function that return a single value of perlin noise for a given 3d position -pub fn (perlin Perlin) perlin3d(x f64, y f64, z f64) f64 { +pub fn (generator Generator) perlin3d(x f64, y f64, z f64) f64 { xi := int(x) & 0xFF yi := int(y) & 0xFF zi := int(z) & 0xFF @@ -12,21 +12,21 @@ pub fn (perlin Perlin) perlin3d(x f64, y f64, z f64) f64 { v := fade(yf) w := fade(zf) - pxi := perlin.perm[xi] - pxi_yi := perlin.perm[pxi + yi] - pxi_yi1 := perlin.perm[pxi + yi + 1] - pxi1 := perlin.perm[xi + 1] - pxi1_yi := perlin.perm[pxi1 + yi] - pxi1_yi1 := perlin.perm[pxi1 + yi + 1] + pxi := generator.perm[xi] + pxi_yi := generator.perm[pxi + yi] + pxi_yi1 := generator.perm[pxi + yi + 1] + pxi1 := generator.perm[xi + 1] + pxi1_yi := generator.perm[pxi1 + yi] + pxi1_yi1 := generator.perm[pxi1 + yi + 1] - aaa := perlin.perm[pxi_yi + zi] - aba := perlin.perm[pxi_yi1 + zi] - aab := perlin.perm[pxi_yi + zi + 1] - abb := perlin.perm[pxi_yi1 + zi + 1] - baa := perlin.perm[pxi1_yi + zi] - bba := perlin.perm[pxi1_yi1 + zi] - bab := perlin.perm[pxi1_yi + zi + 1] - bbb := perlin.perm[pxi1_yi1 + zi + 1] + aaa := generator.perm[pxi_yi + zi] + aba := generator.perm[pxi_yi1 + zi] + aab := generator.perm[pxi_yi + zi + 1] + abb := generator.perm[pxi_yi1 + zi + 1] + baa := generator.perm[pxi1_yi + zi] + bba := generator.perm[pxi1_yi1 + zi] + bab := generator.perm[pxi1_yi + zi + 1] + bbb := generator.perm[pxi1_yi1 + zi + 1] mut x1 := lerp(grad3d(aaa, xf, yf, zf), grad3d(baa, xf - 1, yf, zf), u) mut x2 := lerp(grad3d(aba, xf, yf - 1, zf), grad3d(bba, xf - 1, yf - 1, zf), u) diff --git a/noise/perlin3d_test.v b/noise/perlin3d_test.v index 8edde204c..0f5a67aed 100644 --- a/noise/perlin3d_test.v +++ b/noise/perlin3d_test.v @@ -6,7 +6,7 @@ import vsl.float.float64 fn test_perlin3d() { rand.seed([u32(3155200429), u32(3208395956)]) - mut gen := Perlin.new() + mut gen := Generator.new() gen.randomize() result := gen.perlin3d(0.125, 0.125, 0.125) diff --git a/noise/simplex.v b/noise/simplex.v new file mode 100644 index 000000000..176fe40a6 --- /dev/null +++ b/noise/simplex.v @@ -0,0 +1,389 @@ +module noise + +import math + +// skewing factors for 2D case +const f2 = 0.5 * (math.sqrt(3.0) - 1.0) +const g2 = (3.0 - math.sqrt(3)) / 6.0 +// skewing factors for 3D case +const f3 = f64(1.0 / 3.0) +const g3 = f64(1.0 / 6.0) +// skewing factors for 4D case +const f4 = f64((math.sqrt(5.0) - 1.0) / 4) +const g4 = f64((5.0 - math.sqrt(5.0)) / 20.0) + +// vfmt off +const simplex := [ + [0, 1, 2, 3], [0, 1, 3, 2], [0, 0, 0, 0], [0, 2, 3, 1], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [1, 2, 3, 0], + [0, 2, 1, 3], [0, 0, 0, 0], [0, 3, 1, 2], [0, 3, 2, 1], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [1, 3, 2, 0], + [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], + [1, 2, 0, 3], [0, 0, 0, 0], [1, 3, 0, 2], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [2, 3, 0, 1], [2, 3, 1, 0], + [1, 0, 2, 3], [1, 0, 3, 2], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [2, 0, 3, 1], [0, 0, 0, 0], [2, 1, 3, 0], + [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], + [2, 0, 1, 3], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [3, 0, 1, 2], [3, 0, 2, 1], [0, 0, 0, 0], [3, 1, 2, 0], + [2, 1, 0, 3], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [3, 1, 0, 2], [0, 0, 0, 0], [3, 2, 0, 1], [3, 2, 1, 0] +] +// vfmt on + +// grad_1d returns a gradient value for a given hash and x position +fn grad_1d(hash int, x f64) f64 { + h := hash & 15 + mut grad := 1.0 + f64(h & 7) + grad = if h & 8 == 0 { grad } else { -grad } + return grad * x +} + +// grad_2d returns a gradient value for a given hash and x, y position +fn grad_2d(hash int, x f64, y f64) f64 { + h := hash & 7 + u := if h < 4 { x } else { y } + v := if h < 4 { y } else { x } + return (if h & 1 == 0 { + u + } else { + -u + }) + (if h & 2 == 0 { + 2.0 * v + } else { + -2.0 * v + }) +} + +// grad_3d returns a gradient value for a given hash and x, y, z position +fn grad_3d(hash int, x f64, y f64, z f64) f64 { + h := hash & 15 + u := if h < 8 { x } else { y } + v := if h < 4 { + y + } else { + if h == 12 || h == 14 { x } else { z } + } + return (if h & 1 == 0 { + u + } else { + -u + }) + (if h & 2 == 0 { + v + } else { + -v + }) +} + +// grad_4d returns a gradient value for a given hash and x, y, z, t position +fn grad_4d(hash int, x f64, y f64, z f64, t f64) f64 { + h := hash & 31 + u := if h < 24 { x } else { y } + v := if h < 16 { y } else { z } + w := if h < 8 { z } else { t } + return (if h & 1 == 0 { + u + } else { + -u + }) + (if h & 2 == 0 { + v + } else { + -v + }) + (if h & 4 == 0 { + w + } else { + -w + }) +} + +// simplex_1d returns a simplex noise value for a given x position +pub fn (generator Generator) simplex_1d(x f64) f64 { + i0 := int(x) + i1 := i0 + 1 + x0 := x - i0 + x1 := x0 - 1 + + mut t0 := 1.0 - x0 * x0 + t0 *= t0 + n0 := t0 * t0 * grad_1d(generator.perm[i0 & 0xff], x0) + + mut t1 := 1.0 - x1 * x1 + t1 *= t1 + n1 := t1 * t1 * grad_1d(generator.perm[i1 & 0xff], x1) + + return 0.395 * (n0 + n1) +} + +// simplex_2d returns a simplex noise value for a given x, y position +pub fn (generator Generator) simplex_2d(x f64, y f64) f64 { + s := (x + y) * noise.f2 + i := int(x + s) + j := int(y + s) + + t := f64(i + j) * noise.g2 + x0 := x - (i - t) + y0 := y - (j - t) + + i1, j1 := if x0 > y0 { + 1, 0 + } else { + 0, 1 + } + + x1 := x0 - f64(i1) + noise.g2 + y1 := y0 - f64(j1) + noise.g2 + x2 := x0 - 1.0 + noise.g2 * 2.0 + y2 := y0 - 1.0 + noise.g2 * 2.0 + + ii := i & 0xff + jj := j & 0xff + + mut t0 := 0.5 - x0 * x0 - y0 * y0 + mut n0 := 0.0 + if t0 < 0 { + n0 = 0.0 + } else { + t0 *= t0 + n0 = t0 * t0 * grad_2d(generator.perm[ii + generator.perm[jj]], x0, y0) + } + + mut t1 := 0.5 - x1 * x1 - y1 * y1 + mut n1 := 0.0 + if t1 < 0 { + n1 = 0.0 + } else { + t1 *= t1 + n1 = t1 * t1 * grad_2d(generator.perm[ii + i1 + generator.perm[jj + j1]], x1, + y1) + } + + mut t2 := 0.5 - x2 * x2 - y2 * y2 + mut n2 := 0.0 + if t2 < 0 { + n2 = 0.0 + } else { + t2 *= t2 + n2 = t2 * t2 * grad_2d(generator.perm[ii + 1 + generator.perm[jj + 1]], x2, y2) + } + + return 40.0 * (n0 + n1 + n2) +} + +// simplex_3d returns a simplex noise value for a given x, y, z position +pub fn (generator Generator) simplex_3d(x f64, y f64, z f64) f64 { + s := (x + y + z) * noise.f3 + xs := x + s + ys := y + s + zs := z + s + i := int(xs) + j := int(ys) + k := int(zs) + + t := f64(i + j + k) * noise.g3 + x0 := x - (i - t) + y0 := y - (j - t) + z0 := z - (k - t) + + mut i1 := 0 + mut j1 := 0 + mut k1 := 0 + mut i2 := 0 + mut j2 := 0 + mut k2 := 0 + + if x0 >= y0 { + if y0 >= z0 { + i1 = 1 + i2 = 1 + j2 = 1 + } else if x0 >= z0 { + i1 = 1 + i2 = 1 + k2 = 1 + } else { + k1 = 1 + i2 = 1 + k2 = 1 + } + } else { + if y0 < z0 { + k1 = 1 + j2 = 1 + k2 = 1 + } else if x0 < z0 { + j1 = 1 + j2 = 1 + k2 = 1 + } else { + j1 = 1 + i2 = 1 + j2 = 1 + } + } + + x1 := x0 - i1 + noise.g3 + y1 := y0 - j1 + noise.g3 + z1 := z0 - k1 + noise.g3 + x2 := x0 - i2 + 2.0 * noise.g3 + y2 := y0 - j2 + 2.0 * noise.g3 + z2 := z0 - k2 + 2.0 * noise.g3 + x3 := x0 - 1.0 + 3.0 * noise.g3 + y3 := y0 - 1.0 + 3.0 * noise.g3 + z3 := z0 - 1.0 + 3.0 * noise.g3 + + ii := i & 0xff + jj := j & 0xff + kk := k & 0xff + + mut t0 := 0.6 - x0 * x0 - y0 * y0 - z0 * z0 + mut n0 := 0.0 + if t0 < 0 { + n0 = 0.0 + } else { + t0 *= t0 + n0 = t0 * t0 * grad_3d(generator.perm[ii + generator.perm[jj + generator.perm[kk]]], + x0, y0, z0) + } + + mut t1 := 0.6 - x1 * x1 - y1 * y1 - z1 * z1 + mut n1 := 0.0 + if t1 < 0 { + n1 = 0.0 + } else { + t1 *= t1 + n1 = t1 * t1 * grad_3d(generator.perm[ii + i1 + generator.perm[jj + j1 + generator.perm[kk + + k1]]], x1, y1, z1) + } + + mut t2 := 0.6 - x2 * x2 - y2 * y2 - z2 * z2 + mut n2 := 0.0 + if t2 < 0 { + n2 = 0.0 + } else { + t2 *= t2 + n2 = t2 * t2 * grad_3d(generator.perm[ii + i2 + generator.perm[jj + j2 + generator.perm[kk + + k2]]], x2, y2, z2) + } + + mut t3 := 0.6 - x3 * x3 - y3 * y3 - z3 * z3 + mut n3 := 0.0 + if t3 < 0 { + n3 = 0.0 + } else { + t3 *= t3 + n3 = t3 * t3 * grad_3d(generator.perm[ii + 1 + generator.perm[jj + 1 + generator.perm[kk + + 1]]], x3, y3, z3) + } + + return 32.0 * (n0 + n1 + n2 + n3) +} + +// simplex_4d returns a simplex noise value for a given x, y, z, w position +pub fn (generator Generator) simplex_4d(x f64, y f64, z f64, w f64) f64 { + s := (x + y + z + w) * noise.f4 + xs := x + s + ys := y + s + zs := z + s + ws := w + s + i := int(xs) + j := int(ys) + k := int(zs) + l := int(ws) + + t := f64(i + j + k + l) * noise.g4 + x0 := x - (i - t) + y0 := y - (j - t) + z0 := z - (k - t) + w0 := w - (l - t) + + c1 := if x0 > y0 { 32 } else { 0 } + c2 := if x0 > z0 { 16 } else { 0 } + c3 := if y0 > z0 { 8 } else { 0 } + c4 := if x0 > w0 { 4 } else { 0 } + c5 := if y0 > w0 { 2 } else { 0 } + c6 := if z0 > w0 { 1 } else { 0 } + c := c1 + c2 + c3 + c4 + c5 + c6 + + i1 := if noise.simplex[c][0] >= 3 { 1 } else { 0 } + j1 := if noise.simplex[c][1] >= 3 { 1 } else { 0 } + k1 := if noise.simplex[c][2] >= 3 { 1 } else { 0 } + l1 := if noise.simplex[c][3] >= 3 { 1 } else { 0 } + + i2 := if noise.simplex[c][0] >= 2 { 1 } else { 0 } + j2 := if noise.simplex[c][1] >= 2 { 1 } else { 0 } + k2 := if noise.simplex[c][2] >= 2 { 1 } else { 0 } + l2 := if noise.simplex[c][3] >= 2 { 1 } else { 0 } + + i3 := if noise.simplex[c][0] >= 1 { 1 } else { 0 } + j3 := if noise.simplex[c][1] >= 1 { 1 } else { 0 } + k3 := if noise.simplex[c][2] >= 1 { 1 } else { 0 } + l3 := if noise.simplex[c][3] >= 1 { 1 } else { 0 } + + x1 := x0 - i1 + noise.g4 + y1 := y0 - j1 + noise.g4 + z1 := z0 - k1 + noise.g4 + w1 := w0 - l1 + noise.g4 + x2 := x0 - i2 + 2.0 * noise.g4 + y2 := y0 - j2 + 2.0 * noise.g4 + z2 := z0 - k2 + 2.0 * noise.g4 + w2 := w0 - l2 + 2.0 * noise.g4 + x3 := x0 - i3 + 3.0 * noise.g4 + y3 := y0 - j3 + 3.0 * noise.g4 + z3 := z0 - k3 + 3.0 * noise.g4 + w3 := w0 - l3 + 3.0 * noise.g4 + x4 := x0 - 1.0 + 4.0 * noise.g4 + y4 := y0 - 1.0 + 4.0 * noise.g4 + z4 := z0 - 1.0 + 4.0 * noise.g4 + w4 := w0 - 1.0 + 4.0 * noise.g4 + + ii := i & 0xff + jj := j & 0xff + kk := k & 0xff + ll := l & 0xff + + mut t0 := 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0 + mut n0 := 0.0 + if t0 < 0.0 { + n0 = 0.0 + } else { + t0 *= t0 + n0 = t0 * t0 * grad_4d(generator.perm[ii + generator.perm[jj + generator.perm[kk + + generator.perm[ll]]]], x0, y0, z0, w0) + } + + mut t1 := 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1 + mut n1 := 0.0 + if t1 < 0.0 { + n1 = 0.0 + } else { + t1 *= t1 + n1 = t1 * t1 * grad_4d(generator.perm[ii + i1 + generator.perm[jj + j1 + generator.perm[kk + + k1 + generator.perm[ll + l1]]]], x1, y1, z1, w1) + } + + mut t2 := 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2 + mut n2 := 0.0 + if t2 < 0.0 { + n2 = 0.0 + } else { + t2 *= t2 + n2 = t2 * t2 * grad_4d(generator.perm[ii + i2 + generator.perm[jj + j2 + generator.perm[kk + + k2 + generator.perm[ll + l2]]]], x2, y2, z2, w2) + } + + mut t3 := 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3 + mut n3 := 0.0 + if t3 < 0.0 { + n3 = 0.0 + } else { + t3 *= t3 + n3 = t3 * t3 * grad_4d(generator.perm[ii + i3 + generator.perm[jj + j3 + generator.perm[kk + + k3 + generator.perm[ll + l3]]]], x3, y3, z3, w3) + } + + mut t4 := 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4 + mut n4 := 0.0 + if t4 < 0.0 { + n4 = 0.0 + } else { + t4 *= t4 + n4 = t4 * t4 * grad_4d(generator.perm[ii + 1 + generator.perm[jj + 1 + generator.perm[kk + + 1 + generator.perm[ll + 1]]]], x4, y4, z4, w4) + } + + return 27.0 * (n0 + n1 + n2 + n3 + n4) +} diff --git a/noise/simplex_test.v b/noise/simplex_test.v new file mode 100644 index 000000000..66dd946a7 --- /dev/null +++ b/noise/simplex_test.v @@ -0,0 +1,40 @@ +module noise + +import rand +import vsl.float.float64 + +fn test_simplex_1d() { + rand.seed([u32(3155200429), u32(3208395956)]) + mut generator := Generator.new() + generator.randomize() + result := generator.simplex_1d(0.287) + expected := -0.3544283326507284 + assert float64.tolerance(result, expected, 1.0e-6) +} + +fn test_simplex_2d() { + rand.seed([u32(3075200429), u32(3094395956)]) + mut generator := Generator.new() + generator.randomize() + result := generator.simplex_2d(0.287, 0.475) + expected := -0.09948661872545192 + assert float64.tolerance(result, expected, 1.0e-6) +} + +fn test_simplex_3d() { + rand.seed([u32(3155200429), u32(3208395956)]) + mut generator := Generator.new() + generator.randomize() + result := generator.simplex_3d(0.287, 0.475, 1.917) + expected := -0.06034653476116279 + assert float64.tolerance(result, expected, 1.0e-6) +} + +fn test_simplex_4d() { + rand.seed([u32(3075200429), u32(3094395956)]) + mut generator := Generator.new() + generator.randomize() + result := generator.simplex_4d(0.287, 0.475, 1.917, 0.684) + expected := 0.015098415881100141 + assert float64.tolerance(result, expected, 1.0e-6) +} diff --git a/plot/axis.v b/plot/axis.v index 731d59697..29488a525 100644 --- a/plot/axis.v +++ b/plot/axis.v @@ -5,11 +5,11 @@ pub struct Axis { pub mut: title AxisTitle tickmode string = 'auto' - tick0 f64 @[omitempty] - dtick f64 @[omitempty] - tickvals []f64 @[omitempty] - ticktext []string @[omitempty] - range []f64 @[omitempty] + tick0 f64 @[omitempty] + dtick f64 @[omitempty] + tickvals []f64 @[omitempty] + ticktext []string @[omitempty] + range []f64 @[omitempty] } // AxisTitle handles needed data to render an axis title diff --git a/plot/trace.v b/plot/trace.v index 705ced8df..f1bb44516 100644 --- a/plot/trace.v +++ b/plot/trace.v @@ -26,20 +26,20 @@ pub struct CommonTrace { pub mut: x XType xbins map[string]f32 - y YType @[omitempty] - z ZType @[omitempty] - name string @[omitempty] - mode string @[omitempty] - marker Marker @[omitempty] - line Line @[omitempty] - pull []f64 @[omitempty] - hole f64 @[omitempty] - fill string @[omitempty] - fillcolor string @[omitempty] - customdata [][]string @[omitempty] - colorscale string = 'Viridis' @[omitempty] - textinfo string @[omitempty] - text []string @[omitempty] + y YType @[omitempty] + z ZType @[omitempty] + name string @[omitempty] + mode string @[omitempty] + marker Marker @[omitempty] + line Line @[omitempty] + pull []f64 @[omitempty] + hole f64 @[omitempty] + fill string @[omitempty] + fillcolor string @[omitempty] + customdata [][]string @[omitempty] + colorscale string = 'Viridis' @[omitempty] + textinfo string @[omitempty] + text []string @[omitempty] } // ScatterTrace is a struct for Scatter trace type diff --git a/poly/README.md b/poly/README.md index 2fe043dcb..0e6855662 100644 --- a/poly/README.md +++ b/poly/README.md @@ -1,13 +1,16 @@ # Polynomials -This chapter describes functions for evaluating and solving polynomials. -There are routines for finding real and complex roots of quadratic and -cubic equations using analytic methods. An iterative polynomial solver -is also available for finding the roots of general polynomials with real -coefficients (of any order). The functions are declared in the module `vsl.poly`. +This chapter describes functions for evaluating and solving polynomials. There are routines +for finding real and complex roots of quadratic and cubic equations using analytic methods. +An iterative polynomial solver is also available for finding the roots of general polynomials +with real coefficients (of any order). The functions are declared in the module `vsl.poly`. ## Polynomial Evaluation +```v ignore +fn eval(c []f64, x f64) f64 +``` + The functions described here evaluate the polynomial ```console @@ -16,20 +19,13 @@ P(x) = c[0] + c[1] x + c[2] x^2 + . . . + c[len-1] x^(len-1) using Horner's method for stability. -```v ignore -fn eval(c []f64, x f64) f64 -``` - -This function evaluates a polynomial with real coefficients for the real variable `x`. - ```v ignore fn eval_derivs(c []f64, x f64, lenres u64) []f64 ``` -This function evaluates a polynomial and its derivatives storing the -results in the array `res` of size `lenres`. The output array -contains the values of `d^k P(x)/d x^k` for the specified value of -`x` starting with `k = 0`. +This function evaluates a polynomial and its derivatives, storing the results in the array +`res` of size `lenres`. The output array contains the values of `d^k P(x)/d x^k` for the +specified value of `x`, starting with `k = 0`. ## Quadratic Equations @@ -43,22 +39,17 @@ This function finds the real roots of the quadratic equation, a x^2 + b x + c = 0 ``` -The number of real roots (either zero, one or two) is returned, and -their locations are are returned as `[ x0, x1 ]`. If no real roots -are found then `[]` is returned. If one real root -is found (i.e. if `a=0`) then it is are returned as `[ x0 ]`. When two -real roots are found they are are returned as `[ x0, x1 ]` in -ascending order. The case of coincident roots is not considered -special. For example `(x-1)^2=0` will have two roots, which happen -to have exactly equal values. +The number of real roots (either zero, one or two) is returned, and their locations are +returned as `[ x0, x1 ]`. If no real roots are found then `[]` is returned. If one real root +is found (i.e. if `a=0`) then it is returned as `[ x0 ]`. When two real roots are found they +are returned as `[ x0, x1 ]` in ascending order. The case of coincident roots is not considered +special. For example `(x-1)^2=0` will have two roots, which happen to have exactly equal values. -The number of roots found depends on the sign of the discriminant -`b^2 - 4 a c`. This will be subject to rounding and cancellation -errors when computed in double precision, and will also be subject to -errors if the coefficients of the polynomial are inexact. These errors -may cause a discrete change in the number of roots. However, for -polynomials with small integer coefficients the discriminant can always -be computed exactly. +The number of roots found depends on the sign of the discriminant `b^2 - 4 a c`. This will +be subject to rounding and cancellation errors when computed in double precision, and will +also be subject to errors if the coefficients of the polynomial are inexact. These errors may +cause a discrete change in the number of roots. However, for polynomials with small integer +coefficients the discriminant can always be computed exactly. ## Cubic Equations @@ -72,13 +63,93 @@ This function finds the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0 ``` -with a leading coefficient of unity. The number of real roots (either -one or three) is returned, and their locations are returned as `[ x0, x1, x2 ]`. -If one real root is found then only `[ x0 ]` -is returned. When three real roots are found they are returned as -`[ x0, x1, x2 ]` in ascending order. The case of -coincident roots is not considered special. For example, the equation -`(x-1)^3=0` will have three roots with exactly equal values. As -in the quadratic case, finite precision may cause equal or -closely-spaced real roots to move off the real axis into the complex -plane, leading to a discrete change in the number of real roots. +with a leading coefficient of unity. The number of real roots (either one or three) is +returned, and their locations are returned as `[ x0, x1, x2 ]`. If one real root is found +then only `[ x0 ]` is returned. When three real roots are found they are returned as +`[ x0, x1, x2 ]` in ascending order. The case of coincident roots is not considered special. +For example, the equation `(x-1)^3=0` will have three roots with exactly equal values. As +in the quadratic case, finite precision may cause equal or closely-spaced real roots to move +off the real axis into the complex plane, leading to a discrete change in the number of real roots. + +## Companion Matrix + +```v ignore +fn companion_matrix(a []f64) [][]f64 +``` + +Creates a companion matrix for the polynomial + +```console +P(x) = a_n * x^n + a_{n-1} * x^{n-1} + ... + a_1 * x + a_0 +``` + +The companion matrix `C` is defined as: + +``` +[0 0 0 ... 0 -a_0/a_n] +[1 0 0 ... 0 -a_1/a_n] +[0 1 0 ... 0 -a_2/a_n] +[. . . ... . ........] +[0 0 0 ... 1 -a_{n-1}/a_n] +``` + +## Balanced Companion Matrix + +```v ignore +fn balance_companion_matrix(cm [][]f64) [][]f64 +``` + +Balances a companion matrix `C` to improve numerical stability. It uses an iterative scaling +process to make the row and column norms as close to each other as possible. The output is +a balanced matrix `B` such that `D^(-1)CD = B`, where `D` is a diagonal matrix. + +## Polynomial Operations + +```v ignore +fn add(a []f64, b []f64) []f64 +``` + +Adds two polynomials: + +```console +(a_n * x^n + ... + a_0) + (b_m * x^m + ... + b_0) +``` + +Returns the result as `[a_0 + b_0, a_1 + b_1, ..., a_k + b_k ...]`. + +```v ignore +fn subtract(a []f64, b []f64) []f64 +``` + +Subtracts two polynomials: + +```console +(a_n * x^n + ... + a_0) - (b_m * x^m + ... + b_0) +``` + +Returns the result as `[a_0 - b_0, a_1 - b_1, ..., a_k - b_k, ...]`. + +```v ignore +fn multiply(a []f64, b []f64) []f64 +``` + +Multiplies two polynomials: + +```console +(a_n * x^n + ... + a_0) * (b_m * x^m + ... + b_0) +``` + +Returns the result as `[c_0, c_1, ..., c_{n+m}]` where `c_k = βˆ‘_{i+j=k} a_i * b_j`. + +```v ignore +fn divide(a []f64, b []f64) ([]f64, []f64) +``` + +Divides two polynomials: + +```console +(a_n * x^n + ... + a_0) / (b_m * x^m + ... + b_0) +``` + +Uses polynomial long division algorithm. Returns `(q, r)` where `q` is the quotient and `r` +is the remainder such that `a(x) = b(x) * q(x) + r(x)` and `degree(r) < degree(b)`. diff --git a/poly/poly.v b/poly/poly.v index 31fcce2e3..697c8de86 100644 --- a/poly/poly.v +++ b/poly/poly.v @@ -1,13 +1,18 @@ module poly import math +import vsl.errors const radix = 2 const radix2 = (radix * radix) +// eval is a function that evaluates a polynomial P(x) = a_n * x^n + a_{n-1} * x^{n-1} + ... + a_1 * x + a_0 +// using Horner's method: P(x) = (...((a_n * x + a_{n-1}) * x + a_{n-2}) * x + ... + a_1) * x + a_0 +// Input: c = [a_0, a_1, ..., a_n], x +// Output: P(x) pub fn eval(c []f64, x f64) f64 { if c.len == 0 { - panic('coeficients can not be empty') + errors.vsl_panic('coeficients can not be empty', .efailed) } len := c.len mut ans := 0.0 @@ -19,6 +24,9 @@ pub fn eval(c []f64, x f64) f64 { return ans } +// eval_derivs evaluates a polynomial P(x) and its derivatives P'(x), P''(x), ..., P^(k)(x) +// Input: c = [a_0, a_1, ..., a_n] representing P(x), x, and lenres (k+1) +// Output: [P(x), P'(x), P''(x), ..., P^(k)(x)] pub fn eval_derivs(c []f64, x f64, lenres int) []f64 { mut res := []f64{} lenc := c.len @@ -50,7 +58,11 @@ pub fn eval_derivs(c []f64, x f64, lenres int) []f64 { return res } -pub fn solve_quadratic(a f64, b f64, c f64) []f64 { // Handle linear case +// solve_quadratic solves the quadratic equation ax^2 + bx + c = 0 +// using the quadratic formula: x = (-b Β± √(b^2 - 4ac)) / (2a) +// Input: a, b, c +// Output: Array of real roots (if any) +pub fn solve_quadratic(a f64, b f64, c f64) []f64 { if a == 0 { if b == 0 { return [] @@ -77,6 +89,10 @@ pub fn solve_quadratic(a f64, b f64, c f64) []f64 { // Handle linear case } } +// solve_cubic solves the cubic equation x^3 + ax^2 + bx + c = 0 +// using Cardano's formula and trigonometric solution +// Input: a, b, c +// Output: Array of real roots pub fn solve_cubic(a f64, b f64, c f64) []f64 { q_ := (a * a - 3.0 * b) r_ := (2.0 * a * a * a - 9.0 * a * b + 27.0 * c) @@ -89,15 +105,6 @@ pub fn solve_cubic(a f64, b f64, c f64) []f64 { if r == 0.0 && q == 0.0 { return [-a / 3.0, -a / 3.0, -a / 3.0] } else if cr2 == cq3 { - /* - this test is actually r2 == q3, written in a form suitable - for exact computation with integers - */ - /* - Due to finite precision some double roots may be missed, and - considered to be a pair of complex roots z = x +/- epsilon i - close to the real axis. - */ sqrt_q := math.sqrt(q) if r > 0.0 { return [-2.0 * sqrt_q - a / 3.0, sqrt_q - a / 3.0, sqrt_q - a / 3.0] @@ -122,11 +129,17 @@ pub fn solve_cubic(a f64, b f64, c f64) []f64 { } } +// swap_ swaps two numbers: f(a, b) = (b, a) +// Input: a, b +// Output: (b, a) @[inline] fn swap_(a f64, b f64) (f64, f64) { return b, a } +// sorted_3_ sorts three numbers in ascending order: f(x, y, z) = (min(x,y,z), median(x,y,z), max(x,y,z)) +// Input: x, y, z +// Output: (min, median, max) @[inline] fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) { mut x := x_ @@ -144,6 +157,15 @@ fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) { return x, y, z } +// companion_matrix creates a companion matrix for the polynomial P(x) = a_n * x^n + a_{n-1} * x^{n-1} + ... + a_1 * x + a_0 +// The companion matrix C is defined as: +// [0 0 0 ... 0 -a_0/a_n] +// [1 0 0 ... 0 -a_1/a_n] +// [0 1 0 ... 0 -a_2/a_n] +// [. . . ... . ........] +// [0 0 0 ... 1 -a_{n-1}/a_n] +// Input: a = [a_0, a_1, ..., a_n] +// Output: Companion matrix C pub fn companion_matrix(a []f64) [][]f64 { nc := a.len - 1 mut cm := [][]f64{len: nc, init: []f64{len: nc}} @@ -162,6 +184,10 @@ pub fn companion_matrix(a []f64) [][]f64 { return cm } +// balance_companion_matrix balances a companion matrix C to improve numerical stability +// Uses an iterative scaling process to make the row and column norms as close to each other as possible +// Input: Companion matrix C +// Output: Balanced matrix B such that D^(-1)CD = B, where D is a diagonal matrix pub fn balance_companion_matrix(cm [][]f64) [][]f64 { nc := cm.len mut m := cm.clone() @@ -170,7 +196,7 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 { mut col_norm := 0.0 for not_converged { not_converged = false - for i := 0; i < nc; i++ { // column norm, excluding the diagonal + for i := 0; i < nc; i++ { if i != nc - 1 { col_norm = math.abs(m[i + 1][i]) } else { @@ -178,7 +204,7 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 { for j := 0; j < nc - 1; j++ { col_norm += math.abs(m[j][nc - 1]) } - } // row norm, excluding the diagonal + } if i == 0 { row_norm = math.abs(m[0][nc - 1]) } else if i == nc - 1 { @@ -223,88 +249,68 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 { return m } -// Arithmetic operations on polynomials -// -// In the following descriptions a, b, c are polynomials of degree -// na, nb, nc respectively. -// -// Each polynomial is represented by an array containing its -// coefficients, together with a separately declared integer equal -// to the degree of the polynomial. The coefficients appear in -// ascending order; that is, -// -// a(x) = a[0] + a[1] * x + a[2] * x^2 + ... + a[na] * x^na . -// -// -// -// sum = eval( a, x ) Evaluate polynomial a(t) at t = x. -// c = add( a, b ) c = a + b, nc = max(na, nb) -// c = sub( a, b ) c = a - b, nc = max(na, nb) -// c = mul( a, b ) c = a * b, nc = na+nb -// -// -// Division: -// -// c = div( a, b ) c = a / b, nc = MAXPOL -// -// returns i = the degree of the first nonzero coefficient of a. -// The computed quotient c must be divided by x^i. An error message -// is printed if a is identically zero. -// -// -// Change of variables: -// If a and b are polynomials, and t = a(x), then -// c(t) = b(a(x)) -// is a polynomial found by substituting a(x) for t. The -// subroutine call for this is -// +// add adds two polynomials: (a_n * x^n + ... + a_0) + (b_m * x^m + ... + b_0) +// Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m] +// Output: [a_0 + b_0, a_1 + b_1, ..., max(a_k, b_k), ...] pub fn add(a []f64, b []f64) []f64 { - na := a.len - nb := b.len - nc := int(math.max(na, nb)) - mut c := []f64{len: nc} - for i := 0; i < nc; i++ { - if i > na { - c[i] = b[i] - } else if i > nb { - c[i] = a[i] - } else { - c[i] = a[i] + b[i] - } + mut result := []f64{len: math.max(a.len, b.len)} + for i in 0 .. result.len { + result[i] = if i < a.len { a[i] } else { 0.0 } + if i < b.len { b[i] } else { 0.0 } } - return c + return result } +// subtract subtracts two polynomials: (a_n * x^n + ... + a_0) - (b_m * x^m + ... + b_0) +// Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m] +// Output: [a_0 - b_0, a_1 - b_1, ..., a_k - b_k, ...] pub fn subtract(a []f64, b []f64) []f64 { - na := a.len - nb := b.len - nc := int(math.max(na, nb)) - mut c := []f64{len: nc} - for i := 0; i < nc; i++ { - if i > na { - c[i] = -b[i] - } else if i > nb { - c[i] = a[i] - } else { - c[i] = a[i] - b[i] - } + mut result := []f64{len: math.max(a.len, b.len)} + for i in 0 .. result.len { + result[i] = if i < a.len { a[i] } else { 0.0 } - if i < b.len { b[i] } else { 0.0 } } - return c + return result } +// multiply multiplies two polynomials: (a_n * x^n + ... + a_0) * (b_m * x^m + ... + b_0) +// Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m] +// Output: [c_0, c_1, ..., c_{n+m}] where c_k = βˆ‘_{i+j=k} a_i * b_j pub fn multiply(a []f64, b []f64) []f64 { - na := a.len - nb := b.len - nc := na + nb - mut c := []f64{len: nc} - for i := 0; i < na; i++ { - x := a[i] - for j := 0; j < nb; j++ { - k := i + j - c[k] += x * b[j] + mut result := []f64{len: a.len + b.len - 1} + for i in 0 .. a.len { + for j in 0 .. b.len { + result[i + j] += a[i] * b[j] } } - return c + return result +} + +// divide divides two polynomials: (a_n * x^n + ... + a_0) / (b_m * x^m + ... + b_0) +// Uses polynomial long division algorithm +// Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m] +// Output: (q, r) where q is the quotient and r is the remainder +// such that a(x) = b(x) * q(x) + r(x) and degree(r) < degree(b) +pub fn divide(a []f64, b []f64) ([]f64, []f64) { + mut quotient := []f64{} + mut remainder := a.clone() + b_lead_coef := b[0] + + for remainder.len >= b.len { + lead_coef := remainder[0] / b_lead_coef + quotient << lead_coef + for i in 0 .. b.len { + remainder[i] -= lead_coef * b[i] + } + remainder = unsafe { remainder[1..] } + for remainder.len > 0 && math.abs(remainder[0]) < 1e-10 { + remainder = unsafe { remainder[1..] } + } + } + + if remainder.len == 0 { + remainder = []f64{} + } + + return quotient, remainder } pub fn divide(dividend []f64, divisor []f64) ([]f64, []f64) { diff --git a/poly/poly_test.v b/poly/poly_test.v index 1bde71a5c..13573b3c0 100644 --- a/poly/poly_test.v +++ b/poly/poly_test.v @@ -1,9 +1,8 @@ module poly +import math + fn test_eval() { - // ans = 2 - // ans = 5 + 4 * 2 = 13 - // ans = 4 + 4 * 13 = 56 x := 4 coef := [4.0, 5, 2] assert eval(coef, 4) == 56 @@ -59,26 +58,37 @@ fn test_balance_companion_matrix() { } fn test_add() { - a := [6.0, 777, -3] - b := [1.0, -755, -4] - assert add(a, b) == [7.0, 22, -7] + a := [1.0, 2.0, 3.0] + b := [6.0, 20.0, -10.0] + result := add(a, b) + println('Add result: ${result}') + assert result == [7.0, 22.0, -7.0] } fn test_subtract() { - a := [6.0, 777, -3] - b := [1.0, -755, -4] - assert subtract(a, b) == [5.0, 1532, 1] + a := [6.0, 1532.0, -4.0] + b := [1.0, -1.0, -5.0] + result := subtract(a, b) + println('Subtract result: ${result}') + assert result == [5.0, 1533.0, 1.0] } fn test_multiply() { - a := [9.0, -1, 5] - b := [0.0, -1, 7] - assert multiply(a, b) == [0.0, -9, 64, -12, 35, 0] + // (2+3x+4x^2) * (-3x+2x^2) = (-6x -5x^2 -6x^3 + 8x^4) + a := [2.0, 3.0, 4.0] + b := [0.0, -3.0, 2.0] + result := multiply(a, b) + println('Multiply result: ${result}') + assert result == [0.0, -6.0, -5.0, -6.0, 8.0] } -fn test_division() { - dividend := [6.0, -5.0, 1.0] - divisor := [2.0, 1.0] - quotient, remainder := divide(dividend, divisor) - assert quotient == [-7.0, 1] && remainder == [20.0] +fn test_divide() { + // (x^2 + 2x + 1) / (x + 1) = (x + 1) + a := [1.0, 2.0, 1.0] + b := [1.0, 1.0] + quotient, remainder := divide(a, b) + println('Divide quotient: ${quotient}') + println('Divide remainder: ${remainder}') + assert quotient == [1.0, 1.0] + assert remainder == [] // The empty set indicates that two polynomials divide each other exactly (without remainder). } diff --git a/vcl/image.c.v b/vcl/image.c.v index b2d552d0a..d06fe7cde 100644 --- a/vcl/image.c.v +++ b/vcl/image.c.v @@ -17,7 +17,7 @@ pub interface IImage { width int height int nr_channels int - data voidptr + data &u8 } // Image memory buffer on the device with image data diff --git a/vlas/lapack_common.v b/vlas/lapack_common.v index 8a30e8f5e..578c34c93 100644 --- a/vlas/lapack_common.v +++ b/vlas/lapack_common.v @@ -74,9 +74,13 @@ pub fn dgesv(n int, nrhs int, mut a []f64, lda int, ipiv []int, mut b []f64, ldb // Note that the routine returns V**T, not V. // // NOTE: matrix 'a' will be modified -pub fn dgesvd(jobu &char, jobvt &char, m int, n int, a []f64, lda int, s []f64, u []f64, ldu int, vt []f64, ldvt int, superb []f64) { - info := C.LAPACKE_dgesvd(.row_major, jobu, jobvt, m, n, &a[0], lda, &s[0], &u[0], - ldu, &vt[0], ldvt, &superb[0]) +pub fn dgesvd(jobu rune, jobvt rune, m int, n int, a []f64, lda int, s []f64, u []f64, ldu int, vt []f64, ldvt int, superb []f64) { + // lapack_int LAPACKE_dgesvd( int matrix_order, char jobu, char jobvt, + // lapack_int m, lapack_int n, double* a, + // lapack_int lda, double* s, double* u, lapack_int ldu, + // double* vt, lapack_int ldvt, double* superb ); + info := C.LAPACKE_dgesvd(.row_major, char(jobu), char(jobvt), m, n, &a[0], lda, &s[0], + &u[0], ldu, &vt[0], ldvt, &superb[0]) if info != 0 { errors.vsl_panic('lapack failed', .efailed) } @@ -189,7 +193,7 @@ pub fn dgeev(calc_vl bool, calc_vr bool, n int, mut a []f64, lda int, wr []f64, ldvr = 1 } unsafe { - info := C.LAPACKE_dgeev(.row_major, &char(job_vlr(calc_vl).str().str), &char(job_vlr(calc_vr).str().str), + info := C.LAPACKE_dgeev(.row_major, char(job_vlr(calc_vl)), char(job_vlr(calc_vr)), n, &a[0], lda, &wr[0], &wi[0], &vvl, ldvl, &vvr, ldvr) if info != 0 { errors.vsl_panic('lapack failed', .efailed) diff --git a/vlas/lapack_macos.c.v b/vlas/lapack_macos.c.v index 2c0cbc7c3..1636cc01b 100644 --- a/vlas/lapack_macos.c.v +++ b/vlas/lapack_macos.c.v @@ -1,7 +1,15 @@ module vlas -fn C.LAPACKE_dlange(norm &char, m int, n int, a &f64, lda int, work &f64) f64 +import vsl.vlas.internal.blas + +// double LAPACKE_dlange( int matrix_order, char norm, lapack_int m, +// lapack_int n, const double* a, lapack_int lda ); + +// double LAPACKE_dlange_work( int matrix_order, char norm, lapack_int m, +// lapack_int n, const double* a, lapack_int lda, double* work ); + +fn C.LAPACKE_dlange_work(matrix_order blas.MemoryLayout, norm char, m int, n int, const_a &f64, lda int, work &f64) f64 pub fn dlange(norm rune, m int, n int, a []f64, lda int, work []f64) f64 { - return unsafe { C.LAPACKE_dlange(&char(norm.str().str), m, n, &a[0], lda, &work[0]) } + return unsafe { C.LAPACKE_dlange_work(.row_major, char(norm), m, n, &a[0], lda, &work[0]) } }