Decision tree

"""A decision tree classifier (CART with Gini impurity), standard library only.

Grows a tree on the training days by always making the split that leaves the
two sides purest, then predicts rain for the test days. Prints the tree in
preorder: each node is either a split (feature and threshold) or a leaf
(its prediction).
"""

import csv
import sys

FEATURES = ("humidity", "pressure", "wind")
MAX_DEPTH = 3
MIN_SIZE = 5


def read(path):
    train, test = [], []
    with open(path, newline="") as f:
        for row in csv.DictReader(f):
            day = ([float(row[name]) for name in FEATURES], row["rain"])
            (train if row["split"] == "train" else test).append(day)
    return train, test


def gini(labels):
    # Chance that two days drawn at random (with replacement) disagree.
    p = labels.count("yes") / len(labels)
    return 1 - p * p - (1 - p) * (1 - p)


def majority(labels):
    # An even split predicts "no".
    return "yes" if labels.count("yes") > len(labels) / 2 else "no"


def best_split(days):
    """The (score, feature, threshold) with the lowest weighted Gini impurity.
    Features are tried in column order and thresholds from smallest to largest;
    only a strictly lower score replaces the best so far."""
    best = None
    for f in range(len(FEATURES)):
        values = sorted({x[f] for x, _ in days})
        for a, b in zip(values, values[1:]):
            threshold = (a + b) / 2
            left = [label for x, label in days if x[f] <= threshold]
            right = [label for x, label in days if x[f] > threshold]
            score = (len(left) * gini(left) + len(right) * gini(right)) / len(days)
            if best is None or score < best[0]:
                best = (score, f, threshold)
    return best


def grow(days, depth=0):
    """A leaf is a label; a split is (feature, threshold, left, right)."""
    labels = [label for _, label in days]
    if gini(labels) == 0 or len(days) < MIN_SIZE or depth == MAX_DEPTH:
        return majority(labels)
    split = best_split(days)
    if split is None or split[0] >= gini(labels):
        return majority(labels)  # no split makes the sides purer
    _, f, threshold = split
    left = [d for d in days if d[0][f] <= threshold]
    right = [d for d in days if d[0][f] > threshold]
    return (f, threshold, grow(left, depth + 1), grow(right, depth + 1))


def predict(node, x):
    while isinstance(node, tuple):
        f, threshold, left, right = node
        node = left if x[f] <= threshold else right
    return node


def preorder(node):
    yield node
    if isinstance(node, tuple):
        yield from preorder(node[2])
        yield from preorder(node[3])


def main():
    train, test = read(sys.argv[1])
    tree = grow(train)

    print(f"train_count {len(train)}")
    print(f"test_count {len(test)}")
    leaves = 0
    for k, node in enumerate(preorder(tree), start=1):
        if isinstance(node, tuple):
            print(f"node_{k}_feature {FEATURES[node[0]]}")
            print(f"node_{k}_threshold {node[1]:.6f}")
        else:
            print(f"node_{k}_predict {node}")
            leaves += 1
    print(f"leaves {leaves}")

    train_correct = sum(predict(tree, x) == label for x, label in train)
    correct = sum(predict(tree, x) == label for x, label in test)
    print(f"train_accuracy {train_correct / len(train):.6f}")
    print(f"correct {correct}")
    print(f"accuracy {correct / len(test):.6f}")


if __name__ == "__main__":
    main()

--  A decision tree classifier (CART with Gini impurity), using only Ada's
--  standard library.
--
--  Grows a tree on the training days by always making the split that leaves
--  the two sides purest, then predicts rain for the test days. Prints the
--  tree in preorder: each node is either a split (feature and threshold) or
--  a leaf (its prediction). Ada has no decision tree or classifier to skip.

with Ada.Text_IO;                use Ada.Text_IO;
with Ada.Long_Float_Text_IO;     use Ada.Long_Float_Text_IO;
with Ada.Command_Line;
with Ada.Characters.Handling;    use Ada.Characters.Handling;
with Ada.Strings.Fixed;          use Ada.Strings.Fixed;
with Ada.Containers.Vectors;
with Ada.Containers.Generic_Array_Sort;

procedure Main is
   type Feature is (Humidity, Pressure, Wind);   --  also the column names
   type Features is array (Feature) of Long_Float;

   type Day is record
      X    : Features;
      Rain : Boolean;
   end record;
   type Day_Array is array (Positive range <>) of Day;
   package Day_Vectors is new Ada.Containers.Vectors (Positive, Day);

   function To_Array (V : Day_Vectors.Vector) return Day_Array is
     [for I in 1 .. Natural (V.Length) => V (I)];

   type Float_Array is array (Positive range <>) of Long_Float;
   procedure Sort is new Ada.Containers.Generic_Array_Sort
     (Positive, Long_Float, Float_Array);

   Max_Depth : constant := 3;
   Min_Size  : constant := 5;

   --  A variant record: the discriminant Is_Leaf, fixed when a node is
   --  made, decides which fields the node has. A Node_Access points to a
   --  node, so a split can hold its children. The first line announces
   --  Node (an incomplete declaration) so that Node_Access can name it.
   type Node (Is_Leaf : Boolean);
   type Node_Access is access Node;
   type Node (Is_Leaf : Boolean) is record
      case Is_Leaf is
         when True =>
            Rain : Boolean;
         when False =>
            On          : Feature;
            Threshold   : Long_Float;
            Left, Right : Node_Access;
      end case;
   end record;

   --  The text of field number Col (counting from 1) in a CSV line.
   function Field (Line : String; Col : Positive) return String is
      Start : Positive := Line'First;
   begin
      for I in 1 .. Col - 1 loop
         Start := Index (Line (Start .. Line'Last), ",") + 1;
      end loop;
      declare
         Stop : constant Natural := Index (Line (Start .. Line'Last), ",");
      begin
         return (if Stop = 0 then Line (Start .. Line'Last)
                 else Line (Start .. Stop - 1));
      end;
   end Field;

   function Column (Header, Name : String) return Positive is
   begin
      for C in 1 .. Ada.Strings.Fixed.Count (Header, ",") + 1 loop
         if Field (Header, C) = Name then
            return C;
         end if;
      end loop;
      raise Constraint_Error with "column not found: " & Name;
   end Column;

   procedure Read (Path : String; Train, Test : out Day_Vectors.Vector) is
      File : File_Type;
   begin
      Open (File, In_File, Path);
      declare
         Header : constant String := Get_Line (File);
         --  F'Image is an enumeration value's name in capitals: "WIND".
         Cols : constant array (Feature) of Positive :=
           [for F in Feature => Column (Header, To_Lower (F'Image))];
         Rain_Col  : constant Positive := Column (Header, "rain");
         Split_Col : constant Positive := Column (Header, "split");
      begin
         while not End_Of_File (File) loop
            declare
               Line : constant String := Get_Line (File);
               D    : Day;
            begin
               for F in Feature loop
                  D.X (F) := Long_Float'Value (Field (Line, Cols (F)));
               end loop;
               D.Rain := Field (Line, Rain_Col) = "yes";
               if Field (Line, Split_Col) = "train" then
                  Train.Append (D);
               else
                  Test.Append (D);
               end if;
            end;
         end loop;
      end;
      Close (File);
   end Read;

   function Yes_Count (Days : Day_Array) return Natural is
      N : Natural := 0;
   begin
      for D of Days loop
         if D.Rain then
            N := N + 1;
         end if;
      end loop;
      return N;
   end Yes_Count;

   --  Chance that two days drawn at random (with replacement) disagree.
   function Gini (Days : Day_Array) return Long_Float is
      P : constant Long_Float :=
        Long_Float (Yes_Count (Days)) / Long_Float (Days'Length);
   begin
      return 1.0 - P * P - (1.0 - P) * (1.0 - P);
   end Gini;

   --  An even split predicts "no".
   function Majority (Days : Day_Array) return Boolean is
     (2 * Yes_Count (Days) > Days'Length);

   --  The days on one side of a threshold: Left means X (F) <= T.
   function Side (Days : Day_Array; F : Feature; T : Long_Float;
                  Left : Boolean) return Day_Array is
      Kept : Day_Vectors.Vector;
   begin
      for D of Days loop
         if (D.X (F) <= T) = Left then
            Kept.Append (D);
         end if;
      end loop;
      return To_Array (Kept);
   end Side;

   --  Features are tried in order and thresholds from smallest to largest;
   --  only a strictly lower weighted Gini replaces the best so far.
   procedure Best_Split (Days : Day_Array; Found : out Boolean;
                         Score : out Long_Float; On : out Feature;
                         Threshold : out Long_Float) is
   begin
      Found := False;
      Score := Long_Float'Last;
      On := Feature'First;
      Threshold := 0.0;
      for F in Feature loop
         declare
            Values : Float_Array := [for I in Days'Range => Days (I).X (F)];
         begin
            Sort (Values);
            for I in Values'First .. Values'Last - 1 loop
               if Values (I) /= Values (I + 1) then   --  skip repeats
                  declare
                     T : constant Long_Float :=
                       (Values (I) + Values (I + 1)) / 2.0;
                     L : constant Day_Array := Side (Days, F, T, True);
                     R : constant Day_Array := Side (Days, F, T, False);
                     S : constant Long_Float :=
                       (Long_Float (L'Length) * Gini (L)
                        + Long_Float (R'Length) * Gini (R))
                       / Long_Float (Days'Length);
                  begin
                     if S < Score then
                        Found := True;
                        Score := S;
                        On := F;
                        Threshold := T;
                     end if;
                  end;
               end if;
            end loop;
         end;
      end loop;
   end Best_Split;

   function Grow (Days : Day_Array; Depth : Natural) return Node_Access is
      Found     : Boolean;
      Score, T  : Long_Float;
      F         : Feature;
      Impurity  : constant Long_Float := Gini (Days);
      function Leaf return Node_Access is
        (new Node'(Is_Leaf => True, Rain => Majority (Days)));
   begin
      if Impurity = 0.0 or Days'Length < Min_Size or Depth = Max_Depth then
         return Leaf;
      end if;
      Best_Split (Days, Found, Score, F, T);
      if not Found or Score >= Impurity then
         return Leaf;   --  no split makes the sides purer
      end if;
      return new Node'(Is_Leaf   => False, On => F, Threshold => T,
                       Left  => Grow (Side (Days, F, T, True), Depth + 1),
                       Right => Grow (Side (Days, F, T, False), Depth + 1));
   end Grow;

   function Predict (Tree : Node_Access; X : Features) return Boolean is
      N : Node_Access := Tree;
   begin
      while not N.Is_Leaf loop
         N := (if X (N.On) <= N.Threshold then N.Left else N.Right);
      end loop;
      return N.Rain;
   end Predict;

   procedure Say (Name : String; Value : Long_Float) is
   begin
      Put (Name & " ");
      Put (Value, Fore => 1, Aft => 6, Exp => 0);   --  six decimals, no E
      New_Line;
   end Say;

   K, Leaves : Natural := 0;

   --  Number each node in preorder and print it.
   procedure Print (N : Node_Access) is
   begin
      K := K + 1;
      declare
         Name : constant String := "node_" & Trim (K'Image, Ada.Strings.Left);
      begin
         if N.Is_Leaf then
            Put_Line (Name & "_predict " & (if N.Rain then "yes" else "no"));
            Leaves := Leaves + 1;
         else
            Put_Line (Name & "_feature " & To_Lower (N.On'Image));
            Say (Name & "_threshold", N.Threshold);
            Print (N.Left);
            Print (N.Right);
         end if;
      end;
   end Print;

   function Correct (Tree : Node_Access; Days : Day_Array) return Natural is
      N : Natural := 0;
   begin
      for D of Days loop
         if Predict (Tree, D.X) = D.Rain then
            N := N + 1;
         end if;
      end loop;
      return N;
   end Correct;

   Train_V, Test_V : Day_Vectors.Vector;
begin
   Read (Ada.Command_Line.Argument (1), Train_V, Test_V);
   declare
      Train : constant Day_Array := To_Array (Train_V);
      Test  : constant Day_Array := To_Array (Test_V);
      Tree  : constant Node_Access := Grow (Train, 0);
      Hits  : constant Natural := Correct (Tree, Test);
   begin
      Put_Line ("train_count" & Train'Length'Image);
      Put_Line ("test_count" & Test'Length'Image);
      Print (Tree);
      Put_Line ("leaves" & Leaves'Image);
      Say ("train_accuracy",
           Long_Float (Correct (Tree, Train)) / Long_Float (Train'Length));
      Put_Line ("correct" & Hits'Image);
      Say ("accuracy", Long_Float (Hits) / Long_Float (Test'Length));
   end;
end Main;

⍝ A decision tree classifier (CART with Gini impurity), from scratch in GNU
⍝ APL. Grows a tree on the training days by always making the split that
⍝ leaves the two sides purest, then predicts rain for the test days.
⍝ APL reads right to left, with no operator precedence: 2×3+4 is 2×(3+4).
⍝ ⍝ starts a comment; ← assigns; ⍵ is a function's right argument.

∇NAME SAY N
 ⍝ Print "name value" for a whole number. ⍕ formats a number as text.
 ⎕←NAME,' ',⍕N
∇

∇NAME SAY6 X;T
 ⍝ Print "name value" with six decimal places.
 T←0 6⍕|X                    ⍝ |X is the size of X; 0 6⍕ gives 6 decimals
 T←(∨\T≠' ')/T               ⍝ drop the leading blank ⍕ adds
 T←((T[1]='.')/'0'),T        ⍝ ⍕ writes .94 for 0.94, so put the 0 back
 ⎕←NAME,' ',((X<0)/'-'),T    ⍝ sign as -, since APL's own minus is ¯
∇

FEATURES←'humidity' 'pressure' 'wind'
MAX_DEPTH←3 ⋄ MIN_SIZE←5

⍝ Read the file. ⎕ARG is the command line; ↑¯1↑ takes its last item, the
⍝ data path. ⎕FIO[26] reads the whole file as one string.
TEXT←⎕FIO[26]↑¯1↑⎕ARG
LINES←(TEXT≠⎕UCS 10)⊂TEXT   ⍝ ⊂ splits at the 0s of the mask: one item per line
⍝ Split each (¨) line at commas; ⊃ stacks them into a table.
ROWS←⊃{(⍵≠',')⊂⍵}¨LINES
⍝ A column by header name (⍳), minus the header (1↓).
COLUMN←{1↓ROWS[;ROWS[1;]⍳⊂⍵]}

X←⍉⊃{⍎¨COLUMN ⍵}¨FEATURES    ⍝ one row per day; ⍎ turns text into a number
Y←(COLUMN 'rain')≡¨⊂'yes'    ⍝ 1 for a rainy day (≡ matches whole strings)
TRAIN←(COLUMN 'split')≡¨⊂'train'
TX←TRAIN⌿X ⋄ TY←TRAIN/Y      ⍝ ⌿ keeps the training rows of a table

⍝ Chance that two days drawn at random disagree, from the share ⍵ of rain.
GINI←{1-(⍵*2)+(1-⍵)*2}

∇R←Y SPLIT V;U;T;M;N;NL;YL;S
 ⍝ The best (score, threshold) for one feature's values V at a node.
 U←∪V ⋄ U←U[⍋U]              ⍝ the distinct values, sorted (⍋ is grade up)
 R←1 0 ⋄ →(1=≢U)/0           ⍝ one value: no split (score 1 never wins)
 T←((¯1↓U)+1↓U)÷2            ⍝ midpoints of neighboring values
 ⍝ ∘.≤ makes a table, one row per day and one column per threshold, with
 ⍝ a 1 where the day goes left. +⌿ counts down each column, and Y+.×
 ⍝ counts the rainy days among them, for every threshold at once.
 M←V∘.≤T ⋄ N←≢V ⋄ NL←+⌿M ⋄ YL←Y+.×M
 S←((NL×GINI YL÷NL)+(N-NL)×GINI((+/Y)-YL)÷N-NL)÷N
 R←(⌊/S),T[S⍳⌊/S]            ⍝ ⍳ finds the first minimum: smallest wins
∇

⍝ The tree is four vectors, one item per node in preorder. A split node has
⍝ a feature number (FEAT) and threshold; a leaf has FEAT 0 and a prediction.
⍝ A left child always comes next in preorder; RIGHT records the right one.
FEAT←THR←PRED←RIGHT←⍬

∇DEPTH GROW I;Y;N;G;B;F;K;L
 ⍝ Grow the node holding training rows I, then (recursively) its children.
 Y←TY[I] ⋄ N←≢I ⋄ G←GINI(+/Y)÷N
 ⍝ Add it as a leaf first; an even split predicts "no".
 ⍝ GNU APL has no FEAT,←0 shorthand, so each vector is rebuilt with , (join).
 FEAT←FEAT,0 ⋄ THR←THR,0 ⋄ RIGHT←RIGHT,0 ⋄ PRED←PRED,(+/Y)>N÷2
 K←≢FEAT ⋄ →((G=0)∨(N<MIN_SIZE)∨DEPTH=MAX_DEPTH)/0
 ⍝ ⊂[1] gives each feature's column; ⊃ stacks the results, one row each.
 B←⊃(⊂Y)SPLIT¨⊂[1]TX[I;]
 F←B[;1]⍳⌊/B[;1]             ⍝ the first feature with the lowest score
 →(B[F;1]≥G)/0               ⍝ no split makes the sides purer: stay a leaf
 FEAT[K]←F ⋄ THR[K]←B[F;2]
 L←TX[I;F]≤THR[K]
 (DEPTH+1)GROW L/I
 RIGHT[K]←1+≢FEAT            ⍝ the right subtree starts after the left
 (DEPTH+1)GROW(~L)/I
∇

0 GROW⍳≢TY

∇P←WALK DAY;K
 ⍝ Follow the splits from the root (node 1) down to a leaf.
 K←1
NEXT: →(FEAT[K]=0)/LEAF          ⍝ → jumps to a label when its argument is 1
 ⍝ The two children are K+1 and RIGHT[K]; > picks the right one (index 2).
 K←((K+1),RIGHT[K])[1+DAY[FEAT[K]]>THR[K]]
 →NEXT
LEAF: P←PRED[K]
∇

⍝ How many of the days in rows ⍵ the tree gets right.
CORRECT←{+/Y[⍵]=WALK¨⊂[2]X[⍵;]}
TRAINED←CORRECT TRAIN/⍳≢Y ⋄ TESTED←CORRECT(~TRAIN)/⍳≢Y

∇SHOW K;NODE
 NODE←'node_',(⍕K),'_'
 →(FEAT[K]=0)/LEAF
 ⎕←NODE,'feature ',⊃FEATURES[FEAT[K]]
 (NODE,'threshold')SAY6 THR[K]
 →0
LEAF: ⎕←NODE,'predict ',⊃('no' 'yes')[1+PRED[K]]
∇

'train_count' SAY +/TRAIN
'test_count' SAY +/~TRAIN
SHOW¨⍳≢FEAT
'leaves' SAY +/FEAT=0
'train_accuracy' SAY6 TRAINED÷+/TRAIN
'correct' SAY TESTED
'accuracy' SAY6 TESTED÷+/~TRAIN
)OFF

# A decision tree classifier (CART with Gini impurity), in POSIX awk.
#
# Grows a tree on the training days by always making the split that leaves the
# two sides purest, then predicts rain for the test days. Nodes are numbered in
# preorder as they are grown, so the tree prints in that order.

BEGIN {
    FS = ","
    nf = split("humidity,pressure,wind", feature, ",")
    MAX_DEPTH = 3; MIN_SIZE = 5
}

# The header line: find each column by name.
NR == 1 {
    for (i = 1; i <= NF; i++) {
        for (f = 1; f <= nf; f++)
            if ($i == feature[f]) fcol[f] = i
        if ($i == "rain") c_lab = i
        if ($i == "split") c_split = i
    }
    for (f = 1; f <= nf; f++)
        if (!fcol[f]) missing = missing " " feature[f]
    if (!c_lab) missing = missing " rain"
    if (!c_split) missing = missing " split"
    if (missing != "") {
        print "missing column(s):" missing > "/dev/stderr"
        failed = 1
        exit 1
    }
    next
}

# Every other line goes into the training or the test arrays. x[i, f] is
# feature f of training day i; awk fakes 2-D arrays with a joined key.
$c_split == "train" {
    n++
    for (f = 1; f <= nf; f++) x[n, f] = $fcol[f] + 0
    y[n] = $c_lab
    next
}
$c_split == "test" {
    m++
    for (f = 1; f <= nf; f++) q[m, f] = $fcol[f] + 0
    qy[m] = $c_lab
}

# Chance that two days drawn at random (with replacement) disagree.
function gini(yes, total,    p) {
    p = yes / total
    return 1 - p * p - (1 - p) * (1 - p)
}

# Grow the node holding the training days listed in days[1..count], and return
# its number. awk has no records, so node k lives in parallel arrays: split_f,
# split_t, left, and right for a split, and predict for a leaf. Locals are
# declared as extra parameters, so each recursive call gets its own.
function grow(days, count, depth,
              k, i, j, s, f, v, nv, a, t, nl, yl, yes, score, node_gini,
              best, best_f, best_t, ldays, rdays, nr) {
    k = ++nodes
    yes = 0
    for (i = 1; i <= count; i++) if (y[days[i]] == "yes") yes++
    # An even split predicts "no".
    predict[k] = (yes > count / 2) ? "yes" : "no"
    node_gini = gini(yes, count)
    if (node_gini == 0 || count < MIN_SIZE || depth == MAX_DEPTH) return k

    # Features in column order, thresholds from smallest to largest; only a
    # strictly lower score replaces the best, so the first of a tie wins.
    best = ""
    for (f = 1; f <= nf; f++) {
        # The distinct values, kept sorted by insertion (POSIX awk has no sort).
        nv = 0
        for (i = 1; i <= count; i++) {
            a = x[days[i], f]
            for (j = nv; j > 0 && v[j] > a; j--) ;  # find a's place
            if (j > 0 && v[j] == a) continue         # already listed
            for (s = nv; s > j; s--) v[s + 1] = v[s]
            v[j + 1] = a; nv++
        }
        for (j = 1; j < nv; j++) {
            t = (v[j] + v[j + 1]) / 2
            nl = yl = 0
            for (i = 1; i <= count; i++)
                if (x[days[i], f] <= t) { nl++; if (y[days[i]] == "yes") yl++ }
            score = (nl * gini(yl, nl) + (count - nl) * gini(yes - yl, count - nl)) / count
            if (best == "" || score < best) { best = score; best_f = f; best_t = t }
        }
    }
    if (best == "" || best >= node_gini) return k  # no split makes the sides purer

    nl = nr = 0
    for (i = 1; i <= count; i++)
        if (x[days[i], best_f] <= best_t) ldays[++nl] = days[i]
        else rdays[++nr] = days[i]
    split_f[k] = best_f; split_t[k] = best_t
    left[k] = grow(ldays, nl, depth + 1)
    right[k] = grow(rdays, nr, depth + 1)
    return k
}

# Walk from the root to a leaf; z[f] holds the day's features.
function classify(z,    k) {
    k = 1
    while (k in split_f) k = (z[split_f[k]] <= split_t[k]) ? left[k] : right[k]
    return predict[k]
}

END {
    if (failed) exit 1

    for (i = 1; i <= n; i++) all[i] = i
    grow(all, n, 0)

    printf "train_count %d\n", n
    printf "test_count %d\n", m
    for (k = 1; k <= nodes; k++) {
        if (k in split_f) {
            printf "node_%d_feature %s\n", k, feature[split_f[k]]
            printf "node_%d_threshold %.6f\n", k, split_t[k]
        } else {
            printf "node_%d_predict %s\n", k, predict[k]
            leaves++
        }
    }
    printf "leaves %d\n", leaves

    for (i = 1; i <= n; i++) {
        for (f = 1; f <= nf; f++) z[f] = x[i, f]
        if (classify(z) == y[i]) train_correct++
    }
    for (i = 1; i <= m; i++) {
        for (f = 1; f <= nf; f++) z[f] = q[i, f]
        if (classify(z) == qy[i]) correct++
    }
    printf "train_accuracy %.6f\n", train_correct / n
    printf "correct %d\n", correct
    printf "accuracy %.6f\n", correct / m
}

10 REM A decision tree classifier (CART with Gini impurity), in classic BASIC.
20 REM
30 REM Grows a tree on the training days by always making the split that
40 REM leaves the two sides purest, then predicts rain for the test days.
50 REM GOSUB can call itself, but every variable is global, so a recursive
60 REM grow would lose its place. Instead, nodes wait on a stack: popping
70 REM the left child before the right visits them in preorder, and each
80 REM node is printed as soon as it is decided.
90 REM
100 REM run.sh passes the data path in the DATA environment variable,
110 REM because bwbasic can't read command-line arguments.
120 MAXDEPTH = 3 : MINSIZE = 5
130 DIM NAME$(3), FL$(20)
140 NAME$(1) = "humidity" : NAME$(2) = "pressure" : NAME$(3) = "wind"
150 REM Each day has three features, a label, and DN, the tree node it sits
160 REM in while the tree grows (0 for a test day, so it is never counted).
170 DIM X(1000, 3), RAIN$(1000), ISTRAIN(1000), DN(1000)
180 REM Depth 3 allows at most 15 nodes. A split node has a feature (0 for
190 REM a leaf), a threshold, and two children; a leaf has a prediction.
200 DIM DEPTH(15), FEAT(15), THR(15), LFT(15), RGT(15), PRED$(15), STACK(15)
210 GOSUB 1000 : REM read the days
220 PRINT "train_count"; NTRAIN
230 PRINT "test_count"; NTEST
240 REM Every training day starts in the root, node 1, at depth 0.
250 NODES = 1 : DEPTH(1) = 0
260 FOR I = 1 TO NDAYS
270 IF ISTRAIN(I) THEN DN(I) = 1
280 NEXT I
290 SP = 1 : STACK(1) = 1 : PREORDER = 0 : LEAVES = 0
300 WHILE SP > 0
310 K = STACK(SP) : SP = SP - 1
320 GOSUB 2000 : REM decide node K
330 WEND
340 PRINT "leaves"; LEAVES
350 TRAINOK = 0 : TESTOK = 0
360 FOR I = 1 TO NDAYS
370 GOSUB 8000 : REM predict day I into P$
380 REM A second IF after THEN on the same line is unreliable in bwbasic,
390 REM so both tests go in one IF, joined by AND.
400 IF P$ = RAIN$(I) AND ISTRAIN(I) = 1 THEN TRAINOK = TRAINOK + 1
410 IF P$ = RAIN$(I) AND ISTRAIN(I) = 0 THEN TESTOK = TESTOK + 1
420 NEXT I
430 L$ = "train_accuracy" : Y = TRAINOK / NTRAIN : GOSUB 3000
440 PRINT "correct"; TESTOK
450 L$ = "accuracy" : Y = TESTOK / NTEST : GOSUB 3000
460 SYSTEM : REM without SYSTEM, bwbasic waits at its prompt
1000 REM Read every day into X, RAIN$, and ISTRAIN, finding the columns
1010 REM by their names in the header.
1020 F$ = ENVIRON$("DATA")
1030 OPEN "I", #1, F$
1040 LINE INPUT #1, L$
1050 GOSUB 4000
1060 DIM COL(3)
1070 FOR C = 1 TO NF
1080 FOR F = 1 TO 3
1090 IF FL$(C) = NAME$(F) THEN COL(F) = C
1100 NEXT F
1110 IF FL$(C) = "rain" THEN RAINCOL = C
1120 IF FL$(C) = "split" THEN SPLITCOL = C
1130 NEXT C
1140 NDAYS = 0 : NTRAIN = 0
1150 IF EOF(1) THEN 1250
1160 LINE INPUT #1, L$
1170 GOSUB 4000
1180 NDAYS = NDAYS + 1
1190 FOR F = 1 TO 3
1200 X(NDAYS, F) = VAL(FL$(COL(F)))
1210 NEXT F
1220 RAIN$(NDAYS) = FL$(RAINCOL)
1230 IF FL$(SPLITCOL) = "train" THEN ISTRAIN(NDAYS) = 1 : NTRAIN = NTRAIN + 1
1240 GOTO 1150
1250 CLOSE #1 : NTEST = NDAYS - NTRAIN
1260 RETURN
2000 REM Decide node K: make it a leaf, or split it and push its two
2010 REM children, right first so the left one is popped next.
2020 PREORDER = PREORDER + 1 : ID$ = "node_" + MID$(STR$(PREORDER), 2)
2030 N = 0 : YES = 0
2040 FOR I = 1 TO NDAYS
2050 IF DN(I) = K THEN N = N + 1
2060 IF DN(I) = K AND RAIN$(I) = "yes" THEN YES = YES + 1
2070 NEXT I
2080 P = YES / N : G = 1 - P * P - (1 - P) * (1 - P)
2090 BEST = 9 : REM higher than any Gini score, so "no split found"
2100 IF G > 0 AND N >= MINSIZE AND DEPTH(K) < MAXDEPTH THEN GOSUB 5000
2110 REM No split was tried, or none makes the sides purer: a leaf.
2120 IF BEST < G THEN 2180
2130 LEAVES = LEAVES + 1
2140 REM An even split predicts "no".
2150 IF YES * 2 > N THEN PRED$(K) = "yes" ELSE PRED$(K) = "no"
2160 PRINT ID$; "_predict "; PRED$(K)
2170 RETURN
2180 FEAT(K) = BESTF : THR(K) = BESTT
2190 PRINT ID$; "_feature "; NAME$(BESTF)
2200 L$ = ID$ + "_threshold" : Y = BESTT : GOSUB 3000
2210 LFT(K) = NODES + 1 : RGT(K) = NODES + 2 : NODES = NODES + 2
2220 DEPTH(LFT(K)) = DEPTH(K) + 1 : DEPTH(RGT(K)) = DEPTH(K) + 1
2230 FOR I = 1 TO NDAYS
2240 IF DN(I) <> K THEN 2260 : REM skip days in other nodes
2250 IF X(I, BESTF) <= BESTT THEN DN(I) = LFT(K) ELSE DN(I) = RGT(K)
2260 NEXT I
2270 SP = SP + 1 : STACK(SP) = RGT(K)
2280 SP = SP + 1 : STACK(SP) = LFT(K)
2290 RETURN
3000 REM Print L$ and Y with six decimals.
3010 PRINT L$; " ";
3020 PRINT USING "#.######"; Y
3030 RETURN
4000 REM Split L$ at its commas into FL$(1) to FL$(NF).
4010 NF = 0 : P = 1
4020 Q = INSTR(P, L$, ",")
4030 NF = NF + 1
4040 IF Q = 0 THEN FL$(NF) = MID$(L$, P) : RETURN
4050 FL$(NF) = MID$(L$, P, Q - P) : P = Q + 1
4060 GOTO 4020
5000 REM The lowest weighted Gini impurity over every feature and threshold
5010 REM for node K, into BEST, BESTF, and BESTT. Features go in column
5020 REM order, thresholds from smallest to largest, and only a strictly
5030 REM lower score replaces the best so far.
5040 FOR F = 1 TO 3
5050 REM Start at the node's smallest value, then step through the
5060 REM distinct values in order without sorting them.
5070 LOW = 1E+30
5080 FOR I = 1 TO NDAYS
5090 IF DN(I) = K AND X(I, F) < LOW THEN LOW = X(I, F)
5100 NEXT I
5110 GOSUB 6000 : REM the next value above LOW, into NXT
5120 WHILE HASNEXT
5130 T = (LOW + NXT) / 2
5140 GOSUB 7000 : REM score the split at T, into S
5150 IF S < BEST THEN BEST = S : BESTF = F : BESTT = T
5160 LOW = NXT
5170 GOSUB 6000
5180 WEND
5190 NEXT F
5200 RETURN
6000 REM The smallest value of feature F in node K above LOW, into NXT,
6010 REM with HASNEXT = 0 if there is none.
6020 HASNEXT = 0
6030 FOR I = 1 TO NDAYS
6040 IF DN(I) <> K THEN 6070
6050 IF X(I, F) <= LOW THEN 6070
6060 IF HASNEXT = 0 OR X(I, F) < NXT THEN NXT = X(I, F) : HASNEXT = 1
6070 NEXT I
6080 RETURN
7000 REM (nL * giniL + nR * giniR) / n for feature F at threshold T, into
7010 REM S. Left is value <= T. Both sides hold at least one day, since T
7020 REM lies between two values present in the node.
7030 NL = 0 : YL = 0
7040 FOR I = 1 TO NDAYS
7050 IF DN(I) <> K OR X(I, F) > T THEN 7080
7060 NL = NL + 1
7070 IF RAIN$(I) = "yes" THEN YL = YL + 1
7080 NEXT I
7090 P = YL / NL : GL = 1 - P * P - (1 - P) * (1 - P)
7100 P = (YES - YL) / (N - NL) : GR = 1 - P * P - (1 - P) * (1 - P)
7110 S = (NL * GL + (N - NL) * GR) / N
7120 RETURN
8000 REM Walk day I down from the root to a leaf, into P$.
8010 K = 1
8020 WHILE FEAT(K) > 0
8030 REM X(I, FEAT(K)) in one step reads the wrong element in bwbasic, so
8040 REM copy FEAT(K) into F first.
8050 F = FEAT(K)
8060 IF X(I, F) <= THR(K) THEN K = LFT(K) ELSE K = RGT(K)
8070 WEND
8080 P$ = PRED$(K)
8090 RETURN

/* A decision tree classifier (CART with Gini impurity), using only the C
 * standard library.
 *
 * Grows a tree on the training days by always making the split that leaves
 * the two sides purest, then predicts rain for the test days. Prints the tree
 * in preorder: each node is either a split (feature and threshold) or a leaf
 * (its prediction). */

#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define N_FEATURES 3
#define MAX_DEPTH 3
#define MIN_SIZE 5
#define MAX_LINE 1024
#define MAX_FIELDS 16

static const char *FEATURES[N_FEATURES] = {"humidity", "pressure", "wind"};

struct day {
    double x[N_FEATURES];
    bool rain;
};

/* One struct serves for both kinds of node. A split has two children; a leaf
 * has NULL children and a label. The struct can point to itself because a
 * pointer's size is known before the struct is complete. */
struct node {
    int feature;
    double threshold;
    struct node *left, *right;
    bool label;
};

static bool is_leaf(const struct node *n) { return n->left == NULL; }

/* Split a line at its commas into fields[], returning how many there are.
 * strtok replaces each comma with '\0', so the fields point into line. */
static int split_fields(char *line, char *fields[]) {
    int count = 0;
    for (char *f = strtok(line, ",\r\n"); f && count < MAX_FIELDS; f = strtok(NULL, ",\r\n"))
        fields[count++] = f;
    return count;
}

/* Return the 0-based position of name among the header fields, or exit. */
static int column(char *header[], int count, const char *name) {
    for (int i = 0; i < count; i++)
        if (strcmp(header[i], name) == 0) return i;
    fprintf(stderr, "no column %s\n", name);
    exit(1);
}

/* Append one day to a growable array, doubling its capacity when full. */
static void push(struct day **days, size_t *n, size_t *capacity, struct day d) {
    if (*n == *capacity) {
        *capacity = *capacity ? *capacity * 2 : 64;
        *days = realloc(*days, *capacity * sizeof **days);
    }
    (*days)[(*n)++] = d;
}

static void read_days(const char *path, struct day **train, size_t *n_train,
                      struct day **test, size_t *n_test) {
    FILE *file = fopen(path, "r");
    if (!file) { perror(path); exit(1); }

    char header_line[MAX_LINE], line[MAX_LINE];
    char *fields[MAX_FIELDS];
    if (!fgets(header_line, sizeof header_line, file)) { fprintf(stderr, "empty file\n"); exit(1); }
    int count = split_fields(header_line, fields);
    int feature_column[N_FEATURES];
    for (int f = 0; f < N_FEATURES; f++) feature_column[f] = column(fields, count, FEATURES[f]);
    int rain_column = column(fields, count, "rain");
    int split_column = column(fields, count, "split");

    size_t train_capacity = 0, test_capacity = 0;
    *train = *test = NULL;
    *n_train = *n_test = 0;
    while (fgets(line, sizeof line, file)) {
        if (split_fields(line, fields) < count) continue;
        struct day d;
        for (int f = 0; f < N_FEATURES; f++) d.x[f] = strtod(fields[feature_column[f]], NULL);
        d.rain = strcmp(fields[rain_column], "yes") == 0;
        if (strcmp(fields[split_column], "train") == 0)
            push(train, n_train, &train_capacity, d);
        else
            push(test, n_test, &test_capacity, d);
    }
    fclose(file);
}

static size_t count_yes(const struct day *days, size_t n) {
    size_t yes = 0;
    for (size_t i = 0; i < n; i++) yes += days[i].rain;
    return yes;
}

/* The chance that two days drawn at random (with replacement) disagree. */
static double gini(const struct day *days, size_t n) {
    double p = (double)count_yes(days, n) / n;
    return 1 - p * p - (1 - p) * (1 - p);
}

/* "yes" only for a strict majority; an even split predicts "no". */
static bool majority(const struct day *days, size_t n) {
    return 2 * count_yes(days, n) > n;
}

/* Copy the days at or below the threshold into left, the rest into right.
 * Both arrays are the caller's, each with room for n days. */
static void partition(const struct day *days, size_t n, int f, double threshold,
                      struct day *left, size_t *n_left, struct day *right, size_t *n_right) {
    *n_left = *n_right = 0;
    for (size_t i = 0; i < n; i++) {
        if (days[i].x[f] <= threshold) left[(*n_left)++] = days[i];
        else right[(*n_right)++] = days[i];
    }
}

/* qsort calls this to compare two elements, passed as untyped pointers. */
static int compare_doubles(const void *a, const void *b) {
    double x = *(const double *)a, y = *(const double *)b;
    return (x > y) - (x < y);
}

/* Find the split with the lowest weighted Gini impurity. Features are tried
 * in column order and thresholds from smallest to largest; only a strictly
 * lower score replaces the best so far. Returns false if no feature has two
 * distinct values. */
static bool best_split(const struct day *days, size_t n,
                       double *score, int *feature, double *threshold) {
    bool found = false;
    double *values = malloc(n * sizeof *values);
    struct day *left = malloc(n * sizeof *left), *right = malloc(n * sizeof *right);
    for (int f = 0; f < N_FEATURES; f++) {
        for (size_t i = 0; i < n; i++) values[i] = days[i].x[f];
        qsort(values, n, sizeof *values, compare_doubles);
        /* Drop repeats, which sit next to each other once sorted. */
        size_t distinct = 1;
        for (size_t i = 1; i < n; i++)
            if (values[i] != values[distinct - 1]) values[distinct++] = values[i];

        for (size_t i = 1; i < distinct; i++) {
            double t = (values[i - 1] + values[i]) / 2;
            size_t n_left, n_right;
            partition(days, n, f, t, left, &n_left, right, &n_right);
            double s = (n_left * gini(left, n_left) + n_right * gini(right, n_right)) / n;
            if (!found || s < *score) {
                *score = s; *feature = f; *threshold = t;
                found = true;
            }
        }
    }
    free(values); free(left); free(right);
    return found;
}

static struct node *new_leaf(const struct day *days, size_t n) {
    struct node *leaf = malloc(sizeof *leaf);
    *leaf = (struct node){.left = NULL, .right = NULL, .label = majority(days, n)};
    return leaf;
}

/* Grow a tree on days, calling itself on each side of the best split. */
static struct node *grow(const struct day *days, size_t n, int depth) {
    double g = gini(days, n);
    if (g == 0 || n < MIN_SIZE || depth == MAX_DEPTH) return new_leaf(days, n);

    double score, threshold;
    int f;
    if (!best_split(days, n, &score, &f, &threshold) || score >= g)
        return new_leaf(days, n); /* no split makes the sides purer */

    struct day *left = malloc(n * sizeof *left), *right = malloc(n * sizeof *right);
    size_t n_left, n_right;
    partition(days, n, f, threshold, left, &n_left, right, &n_right);

    struct node *split = malloc(sizeof *split);
    *split = (struct node){.feature = f, .threshold = threshold,
                           .left = grow(left, n_left, depth + 1),
                           .right = grow(right, n_right, depth + 1)};
    free(left); free(right); /* the children keep no pointers into these */
    return split;
}

static bool predict(const struct node *n, const double x[]) {
    while (!is_leaf(n)) n = x[n->feature] <= n->threshold ? n->left : n->right;
    return n->label;
}

/* Print every node, parents before children, left before right. *k is the
 * running node number and *leaves the running leaf count. */
static void print_preorder(const struct node *n, int *k, int *leaves) {
    int id = ++*k;
    if (is_leaf(n)) {
        printf("node_%d_predict %s\n", id, n->label ? "yes" : "no");
        ++*leaves;
        return;
    }
    printf("node_%d_feature %s\n", id, FEATURES[n->feature]);
    printf("node_%d_threshold %.6f\n", id, n->threshold);
    print_preorder(n->left, k, leaves);
    print_preorder(n->right, k, leaves);
}

static size_t count_correct(const struct node *tree, const struct day *days, size_t n) {
    size_t correct = 0;
    for (size_t i = 0; i < n; i++) correct += predict(tree, days[i].x) == days[i].rain;
    return correct;
}

/* Nothing frees memory for you: release the children, then the node. */
static void free_tree(struct node *n) {
    if (!is_leaf(n)) { free_tree(n->left); free_tree(n->right); }
    free(n);
}

int main(int argc, char **argv) {
    if (argc < 2) { fprintf(stderr, "usage: %s data.csv\n", argv[0]); return 1; }

    struct day *train, *test;
    size_t n_train, n_test;
    read_days(argv[1], &train, &n_train, &test, &n_test);
    struct node *tree = grow(train, n_train, 0);

    printf("train_count %zu\n", n_train);
    printf("test_count %zu\n", n_test);
    int k = 0, leaves = 0;
    print_preorder(tree, &k, &leaves);
    printf("leaves %d\n", leaves);

    size_t correct = count_correct(tree, test, n_test);
    printf("train_accuracy %.6f\n", (double)count_correct(tree, train, n_train) / n_train);
    printf("correct %zu\n", correct);
    printf("accuracy %.6f\n", (double)correct / n_test);

    free_tree(tree);
    free(train);
    free(test);
    return 0;
}

// A decision tree classifier (CART with Gini impurity), using only the C++
// standard library.
//
// Grows a tree on the training days by always making the split that leaves
// the two sides purest, then predicts rain for the test days. Prints the tree
// in preorder: each node is either a split (feature and threshold) or a leaf
// (its prediction).

#include <algorithm>
#include <array>
#include <fstream>
#include <memory>
#include <optional>
#include <print>
#include <set>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>

constexpr std::array<const char*, 3> FEATURES = {"humidity", "pressure", "wind"};
constexpr int MAX_DEPTH = 3;
constexpr std::size_t MIN_SIZE = 5;

struct Day {
    std::array<double, FEATURES.size()> x;
    bool rain;
};

// One struct serves for both kinds of node; a leaf has no children. A
// std::unique_ptr owns the node it points to and deletes it when it goes out
// of scope, so dropping the root frees the whole tree with no free_tree.
// It can't be copied, only moved, so each node has exactly one owner.
struct Node {
    int feature = 0;
    double threshold = 0;
    std::unique_ptr<Node> left, right;
    bool rain = false;

    bool is_leaf() const { return left == nullptr; }
};

// A candidate split and its weighted impurity, before any children exist.
struct Candidate {
    double score;
    int feature;
    double threshold;
};

// Split one line at every comma.
std::vector<std::string> split(const std::string& line) {
    std::vector<std::string> fields;
    std::stringstream stream(line);
    std::string field;
    while (std::getline(stream, field, ',')) {
        fields.push_back(field);
    }
    return fields;
}

// The named column of a CSV file, as strings; the caller parses numbers.
std::vector<std::string> read_column(const std::string& path, const std::string& name) {
    std::ifstream file(path);
    std::string line;
    std::getline(file, line);
    std::vector<std::string> header = split(line);
    auto found = std::find(header.begin(), header.end(), name);
    if (found == header.end()) {
        throw std::runtime_error("no column " + name);
    }
    std::size_t index = found - header.begin();

    std::vector<std::string> values;
    while (std::getline(file, line)) {
        if (!line.empty()) {
            values.push_back(split(line).at(index));
        }
    }
    return values;
}

std::size_t count_yes(const std::vector<Day>& days) {
    // count_if counts the elements the lambda (an inline function) accepts.
    return std::ranges::count_if(days, [](const Day& d) { return d.rain; });
}

// The chance that two days drawn at random (with replacement) disagree.
double gini(const std::vector<Day>& days) {
    double p = static_cast<double>(count_yes(days)) / days.size();
    return 1 - p * p - (1 - p) * (1 - p);
}

// Rain only on a strict majority; an even split predicts no rain.
bool majority(const std::vector<Day>& days) {
    return 2 * count_yes(days) > days.size();
}

// The days at or below the threshold (left) or above it (right).
std::vector<Day> side(const std::vector<Day>& days, int f, double threshold, bool left) {
    std::vector<Day> out;
    for (const Day& d : days) {
        if ((d.x[f] <= threshold) == left) {
            out.push_back(d);
        }
    }
    return out;
}

// The candidate with the lowest weighted Gini impurity, or nothing
// (std::nullopt) if no feature has two distinct values. Features are tried in
// column order and thresholds from smallest to largest; only a strictly lower
// score replaces the best so far.
std::optional<Candidate> best_split(const std::vector<Day>& days) {
    std::optional<Candidate> best;
    for (int f = 0; f < static_cast<int>(FEATURES.size()); f++) {
        // A std::set keeps its values sorted and drops repeats.
        std::set<double> values;
        for (const Day& d : days) {
            values.insert(d.x[f]);
        }
        std::optional<double> previous;
        for (double value : values) {
            if (previous) {
                double t = (*previous + value) / 2;
                std::vector<Day> left = side(days, f, t, true);
                std::vector<Day> right = side(days, f, t, false);
                double score = (left.size() * gini(left) + right.size() * gini(right))
                    / days.size();
                if (!best || score < best->score) {
                    best = Candidate{score, f, t};
                }
            }
            previous = value;
        }
    }
    return best;
}

std::unique_ptr<Node> leaf(const std::vector<Day>& days) {
    auto node = std::make_unique<Node>();
    node->rain = majority(days);
    return node;
}

// Grow a tree on days, calling itself on each side of the best split.
std::unique_ptr<Node> grow(const std::vector<Day>& days, int depth) {
    double impurity = gini(days);
    if (impurity == 0 || days.size() < MIN_SIZE || depth == MAX_DEPTH) {
        return leaf(days);
    }
    std::optional<Candidate> best = best_split(days);
    if (!best || best->score >= impurity) {
        return leaf(days); // no split makes the sides purer
    }
    auto node = std::make_unique<Node>();
    node->feature = best->feature;
    node->threshold = best->threshold;
    node->left = grow(side(days, best->feature, best->threshold, true), depth + 1);
    node->right = grow(side(days, best->feature, best->threshold, false), depth + 1);
    return node;
}

bool predict(const Node& node, const Day& day) {
    const Node* n = &node;
    while (!n->is_leaf()) {
        n = day.x[n->feature] <= n->threshold ? n->left.get() : n->right.get();
    }
    return n->rain;
}

// Print every node, parents before children, left before right. k is the
// running node number and leaves the running leaf count, both passed by
// reference (&) so the calls below share them.
void print_preorder(const Node& node, int& k, int& leaves) {
    int id = ++k;
    if (node.is_leaf()) {
        std::println("node_{}_predict {}", id, node.rain ? "yes" : "no");
        leaves++;
        return;
    }
    std::println("node_{}_feature {}", id, FEATURES[node.feature]);
    std::println("node_{}_threshold {:.6f}", id, node.threshold);
    print_preorder(*node.left, k, leaves);
    print_preorder(*node.right, k, leaves);
}

std::size_t count_correct(const Node& tree, const std::vector<Day>& days) {
    return std::ranges::count_if(days, [&](const Day& d) { return predict(tree, d) == d.rain; });
}

int main(int argc, char* argv[]) {
    if (argc != 2) {
        std::println(stderr, "usage: program data.csv");
        return 1;
    }
    std::array<std::vector<std::string>, FEATURES.size()> columns;
    for (std::size_t f = 0; f < FEATURES.size(); f++) {
        columns[f] = read_column(argv[1], FEATURES[f]);
    }
    std::vector<std::string> rain = read_column(argv[1], "rain");
    std::vector<std::string> split_name = read_column(argv[1], "split");

    std::vector<Day> train, test;
    for (std::size_t row = 0; row < rain.size(); row++) {
        Day day{};
        for (std::size_t f = 0; f < FEATURES.size(); f++) {
            day.x[f] = std::stod(columns[f][row]);
        }
        day.rain = rain[row] == "yes"; // std::string compares by content
        (split_name[row] == "train" ? train : test).push_back(day);
    }
    std::unique_ptr<Node> tree = grow(train, 0);

    std::println("train_count {}", train.size());
    std::println("test_count {}", test.size());
    int k = 0, leaves = 0;
    print_preorder(*tree, k, leaves);
    std::println("leaves {}", leaves);

    std::size_t correct = count_correct(*tree, test);
    std::println("train_accuracy {:.6f}",
                 static_cast<double>(count_correct(*tree, train)) / train.size());
    std::println("correct {}", correct);
    std::println("accuracy {:.6f}", static_cast<double>(correct) / test.size());
    // No cleanup: when tree goes out of scope, every node is deleted.
}