diff --git a/pyscriptjs/examples/fractals.py b/pyscriptjs/examples/fractals.py
index 8d4ad78e..65e0f31a 100644
--- a/pyscriptjs/examples/fractals.py
+++ b/pyscriptjs/examples/fractals.py
@@ -1,4 +1,6 @@
+from typing import Tuple
import numpy as np
+from numpy.polynomial import Polynomial
def mandelbrot(width: int, height: int, *,
x: float = -0.5, y: float = 0, zoom: int = 1, max_iterations: int = 100) -> np.array:
@@ -60,3 +62,48 @@ def julia(width: int, height: int, *,
div_time[m] = i
return div_time
+
+Range = Tuple[float, float]
+
+def newton(width: int, height: int, *,
+ p: Polynomial, a: complex, xr: Range = (-2.5, 1), yr: Range = (-1, 1), max_iterations: int = 100) -> (np.array, np.array):
+ """ """
+ # To make navigation easier we calculate these values
+ x_from, x_to = xr
+ y_from, y_to = yr
+
+ # Here the actual algorithm starts
+ x = np.linspace(x_from, x_to, width).reshape((1, width))
+ y = np.linspace(y_from, y_to, height).reshape((height, 1))
+ z = x + 1j*y
+
+ # Compute the derivative
+ dp = p.deriv()
+
+ # Compute roots
+ roots = p.roots()
+ epsilon = 1e-5
+
+ # Set the initial conditions
+ a = np.full(z.shape, a)
+
+ # To keep track in which iteration the point diverged
+ div_time = np.zeros(z.shape, dtype=int)
+
+ # To keep track on which points did not converge so far
+ m = np.full(a.shape, True, dtype=bool)
+
+ # To keep track which root each point converged to
+ r = np.full(a.shape, 0, dtype=int)
+
+ for i in range(max_iterations):
+ z[m] = z[m] - a[m]*p(z[m])/dp(z[m])
+
+ for j, root in enumerate(roots):
+ converged = (np.abs(z.real - root.real) < epsilon) & (np.abs(z.imag - root.imag) < epsilon)
+ m[converged] = False
+ r[converged] = j + 1
+
+ div_time[m] = i
+
+ return div_time, r
diff --git a/pyscriptjs/examples/numpy_canvas_fractals.html b/pyscriptjs/examples/numpy_canvas_fractals.html
index cbee92fb..2247747d 100644
--- a/pyscriptjs/examples/numpy_canvas_fractals.html
+++ b/pyscriptjs/examples/numpy_canvas_fractals.html
@@ -1,6 +1,6 @@
- Visualization of Mandelbrot and Julia sets with NumPy and HTML5 canvas
+ Visualization of Mandelbrot, Julia and Newton sets with NumPy and HTML5 canvas
@@ -17,6 +17,10 @@
animation: spin 1s ease-in-out infinite;
}
+ canvas {
+ display: none;
+ }
+
@keyframes spin {
to {
transform: rotate(360deg);
@@ -28,35 +32,68 @@
-
-
+
+
Mandelbrot set
-
-
+
- numpy
+- sympy
- paths:
- /palettes.py
- /fractals.py
-from pyodide import to_js
+from pyodide import to_js, create_proxy
import numpy as np
+import sympy
from palettes import Magma256
-from fractals import mandelbrot, julia
+from fractals import mandelbrot, julia, newton
from js import (
console,
@@ -65,28 +102,20 @@ from js import (
ImageData,
Uint8ClampedArray,
CanvasRenderingContext2D as Context2d,
+ requestAnimationFrame,
)
-def create_canvas(width: int, height: int, target: str) -> Context2d:
- pixel_ratio = devicePixelRatio
-
- canvas = document.createElement("canvas")
+def prepare_canvas(width: int, height: int, canvas: Element) -> Context2d:
ctx = canvas.getContext("2d")
canvas.style.width = f"{width}px"
canvas.style.height = f"{height}px"
- canvas.width = width*pixel_ratio
- canvas.height = height*pixel_ratio
-
- ctx.scale(pixel_ratio, pixel_ratio)
- ctx.translate(0.5, 0.5)
+ canvas.width = width
+ canvas.height = height
ctx.clearRect(0, 0, width, height)
- el = document.querySelector(target)
- el.replaceChildren(canvas)
-
return ctx
def color_map(array: np.array, palette: np.array) -> np.array:
@@ -105,30 +134,184 @@ def draw_image(ctx: Context2d, image: np.array) -> None:
image_data = ImageData.new(data, width, height)
ctx.putImageData(image_data, 0, 0)
-def draw_mandelbrot(width: int, height: int) -> None:
- ctx = create_canvas(width, height, "#mandelbrot")
+width, height = 600, 600
+
+async def draw_mandelbrot() -> None:
+ spinner = document.querySelector("#mandelbrot .loading")
+ canvas = document.querySelector("#mandelbrot canvas")
+
+ spinner.style.display = ""
+ canvas.style.display = "none"
+
+ ctx = prepare_canvas(width, height, canvas)
console.log("Computing Mandelbrot set ...")
console.time("mandelbrot")
- array = mandelbrot(width, height)
+ iters = mandelbrot(width, height)
console.timeEnd("mandelbrot")
- image = color_map(array, Magma256)
+ image = color_map(iters, Magma256)
draw_image(ctx, image)
-def draw_julia(width: int, height: int) -> None:
- ctx = create_canvas(width, height, "#julia")
+ spinner.style.display = "none"
+ canvas.style.display = "block"
+
+async def draw_julia() -> None:
+ spinner = document.querySelector("#julia .loading")
+ canvas = document.querySelector("#julia canvas")
+
+ spinner.style.display = ""
+ canvas.style.display = "none"
+
+ ctx = prepare_canvas(width, height, canvas)
console.log("Computing Julia set ...")
console.time("julia")
- array = julia(width, height)
+ iters = julia(width, height)
console.timeEnd("julia")
- image = color_map(array, Magma256)
+ image = color_map(iters, Magma256)
draw_image(ctx, image)
-draw_mandelbrot(600, 600)
-draw_julia(600, 600)
+ spinner.style.display = "none"
+ canvas.style.display = "block"
+
+def ranges():
+ x0_in = document.querySelector("#x0")
+ x1_in = document.querySelector("#x1")
+ y0_in = document.querySelector("#y0")
+ y1_in = document.querySelector("#y1")
+
+ xr = (float(x0_in.value), float(x1_in.value))
+ yr = (float(y0_in.value), float(y1_in.value))
+
+ return xr, yr
+
+current_image = None
+async def draw_newton() -> None:
+ spinner = document.querySelector("#newton .loading")
+ canvas = document.querySelector("#newton canvas")
+
+ spinner.style.display = ""
+ canvas.style.display = "none"
+
+ ctx = prepare_canvas(width, height, canvas)
+
+ console.log("Computing Newton set ...")
+
+ poly_in = document.querySelector("#poly")
+ coef_in = document.querySelector("#coef")
+ conv_in = document.querySelector("#conv")
+ iter_in = document.querySelector("#iter")
+
+ xr, yr = ranges()
+
+ # z**3 - 1
+ # z**8 + 15*z**4 - 16
+ # z**3 - 2*z + 2
+
+ expr = sympy.parse_expr(poly_in.value)
+ coeffs = [ complex(c) for c in reversed(sympy.Poly(expr, sympy.Symbol("z")).all_coeffs()) ]
+ poly = np.polynomial.Polynomial(coeffs)
+
+ coef = complex(sympy.parse_expr(coef_in.value))
+
+ console.time("newton")
+ iters, roots = newton(width, height, p=poly, a=coef, xr=xr, yr=yr)
+ console.timeEnd("newton")
+
+ if conv_in.checked:
+ n = poly.degree() + 1
+ k = int(len(Magma256)/n)
+
+ colors = Magma256[::k, :][:n]
+ colors[0, :] = [255, 0, 0] # red: no convergence
+
+ image = color_map(roots, colors)
+ else:
+ image = color_map(iters, Magma256)
+
+ global current_image
+ current_image = image
+ draw_image(ctx, image)
+
+ spinner.style.display = "none"
+ canvas.style.display = "block"
+
+handler = create_proxy(lambda _event: draw_newton())
+document.querySelector("#newton fieldset").addEventListener("change", handler)
+
+canvas = document.querySelector("#newton canvas")
+
+is_selecting = False
+init_sx, init_sy = None, None
+sx, sy = None, None
+async def mousemove(event):
+ global is_selecting
+ global init_sx
+ global init_sy
+ global sx
+ global sy
+
+ def invert(sx, source_range, target_range):
+ source_start, source_end = source_range
+ target_start, target_end = target_range
+ factor = (target_end - target_start)/(source_end - source_start)
+ offset = -(factor * source_start) + target_start
+ return (sx - offset) / factor
+
+ bds = canvas.getBoundingClientRect()
+ event_sx, event_sy = event.clientX - bds.x, event.clientY - bds.y
+
+ ctx = canvas.getContext("2d")
+
+ pressed = event.buttons == 1
+ if is_selecting:
+ if not pressed:
+ xr, yr = ranges()
+
+ x0 = invert(init_sx, xr, (0, width))
+ x1 = invert(sx, xr, (0, width))
+ y0 = invert(init_sy, yr, (0, height))
+ y1 = invert(sy, yr, (0, height))
+
+ document.querySelector("#x0").value = x0
+ document.querySelector("#x1").value = x1
+ document.querySelector("#y0").value = y0
+ document.querySelector("#y1").value = y1
+
+ is_selecting = False
+ init_sx, init_sy = None, None
+ sx, sy = init_sx, init_sy
+
+ await draw_newton()
+ else:
+ ctx.save()
+ ctx.clearRect(0, 0, width, height)
+ draw_image(ctx, current_image)
+ sx, sy = event_sx, event_sy
+ ctx.beginPath()
+ ctx.rect(init_sx, init_sy, sx - init_sx, sy - init_sy)
+ ctx.fillStyle = "rgba(255, 255, 255, 0.4)"
+ ctx.strokeStyle = "rgba(255, 255, 255, 1.0)"
+ ctx.fill()
+ ctx.stroke()
+ ctx.restore()
+ else:
+ if pressed:
+ is_selecting = True
+ init_sx, init_sy = event_sx, event_sy
+ sx, sy = init_sx, init_sy
+
+canvas.addEventListener("mousemove", create_proxy(mousemove))
+
+import asyncio
+
+_ = await asyncio.gather(
+ draw_mandelbrot(),
+ draw_julia(),
+ draw_newton(),
+)