#! /usr/bin/gawk -f
# Last edited on 2019-11-06 19:48:57 by jstolfi

BEGIN {
  # Generates the data for the plot of the activity evolution function {F}
  # and its iterates, showing the fixpoints, for a very large network of 
  # GL neurons under the mean-field hypothesis.
  
  # Requires library "-f neuromat_mean_field_lib.gawk"
  
  # Caller must define (with "-v") parameters {kappa,lambda}, the GL dead-step option
  # {dead}, the number of iterations {per}, and the output file name {fname}
  # (without directory or extension.
  # 
  # Optionally the caller may also define a string {var} in {"mm", "mo", "mp", ... "pp"}
  # and a value {dv} to apply a perturbation {+dv} or {-dv} to {kappa} and/or {lambda}.
  #
  # The plot data is written to file "out/{fname}.txt".  
  # The caption in LaTeX format is written to "out/{fname}.tex"
  
  UNDEF = 99999999; # Substitute for {NaN}.
  
  # Ensure parameters are numeric:
  kappa = mknum(kappa) 
  lambda = mknum(lambda)
  dead = mknum(dead)
  per=mknum(per)  
  
  pfile = ("out/" fname ".txt") # File with plot data.
  tfile = ("out/" fname ".tex") # File with caption in LaTeX format.
  
  printf "----------------------------------------------------------------------\n" > "/dev/stderr";
  printf "neuron equations %s dead step\n", (dead == 0 ? "without" : "with") > "/dev/stderr";
  
  # Apply variation if specified:
  if (var != "") {
    # Each letter of {var} is "o" for no change, "p" for increment, "m" for decrement:
    dv = mknum(dv);
    split("", pert);
    for (k = 1; k <= 2; k++) { 
      vch = substr(var, k, 1);
      pert[k] = 0.0
      if (vch == "p") { pert[k]  += dv; } else if (vch == "m") { pert[k] += -dv; }
    }
    printf "perturbation kappa = %+12.8f  lambda = %+12.8f\n", pert[0], pert[1] > "/dev/stderr";
    kappa  += pert[0]; lambda += pert[1]
  }

  printf "parameters kappa = %15.8f  lambda = %15.8f\n", kappa,lambda > "/dev/stderr";
  printf "iterations = %d\n", per > "/dev/stderr";
  if (per < 1) { error(("bad per = " per)); }

  printf "$\\kappa=%+.8f$, $\lambda=%+.8f$;\\\\\n", kappa, lambda > tfile

  # Compute and write the breakpoints {zA,zB}:
  if (kappa != 0.0) {
    # Breakpoints:
    zA = -lambda/kappa; zB = (1-lambda)/kappa
    printf "breakpoints zA = %15.6f  zB = %15.6f\n", zA, zB > "/dev/stderr";
    printf "$\\rhoA=%+.8f$, $\\rhoB=%+.8f$;\\\\\n", zA, zB >> tfile
  } else {
    # Just in case:
    zA = (lambda <= 0 ? +101.0 : -101.0)
    zB = (lambda >= 1 ? -100.0 : +102.0)
    printf "$\\phantom{\\rhoA=0.0,;}$\\\\\n" >> tfile
  }

  # Compute and write the fixpoint at zero, if any:
  if (F(0,kappa,lambda,dead,per) == 0) { zD = 0.0; } else { zD = UNDEF; }
  write_fixpoint("quiescent fixpoint","D","D",zD,zA,zB,0,kappa,lambda,dead,per)

  # Compute non-trivial fixpoints and notable points (only for {per == 1})
  if (dead == 0) {
    # Single fixpoint {zL} from affine section:
    zL = UNDEF; # Unless defined
    if (per == 1) {
      # Fixpoint of the linear section:
      if (kappa != 1.0) {
        # Linear root:
        zL = -lambda/(kappa - 1)
      } else if ((0 < lambda) && (lambda < 1)) {
        # Degenerate linear root
        zL = lambda 
      }
    }
    write_fixpoint("affine fixpoint","L","L",zL,zA,zB,1,kappa,lambda,dead,per)
  } else {
    # Two potential fixpoints {zQm,zQp} from quadratic, plus maximum point {zM}:
    zQm = UNDEF
    zQp = UNDEF
    if {per == 1) {
      if (kappa != 0.0) {
        # Fixpoints of the quadratic section:
        A = kappa
        B = 1 + lambda - kappa; 
        C = -lambda
        Delta = B*B - 4*A*C;
        printf "quadratic B = %15.6f  Delta = %+17.8f\n", B, Delta > "/dev/stderr";
        # Quadratic maximum and roots:
        zM = (kappa - lambda)/(2*kappa)
        if (Delta >= 0) {
          zQm = (- B - sqrt(Delta))/(2*kappa)
          zQp = (- B + sqrt(Delta))/(2*kappa)
        }
      } else if (lambda > 0) {
        # Quadratic is actually linear:
        zQm = lambda/(1 + lambda)
      } 
    }
    write_fixpoint("quadratic maximum","M","M",zM,zA,zB,1,kappa,lambda,dead,per)
    write_fixpoint("quadratic low fixpoint","Qp","p",zQp,zA,zB,1,kappa,lambda,dead,per)
    write_fixpoint("quadratic high fixpoint","Qp","p",zQp,zA,zB,1,kappa,lambda,dead,per)
  }

  # Compute and write the saturated fixpoint, if any:
  if (F(1,kappa,lambda,dead,per) == 1) { zS = 1.0; } else { zS = UNDEF; } 
  write_fixpoint("saturated fixpoint","S","S",zS,zA,zB,0,kappa,lambda,dead,per)

  if (per > 1) {
     printf "Iteration %d;\\\\\n", per >> tfile
  }
  plot_zu = 0; # Set to 1 to plot the point where the quadratic slope is 1.
  
  # List of roots of {z = F(z)} found.
  split("", root); nroot = 0;
  
  # Plot the {per} iterate of {F}:
  np = 500; # Plot intervals.
  for (k = 0; k <= np; k++) { 
    rho = (k+0.0)/np;
    Frho = F(rho,kappa,lambda,dead,per); 
    printf "%9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", rho, Frho, -1, -1, -1, -1;
    if (k > 0) { check_for_root(ant,rho,kappa,lambda,dead,per); }
    # Save for next loop:
    ant = rho;
    Fant = Frho;    
  }
  if (per == 1) { 
    # Plot the characteristic points of {F}:
    if ((kappa) != 0.0) { 
      # Value of {rho} where middle section of {F} starts (may be out of {[0_1]}):
      za = -lambda/kappa; Fa = F(za,kappa,lambda,dead,1); ea = Fa - za;
      printf "%9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", -1, -1, za, Fa, -1, -1;
      printf "za = %9.6f  F(za) = %9.6f  ea = %9.6f\n", za, Fa, ea > "/dev/stderr";
      
      # Value of {rho} where middle section of {F} ends (may be out of {[0_1]}):
      zb = (1 - lambda)/kappa; Fb = F(zb,kappa,lambda,dead,1); eb = Fb - zb;
      printf "%9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", -1, -1, zb, Fb, -1, -1;
      printf "zb = %9.6f  F(zb) = %9.6f  eb = %9.6f\n", zb, Fb, eb > "/dev/stderr";
      
      if (dead != 0) {
        # Value of {rho} where the middle part has slope 1 (only if in {[za _ zb]}):
        zu = -lambda/(2*kappa);  
        if (((zu-za)*(zu-zb) <= 0) && plot_zu) { 
          Fu = F(zu,kappa,lambda,dead,1); eu = Fu - zu;
          printf "%9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", -1, -1, zu, Fu, -1, -1;
          printf "zu = %9.6f  F(zu) = %9.6f  eu = %9.6f\n", zu, Fu, eu > "/dev/stderr";
        }
      }
    }
  }
  # Plot the roots of {z = F(z)}:
  condense_roots();
  printf "plotting %d apparent roots\n", nroot > "/dev/stderr"; 
  for (r = 0; r < nroot; r++) {
    zr = root[r]; Fr = F(zr,kappa,lambda,dead,per); er = Fr - zr;
    DFr = DF(zr,kappa,lambda,dead,per);
    printf "%9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", -1, -1, -1, -1, zr, Fr;
    printf "root zr = %9.6f  Fr = %9.6f  DFr = %9.6f  er = %9.6f\n", zr, Fr, DFr, er > "/dev/stderr";
  }

  printf "----------------------------------------------------------------------\n" > "/dev/stderr";
  exit(0);
}

function write_fixpoint(name,ptag,ttag,z,zA,zB,strict,kappa,lambda,dead,per,  tmp,ze,Fz,DFz,ok) {
  # Writes the argument {z} to {stderr}, unless it is {UNDEF}.  
  # If {z} is in {[0 _ 1]} and in {[zA _ zB]}, also writes it to {tfile}, 
  # otherwise writes a phantom placeholder.  If {strict} is true then requires 
  # strict inclusion in those ranges. Works even if {zA > zB}
  # (case {kappa < 0}).
  #
  # Also computes the function value, and writes it if significantly
  # different from {z}.  Also computes the derivative. If {z} is 
  # equal to {zA} or {zB} computes the derivative 
  # slightly offset from {z} inside the range.
  # 
  # The symbol of the fixpoint is "z{ptag}" for {stderr}, and "\\rho{ttag}" for {tfile}.
  
  # Sort the range endpoints:
  if (zA > zB) { tmp = zA; zA = zB; zB = tmp }
  
  if (z != UNDEF) {
    # Compute the function's value:
    Fz = F(z,kappa,lambda,dead,per)
    
    # Compute the derivative:
    if (z == zA) {
      ze = z + 1.0e-10
    } else if (z == zB)
      ze = z - 1.0e-10
    } else {
      ze = z
    }
    DFz = DF(ze,kappa,lambda,dead,per)
    
    printf "%s fixpoint z%s = %15.6f  F(z%s) = %15.6f  F'(z%s) = %+15.6f\n", \
       name, ptag, z, ptag, Fz, ptag, DFz > "/dev/stderr";
       
    if (strict) {
      ok = ((z > zA) && (z < zB) && (z > 0.0) && (z < 1.0))
    } else {
      ok = ((z >= zA) && (z <= zB) && (z >= 0.0) && (z <= 1.0))
    }
    if (ok) {
      printf "$\\rho%s=%.8f$", ttag, z >> tfile
      if (abs(Fz - z) > 1.0e-8) { printf ", $F = %.8f$", Fz >> tfile }
      printf ", $F'=%+.8f$;\\\\\n", DFz >> tfile
    } else {
      printf "$\\phantom{\\rho%s=0,;F'}$;\\\\\n", ttag >> tfile
    }   
  } else {
    printf "$\\phantom{\\rho%d,;F'}$;\\\\\n", ttag >> tfile
  }
  

function check_for_root(za,zb,kappa,lambda,dead,per,  zm,Fm,em,Fa,ea,Fb,eb,zo,eo,Fo,nr,r,ztol,etol) {
  # Check for one root of {z = F(z)} in the interval {[za _ zb]},
  # append to global list {root[0..nroot-1], ???[0..nroot-1]}:
  ztol = 0.1/np;
  etol = 0.0000001;
  Fa = F(za,kappa,lambda,dead,per); ea = Fa-za;
  Fb = F(zb,kappa,lambda,dead,per); eb = Fb-zb;
  # printf "za = %9.6f zb = %9.6f ea*eb = %15.12f\n", za, zb, ea*eb > "/dev/stderr";
  if (ea*eb <= 0) {
    # Bisection search:
    printf "using bisection on za = %9.6f zb = %9.6f ea*eb = %15.12f\n", za, zb, ea*eb > "/dev/stderr";
    while (1) {
      zm = (za+zb)/2;
      Fm = F(zm,kappa,lambda,dead,per); em = Fm-zm;
      if (((zb-za) > ztol) || (abs(em) < etol)) { zo = zm; Fo = Fm; eo = em; break; }
      if (ea*em <= 0) 
        { zb = zm; Fb = Fm; eb = em; }
      else
        { za = zm; Fa = Fm; ea = em; }
    }
  } else if ( (abs(ea) < 0.0001) || (abs(eb) < 0.0001) ) {
    # Linear search:
    printf "using linear search on za = %9.6f zb = %9.6f ea*eb = %15.12f\n", za, zb, ea*eb > "/dev/stderr";
    zo = -1; Fo = UNDEF.99; eo = UNDEF.999;
    nr = 20; # Sub-steps.
    for (r = 0; r < nr; r++) {
      zm = za + (zb-za)*(r +0.5)/nr;
      Fm = F(zm,kappa,lambda,dead,per); em = Fm-zm;
      if (abs(em) < abs(eo)) { zo = zm; Fo = Fm; eo = em; }
    }
    if (abs(eo) > etol){ return; }
  } else {
    # Interval does not look promising:
    return;
  }
  # Append {zo} to global root list:
  printf "found apparent root zo = %9.6f  Fo = %9.6f  eo = %9.6f\n", zo, Fo, eo > "/dev/stderr";
  root[nroot] = zo; nroot++;
}

function condense_roots(   r,rootr,rini,rfin,ncond){
  # Condense the roots of {z = F(z)} so that they are about {1/np} apart:
  printf "found %d apparent roots\n", nroot > "/dev/stderr"; 
  ncond = 0; # Number of condensed roots.
  rini = 0; # First non-condensed root.
  for (r = 1; r <= nroot; r++) {
    rootr = (r >= nroot ? UNDEF : root[r]);
    if (abs(rootr-root[rini]) >= 1.0/np) {
      # Merge the bunch of roots:
      rfin = r - 1;
      if (rini < rfin) {
        root[ncond] = (root[rini]+root[rfin])/2;
      } else if (rini == rfin) {
        root[ncond] = root[rini];
      } else {
        error("duh?");
      }
      ncond++;
      rini = r;
    }
  }
  if (ncond < nroot) {
    nroot = ncond;
    printf "condensed to %d apparent roots\n", nroot > "/dev/stderr";
  }
}