Perceptron Polynomial Regression

perceptron-regression.red

Red [
    Title:  "Perceptron Polynomial Regression Engine"
    Author: "hinjolicious"
    File:   %polynomial-regression.red
    Needs:  'view
    Notes:  {
        A pure, single-pass Stochastic Gradient Descent (SGD) engine 
        that discovers the hidden coefficients of polynomial data 
        using dynamic feature mapping and real-time convergence tracking.
    }
	Resources: "Gemini AI, P5js Code from Open-Processing, Nature of Code, etc."
]

; ==============================================================================
; == 1. CONFIGURATION & PARAMETERS
; ==============================================================================

width:  400
height: 400

random/seed now/time/precise

; Hyperparameters
learning-rate: 0.075          ; Optimal zone for stable gradient descent
threshold:     0.000001      ; Precision limit to determine mathematical convergence

data-size: 20 ; data points: random from 3 to data-size!

; ==============================================================================
; == 2. CORE MATHEMATICAL CORE (METAPROGRAMMING & UTILITIES)
; ==============================================================================

; Helper to dynamically compile a custom quadratic polynomial function block at runtime
make-poly: func [a b c /local bd][
    bd: [a * (x ** 2) + (b * x) + c]
    bd/1: a
    bd/5/1: b
    bd/7: c
    func [x] bd
]

; The underlying hidden target formula we want the model to discover
f: make-poly (2 - random 4.0) (2 - random 4.0) (0.5 - random 1.0)

; Linear interpolation mapper to scale values between coordinate spaces
map: func [v a1 z1 a2 z2 /local r][
    r: (v - a1) / (z1 - a1)
    a2 + (r * (z2 - a2))
]

; Object factory helper that instantiates a block and executes its init constructor
new: func [obj args][
    inst: make obj []
    if in inst 'init [apply :inst/init args]
    inst
]

; ==============================================================================
; == 3. MACHINE LEARNING OBJECTS
; ==============================================================================

REGRESSOR: make object! [
    weights: copy []
    
    ; Initializes weight vectors randomly between -1.0 and 1.0
    init: func [n][
        weights: collect [loop n [keep (random 2.0) - 1.0]]
    ]
    
    ; Computes the dot product of given feature inputs and internal weights
    guess: func [inputs [block!] /local sum][
        sum: 0.0
        repeat i (length? inputs) [
            sum: sum + (inputs/:i * weights/:i)
        ]
        sum
    ]
    
    ; Evaluates a single spatial X coordinate using current discovered coefficients
    ; Mapping convention: weights/1 = b (x), weights/2 = a (x²), weights/3 = c (bias)
    guess-y: func [x][
        (weights/1 * x) + (weights/2 * x * x) + weights/3
    ]
    
    ; Adjusts internal weights using Delta Rule / Stochastic Gradient Descent
    train: func [inputs target /local error][
        error: target - guess inputs
        repeat i length? weights [
            weights/:i: weights/:i + (error * inputs/:i * learning-rate)
        ]
    ]
]

POINT: make object! [
    x: 0.0  y: 0.0  target: 0.0  bias: 1.0
    
    init: func [_x _y][
        x: _x  y: _y  
        target: y  ; Target target is the true vertical coordinate
    ]
    
    ; Map internal normalized coordinates (-1.0 to 1.0) to actual visual UI pixels
    pixel-x: func [][map x -1.0 1.0 0.0 width]
    pixel-y: func [][map y -1.0 1.0 height 0.0]
    
    ; Append visual representation of the point to the drawing block
    show: func [/local pxy][
        pxy: as-pair pixel-x pixel-y
        append blk compose [fill-pen 255.0.0.100 circle (pxy) 5]
    ]
]

; == other stuff

regression: function [
    "Fit a quadratic y = a + bx + cx^2 to data points via least squares. Prints coefficients and residuals."
    xa [block! vector!] "Block of x values"
    ya [block! vector!] "Block of y values; must match length of xa"
][
    n: length? xa
    ;; accumulate raw moment sums
    xm: ym: x2m: x3m: x4m: xym: x2ym: 0.0

    repeat i n [
        xi: xa/:i
        yi: ya/:i
        xm:   xm   +  xi
        ym:   ym   +  yi
        x2m:  x2m  + (xi * xi)
        x3m:  x3m  + (xi * xi * xi)
        x4m:  x4m  + (xi * xi * xi * xi)
        xym:  xym  + (xi * yi)
        x2ym: x2ym + (xi * xi * yi)
    ]

    ;; convert sums to means
    xm:   xm   / n
    ym:   ym   / n
    x2m:  x2m  / n
    x3m:  x3m  / n
    x4m:  x4m  / n
    xym:  xym  / n
    x2ym: x2ym / n

    ;; central moments (variance/covariance terms)
    sxx:   x2m  - (xm  * xm)
    sxy:   xym  - (xm  * ym)
    sxx2:  x3m  - (xm  * x2m)
    sx2x2: x4m  - (x2m * x2m)
    sx2y:  x2ym - (x2m * ym)

    ;; solve 3x3 normal equations for a, b, c
    denom:  sxx  * sx2x2 - (sxx2 * sxx2)
        b: (sxy  * sx2x2 - (sx2y * sxx2)) / denom
        a: (sx2y * sxx   - (sxy  * sxx2)) / denom
        c: ym - (b * xm) - (a * x2m)

	reduce [a b c] ; return polynomial coefficients
]

; ==============================================================================
; == 4. DATASET GENERATION & INITIALIZATION
; ==============================================================================

; Generate scattered data points anchored to the hidden formula with artificial noise
data-size: (random 17) + 3
points: collect [
	repeat step data-size [
		x: map (step - 1) 0.0 (data-size - 1) -1.0 1.0
		noise: ((random 30.0) / 100.0) - 0.15
		keep new POINT [x (f x) + noise]
	]
]    
        
; Instantiate the Perceptron Brain with 3 weight nodes: [x, x², bias]
brain: new REGRESSOR [3]

; Global drawing canvas script block
blk: copy [] 

; Convergence status trackers
is-converged?:        false
previous-epoch-error: 0.0

; == compare with mathematical polynomial regression

xa: copy []  ya: copy []
foreach p points [ append xa p/x  append ya p/y ]
coeff: regression xa ya
reg: make-poly coeff/1 coeff/2 coeff/3

; ==============================================================================
; == 5. RENDERING & OPTIMIZATION PIPELINE
; ==============================================================================

update: does [
    clear blk
	append blk [pen off]

    ; Step 1: Render historical scatter plot data
    repeat i length? points [points/(i)/show 
		;print [i ":" points/(i)/x " " points/(i)/y]
		]
    
    ; Step 2: Render the true underlying mathematical distribution (Light Gray)
    ideal-curve: collect [
        repeat step 21 [
            nx: map (step - 1) 0.0 20.0 -1.0 1.0
            ny: f nx
            keep as-pair map nx -1.0 1.0 0.0 400.0 map ny -1.0 1.0 400.0 0.0
        ]
    ]
    append blk compose [fill-pen off line-width 12 pen (200.200.200.150) line (ideal-curve)]
	
	; overlay with mathematical polynomial regression result
	
    reg-curve: collect [
        repeat step 21 [
            nx: map (step - 1) 0.0 20.0 -1.0 1.0
            ny: reg nx
            keep as-pair map nx -1.0 1.0 0.0 400.0 map ny -1.0 1.0 400.0 0.0
        ]
    ]
    append blk compose [fill-pen off line-width 6 pen (0.200.200.150) line (reg-curve)]	

    ; Step 3: Combined Analysis Pass (Calculates Error Matrix & Trains Simultaneously)
    current-epoch-error: 0.0
    foreach pt points [
        guess: brain/guess reduce [pt/x (pt/x ** 2.0) pt/bias]
        error: absolute (pt/target - guess)
        current-epoch-error: current-epoch-error + error
        
        if not is-converged? [
            brain/train reduce [pt/x (pt/x ** 2.0) pt/bias] pt/target
        ]
    ]
    current-epoch-error: current-epoch-error / (length? points)

    ; Step 4: Evaluate Convergence Velocity against threshold
    error-delta: absolute (current-epoch-error - previous-epoch-error)
    if (error-delta < threshold) [ is-converged?: true ]
    previous-epoch-error: current-epoch-error    
    
    ; Step 5: Render Brain's current optimized predictive path (Blue Line)
    fitted-curve: collect [
        repeat step 21 [
            nx: map (step - 1) 0.0 20.0 -1.0 1.0
            ny: brain/guess-y nx
            keep as-pair map nx -1.0 1.0 0.0 400.0 map ny -1.0 1.0 400.0 0.0
        ]
    ]
    append blk compose [fill-pen off line-width 1 pen blue line (fitted-curve)]
    
    ; Step 6: Formulate and display algebraic expression on the GUI view screen
    pa: round/to brain/weights/2 0.01    ; Quad coefficient (a)
    pb: round/to brain/weights/1 0.01    ; Linear coefficient (b)   
    pc: round/to brain/weights/3 0.01    ; Constant intercept (c)
    
    poly: rejoin [
        "Regression: y = " pa "x^^2" 
        either pb > 0 [rejoin [" + " pb]][rejoin [" - " negate pb]] "x"
        either pc > 0 [rejoin [" + " pc]][rejoin [" - " negate pc]]
    ]
	
	; High average distance = low score, zero average distance = 100% score
	max-expected-error: 0.5 ; The maximum messy deviation you'd expect in your coordinate space
	fit-confidence: 1.0 - (current-epoch-error / max-expected-error)

	; Convert to a clean percentage
	confidence-pct: rejoin [round/to (max 0.0 (fit-confidence * 100.0)) 0.1 "%"]	
	
	bot: height - 20
	acc: rejoin ["Confidence: " confidence-pct]
	datsiz: rejoin ["Data size: " data-size]
    append blk compose [ 
		pen black text 10x10 (poly) 
		text (as-pair 10 bot) "Click to restart"
		text (as-pair 10 bot - 15) (acc)
		text (as-pair 10 bot - 30) (datsiz)
	]
    
    ; Step 7: Freeze engine when stable convergence state is unlocked
    if is-converged? [
        append blk compose [ pen black text 10x25 "Fitting done!" ]
        canv/rate: none
        show canv
    ]
]

reset: does [
    ; 1. Generate a brand new true underlying polynomial function
    f: make-poly (2 - random 4.0) (2 - random 4.0) (0.5 - random 1.0)
    
    ; 2. Generate a fresh block of noisy scatter data points
	data-size: (random 17) + 3
	points: collect [
		repeat step data-size [
			x: map (step - 1) 0.0 (data-size - 1) -1.0 1.0
			noise: ((random 30.0) / 100.0) - 0.15
			keep new POINT [x (f x) + noise]
		]
	]    
	
    ; 3. Re-initialize the Perceptron Brain weights randomly
    brain/init 3
    
    ; 4. Reset tracking variables
    is-converged?: false
    previous-epoch-error: 0.0
	
	xa: copy []  ya: copy []
	foreach p points [ append xa p/x  append ya p/y ]
	coeff: regression xa ya
	reg: make-poly coeff/1 coeff/2 coeff/3	
    
    ; 5. Wake up the canvas timer loop (set back to 60 FPS)
    canv/rate: 60
]

; ==============================================================================
; == 6. GRAPHICAL USER INTERFACE LAYOUT
; ==============================================================================

view/tight compose [
    title "Perceptron Polynomial Data Fitting"
    canv: base (as-pair width height) white draw blk 
	on-down [reset]
    rate 60 on-time [update]
	;do [reset]
]

Updated!
Now data points will have an even spacing across the x axis.
Number of data vary randomly from 3 to 20 items.