Ensino (Teaching)

// -----------------------------------------------------------------------------
// rungek4: resolve a equação diferencial 
// dy/dx + y/x = sen(x)
// usando o método de Runge-Kutta de ordem 4
// -----------------------------------------------------------------------------
use Math only sin, cos;
use IO only openWriter;
config const h = 0.1;         // passo em x
const n = round(50/h):int;    // número de passos
var
   x,                         // variável independente
   y:                         // variável dependente
   [0..n] real;
x[0] = 0.0;                   // x inicial
y[0] = 0.0;                   // y inicial
// -----------------------------------------------------------------------------
// função que define a EDO dy/dx = sen(x) - y/x
// -----------------------------------------------------------------------------
proc ff(
   const in x: real,
   const in y: real
   ): real {
   if x == 0.0 then {
      return 0.0 ;
   }
   else {
      return sin(x) - y/x ;
   }
}
// -----------------------------------------------------------------------------
// rk4 implementa um passo do método de Runge-Kutta de ordem 4
// -----------------------------------------------------------------------------
proc rk4(
   const in x: real,
   const in y: real,
   const in h: real,
   const ref af: proc(
      const in ax: real,
      const in ay: real
      ): real
   ): real {
   var k1 = h*af(x,y);
   var k2 = h*af(x+h/2,y+k1/2);
   var k3 = h*af(x+h/2,y+k2/2);
   var k4 = h*af(x+h,y+k3);
   var yn = y + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0;
   return yn;
}
for i in 0..n-1 do {          // loop da solução numérica
   var xn1 = (i+1)*h;
   var yn1 = rk4(x[i],y[i],h,ff);
   x[i+1] = xn1;
   y[i+1] = yn1;
}
var erro = 0.0;               // calcula o erro relativo médio
   for i in 1..n do {
   var yana = sin(x[i])/x[i] - cos(x[i]);
   erro += abs( (y[i] - yana)/yana );
}
erro /= n ;
writef("erro relativo médio = %10.5dr",erro);
writeln();
const fou = openWriter("rungek4.out",locking=false);
for i in 0..n do {            // imprime o arquivo de saída
   fou.writef("%12.6dr %12.6dr\n",x[i],y[i]);
}
fou.close();

C ====================================================================
C ==> forroots : a Fortran 66 program to calculate the roots of a
C 2 nd - degree equation
C ====================================================================
00001 A = 1.0
00002 B = 2.0
00003 C = 0.5
00004 X1 = (-B-SQRT(B*B-4*A*C))/(2*A)
00005 X2 = (-B+SQRT(B*b+4*A*C))/(2*A)
00006 WRITE (6,*) X1
00007 WRITE (6,*) X2
00008 END

(* ================================================================
   ==> pasroots: a Pascal program: roots of a 2nd-degree equation
   ================================================================ *)
PROGRAM ROOTS;
CONST
   A = 1.0;
   B = 2.0;
   C = 0.5;
VAR X1, X2: REAL;
BEGIN
   X1 := (-B - SQRT(B*B - 4*A*C))/(2*A);
   X2 := (-B + SQRT(B*b + 4*A*C))/(2*A);
   WRITELN(X1);
   WRITELN(X2);
END.

// ===================================================================
// ==> croots: roots of the 2nd degree equation in C
// ===================================================================
#include <stdio.h>
#include <math.h>
int main(void) {
#define A 1.0
#define B 2.0
#define C 0.5
float x1,x2;
x1 = (-B - sqrt(B*B - 4*A*C))/(2*A);
x2 = (-B + sqrt(B*B + 4*A*C))/(2*A);
printf("x1 = %g\n",x1);
printf("x2 = %g\n",x2);
}

// ===================================================================
// ==> chplroots: roots of the 2nd degree equation in Chapel
// ===================================================================
const
   a = 1.0,
   b = 2.0,
   c = 0.5;
var x1 = (-b - sqrt(b*b - 4*a*c))/(2*a);
var x2 = (-b + sqrt(b*b + 4*a*c))/(2*a);
writef("x1 = %r\n",x1);
writef("x2 = %r\n",x2);

// ===================================================================
// ==> filewriter: separate use of file and channel
//
// note that enum ioMode { r = 1, cw = 2, rw = 3, cwr = 4, a = 5 }
// ===================================================================
use IO only ioMode, open;
var fou = open("filewriter.out", ioMode.cw);      // the file
var wou = fou.writer();                           // the file writer
wou.writeln("test of file and writer concepts");  // write w/ writer
wou.close();                                      // close writer
fou.close();                                      // close file

// ===================================================================
// ==> lower: read from stdin and print to stdout with all characters
// converted to lowercase
// ===================================================================
use IO, nnstrings;
var xc, yc: int;
while stdin.readCodepoint(xc) do { // read from stdin
   var sxc = chr(xc:int(32));      // convert xc to a 1-char string
   sxc = sxc.toLower();            // convert to lowercase
   write(sxc);                     // write single-character string
}

chpl lower.chpl
./lower < upper.txt 

// ===================================================================
// ==> wmixbin: writes a binary file containing different types of
// variables
// ===================================================================
use IO;
use nnstrings;
var fou = openWriter("mixbin.dat",serializer=new binarySerializer());
var i8: int(8) = -99:int(8);
var xx: real(64) = 17981.347;
var st = "Somewhere";
fou.write(i8);
fou.write(xx);
fou.write(st);
fou.close();

// ===================================================================
// ==> rmixbin: reads a binary file containing different types of
// variables
// ===================================================================
use IO;
var fin = openReader("mixbin.dat",
                     deserializer=new binaryDeserializer());
var i8: int(8);
var xx: real(64);
var st: string;
fin.read(i8);
fin.read(xx);
fin.readBinary(st,9);
writef("%+04i\n",i8);
writef("%+12.4dr\n",xx);
writef("%s\n",st);
fin.close();

// ===================================================================
// ==> exfra1: forall does not always work
// ===================================================================
use Random;
config const n = 1_600_000;
var A: [1..n] real;
fillRandom(A);
var sum = 0.0;
for i in 1..n do {
   sum += A[i];
}
writeln("sum from for    = ", sum);
sum = 0.0;
forall i in 1..n do {
   sum += A[i];          // sum acts as a constant inside the forall    @\label{lin:exfra1-ctsum}@
}
writeln("sum from forall = ", sum);

// ===================================================================
// ==> exfra2: forall does not always work either
// ===================================================================
use Random;
config const n = 1_600_000;
var A: [1..n] real;
fillRandom(A,0);
var sum = 0.0;
for i in 1..n do {
   sum += A[i];
}
writeln("sum from for    = ", sum);
sum = 0.0;
forall i in 1..n with (ref sum) do {
   sum += A[i];
}
writeln("sum from forall = ", sum);

// ===================================================================
// ==> exfra3: forall does not always work; time that it takes to run
// ===================================================================
use Random;
config const n = 1_600_000;
var A: [1..n] real;
fillRandom(A,0);
use Time;                // measure runtime
var runtime: stopwatch;  // with variable by the same name
var sum = 0.0;
runtime.start();
for i in 1..n do {
   sum += A[i];
}
runtime.stop();
writef("sum from for    = %12.4dr  took %12.8dr s\n",
       sum, runtime.elapsed());
sum = 0.0;
writeln("We have ",                          // run at Locales[0]       
   Locales[0].numPUs(logical=true)," cores");// with numPUs cores       
runtime.clear();
runtime.start();
forall i in 1..n with (ref sum) do {
   sum += A[i];
}
runtime.stop();
writef("sum from forall = %12.4dr  took %12.8dr s\n",
       sum, runtime.elapsed());

chpl --fast exfra3.chpl
./exfra3

// ===================================================================
// ==> exfra4: forall with the right "with intent" can be made to work
// easily
// ===================================================================
use Random;
config const n = 1_600_000;
var A: [1..n] real;
fillRandom(A,0);
use Time;                // measure runtime
var runtime: stopwatch;  // with variable by the same name
var sum = 0.0;
runtime.start();
for i in 1..n do {
   sum += A[i];
}
runtime.stop();
writef("sum from for    = %12.4dr  took %12.8dr s\n",
       sum, runtime.elapsed());
sum = 0.0;
writeln("We have ",                          // run at Locales[0]       
   Locales[0].numPUs(logical=true)," cores");// with numPUs cores       
runtime.clear();
runtime.start();
forall i in 1..n with (+ reduce sum) do {
   sum += A[i];
}
runtime.stop();
writef("sum from forall = %12.4dr  took %12.8dr s\n",
       sum, runtime.elapsed());

// ===================================================================
// ==> mansum-s: manual distribution of the sum among available
// processing units
// ===================================================================
use qdran;
use defNn;
use Time;
var A: [1..Nn] real;
seeqd(0);
for i in 1..Nn do A[i] = uranqd();
var npu = Locales[0].numPUs(logical=true);   // # of processing units
assert( Nn % npu == 0);   // make sure we can distribute the array 
var m = Nn/npu;      // # of elements to be summed by each proc unit
var psum: [0..npu-1] real = 0.0;
var runtime: stopwatch;
runtime.start();
for pu in 0..npu-1 do {
   var beg = pu*m + 1;
   var end = (pu+1)*m;
   for i in beg..end do {
      psum[pu] += A[i];
   }
}
var sum = 0.0;
for pu in 0..npu-1 do {
   sum += psum[pu];
}
runtime.stop();
writef("sum = %12.4dr  took %12.8dr s\n", sum,runtime.elapsed());

// ===================================================================
// ==> mansum-p: manual distribution of the sum among available
// processing units
// ===================================================================
use qdran;
use defNn;
use Time;
var A: [1..Nn] real;
seeqd(0);
for i in 1..Nn do A[i] = uranqd();
var npu = Locales[0].numPUs(logical=true);
assert( Nn % npu == 0);   // make sure we can distribute the array 
var m = Nn/npu;
var psum: [0..npu-1] real = 0.0; 
var runtime: stopwatch;
runtime.start();
coforall pu in 0..npu-1 do {
   var beg = pu*m + 1;
   var end = (pu+1)*m;
   for i in beg..end do {
      psum[pu] += A[i];
   }
}
var sum = 0.0;
for pu in 0..npu-1 do {
   sum += psum[pu];
}
runtime.stop();
writef("sum = %12.4dr  took %12.8dr s\n", sum,runtime.elapsed());

// ===================================================================
// ==> redsum: automatic distribution of the sum among available
// processing units with reduce
// ===================================================================
use qdran;          // quick and dirty random number generator from
                    // numerical recipes
use defNn;          // defines Nn
use Time;
// -------------------------------------------------------------------
// fill A with a simple random (serial) generator
// -------------------------------------------------------------------
var A: [1..Nn] real;
seeqd(0);
for i in 1..Nn do A[i] = uranqd();
var runtime: stopwatch;
runtime.start();
var sum = + reduce A;
runtime.stop();
writef("sum = %12.4dr  took %12.8dr s\n",sum, runtime.elapsed());

// ===================================================================
// ==> mansum-g: manual distribution of the sum among available gpu
// cores. uses Time.stopwatch
// ===================================================================
use qdran;                    // quick and dirty random number generator
                              // from Numerical Recipes
use defNn;                    // defines Nn
use GPU;
use Time;
var A: [1..Nn] real;
seeqd(0);
for i in 1..Nn do A[i] = uranqd();
on here.gpus[0] {
   var ah = A;                // arrays need to be copied to the GPU
   const npu = 3200;          // the number of GPU cores needs to be supplied
   assert( Nn % npu == 0);    // make sure we can distribute the array 
   var m = Nn/npu;            // # of elements to be summed by each core
   var psum: [0..npu-1] real = 0.0;
   var rt: stopwatch;
   rt.start();                // will run on the CPU!
   @assertOnGpu
   forall pu in 0..npu-1 do {
      var beg = pu*m + 1;
      var end = (pu+1)*m;
      for i in beg..end do {
         psum[pu] += ah[i];
      }
   }
   var sum = 0.0;
   for pu in 0..npu-1 do {
      sum += psum[pu];
   }
   rt.stop();
   writef("sum = %12.4dr  took %12.8dr s\n", sum,rt.elapsed());
}

// ===================================================================
// ==> abc-s very simple floating point operation on arrays; serial
// ===================================================================
use Random;
use Time;
use defNn;
const alpha = 3.0;
var runtime: stopwatch;
runtime.start();
var A,B,C: [1..Nn] real;
for i in 1..Nn do {
   B[i] = 0.5;
   C[i] = 0.5;
}
for i in 1..Nn do {
   A[i] = B[i] + alpha * C[i];
}
runtime.stop();
writef("took %12.8dr s\n", runtime.elapsed());

// ===================================================================
// ==> abc-p very simple floating point operation on arrays;
// parallel in cpu
// ===================================================================
use Random;
use Time;
use defNn;
const alpha = 3.0;
var runtime: stopwatch;
runtime.start();
var A,B,C: [1..Nn] real;
forall i in 1..Nn do {
   B[i] = 0.5;
   C[i] = 0.5;
}
forall i in 1..Nn do {
   A[i] = B[i] + alpha * C[i];
}
runtime.stop();
writef("took %12.8dr s\n", runtime.elapsed());

// ===================================================================
// ==> abc-g very simple floating point operation on arrays;
//     parallel in gpu
// ===================================================================
use Random;
use Time;
use defNn;
const alpha = 3.0;
var runtime: stopwatch;
runtime.start();
on here.gpus[0] {
   var A,B,C: [1..Nn] real;
   forall i in 1..Nn do {
      B[i] = 0.5;
      C[i] = 0.5;
   }
   forall i in 1..Nn do {
      A[i] = B[i] + alpha * C[i];
   }
}
runtime.stop();
writef("took %12.8dr s\n", runtime.elapsed());

// ===================================================================
// ==> mmv-s: test of sequential matrix-vector multiplication
// ===================================================================
use defmmv;
use Random;
use Time;
var runtime: stopwatch;
var totalsum = 0.0;
for r in 1..nrep do {
   serial {
      fillRandom(A,r);
      fillRandom(x,nrep+r);
   }
   runtime.start();
   mmv(A,x,y);
   runtime.stop();
   totalsum += (+ reduce y);
}
writeln(runtime.elapsed());
writeln("totalsum = ",totalsum);
// -------------------------------------------------------------------
// --> mmv: a serial implementation of matrix-vector multiplication
// -------------------------------------------------------------------
proc mmv(ref aA: [] real, ref ax: [] real, ref ay: [] real) {
   const (m,n) = aA.shape;
   assert (ax.shape == (n,));
   assert (ay.shape == (m,));
   ref A = aA.reindex({1..m,1..n});
   ref x = ax.reindex({1..n});
   ref y = ay.reindex({1..m});
   for i in 1..m do {
      var sum = 0.0;
      for j in 1..n do {
         sum += A[i,j]*x[j];
      }
      y[i] = sum;
   }
}

// ===================================================================
// ==> mmv-p: test of parallel matrix-vector multiplication
// ===================================================================
config const nrep = 1000;
use defmmv;
use Random;
use Time;
var runtime: stopwatch;
for r in 1..nrep do {
   fillRandom(A,r);
   fillRandom(x,nrep+r);
   runtime.start();
   mmv(A,x,y);
   runtime.stop();
}
writeln(runtime.elapsed());
// -------------------------------------------------------------------
// --> mmv: a parallel implementation of matrix-vector multiplication
// -------------------------------------------------------------------
proc mmv(
   const ref aA: [] real,
   const ref ax: [] real,
   ref ay: [] real
   ) where (aA.rank == 2) && (ax.rank == 1) && (ay.rank == 1) {
   const (m, n) = aA.shape;
   assert (ax.shape == (n,));
   assert (ay.shape == (m,));
   ref A = aA.reindex({1..m,1..n});
   ref x = ax.reindex({1..n});
   ref y = ay.reindex({1..m});
   forall i in 1..m do {
      var sum = 0.0;
      for j in 1..n do {
         sum += A[i,j]*x[j];
      }
      y[i] = sum;
   }
}

// ===================================================================
// ==> mmv-r: test of parallel matrix-vector multiplication using
// reduce
// ===================================================================
use defmmv;
use Random;
use Time;
var runtime: stopwatch;
for r in 1..nrep do {
   fillRandom(A,r);
   fillRandom(x,nrep+r);
   runtime.start();
   mmv(A,x,y);
   runtime.stop();
}
writeln(runtime.elapsed());
// -------------------------------------------------------------------
// --> mmv: a parallel implementation of matrix-vector multiplication
// -------------------------------------------------------------------
proc mmv(
   const ref A: [] real,
   const ref x: [] real,
   ref y: [] real
   ) where (A.rank == 2) && (x.rank == 1) && (y.rank==1) {
   const (m,n) = A.shape;
   assert (x.shape == (n,));
   assert (y.shape == (m,));
   y = [i in A.dim(0)] + reduce (A[i,..]*x);
}

// ===================================================================
// ==> morered: more reduction examples
// ===================================================================
const n = 100;
var s1 = + reduce ([k in 1..n] k);
var s2 = + reduce ([k in 1..n] k**2);
var s3 = + reduce ([k in 1..n] k**3);
writeln(s1);
writeln(s2);
writeln(s3);
var A: [1..n,1..n] real;
use Random;
fillRandom(A);
var amax = max reduce A;
writeln("amax = ",amax);

// ===================================================================
// ==> secder-c: a classical approach, using forall, to calculate the
// 2nd derivative of sin(x)
// ===================================================================
use Math only pi, sin;
const a = 0.0;
const b = 2*pi;
const n = 200;
const dx = (b-a)/n;
var x,s: [0..n] real;
var s2: [1..n-1] real;
forall i in 0..n do {
   x[i] = i*dx;
}
forall i in 0..n do {
   s[i] = sin(x[i]);
}
forall i in 1..n-1 do {
   s2[i] = (s[i+1] - 2*s[i] + s[i-1])/(dx**2);
}	
use IO only openWriter;
var fou = openWriter("secder-c.out");
for i in 1..n-1 do {
   fou.writef("%8.3dr %8.3dr %8.3dr\n",x[i],s2[i],-sin(x[i]));
}	
fou.close();

// ===================================================================
// ==> secder-m: a parallel program to calculate the 2nd derivative of
// sin(x)
// ===================================================================
use Math;
const a = 0.0;
const b = 2*pi;
const n = 200;
const dx = (b-a)/n;
const x = [i in 0..n] i*dx;   // expression
const s = sin(x);             // promotion
const s2: [1..n-1] real =                      // declare s2 and ...    
   (s[2..n] - 2*s[1..n-1] + s[0..n-2])/(dx**2);// fill with slicing     
use IO only openWriter;
var fou = openWriter("secder-m.out");
for i in 1..n-1 do {
   fou.writef("%8.3dr %8.3dr %8.3dr\n",x[i],s2[i],-sin(x[i]));
}  
fou.close();

// ===================================================================
// ==> promsinc: promotion
// ===================================================================
use Math;
const a = -pi;
const b = +pi;
const n = 100;
const dx = (b-a)/(2*n);
const x = [i in -n..n] i*dx;       // expression
const s = sinc(x);                 // promotion
use IO only openWriter;
var fou = openWriter("promsinc.out");
for i in -n..n do {
   fou.writef("%8.5dr %8.5dr\n",x[i],s[i]);
}  
fou.close();
// -------------------------------------------------------------------
// sinc is declared as a function of a scalar returning a scalar, but
// it operates seamlessly on an array with Chapel's promotion
// mechanism
// -------------------------------------------------------------------
proc sinc(const in x: real): real {
   if x == 0 then {
      return 1.0;
   }
   else {
      return sin(x)/x;
   }
}

// ===================================================================
// ==> strlin: counts occurrences of substr in standard input
// ===================================================================
use IO;
use nnstrings only ifind;
config const sbstr: string;
writeln("module strlin drives main()"); // I run first
// -------------------------------------------------------------------
// --> main: if it exists, main() always runs immediately after the
// module's executable statements
// -------------------------------------------------------------------
proc main(): int {
   var line: string;   // one line
   var found = 0;      // how many instances of sbstr were found?
   while stdin.readLine(line) do {
      var ret = line.ifind(sbstr);
      if ret > -1 then {
         write(line);
         found += 1;
      }
   }
   writeln();
   return found;
}

// ===================================================================
// ==> chintents: test chapel intents
// ===================================================================
for n in (0,1,4,6,-2) do {    // any iterable obj is game for a for   
   writefact(n);
}
// -------------------------------------------------------------------
// writefact: print the factorial of n
// -------------------------------------------------------------------
proc writefact(const in n: int) { 
   var thefact: int = max(int);
   writef("factorial of %i:\n",n);
   writeln(recfact(n));
   writeln(infact(n));
   outfact(n,thefact);
   writeln(thefact);
}
// -------------------------------------------------------------------
// --> recfact: calculate the factorial of n recursively
// -------------------------------------------------------------------
proc recfact(const in n: int): int {
   if n < 0 then {
      writeln("recfact: n cannot be negative");                         
      return(min(int));
   }
   else if n == 0 then {
      return 1;
   }
   else {
      return n*recfact(n-1);
   }
}
// -------------------------------------------------------------------
// --> infact: calculate the factorial of n non-recursively, using up
// n during the calculation. return as function
// -------------------------------------------------------------------
proc infact(in n: int): int {
   if n < 0 then {       // deal with negative n                      
      writeln("infact: n cannot be negative");                                                       
      return(min(int));                                                                              
   }                     // dealt                                    
   var fact = 1;
   while n > 1 do {
      fact *= n;
      n = n-1;
   }
   return fact;
}
// -------------------------------------------------------------------
// --> outfact: calculate the factorial of n non-recursively, using up
// n during the calculation. return as 'out' argument
// -------------------------------------------------------------------
proc outfact(in n: int, out fact: int) {
   if n < 0 then {
      writeln("outfact: n cannot be negative");
      fact = min(int);
      return;
   }
   fact = 1;
   var count=0;
   while n > 1 do {
      fact *= n;
      n = n-1;
   }
}

// ===================================================================
// ==> refino: inout and ref intents used to sort an array
// ===================================================================
use Random;
var rs = new randomStream(real,0); // a stream of rans seeded with 0    
const iras = 15;        	   // the size of the random arrays
var
   ira,                            // an array of random ints in [0,99]
   jra: [1..iras] int;	 	   // a copy of ira
for e in ira do {
   var r = rs.next();		   // get a real between 0.0 and 1.0
   e = (r*100):int; 		   // get an int between 0 and 99
   writef("%4i",e); 		   // print it
}
writeln();                         // newline
jra = ira;                         // copy ira to jra
writeln('-'*(4*iras));             // print a separator
binosort(ira);                     // sort with inout intent
for e in ira do {                  // print sorted ira
   writef("%4i",e);
}
writeln();                         // newline
writeln('-'*(4*iras));             // print a separator
brefsort(jra);                     // sort copy with ref intent
for e in jra do {                  // print sorte jra
   writef("%4i",e);
}
writeln();                         // newline (last)
// -------------------------------------------------------------------
// --> binosort: good old (and terribly inefficient) bubble sort: sort
// x inplace with inout intent: x is copied twice from and to the
// argument, on call and return
// -------------------------------------------------------------------
proc binosort(
   inout x: [] int
   ) where (x.rank == 1) {   // x size is any but only 1D arrays allowed   		
   for e in x do {           // no need for indices               
      for f in x do {
         if e < f then {
            e <=> f;
         }
      }
   }
}
// -------------------------------------------------------------------
// --> brefsort: good old (and terribly inefficient) bubble sort: sort
// x inplace with ref intent: all changes apply to the calling
// argument
// -------------------------------------------------------------------
proc brefsort(
   ref x: [] int
   ) where (x.rank == 1) {
   for e in x do {
      for f in x do {
         if e < f then {
            e <=> f;
         }
      }
   }
}

// ===================================================================
// ==> badintent: this program will not compile, because you cannot
// take the address of an expression 3*a + 4, etc.: it works as
// ------------------
// ref x = 3*a + 4;
// writeln(x);
// ------------------
// in an inlined version of badint, and this is an error
// ===================================================================
var a: int = 1;
var b: [1..4] int = 1;
badint(3*a + 4);         // cannot pass 3*a + 4 to ref
badarray(3*b + 4);	 // cannot pass 3*b + 4 to ref
// -------------------------------------------------------------------
// --> badint: cannot be used with an expression as its argument
// -------------------------------------------------------------------
proc badint(ref x: int) {
   writeln(x);
}
// -------------------------------------------------------------------
// --> badarray: cannot be used with an expression as its argument
// -------------------------------------------------------------------
proc badarray(ref x: [] int) {
   writeln(x);
}

// ===================================================================
// ==> goodintent: this program will compile, because you the formal
// argument has "in" intent: it works as
// ------------------
// var x = 3*a + 4;
// writeln(x);
// ------------------
// in an inlined version of goodint (same thing for goodarray)
// ===================================================================
var a: int = 1;
var b: [1..4] int = 1;
goodint(3*a + 4);        // can copy expression to x
goodarray(3*b + 4);	 // can copy expression to x
// -------------------------------------------------------------------
// --> goodint: can be used with an expression as its argument
// -------------------------------------------------------------------
proc goodint(in x: int) {
   writeln(x);
}
// -------------------------------------------------------------------
// --> goodarray: can be used with an expression as its argument
// -------------------------------------------------------------------
proc goodarray(in x: [] int) {
   writeln(x);
}

// ===================================================================
// ==> kwave1d-lax-s: solve a 1d kinematic kwave equation with an
// explicit method (Lax) serially. Uses memory efficiently.
// ===================================================================
param Nt = 20_000;                 // number of points in t             
param Nx = 10_000;                 // number of points in x             
param dt = 1.0/Nt;                 // delta t                           
param dx = 10.0/Nx;                // delta x                           
writef("# Nt = %9i\n",Nt);         // write Nt                          
writef("# Nx = %9i\n",Nx);         // write Nx                          
writef("# dt = %9.4dr\n",dt);      // write dt                          
writef("# dx = %9.4dr\n",dx);      // write dx                          
var u: [0..1,0..Nx] real = 0.0;    // initialize u to 0 for all values  
var un: [0..3,0..Nx] real = 0.0;   // store t=0.00, 0.25, 0.50, 0.75        
const n0 = 0;                      // prepare to write 4 times: 0,      
const n1 = round(0.25/dt):int;     // 0.25,                             
const n2 = round(0.50/dt):int;     // 0.50,                             
const n3 = round(0.75/dt):int;     // 0.75.
writeln("n1, n2, n3: ",n1," ",n2," ",n3);
for i in 0..Nx do {                // the initial condition loop        
   var xi = i*dx;                  // xi is local to the for            
   u[0,i] = IC(xi);                // fill the IC values                
}                                  // end loop                          
un[0,..] = u[0,..];                // store IC                          
config const cel = 2.0;  // wave celerity: can be changed at runtime    
const cour = cel*dt/dx;            // the Courant number                
writef("# Co = %9.2dr\n",cour);    // print it                          
assert (cour < 1.0);               // is the method stable?             
var nold = 0;                      // (n-1)th time step
var nnew = 1;                      // nth time step
var k = 1;                         // to store u[n0], u[n1], ...
use Time;
var runtime: stopwatch;
runtime.start();
for n in 1..Nt do {                // loop in time                      
   if n % 100 == 0 then writef("t = %8.4dr\n",n*dt);
   for i in 1..Nx-1 do {           // serial updating                   
      u[nnew,i] = 0.5*((u[nold,i+1] + u[nold,i-1]) -
                  cour*(u[nold,i+1] - u[nold,i-1]));
   }                               // end of serial updating            
   u[nnew,0] = 0.0;                // boundary conditions
   u[nnew,Nx] = 0.0;               // boundary conditions
   select n {                      // for t corresponding to        
      when n1,n2,n3 do {           // 0.25, 0.50, 0.75,
         un[k,..] = u[nnew,..];    // store u[n1], u[n2], u[n3]
         writef("k = %i\n",k);     // is k all right?
         k += 1;
      }
   }
   nnew <=> nold;                  // swap nnew and nold                                   
}                                  // end loop in time                  
runtime.stop();
writef("Elapsed time = %12.6dr  s\n",runtime.elapsed());
use IO;                            // prepare the output (text) file
const fou = openWriter("kwave1d-lax-s.out"); // declare fou channel     
fou.writef("#n=        " + "%10i "*4 + "\n",n0,n1,n2,n3);   // header   
for i in 0..Nx do {                // write n0, n1, n2, n3 timesteps    
   fou.writef("%10.4dr "*5 + "\n", // repeat string 5 times             
              i*dx,un[0,i],un[1,i],un[2,i],un[3,i]);
}
fou.close();
// -------------------------------------------------------------------
// --> IC: calculate the initial condition
// -------------------------------------------------------------------
proc IC(x: real): real {
   if (0 <= x) && (x <= 1.0) then {
      return 2.0*x*(1.0-x);
   }
   else {
      return 0.0;
   }
}

// ===================================================================
// ==> wave1d-lax-p1: solve a 1d kinematic wave equation with an
// explicit method (Lax) in parallel (forall). Uses memory
// efficiently.
// ===================================================================
param Nt = 20000;                  // number of points in t             
param Nx = 10000;                  // number of points in x             
param dt = 1.0/Nt;                 // delta t                           
param dx = 10.0/Nx;                // delta x                           
writef("# Nt = %9i\n",Nt);         // write Nt                          
writef("# Nx = %9i\n",Nx);         // write Nx                          
writef("# dt = %9.4dr\n",dt);      // write dt                          
writef("# dx = %9.4dr\n",dx);      // write dx                          
var u: [0..1,0..Nx] real = 0.0;    // initialize u to 0 for all values  
var un: [0..3,0..Nx] real = 0.0;   // store t=0, 0.25, 0.5, 0.75        
const n0 = 0;                      // prepare to write 4 times: 0,      
const n1 = round(0.25/dt):int;     // 0.25,                             
const n2 = round(0.50/dt):int;     // 0.50,                             
const n3 = round(0.75/dt):int;     // 0.75.                             
writeln("n1, n2, n3: ",n1," ",n2," ",n3);
for i in 0..Nx do {                // the initial condition loop        
   var xi = i*dx;                  // xi is local to the for            
   u[0,i] = IC(xi);                // fill the IC values                
}                                  // end loop                          
un[0,..] = u[0,..];                // store IC                          
config const cel = 2.0;  // wave celerity: can be changed at runtime    
const cou = cel*dt/dx;             // the Courant number                
writef("# Co = %9.2dr\n",cou);     // print it                          
assert (cou < 1.0);                // is the method stable?             
var nold = 0;                      // n-1 time
var nnew = 1;                       // n time
var k = 1;                         // to store u[n0], u[n1], ...
use Time;
var runtime: stopwatch;
runtime.start();
for n in 1..Nt do {                // loop in time                      
   if n % 100 == 0 then writef("t = %8.4dr\n",n*dt);
   forall i in 1..Nx-1 do {        // parallel updating                 
      u[nnew,i] = 0.5*((u[nold,i+1] + u[nold,i-1]) -
                       cou*(u[nold,i+1] - u[nold,i-1]));
   }                               // end of parallel updating          
   u[nnew,0] = 0.0;
   u[nnew,Nx] = 0.0;
   select n {                      // for t corresponding to        
      when n1,n2,n3 do {           // 0.25, 0.50, 0.75,
         un[k,..] = u[nnew,..];    // store u[n1], u[n2], u[n3]
         writef("k = %i\n",k);     // is k all right?
         k += 1;
      }
   }
   nnew <=> nold;                  // swap nnew and nold                                   
}                                  // end loop in time                  
runtime.stop();
writef("Elapsed time = %12.6dr  s\n",runtime.elapsed());
use IO;                            // prepare the output (text) file
const fou = openWriter("kwave1d-lax-p1.out"); // declare fou channel    
fou.writef("#n=        " + "%10i "*4 + "\n",n0,n1,n2,n3);   // header   
for i in 0..Nx do {                // write n0, n1, n2, n3 timesteps    
   fou.writef("%10.4dr "*5 + "\n", // repeat string 5 times             
         i*dx,un[0,i],un[1,i],un[2,i],un[3,i]);
}
fou.close();
// -------------------------------------------------------------------
// --> IC: calculate the initial condition
// -------------------------------------------------------------------
proc IC(x: real): real {
   if (0 <= x) && (x <= 1.0) then {
      return 2.0*x*(1.0-x);
   }
   else {
      return 0.0;
   }
}

// ===================================================================
// ==> wave1d-lax-p2: solve a 1d kinematic wave equation with an
// explicit method (Lax) in parallel (slicing and array
// operators). Uses memory efficiently.
// ===================================================================
param Nt = 20000;                  // number of points in t             
param Nx = 10000;                  // number of points in x             
param dt = 1.0/Nt;                 // delta t                           
param dx = 10.0/Nx;                // delta x                           
writef("# Nt = %9i\n",Nt);         // write Nt                          
writef("# Nx = %9i\n",Nx);         // write Nx                          
writef("# dt = %9.4dr\n",dt);      // write dt                          
writef("# dx = %9.4dr\n",dx);      // write dx                          
var u: [0..1,0..Nx] real = 0.0;    // initialize u to 0 for all values  
var un: [0..3,0..Nx] real = 0.0;   // store t=0, 0.25, 0.5, 0.75        
const n0 = 0;                      // prepare to write 4 times: 0,      
const n1 = round(0.25/dt):int;     // 0.25,                             
const n2 = round(0.50/dt):int;     // 0.50,                             
const n3 = round(0.75/dt):int;     // 0.75.                             
writeln("n1, n2, n3: ",n1," ",n2," ",n3);
for i in 0..Nx do {                // the initial condition loop        
   var xi = i*dx;                  // xi is local to the for            
   u[0,i] = IC(xi);                // fill the IC values                
}                                  // end loop                          
un[0,..] = u[0,..];                // store IC                          
config const cel = 2.0;  // wave celerity: can be changed at runtime    
const cou = cel*dt/dx;             // the Courant number                
writef("# Co = %9.2dr\n",cou);     // print it                          
assert (cou < 1.0);                // is the method stable?             
var nold = 0;                      // n-1 time
var nnew = 1;                      // n time
var k = 1;                         // to store u[n0], u[n1], ...
use Time;
var runtime: stopwatch;
runtime.start();
for n in 1..Nt do {                // loop in time                      
   if n % 100 == 0 then writef("t = %8.4dr\n",n*dt);
   u[nnew,1..Nx-1] =               // operate on ...                    
      0.5*((u[nold,2..Nx] + u[nold,0..Nx-2]) -     // array ...
           cou*(u[nold,2..Nx] - u[nold,0..Nx-2])); // slices            
   u[nnew,0] = 0.0;
   u[nnew,Nx] = 0.0;
   select n {                      // for t corresponding to        
      when n1,n2,n3 do {           // 0.25, 0.50, 0.75,
         un[k,..] = u[nnew,..];    // store u[n1], u[n2], u[n3]
         writef("k = %i\n",k);     // is k all right?
         k += 1;
      }
   }
   nnew <=> nold;                  // swap nnew and nold                                   
}                                  // end loop in time                  
runtime.stop();
writef("Elapsed time = %12.6dr  s\n",runtime.elapsed());
use IO;                            // prepare the output (text) file
const fou = openWriter("kwave1d-lax-p2.out");  // declare fou channel   
fou.writef("#n=        " + "%10i "*4 + "\n",n0,n1,n2,n3);   // header   
for i in 0..Nx do {                // write n0, n1, n2, n3 timesteps    
        fou.writef("%10.4dr "*5 + "\n", // repeat string 5 times        
                   i*dx,un[0,i],un[1,i],un[2,i],un[3,i]);
}
fou.close();
// -------------------------------------------------------------------
// --> IC: calculate the initial condition
// -------------------------------------------------------------------
proc IC(x: real): real {
   if (0 <= x) && (x <= 1.0) then {
      return 2.0*x*(1.0-x);
   }
   else {
      return 0.0;
   }
}

// =============================================================================
// ==> difadi2d-s: solve the 2-dimensional diffusion equation using the
// alternating direction implicit method with serial processing
//
// dphi/dt = D(d^2phi/dx^2 + d^2phi/dy^2),
//
// phi(0,x,y) = 1,
// phi(t,0,y) = phi(t,L,y) = 0,
// phi(t,x,0) = phi(t,x,M) = 0.
// =============================================================================
use Time only stopwatch;
var runtime: stopwatch;
runtime.start();
use difgrid2d;
const Fon = Dif*dt/((dl)**2);           // Fourier number
writef("Fo = %10.6dr\n",Fon);
use IO only openWriter, binarySerializer;
const fou = openWriter("difadi2d-s.dat",
               serializer = new binarySerializer(),locking=false);
var phi: [0..1,0..Nx,0..Ny] real = 0.0; // only two timesteps are necessary
for i in 0..Nx do {                     // sets up the initial condition
   for j in 0..Ny do {           
      phi[0,i,j] = IC(i,j);
   }
}
fou.write(phi[0,0..Nx,0..Ny]);          // imprime a condição inicial
var
   A,                                   // system matrix lower diagonal
   B,                                   // system matrix main diagonal
   C,                                   // system matrix upper diagonal
   D                                    // forcing vector
   : [1..Nn-1] real = 0.0;              // I am assuming Nx == Ny
// -----------------------------------------------------------------------------
// monta a matriz do sistema
// -----------------------------------------------------------------------------
A[1]      = 0.0;                        // just a marker
A[2..Nn-1] = -Fon;                      // fills lower diag
B[1..Nn-1] = 1.0 + 2*Fon;               // fills main diag
C[1..Nn-2] = -Fon;                      // fills upper diag
C[Nn-1]   = 0.0;                        // just a marker
// -----------------------------------------------------------------------------
// importa tridiag
// -----------------------------------------------------------------------------
use tridiag;
var iold = 0;                           // alternates memory positions
var inew = 1;
var n = 0;
while n < Nt do {                       // loop in time
   n += 1;                              // increment n
   // --------------------------------------------------------------------------
   // sweep in x direction
   // --------------------------------------------------------------------------
   for j in 0..Ny do {                  // BCs for x = 0, x = L
      phi[inew,0,j]  = BCy(n,j);
      phi[inew,Nx,j] = BCy(n,j);
   }
   for j in 1..Ny-1 do {                // Ny-1 sweeps in x
      D[1..Nx-1] = Fon*phi[iold,1..Nx-1,j-1]
         + (1.0 - 2*Fon)*phi[iold,1..Nx-1,j]
         + Fon*phi[iold,1..Nx-1,j+1];  
      D[1] += Fon*phi[inew,0,j];        // left BC
      D[Nx-1] += Fon*phi[inew,Nx,j];    // right BC 
      tridiag(A,B,C,D,phi[inew,1..Nx-1,j]);  // solve system in x
   }                                    
   if n % deln == 0 then {
      writeln(n);
      fou.write(phi[inew,0..Nx,0..Ny]);
   }
  // -----------------------------------------------------------------------------
  // sweep in y direction
  // -----------------------------------------------------------------------------
   inew <=> iold;
   n += 1;
   for i in 0..Nx do {                  // BCs for y = 0, y = L
      phi[inew,i,0]  = BCx(n,i);
      phi[inew,i,Ny] = BCx(n,i);
   }
   for i in 1..Nx-1 do {                // Nx-1 sweeps in y 
      D[1..Ny-1] = Fon*phi[iold,i-1,1..Ny-1]
         + (1.0 - 2*Fon)*phi[iold,i,1..Ny-1]
         + Fon*phi[iold,i+1,1..Ny-1];
      D[1] += Fon*phi[inew,i,0];
      D[Ny-1] += Fon*phi[inew,i,Ny];
      tridiag(A,B,C,D,phi[inew,i,1..Ny-1]);  // solve system in y
   }                                    
   if n % deln == 0 then {
      fou.write(phi[inew,0..Nx,0..Ny]);
      writeln(n);
   }
   inew <=> iold;
}
fou.close();                            // close output file
runtime.stop();
writeln("rtime = ",runtime.elapsed(), " ");
// -----------------------------------------------------------------------------
// initial condition
// -----------------------------------------------------------------------------
inline proc IC(
   const in i: int,
   const in j: int
   ): real {
   if i == 0 || i == Nx || j == 0 || j == Ny then {
      return 0.0;
   }
   return 1.0;
}
// -------------------------------------------------------------------
// boundary conditions
// -------------------------------------------------------------------                                                         
inline proc BCy(
   const in n: int,
   const in j: int
   ): real {
   return 0.0;
}
inline proc BCx(
   const in n: int,
   const in i: int
   ): real {
   return 0.0;
}

// =============================================================================
// ==> difadi2d-p: solve the 2-dimensional diffusion equation using the
// alternating direction implicit method with parallel processing
//
// dphi/dt = D(d^2phi/dx^2 + d^2phi/dy^2),
//
// phi(0,x,y) = 1,
// phi(t,0,y) = phi(t,L,y) = 0,
// phi(t,x,0) = phi(t,x,M) = 0.
// =============================================================================
use Time only stopwatch;
var runtime: stopwatch;
runtime.start();
use difgrid2d;
const Fon = Dif*dt/((dl)**2);           // Fourier number
writef("Fo = %10.6dr\n",Fon);
use IO only openWriter, binarySerializer;
const fou = openWriter("difadi2d-p.dat",
               serializer = new binarySerializer(),locking=false);
var phi: [0..1,0..Nx,0..Ny] real = 0.0; // only two timesteps are necessary
forall i in 0..Nx do {                  // sets up the initial condition
   forall j in 0..Ny do {           
      phi[0,i,j] = IC(i,j);
   }
}
fou.write(phi[0,0..Nx,0..Ny]);          // imprime a condição inicial
var
   A,                                   // system matrix lower diagonal
   B,                                   // system matrix main diagonal
   C                                    // system matrix upper diagonal
   : [1..Nn-1] real = 0.0;              // I am assuming Nx == Ny
// -----------------------------------------------------------------------------
// monta a matriz do sistema
// -----------------------------------------------------------------------------
A[1]      = 0.0;                        // just a marker
A[2..Nn-1] = -Fon;                      // fills lower diag
B[1..Nn-1] = 1.0 + 2*Fon;               // fills main diag
C[1..Nn-2] = -Fon;                      // fills upper diag
C[Nn-1]   = 0.0;                        // just a marker
// -----------------------------------------------------------------------------
// importa tridiag
// -----------------------------------------------------------------------------
use tridiag;
var iold = 0;                           // alternates memory positions
var inew = 1;
var n = 0;
while n < Nt do {                       // loop in time
   n += 1;                              // increment n
   // --------------------------------------------------------------------------
   // sweep in x direction
   // --------------------------------------------------------------------------
   forall j in 0..Ny do {               // BCs for x = 0, x = L
      phi[inew,0,j]  = BCy(n,j);
      phi[inew,Nx,j] = BCy(n,j);
   }
   forall j in 1..Ny-1 do {             // Ny-1 sweeps in x
      var D: [1..Nx-1] real;
      D[1..Nx-1] = Fon*phi[iold,1..Nx-1,j-1]
         + (1.0 - 2*Fon)*phi[iold,1..Nx-1,j]
         + Fon*phi[iold,1..Nx-1,j+1];  
      D[1] += Fon*phi[inew,0,j];        // left BC
      D[Nx-1] += Fon*phi[inew,Nx,j];    // right BC 
      tridiag(A,B,C,D,phi[inew,1..Nx-1,j]);  // solve system in x
   }                                    
   if n % deln == 0 then {
      writeln(n);
      fou.write(phi[inew,0..Nx,0..Ny]);
   }
  // -----------------------------------------------------------------------------
  // sweep in y direction
  // -----------------------------------------------------------------------------
   inew <=> iold;
   n += 1;
   forall i in 0..Nx do {               // BCs for y = 0, y = L
      phi[inew,i,0]  = BCx(n,i);
      phi[inew,i,Ny] = BCx(n,i);
   }
   forall i in 1..Nx-1 do {             // Nx-1 sweeps in y 
      var D: [1..Ny-1] real;
      D[1..Ny-1] = Fon*phi[iold,i-1,1..Ny-1]
         + (1.0 - 2*Fon)*phi[iold,i,1..Ny-1]
         + Fon*phi[iold,i+1,1..Ny-1];
      D[1] += Fon*phi[inew,i,0];
      D[Ny-1] += Fon*phi[inew,i,Ny];
      tridiag(A,B,C,D,phi[inew,i,1..Ny-1]);  // solve system in y
   }                                    
   if n % deln == 0 then {
      fou.write(phi[inew,0..Nx,0..Ny]);
      writeln(n);
   }
   inew <=> iold;
}
fou.close();                            // close output file
runtime.stop();
writeln("rtime = ",runtime.elapsed(), " ");
// -----------------------------------------------------------------------------
// initial condition
// -----------------------------------------------------------------------------
inline proc IC(
   const in i: int,
   const in j: int
   ): real {
   if i == 0 || i == Nx || j == 0 || j == Ny then {
      return 0.0;
   }
   return 1.0;
}
// -------------------------------------------------------------------
// boundary conditions
// -------------------------------------------------------------------                                                         
inline proc BCy(
   const in n: int,
   const in j: int
   ): real {
   return 0.0;
}
inline proc BCx(
   const in n: int,
   const in i: int
   ): real {
   return 0.0;
}

// =============================================================================
// ==> nlowess: my implementation of the lowess algorithm
//
// with most of the inspiration from
// https://medium.com/data-science-collective/loess-373d43b03564
//
// 2022-04-22T12:11:44 a new star is born
// 2022-04-22T17:34:40 essentially done ... hopefully
// 2022-04-23T18:54:06 moved into nstat
// 2025-02-19T08:33:44 is n_xx a single variable for all lowess_estimate calls?
//    yes, because n_xx is only used by lowess_estimate to calculate distances
//    locally
//
// 2025-05-05T13:54:12 using only the diagonal of W???
// =============================================================================
// -----------------------------------------------------------------------------
// I use tricubic as a scalar function only
// -----------------------------------------------------------------------------
private proc tricubic(const in x: real): real {
   var ax = abs(x) ;
   return if ax < 1.0 then (1-ax**3)**3 else 0.0 ;
}
// -----------------------------------------------------------------------------
// --> normalize_array: gess what: normalizes data
// -----------------------------------------------------------------------------
private proc normalize_array(
   const ref a: [] real,            // raw data
   ref n_a: [] real           // normalized data
   ): (real,real)
   where (a.rank == 1) {
   assert (a.shape == n_a.shape);
   var a_min = amin(a);
   var a_max = amax(a);
   n_a = (a - a_min) / (a_max - a_min);
   return (a_min, a_max);
}
// -----------------------------------------------------------------------------
// private variables!
// -----------------------------------------------------------------------------
private var
   degree = 1;      // degree of linear regression of this lowess
private var
   window: int;     // window of this lowess
private var
   xx_min,          // minimum of xx data
   xx_max,          // maximum of xx data
   yy_min,          // minimum of yy data
   yy_max:          // maximum of yy data
   real;
private var
   ldom = {0..-1};  // data domain initially empty
private var
   n_xx,            // normalized x data
   n_yy:            // normalized y data
   [ldom] real;

// -----------------------------------------------------------------------------
// --> init_lowess: normalize data arrays making local copies to n_xx and n_yy
// via normalize_array. For clarity, nlowess will use 1-based arrays throughout
// -----------------------------------------------------------------------------
proc init_lowess(
   ref xx: [] real,           // the raw x data
   ref yy: [] real,           // the raw y data
   const in win: int,         // the window, once and for all, of this loess
   const in deg: int = 1      // the degree of this loess
)  where (xx.rank == 1) && (yy.rank == 1) {
   assert (xx.shape == yy.shape);
   // --------------------------------------------------------------------------
   // set these global variables once and for all
   // --------------------------------------------------------------------------
   window = win;
   degree = deg;
   var nd: int = xx.size;
   ldom = {1..nd};            // allocate mem for n_xx and n_yy
   (xx_min,xx_max) = normalize_array(xx,n_xx);
   (yy_min,yy_max) = normalize_array(yy,n_yy);
   return;
}

// -----------------------------------------------------------------------------
// --> get_min_range: the indices of the window smallest
// -----------------------------------------------------------------------------
private proc get_min_range(
   ref distances: [] real,   // the local distances
   ref min_range: [] int,
   out maxdist: real
   ) where distances.rank == 1 {      
   var dinds = distances.indices;
   var indx = [i in dinds] i ;
   indxquickselect(indx,distances,dinds.first,dinds.last,window);
   min_range =  [k in 1..window] indx[k];
   maxdist = distances[indx[window]];
}
// -----------------------------------------------------------------------------
// --> get_weights: I suspect that max_distance == distances[indx[wind]]
// because of partial sorting
// -----------------------------------------------------------------------------
private proc get_weights(
   in maxdist: real,
   ref distances: [] real,
   ref min_range: [] int,
   ref weights: [] real
   ) where distances.rank == 1 && min_range.rank == 1 {
   assert(min_range.shape == weights.shape);
   var wind = min_range.size;
// -----------------------------------------------------------------------------
// normalize with weights
// -----------------------------------------------------------------------------
   for k in 1..wind do {
      weights[k] = tricubic(distances[min_range[k]]/maxdist);
   }
}
// -----------------------------------------------------------------------------
// we will need to come back to normalize
// -----------------------------------------------------------------------------
private inline proc normalize_x(
   const in x: real        // the data array?
   ): real {
   return (x - xx_min)/(xx_max - xx_min);
}

private inline proc denormalize_y(
   const in y: real) : real {
   return y * (yy_max - yy_min) + yy_min;
}
// -----------------------------------------------------------------------------
// --> estimate is the core loess algorithm
// -----------------------------------------------------------------------------
proc lowess_estimate(
      in x: real                 // point around which to do linear regression
): real {
   const m: int = window;     // size of neighborhood == window
   const n: int = degree;     // degree of linear regression == degree
   var n_x = normalize_x(x);  // normalize independent variable
   // --------------------------------------------------------------------------
   // do I really need to calculate all distances for each point?
   // --------------------------------------------------------------------------
   var nd = n_xx.size;
   var maxdist: real;
   var distances: [1..nd] real = abs(n_x - n_xx);
   var min_range: [1..window] int;
   var weights: [1..window] real;
   get_min_range(distances, min_range, maxdist);
   get_weights(maxdist, distances, min_range, weights);
   // --------------------------------------------------------------------------
   // here I use m,n instead of window, degree
   // estimate beta = [Xtr(n+1,m)W(m,m)X(m,n+1)]^(-1)Xtr(n+1,m)W(m,m)n_yy(m,1)
   // --------------------------------------------------------------------------
   var W: [1..m] real = 0.0;
   var X: [1..m,0..n] real;
   var Y: [1..m] real;
   // --------------------------------------------------------------------------
   // fill diagonal of W with weights
   // --------------------------------------------------------------------------
   // foreach i in 1..m do {
   //    W[i,i] = weights[i];
   // }
   foreach i in 1..m do {
      W[i] = weights[i];
   }   
   // -----------------------------------------------------------------------------
   // trickiest part is to build X
   // -----------------------------------------------------------------------------
   for i in 1..m do {
      X[i,0] = 1.0 ;
      X[i,1] = n_xx[min_range[i]];
      for j in 2..n do {
         X[i,j] = (n_xx[min_range[i]])**j;
      }
   }
   // -----------------------------------------------------------------------------
   // Y also needs to be built
   // -----------------------------------------------------------------------------
   foreach i in 1..m do {
      Y[i] = n_yy[min_range[i]];
   }
   // -----------------------------------------------------------------------------
   // ready for linear regression
   // -----------------------------------------------------------------------------
   var A: [0..n,1..m] real;
   var B: [0..n,0..n] real;
   var C: [0..n,1..m] real;
   var beta: [0..n] real;
   // -----------------------------------------------------------------------------
   // lots of matrix multiplications!
   // -----------------------------------------------------------------------------
   dot_mt_diagm(X,W,A);
   dot_mm(A,X,B);
   minvgj(B);
   dot_mmt(B,X,A);            // re-using A !!!
   dot_m_diagm(A,W,C);
   dot_mv(C,Y,beta);          // finally the parameters of the LLR
   // --------------------------------------------------------------------------
   // now we can estimate!
   // --------------------------------------------------------------------------
   var n_xp: [0..n] real ;
   n_xp[0] = 1.0;
   n_xp[1] = n_x ;
   foreach j in 2..n do {
      n_xp[j] = n_x**j;
   }      
   var y = dot_vtv(beta,n_xp);  // and the local estimate!
   return denormalize_y(y);
}

use nnstat;
use IO only openWriter;
use nprob only fillNormal;
use Math;
const n = 200;
var dx = 2*pi/n ;
var xx = [ i in 0..n-1] i*dx;           // data points in x
var ya = sin(xx);                       // the "real" data
var zz: [0..n-1] real;                  // noise around it
fillNormal(zz);                         // gaussian noise
var yy = ya + 0.25*zz;                  // the noisy data
var ye: [0..n-1] real;                  // the estimates of the real data
const fou = openWriter("exlow.out");    // open the output file
init_lowess(xx,yy,n/8,1);               // init the lowess filter
forall i in 0..n-1 do {                 // apply the lowess filter
   ye[i] = lowess_estimate(xx[i]);
}
for i in 0..n-1 do {                    // write the output file
   fou.writef("%12.8dr %12.8dr %12.8dr %12.8dr\n",xx[i],yy[i],ya[i],ye[i]);
}
fou.close();                            // close the output

config const describe = false;
const doc = "\
=================================================================================\
==> makeRaRn: after bestrad2 and polyalb, calculates Ra and Rn from CASTm.dat    \
                                                                                 \
./makeRaRn                                                                       \
=================================================================================\
";
if describe then {
   writeln(doc);
   exit(0);
}
// -----------------------------------------------------------------------------
// 2024-03-16T10:44:47 start over again!
// 2024-03-19T12:27:43 low-pass (lowess) filter
// -----------------------------------------------------------------------------
// =============================================================================
// modules 
// =============================================================================
use nangles only dec2rad, todec;
use IO;
use ndgrow;
use nevap only IniPar,Prescott, Brutsaert, Bolz, Radiation, SPrescott, ebrac, drac;
use nsunearth only ddse, rsds;
use nnstat only init_lowess, lowess_estimate, performance;
use Time only dateTime;
// =============================================================================
// declare variables that will be used (mostly) throughout: there are *14*
// fields in the datafile
// =============================================================================
const nf: int = 13;
// -----------------------------------------------------------------------------
// latitude is -30 S 27 minutes --> decimal --> radians
// -----------------------------------------------------------------------------
const lat = dec2rad(todec(-30.27));
IniPar(40.5,0.97);
// -----------------------------------------------------------------------------
// best radiation parameters
// -----------------------------------------------------------------------------
Prescott(0.2,0.575);
Brutsaert(0.68537,0.0803601);
Bolz(0.204501,0.444683);
drac(ebrac.bolz);
const p = [0.198969,0.000532508,5.60035e-06,-5.74951e-08,1.04339e-10];
inline proc albedo(nn: real): real {
   return p[0] + p[1]*nn + p[2]*nn**2 + p[3]*nn**3 + p[4]*nn**4 ;
}
// -----------------------------------------------------------------------------
// arbitrarily set the number of days at 10; will grow as needed
// -----------------------------------------------------------------------------
var xdat = {1..10};
var line: string;
var sdate: [xdat] string;     // the table of dates
var
   Ra4,                       // Ra measured by the CNR4
   Ra,                        // Ra estimated by evap.Radiation
   Rn4,                       // Rn measured by the CNR4
   Rn                         // Rn estimated by evap.Radiation
   : [xdat] real;
var n: int = 0;               // count the days
// =============================================================================
// read the data and organize the arrays
// =============================================================================
var cdat = openReader("CASTm.dat",locking=false);
while cdat.readLine(line) do {
   if line[0] == '#' || line[0] == '\n' then {
      continue;
   }
// -----------------------------------------------------------------------------
// read data from cdat
// -----------------------------------------------------------------------------
   line = line.strip();
   var field = line.split();
   assert(field.size == (nf + 1));
   n += 1;
   dgrow(n,xdat);
// -----------------------------------------------------------------------------
// sequence to calculate day of year
// -----------------------------------------------------------------------------   
   sdate[n] = field[0];
   var now = dateTime.strptime(sdate[n],"%Y-%m-%d");
   var yy = now.year;
   var mm = now.month;
   var dd = now.day;
   var janfirst = new dateTime(yy,1,1);
   var timediff = now - janfirst;
   var dayofyear = timediff.days + 1;
// -----------------------------------------------------------------------------
// calculate the albedo
// -----------------------------------------------------------------------------
   var alb = albedo(dayofyear:real);
   var Ta = field[1]:real;
   var ea = field[2]:real;
   var Rs4 = field[6]:real;
   var Rsr = -(field[7]:real);
   Ra4[n] = field[8]:real;
   var Lu4 = -(field[9]:real);
   Rn4[n] = Rs4 - Rsr + Ra4[n] - Lu4;
   var T0m = field[13]:real;
// -----------------------------------------------------------------------------
// now estimate all radiations needed
// -----------------------------------------------------------------------------
   var (delta,rr) = ddse(yy,mm,dd);
   var (Rsea,dsmax) = rsds(lat,rr,delta);
   var S = SPrescott(Rsea,Rs4);
   var Re: real;
   Radiation(alb,ea,Ta+273.15,T0m+273.15,S,Rs4,Ra[n],Re,Rn[n]);
   writef("%8.4dr %8.2dr %8.2dr %8.2dr\n",Ra4[n],Ra[n],Rn4[n],Rn[n]);
   if (ea == 0.0) then {
      writef("%s",sdate[n]);
      halt("perilous bend");
   }
}
writeln(n," days");
xdat = {1..n};
// // =============================================================================
// // now lowess comes into play
// // =============================================================================
// const win: int = 21;
// init_lowess(xx,Ts,win);
// forall i in 1..n do {
//    lowTs[i] = lowess_estimate(xx[i]);
// }
// =============================================================================
// print results for daily data
// =============================================================================
const dheader: string = "#YYY-mm-dd      Ra4     Ra^     Rn4     Rn^\n";
var cou = openWriter("RaRn.dat",locking=false);
cou.writef(dheader);
for i in 1..n do {
   cou.writef("%s ",sdate[i]);
   cou.writef(" %7.2dr %7.2dr %7.2dr %7.2dr", Ra4[i], Ra[i], Rn4[i], Rn[i]);
   cou.write("\n");
}
cou.close();
// =============================================================================
// performance, for Ra and Rn
// =============================================================================
var ldat = {1..n};
var x,y: [ldat] real;
var m = 0;
// -----------------------------------------------------------------------------
// first atmospheric radiation
// -----------------------------------------------------------------------------
for i in 1..n do {
   if (!isNan(Ra4[i])) && (!isNan(Ra[i])) then {
      m += 1;
      x[m] = Ra4[i];
      y[m] = Ra[i];
   }
}
ldat = {1..m};
var (r,Cd,bias,mae,rmse,dr) = performance(x,y);
writeln("For Ra:");
writef("r    = %8.4dr\n",r);
writef("Cd   = %8.4dr\n",Cd);
writef("bias = %8.2dr\n",bias);
writef("mae  = %8.2dr\n",mae);
writef("rmse = %8.2dr\n",rmse);
writef("dr   = %8.4dr\n",dr);
// -----------------------------------------------------------------------------
// now net radiation
// -----------------------------------------------------------------------------
ldat = {1..n};
x = 0.0;
y = 0.0;
m = 0;
for i in 1..n do {
   if (!isNan(Rn4[i])) && (!isNan(Rn[i])) then {
      m += 1;
      x[m] = Rn4[i];
      y[m] = Rn[i];
   }
}
ldat = {1..m};
(r,Cd,bias,mae,rmse,dr) = performance(x,y);
writeln("For Rn:");
writef("r    = %8.4dr\n",r);
writef("Cd   = %8.4dr\n",Cd);
writef("bias = %8.2dr\n",bias);
writef("mae  = %8.2dr\n",mae);
writef("rmse = %8.2dr\n",rmse);
writef("dr   = %8.4dr\n",dr);

C ====================================================================
C ==> forroots : a Fortran 66 program to calculate the roots of a
C 2 nd - degree equation
C ====================================================================
00001 A = 1.0
00002 B = 2.0
00003 C = 0.5
00004 X1 = (-B-SQRT(B*B-4*A*C))/(2*A)
00005 X2 = (-B+SQRT(B*b+4*A*C))/(2*A)
00006 WRITE (6,*) X1
00007 WRITE (6,*) X2
00008 END

(* ================================================================
   ==> pasroots: a Pascal program: roots of a 2nd-degree equation
   ================================================================ *)
PROGRAM ROOTS;
CONST
   A = 1.0;
   B = 2.0;
   C = 0.5;
VAR X1, X2: REAL;
BEGIN
   X1 := (-B - SQRT(B*B - 4*A*C))/(2*A);
   X2 := (-B + SQRT(B*b + 4*A*C))/(2*A);
   WRITELN(X1);
   WRITELN(X2);
END.

// ===================================================================
// ==> croots: roots of the 2nd degree equation in C
// ===================================================================
#include <stdio.h>
#include <math.h>
int main(void) {
#define A 1.0
#define B 2.0
#define C 0.5
float x1,x2;
x1 = (-B - sqrt(B*B - 4*A*C))/(2*A);
x2 = (-B + sqrt(B*B + 4*A*C))/(2*A);
printf("x1 = %g\n",x1);
printf("x2 = %g\n",x2);
}

// ===================================================================
// ==> chplroots: roots of the 2nd degree equation in Chapel
// ===================================================================
const
   a = 1.0,
   b = 2.0,
   c = 0.5;
var x1 = (-b - sqrt(b*b - 4*a*c))/(2*a);
var x2 = (-b + sqrt(b*b + 4*a*c))/(2*a);
writef("x1 = %r\n",x1);
writef("x2 = %r\n",x2);