{"id":11220,"date":"2025-07-18T01:03:00","date_gmt":"2025-07-18T05:03:00","guid":{"rendered":"https:\/\/www.both.org\/?p=11220"},"modified":"2025-07-16T12:41:42","modified_gmt":"2025-07-16T16:41:42","slug":"a-small-algol-68-project-part-3","status":"publish","type":"post","link":"https:\/\/www.both.org\/?p=11220","title":{"rendered":"A Small Algol 68 Project, Part 3"},"content":{"rendered":"
\r\n
\r\n \r\n <\/i>\r\n <\/a>\r\n <\/span>\r\n<\/div><\/div>\n


In memory of J. Kevin Douglas, a good friend and fellow fan of Algol 68<\/em><\/p>\n\n\n\n


In the previous article in this series<\/u><\/a>, we diverged from the main goal of this three article mini-series – determining a best fit of a line to a data set consisting of pairs of (x,y) values. We did so in order to have some utility routines for reading delimited text files (like CSV files), and to learn a bit more about Algol 68.<\/p>\n\n\n\n


In this article, we\u2019re going to program the calculations necessary. It could be a good moment to review the first article in the series in order to remember that we are solving a matrix equation<\/p>\n\n\n\n


X<\/em>T<\/em><\/sup> X c<\/u> = X<\/em>T<\/em><\/sup>y<\/u><\/em><\/p>\n\n\n\n


where the matrix X<\/em> has n<\/em> rows (n<\/em> being the number of data points) and two columns (the first being the x<\/em> data points, the second being all 1), and the vector y<\/u><\/em> has n<\/em> rows (containing all the y<\/em> data points).<\/p>\n\n\n\n


For those of you struggling to remember how to come up with the matrix X<\/em>T<\/em><\/sup>, you can think of spinning X <\/em>diagonally along the upper left – lower right dimension so that the first row is all the x<\/em> values (the first column of X<\/em>) and the second row is all the 1s (the second column of X<\/em>).<\/p>\n\n\n\n


Then we calculate the following dot products and assign them to the elements of X<\/em>T<\/em><\/sup> X<\/em><\/p>\n\n\n\n

<\/td>		(row 1 of X<\/em>T<\/em><\/sup>) dot<\/strong> (column 1 of X<\/em>)<\/td>(row 1 of X<\/em>T<\/em><\/sup>) dot<\/strong> (column 2 of X<\/em>)<\/td>
<\/td><\/tr>

<\/td>		(row 2 of X<\/em>T<\/em><\/sup>) dot<\/strong> (column 1 of X<\/em>)<\/td>(row 2 of X<\/em>T<\/em><\/sup>) dot<\/strong> (column 2 of X<\/em>)<\/td>
<\/td><\/tr><\/tbody><\/table>\n\n\n\n


We can express these dot products as sums of multiplied terms as follows<\/p>\n\n\n\n\n\n\n \n

<\/td>\n 		x1<\/sub>2<\/sup> + x2<\/sub>2<\/sup> + … + xn<\/sub>2<\/sup><\/em><\/td>\n x1<\/sub> + x2<\/sub> + … + xn<\/sub><\/em><\/td>\n
<\/td> \n<\/tr>\n

<\/td>\n 		x1<\/sub> + x2<\/sub> + … + xn<\/sub><\/em><\/td>\n 1 + 1 + … + 1<\/em><\/td>\n
<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n\n\n\n

On the right hand side, we proceed in a similar fashion:<\/p>\n\n\n\n

<\/td>		(row 1 of X<\/em>T<\/em><\/sup>) dot<\/strong> y<\/u><\/em><\/td>
<\/td><\/tr>

<\/td>		(row 2 of X<\/em>T<\/em><\/sup>) dot<\/strong> y<\/u><\/em><\/td>
<\/td><\/tr><\/tbody><\/table>\n\n\n\n

to calculate X<\/em>T<\/em><\/sup>y<\/u><\/em>:<\/p>\n\n\n\n\n\n\n \n

<\/td>\n 		x1<\/sub> \u00b7 y1<\/sub> + x2<\/sub> \u00b7 y2<\/sub> + … + xn<\/sub> \u00b7 yn<\/sub><\/em><\/td>\n
<\/td> \n<\/tr>\n

<\/td>\n 		y1<\/sub> + y2<\/sub> + … + yn<\/sub><\/em><\/td>\n
<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n\n\n\n

<\/a> Assembling the matrix equation from the data<\/h2>\n\n\n\n

Inspecting the above sums, we see that we can accumulate them “on the fly”, as we read the input data file; so no need to create an n<\/em> rows by 2 columns matrix, nor its transpose, nor an n<\/em> row vector.<\/p>\n\n\n\n

Note as well that the matrix X<\/em>T<\/em><\/sup> X <\/em>is symmetric, that is the value in row 1 column 2 is the same as the value in row 2 column 1; so we only need to accumulate one of those; we’ll do that into row 1 column 2.<\/p>\n\n\n\n

So now we know how to accumulate the values in the matrix equation as we read them; but how do we solve the matrix equation to get the slope m<\/em> and the y-intercept b<\/em> of the line?<\/p>\n\n\n\n

<\/a>Solving the matrix equation<\/h2>\n\n\n\n

In a simple case like this, we can use Gaussian elimination<\/u><\/a>, followed by back substitution. That is, we determine a factor a<\/em> that when row 1 values of X<\/em>T<\/em><\/sup> X<\/em> and X<\/em>T<\/em><\/sup>y<\/u><\/em> are muliplied by a <\/em>and added to the values of row 2, <\/em>the new value of X<\/em>T<\/em><\/sup> X<\/em>1,1<\/em><\/sub> is zero.<\/p>\n\n\n\n

Therefore, the value of a<\/em> needs to be calculated as<\/p>\n\n\n\n

a = -X<\/em>T<\/em><\/sup> X<\/em>2,1<\/em><\/sub> \u00f7 X<\/em>T<\/em><\/sup> X<\/em>1,1<\/em><\/sub><\/p>\n\n\n\n

and we update the second row of X<\/em>T<\/em><\/sup> X<\/em> as follows<\/p>\n\n\n\n

X<\/em>T<\/em><\/sup> X<\/em>2,1 <\/em><\/sub>\u2190 X<\/em>T<\/em><\/sup> X<\/em>2,1 <\/em><\/sub>+ a \u00b7 X<\/em>T<\/em><\/sup> X<\/em>1,1<\/em><\/sub> = X<\/em>T<\/em><\/sup> X<\/em>2,1 <\/em><\/sub>\u2013 X<\/em>T<\/em><\/sup> X<\/em>2,1 <\/em><\/sub>= 0<\/em><\/p>\n\n\n\n

X<\/em>T<\/em><\/sup> X<\/em>2,2 <\/em><\/sub>\u2190 X<\/em>T<\/em><\/sup> X<\/em>2,2 <\/em><\/sub>+ a \u00b7 X<\/em>T<\/em><\/sup> X<\/em>1,2<\/em><\/sub><\/p>\n\n\n\n

and the second element of X<\/em>T<\/em><\/sup> y<\/u><\/em><\/em> as follows<\/p>\n\n\n\n

X<\/em>T<\/em><\/sup> y<\/u><\/em>2<\/em><\/sub><\/u><\/em><\/em> <\/em><\/sub>\u2190 X<\/em>T<\/em><\/sup><\/em> y<\/u><\/em>2<\/em><\/sub><\/u><\/em><\/em>+ a \u00b7 X<\/em>T<\/em><\/sup><\/em> y<\/u><\/em><\/u><\/em><\/em><\/em><\/em>1<\/em><\/sub><\/p>\n\n\n\n

At this point we can back substitute for the value of b<\/em> by taking the dot product of the second row of X<\/em>T<\/em><\/sup> X <\/em>with the vector [m,b] <\/em>to solve for b<\/em><\/p>\n\n\n\n

0 \u00b7 m + X<\/em>T<\/em><\/sup> X<\/em>2,2 <\/em><\/sub>\u00b7 b = X<\/em>T<\/em><\/sup><\/em> y<\/u><\/em>2<\/em><\/sub><\/p>\n\n\n\n

and so<\/p>\n\n\n\n

b = X<\/em>T<\/em><\/sup><\/em> y<\/u><\/em>2<\/em><\/sub><\/em> <\/em><\/sub>\u00f7 X<\/em>T<\/em><\/sup> X<\/em>2,2<\/em><\/sub><\/p>\n\n\n\n

And then, knowing b<\/em>, similarly with the dot product of the first row of X<\/em>T<\/em><\/sup> X <\/em>for m<\/em><\/p>\n\n\n\n

X<\/em>T<\/em><\/sup> X<\/em>1,1<\/em><\/sub> \u00b7 m + X<\/em>T<\/em><\/sup> X<\/em>2,2<\/em><\/sub> \u00b7 b = X<\/em>T<\/em><\/sup><\/em> y<\/u><\/em><\/em><\/em>1<\/em><\/sub><\/p>\n\n\n\n

or<\/p>\n\n\n\n

m = <\/em>(X<\/em>T<\/em><\/sup> y<\/u><\/em><\/em><\/em>1<\/em><\/sub> \u2013 b \u00b7 X<\/em>T<\/em><\/sup> X<\/em>2,2 <\/em><\/sub>)\u00f7 X<\/em>T<\/em><\/sup> X<\/em>1,1<\/em><\/sub>
<\/p>\n\n\n\n

<\/a> Never mind the formulas, let\u2019s see the code<\/h2>\n\n\n\n

As I began this article, I realized I forgot one small but important utility routine, in order to convert the STRING fields on the data records to REAL. Here it is:<\/p>\n\n\n\n

     1\tPROC str to real = (STRING s) REAL:\n     2\tBEGIN\n     3\t    REAL v := 0.0, pt := 1.0, scale := 1.0;\n     4\t    BOOL negative := FALSE;\n     5\t    INT a0 = ABS \"0\";\n     6\t    FOR d FROM UPB s BY -1 TO LWB s DO\n     7\t        CHAR sd = s[d];\n     8\t        IF is digit(sd) THEN\n     9\t            v +:= (ABS sd - a0) * pt;\n    10\t            pt *:= 10.0\n    11\t        ELIF sd = \"-\" THEN\n    12\t            negative := TRUE\n    13\t        ELIF sd = \".\" THEN\n    14\t            scale := pt\n    15\t        FI\n    16\t    OD;\n    17\t    (negative | -v | v) \/ scale\n    18\tEND;<\/pre>\n\n\n\n
Nothing too crazy here. Worth emphasizing:<\/p>\n\n\n\n
This procedure returns a real value; hence the final expression – a conditional-clause in Algol 68 terms – yields a REAL<\/code> value<\/p>\n\n\n\n
Algol 68 Genie only uses ASCII characters in names, so the numeric characters are therefore known to be contiguous and the corresponding integer value can be calculated by taking the ABS of the character minus ABS \u201c0\u201d<\/code><\/p>\n\n\n\n
And finally, the least squares fitting code:<\/p>\n\n\n\n
     1\t#\n     2\t    Utility routines\n     3\t#\n       \n     4\tPR read \"str_to_real.a68\" PR\n     5\tPR read \"split.a68\" PR\n     6\tPR read \"each_delimited_line.a68\" PR\n       \n     7\t#\n     8\t    Perform a least-squares fit of a line to a data file containing (x,y) values\n     9\t#\n       \n    10\tBEGIN\n       \n    11\t    [1:2,1:2] REAL xtx := ((0.0, 0.0),(0.0,0.0));\n    12\t    [1:2] REAL xty := (0.0, 0.0);\n       \n    13\t    # Accumulate XTX and XTy #\n       \n    14\t    each delimited line(\"test_data.txt\", \",\", (STRING line, []STRING fields, INT line count) VOID: BEGIN\n    15\t        # skip the header line #\n    16\t        IF line count > 1 THEN\n    17\t            REAL xi = str to real(fields[1]), yi = str to real(fields[2]);\n    18\t            xtx[1,1] +:= xi * xi;\n    19\t            xtx[1,2] +:= xi;\n    20\t            xtx[2,2] +:= 1.0;\n    21\t            xty[1] +:= xi * yi;\n    22\t            xty[2] +:= yi\n    23\t        FI\n    24\t    END);\n    25\t    xtx[2,1] := xtx[1,2];\n       \n    26\t    # let's see the matrix and vector before Gauss #\n       \n    27\t    print((\"before Gauss:\",new line));\n    28\t    print((\"\u23a1\",xtx[1,1],\"   \",xtx[1,2],\"\u23a4  =  \u23a1\",xty[1],\"\u23a4\",new line));\n    29\t    print((\"\u23a3\",xtx[2,1],\"   \",xtx[2,2],\"\u23a6     \u23a3\",xty[2],\"\u23a6\",new line));\n       \n    30\t    # Gaussian elimination #\n       \n    31\t    REAL a = -xtx[2,1] \/ xtx[1,1];\n    32\t    xtx[2,1] +:= a * xtx[1,1];\n    33\t    xtx[2,2] +:= a * xtx[1,2];\n    34\t    xty[2] +:= a * xty[1];\n       \n    35\t    # let's see the matrix and vector after Gauss #\n       \n    36\t    print((\"after Gauss:\",new line));\n    37\t    print((\"\u23a1\",xtx[1,1],\"   \",xtx[1,2],\"\u23a4  =  \u23a1\",xty[1],\"\u23a4\",new line));\n    38\t    print((\"\u23a3\",xtx[2,1],\"   \",xtx[2,2],\"\u23a6     \u23a3\",xty[2],\"\u23a6\",new line));\n       \n    39\t    # back substitution #\n       \n    40\t    REAL b = xty[2] \/ xtx[2,2];\n    41\t    REAL m = (xty[1] - xtx[1,2] * b) \/ xtx[1,1];\n       \n    42\t    # let's see the equation #\n       \n    43\t    print((\"equation:\",new line));\n    44\t    print((\"y = \",fixed(m,0,7),\" * x + \",fixed(b,0,7),new line))\n       \n    45\tEND<\/pre>\n\n\n\n
Not too much new here, either, but worth mentioning:<\/p>\n\n\n\n