Bugs in "K-Means" BASIC code

[thecloudexpanse]

@mdcramer
I attempted to follow along with the "K-Means" program described in https://mdcramer.github.io/apple-2-blog/k-means/, once I had the previous "Refactored" program running successfully. There appears to be an omission in the published BASIC source and it is unable to successfully converge.

The "Draw Decision Boundary" code calls a non-existent subroutine at line 3060:
3060 GOSUB 3500

It's not clear what this subroutine should do or whether some previously presented code has been renumbered. If I just remove it, e.g.:
3060 REM GOSUB 3500 -- WHAT DOES THIS DO?

then the program will execute but it never converges, producing ACCURACY = 50% output for multiple iterations until it eventually hits a bad subscript error.

Also, I read that you are using ChatGPT to scan photos of the screen. Looks like it has trouble distinguishing 0 (zero) from O (upper-case O) on these lines:

- 2620   K0(J,0) = KM(J,0)
- 2630   K0(J,1) = KM(J,1)
+ 2620   KO(J,0) = KM(J,0)
+ 2630   KO(J,1) = KM(J,1)

UPDATE: After correcting the additional typos found on this line, the k-means program now successfully converges, though it is still reporting an ACCURACY = 50% so not really trusting the results yet.

- 2830   DI = DI + (KM%(I,0) - KO%(I,0)) ^ 2 + (KM%(I,1) - KO%(I,1)) ^ 2
+ 2830   DI = DI + (KM(I,0) - KO(I,0)) ^ 2 + (KM(I,1) - KO(I,1)) ^ 2

UPDATE 2: I think it's fixed now. Those were not the only references to a KM% array (instead of just KM). Once I changed all the references, it seems to be computing convergence correctly!

2690   KM(Y,0) = KM(Y,0) + DS%(I,0)
2700   KM(Y,1) = KM(Y,1) + DS%(I,1)
2710   KM(Y,2) = KM(Y,2) + 1

2740   KM(I,0) = KM(I,0) / KM(I,2)
2750   KM(I,1) = KM(I,1) / KM(I,2)

3030   IF KM(I - 1,1) - KM(I,1) <> 0 THEN M = -1 * (KM(I - 1,0) - KM(I,0)) / (KM(I - 1,1) - KM(I,1))
3040   P%(0,0) = (KM(I,0) - KM(I - 1,0)) / 2 + KM(I - 1,0)
3050   P%(0,1) = (KM(I,1) - KM(I - 1,1)) / 2 + KM(I - 1,1)

[thecloudexpanse]

Here is the full source of the "K-Means" program that I've been able to get running so far.

10  HOME : VTAB 21
20  PI = 3.14159265
30  GOSUB 1000 : REM  DRAW AXIS
40  GOSUB 100 : REM  GENERATE DATA
50  GOSUB 900 : REM  WAIT FOR KEY
60  GOSUB 2000 : REM  RUN K-MEANS
70  END
100 REM == HYPERPARAMETERS ==
110 KN = 2: REM # OF CLASSES
120 DIM KN%(KN - 1): REM # OF SAMPLES FROM EACH CLASS
130 KN%(0) = 5
140 KN%(1) = 5
150 REM -- GENERATORS --
160 REM -- CLASS - MU-X0, MU-X1, SIG-X0, SIG-X1, RHO --
170 DIM GE(KN - 1,5)
180 REM -- CLASS 0 --
190 GE(0,0) = 100: REM MU-X0
200 GE(0,1) = 50: REM MU-X1
210 GE(0,2) = 30: REM SIG-X0
220 GE(0,3) = 10: REM SIG-X1
230 GE(0,4) = 0.5: REM RHO
240 GE(0,5) = SQR(1 - GE(0,4) ^ 2)
250 REM -- CLASS 1 --
260 GE(1,0) = 200: REM MU-X0
270 GE(1,1) = 100: REM MU-X1
280 GE(1,2) = 50: REM SIG-X0
290 GE(1,3) = 25: REM SIG-X1
300 GE(1,4) = -0.5: REM RHO
310 GE(1,5) = SQR(1 - GE(1,4) ^ 2)
320 REM == END HYPERPARAMETERS ==
330 NS = 0: REM TOTAL # SAMPLES
340 FOR I = 0 TO KN - 1
350   NS = NS + KN%(I)
360 NEXT
370 REM -- DATA TABLE --
380 REM -- # - X0, X1, Y-HAT, Y --
390 DIM DS%(NS - 1,3)
400 REM -- REMAINING SAMPLES COUNT --
410 DIM RK%(KN - 1)
420 FOR I = 0 TO KN - 1
430   RK%(I) = KN%(I)
440 NEXT
450 DIM P%(2,1) : REM  RANDOM POINTS
460 REM  == K-MEANS DATA TABLES ==
470 DIM DI(NS - 1,KN - 1)
480 REM  -- K - MU-XO, MU-X1, N-K --
490 DIM KM(KN - 1,2)
500 REM  -- K - OLD MU-X0, OLD MU-X1 --
510 DIM KO(KN - 1,1)
520 REM == GENERATE SAMPLES ==
530 FOR S = 0 TO NS - 1
540   REM -- ADAPTIVE QUOTA SAMPLING --
550   REM -- PICK CLASS BASED ON REMAINING SAMPLES TO GENERATE --
560   PC% = RND(1) * (NS - S)
570   CS = 0
580   J = 0
590   CS = CS + RK%(J)
600   K = J
610   J = J + 1
620   IF PC% >= CS GOTO 590
630   RK%(K) = RK%(K) - 1
640   REM -- GENERATE RANDOM NORMAL VARIABLES Z0, Z1 --
650   GOSUB 1400
660   REM -- GENERATE SAMPLE FROM Z0, Z1 --
670   REM X0 = MU-X0 + SIG-X0 * Z0
680   X0% = GE(K,0) + GE(K,2) * Z0
690   REM X1 = MU-X1 + RHO * SIG-X1 * Z0 + SQR(1 - RHO^2) * SIG-X1 * Z1
700   X1% = GE(K,1) + GE(K,4) * GE(K,3) * Z0 + GE(K,5) * GE(K,3) * Z1
710   DS%(S,0) = X0%
720   DS%(S,1) = X1%
730   DS%(S,2) = K
740   PRINT S + 1; TAB(5); "SAMPLE "; KN%(K) - RK%(K);
750   PRINT TAB(16); "IN CLASS "; K;
760   PRINT " AT ("; X0%; ","; X1%; ")"
770   ON K + 1 GOSUB 1200,1300: REM PLOT SAMPLE
780 NEXT
790 REM == COUNT SAMPLES OF EACH CLASS TO VERIFY RESULT ==
800 DIM CT%(KN - 1)
810 FOR I = 0 TO NS - 1
820   K = DS%(I,2)
830   CT%(K) = CT%(K) + 1
840 NEXT
850 PRINT "TOTAL SAMPLES FOR EACH CLASS:"
860 FOR I = 0 TO KN - 1
870   PRINT CT%(I); TAB(5); "SAMPLES IN CLASS "; I
880 NEXT
890 RETURN
900 REM  == WAIT FOR KEYSTROKE ==
905 PRINT "TYPE ANY KEY TO CONTINUE..."
910 POKE 49168,0 : REM  CLEAR BUFFER
920 IF PEEK(49152) < 128 GOTO 920
930 POKE 49168,0
940 RETURN
1000 REM  == SET GRAPHICS AND DRAW Y AXIS WITH TICK MARKS ==
1010 HGR : HCOLOR=7
1020 HPLOT 0,159 TO 279,159
1030 HPLOT 0,159 TO 0,0
1040 J = 0
1050 FOR I = 159 TO 0 STEP -10
1060   HPLOT 0,I TO 1,I
1070   IF J / 5 = INT(J / 5) THEN HPLOT 0,I TO 2,I
1080   IF J / 10 = INT(J / 10) THEN HPLOT 0,I TO 3,I
1090   J = J + 1
1100 NEXT I
1110 J = 0
1120 FOR I = 0 TO 279 STEP 10
1130   HPLOT I,159 TO I,158
1140   IF J / 5 = INT(J / 5) THEN HPLOT I,159 TO I,157
1150   IF J / 10 = INT(J / 10) THEN HPLOT I,159 TO I,156
1160   J = J + 1
1170 NEXT I
1180 RETURN
1200 REM == DRAW + AT X0, X1 ==
1210 X1% = 159 - X1%
1220 IF X0% < 1 OR X0% > 270 THEN RETURN
1230 IF X1% < 1 OR X1% > 150 THEN RETURN
1240 HPLOT X0% - 1,X1% TO X0% + 1,X1%
1250 HPLOT X0%,X1% - 1 TO X0%,X1% + 1
1260 RETURN
1300 REM == DRAW BOX AT X0, X1 ==
1310 X1% = 159 - X1%
1320 IF X0% < 1 OR X0% > 270 THEN RETURN
1330 IF X1% < 1 OR X1% > 150 THEN RETURN
1340 HPLOT X0% - 1,X1% - 1 TO X0% + 1,X1% - 1
1350 HPLOT X0% - 1,X1%: HPLOT X0% + 1,X1%
1360 HPLOT X0% - 1,X1% + 1 TO X0% + 1,X1% + 1
1370 RETURN
1400 REM == IRWIN-HALL DISTRIBUTION ==
1490 Z0 = -6:Z1 = -6
1500 FOR IH = 1 TO 12
1510   Z0 = Z0 + RND(1)
1520   Z1 = Z1 + RND(1)
1530 NEXT
1540 RETURN
1550 REM == RANDOM GEN SPEED TEST ==
1560 FOR RG = 1 TO 500
1570   GOSUB 1400
1580 NEXT
2000 REM  == K-MEANS ==
2010 PRINT "RUN K-MEANS"
2020 REM  -- CLEAR PREDICTIONS --
2030 FOR I = 0 TO NS - 1
2040   DS%(I,3) = 0
2050 NEXT I
2100 REM  -- INITIALIZE CENTROIDS --
2110 FOR I = 0 TO KN - 1
2120   J = INT(RND(1) * NS)
2130   IF DS%(J,3) = 1 GOTO 2120
2140   KM(I,1) = DS%(J,1)
2150   KM(I,2) = DS%(J,2)
2160   DS%(J,3) = 1
2170 NEXT I
2200 REM  -- DRAW LINES BETWEEN CENTROIDS --
2210 FOR I = 1 TO KN - 1
2220   HPLOT KM(I-1,0), 159-KM(I-1,1) TO KM(I,0), 159-KM(I,1)
2230 NEXT I
2240 GOSUB 3000: REM  DRAW DECISION BOUNDARY
2300 REM  -- COMPUTE ASSIGNMENTS --
2310 FOR I = 0 TO NS - 1
2320   PRINT "POINT ";I;" AT ";DS%(I,0);",";DS%(I,1);
2330   DS%(I,3) = 0
2340   FOR J = 0 TO KN - 1
2350     DI(I,J) = (DS%(I,0)-KM(J,0))^2 + (DS%(I,1)-KM(J,1))^2
2360     IF J >0 AND (DI(I,J) < DI(I,DS%(I,3))) THEN DS%(I,3) = J
2370   NEXT J
2380   PRINT " -> ";DS%(I,3);" Y^=";DS%(I,2)
2390 NEXT I
2500 REM  -- COMPUTE ACCURACY --
2510 CT = 0
2520 FOR I = 0 TO NS - 1
2530   IF DS%(I,2) = DS%(I,3) THEN CT = CT + 1
2540 NEXT I
2550 A = CT / NS
2560 IF A < 0.5 THEN A = 1 - A
2570 PRINT "ACCURACY = "; INT(A*10000+0.5)/100;"%"
2600 REM  -- COMPUTE CENTROIDS --
2610 FOR J = 0 TO KN - 1
2620   KO(J,0) = KM(J,0)
2630   KO(J,1) = KM(J,1)
2640   KM(J,0) = 0: KM(J,1) = 0
2650   KM(J,2) = 0
2660 NEXT
2670 FOR I = 0 TO NS - 1
2680   Y = DS%(I,3)
2690   KM(Y,0) = KM(Y,0) + DS%(I,0)
2700   KM(Y,1) = KM(Y,1) + DS%(I,1)
2710   KM(Y,2) = KM(Y,2) + 1
2720 NEXT
2730 FOR I = 0 TO KN - 1
2740   KM(I,0) = KM(I,0) / KM(I,2)
2750   KM(I,1) = KM(I,1) / KM(I,2)
2760 NEXT
2800 REM  -- DETERMINE CONVERGENCE --
2810 DI = 0
2820 FOR I = 0 TO KN - 1
2830   DI = DI + (KM(I,0) - KO(I,0)) ^ 2 + (KM(I,1) - KO(I,1)) ^ 2
2840 NEXT
2850 IF DI > 0.01 THEN GOTO 2200
2860 PRINT "K-MEANS CONVERGED"
2900 REM  -- CLEAR GRAPHICS AND REDRAW WITH DECISION BOUNDARY --
2910 GOSUB 1000
2920 FOR IR = 0 TO NS - 1
2930   X0% = DS%(IR,0)
2940   X1% = DS%(IR,1)
2950   K = DS%(IR,2)
2960   ON K + 1 GOSUB 1200,1300
2970 NEXT
2980 GOSUB 3000
2990 RETURN
3000 REM  -- DRAW DECISION BOUNDARY --
3010 FOR I = 1 TO KN - 1
3020   M = 1E6
3030   IF KM(I - 1,1) - KM(I,1) <> 0 THEN M = -1 * (KM(I - 1,0) - KM(I,0)) / (KM(I - 1,1) - KM(I,1))
3040   P%(0,0) = (KM(I,0) - KM(I - 1,0)) / 2 + KM(I - 1,0)
3050   P%(0,1) = (KM(I,1) - KM(I - 1,1)) / 2 + KM(I - 1,1)
3060   REM GOSUB 3500 -- WHAT DOES THIS DO?
3070 NEXT
3080 REM  -- DRAW LINE FROM SLOPE AND POINT --
3090 NX = 1 : REM  -- REM NUMBER OF INTERSECTIONS --
3100 IF ABS(M) > 1E5 THEN GOSUB 3240 : GOTO 3210 : REM  VERTICAL LINE
3110 P%(NX,1) = M * (10 - P%(0,0)) + P%(0,1)
3120 IF P%(NX,1) > 10 AND P%(NX,1) < 149 THEN P%(NX,0) = 10 : NX = NX + 1
3130 P%(NX,1) = M * (269 - P%(0,0)) + P%(0,1)
3140 IF P%(NX,1) > 10 AND P%(NX,1) < 149 THEN P%(NX,0) = 269 : NX = NX + 1
3150 IF NX = 3 THEN GOTO 3210
3160 IF M <> 0 THEN P%(NX,0) = (10 - P%(0,1)) / M + P%(0,0)
3170 IF M <> 0 AND P%(NX,0) > 10 AND P%(NX,0) < 269 THEN P%(NX,1) = 10 : NX = NX + 1
3180 IF NX = 3 THEN GOTO 3210
3190 IF M <> 0 THEN P%(NX,0) = (149 - P%(0,1)) / M + P%(0,0)
3200 IF M <> 0 AND P%(NX,0) > 10 AND P%(NX,0) < 269 THEN P%(NX,1) = 149 : NX = NX + 1
3210 REM  -- DRAW LINE --
3220 IF NX = 3 THEN HPLOT P%(1,0),159 - P%(1,1) TO P%(2,0),159 - P%(2,1)
3230 RETURN
3240 REM  -- VERTICAL LINE --
3250 P%(1,0) = P%(0,0)
3260 P%(2,0) = P%(0,0)
3270 P%(1,1) = 10
3280 P%(2,1) = 269
3290 RETURN

[mdcramer]

@thecloudexpanse , thank you again. I think I'm seeing an issue with not running the code pulled from the chatGPT scanned photos. There is no line 3500, as you saw, but it's running well on my machine, so I'll boot it up tomorrow and check out that line. I'll report back... thank you for doing this!

[thecloudexpanse]

@mdcramer Just updated the source, changing all KM% array references to KM. That seems to be consistent and allows convergence to occur with >50% accuracy.