#! /usr/bin/gawk -f
# Last edited on 2020-01-17 16:47:27 by jstolfi

BEGIN {
  # Writes to standard out the data for the plot of the phase diagram 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}_{min,max}},
  # the GL dead-step option {dead}, the {dense} option to explore the 
  # diagram systemaptically, and the desired {region} (-1 for all).
  
  # Ensure parameters are numeric:
  kappa_min = mknum(kappa_min) 
  kappa_max = mknum(kappa_max) 
  lambda_min = mknum(lambda_min)
  lambda_max = mknum(lambda_max)
  dead = mknum(dead)
  dense = mknum(dense)
  region = mknum(region)
  
  printf "neuron equations %s dead step\n", (dead == 0 ? "without" : "with") > "/dev/stderr";
  printf "parameters ranges:\n" > "/dev/stderr";
  printf "  kappa = [ %.8f _ %.8f ]\n", kappa_min, kappa_max  > "/dev/stderr";
  printf "  lambda = [ %.8f _ %.8f ]\n", lambda_min, lambda_max  > "/dev/stderr";
  if (region >= 0) { printf "  region = %d\n", region  > "/dev/stderr"; }
  
  # Signals "no number":
  NIL = -999999999
  
  # Shorter names:
  xlo = kappa_min; xhi = kappa_max
  ylo = lambda_min; yhi = lambda_max
  
  # Number of points in plots:
  NP = 100
  
  # Largest abs coordinate in window:
  M = 2*max(max(abs(xlo), abs(xhi)), max(abs(ylo), abs(yhi)))
  
  # RF = 0.028*(xhi-xlo)  # Half-side of plot.
  RF = 0.032*(xhi-xlo)  # Half-side of plot.

  if (dead == 0) {
    draw_phases_no_dead(NP,M,RF)
  } else {
    draw_phases_dead(NP,M,RF,xlo,xhi,ylo,yhi)
  }
  exit(0);
}
    
function draw_phases_no_dead(NP,M,RF,  r) {
  
    # Model without dead step:
    
    # Grid lines:
    write_grid_lines(M)

    # Symmetry line:
    write_segment(-M, (1+M)/2,  +M, (1-M)/2,   6)
    
    if ((region < 0) || (region == 1)) { 
      # Fill of region 1:
      write_fill_corner(+1-M,   0,  11)
      write_fill_corner(  +1,   0,  11)
      write_fill_corner(+1+M,  -M,  11)
      write_fill_corner(+1-M,  -M,  11)
      write_fill_corner(+1-M,   0,  11)
      printf "\n"
    }
    
    if ((region < 0) || (region == 2)) { 
      # Fill of region 2:
      write_fill_corner(+1-M, +M,  12)
      write_fill_corner(  +1,  0,  12)
      write_fill_corner(+1+M,  0,  12)
      write_fill_corner(+1+M, +M,  12)
      write_fill_corner(+1-M, +M,  12)
      printf "\n"
    }
    
    if ((region < 0) || (region == 3)) { 
      # Fill of region 3:
      write_fill_corner(-1,  0,  13)
      write_fill_corner(+1,  0,  13)
      write_fill_corner(-1, +2,  13)
      write_fill_corner(-1,  0,  13)
      printf "\n"
    }
    
    if ((region < 0) || (region == 4)) { 
      # Fill of region 4:
      write_fill_corner(+1,    0,  15)
      write_fill_corner(+1+M,  0,  15)
      write_fill_corner(+1+M, -M,  15)
      write_fill_corner(+1,    0,  15)
      printf "\n"
    }
    
    if ((region < 0) || (region == 5)) {
      # Fill of region 5:
      write_fill_corner(+1-M,  0,  14)
      write_fill_corner(-1,    0,  14)
      write_fill_corner(-1,   +2,  14)
      write_fill_corner(+1-M, +M,  14)
      write_fill_corner(+1-M,  0,  14)
      printf "\n"
    }
    
    # Boundaries between major phases:
    write_segment(-M,   0,   +M,    0,  0)  # 1-5, 1-3, 2-4
    write_segment(-1,   0,   -1,   +2,  0)  # 3-5
    write_segment(-M, +M+1,  +M, -M+1,  0)  # 2-5, 2-3, 1-4
    
    # Boundaries between sub-phases:
    write_segment(-M,  1,    0,  1,  1) # 5a-5c, 3b-3d
    write_segment( 0,  0,    0, +1,  1) # 3c-3d     
    write_segment(-M, +M,    0,  0,  1) # 5b-5c, 3a-3d

    for (r = 1; r <= 4; r++) {
      write_ratbound( 0, +M, r, 0, 1) # 1{r}-1{r+1}
      write_ratbound( 0, +M, r, 1, 1) # 2{r}-2{r+1}
    }

    # Region labels:
    if ((region < 0) || (region == 1)) { write_label(-0.2, -0.6, 2,  NIL,  NIL,  16, "1") }
    if ((region < 0) || (region == 2)) { write_label(+1.5, +1.5, 2,  NIL,  NIL,  16, "2") }
    if ((region < 0) || (region == 3)) { write_label(-0.3, +0.4, 2,  NIL,  NIL,  16, "3") }
    if ((region < 0) || (region == 4)) { write_label(+2.4, -0.6, 2,  NIL,  NIL,  16, "4") }
    if ((region < 0) || (region == 5)) { write_label(-1.7, +1.3, 2,  NIL,  NIL,  16, "5") }

    if (dense) {
      write_dense_F_plots(NP,M,RF,xlo,xhi,ylo,yhi)
    } else {
      if ((region < 0) || (region == 1)) { 
        # Sample F-plots in phase 1:
        write_F_plot("xxx",  NIL,  NIL, -1.0, -1.2,  RF)
        write_F_plot("xxx",  NIL,  NIL, +1.5, -1.0,  RF)
      }

      if ((region < 0) || (region == 2)) {
        # Sample F-plots in phase 2:
        write_F_plot("xxx",  NIL,  NIL, +1.5, +0.5,  RF)
        write_F_plot("xxx",  NIL,  NIL, +0.5, +2.0,  RF)
      }

      if ((region < 0) || (region == 3)) {
        # Sample F-plots in phase 3: 
        write_F_plot("xxx",  NIL,  NIL, +0.3, +0.3,  RF)
        write_F_plot("xxx",  NIL,  NIL, -0.7, +1.3,  RF)
        write_F_plot("xxx",  NIL,  NIL, -0.7, +0.3,  RF)
        write_F_plot("xxx",  NIL,  NIL, -0.3, +0.8,  RF)
      }

      if ((region < 0) || (region == 4)) {
        # Sample F-plots in phase 4:
        write_F_plot("xxx",  NIL,  NIL, +1.9, -0.3,  RF)
      }

      if ((region < 0) || (region == 5)) {
        # Sample F-plots in phase 5:
        write_F_plot("xxx",  NIL,  NIL, -1.8, +2.3,  RF)
        write_F_plot("xxx",  NIL,  NIL, -2.2, +1.6,  RF)
        write_F_plot("xxx",  NIL,  NIL, -1.7, +0.5,  RF)
      }
    }
}

function draw_phases_dead(NP,M,RF,xlo,xhi,ylo,yhi,  r,xq,yq,dxq,RQ) {
  
    # Model WITH dead step:
    
    # Grid lines:
    write_grid_lines(M)
    
    # Boundaries between major phases:
    write_segment( -M,     0,  +M,     0,  0)  # {\lam = 0} 
    write_segment( -M, 1+M/2,  +M, 1-M/2,  0)  # {\lam = 1-\kap/2}
    write_segment(  0,     1, 0.5,     0,  0)  # {\lam = 1-2\kap}
    write_segment(1-M,    +M,   1,     0,  0)  # {\lam = 1-\kap}
    write_segment(  0,     1,  +M,     1,  0)  # {\lam = 1}

    write_parabola(1,      4,   4,     0,  +1, 0)  # {(1 + \lam + \kap)^2 = 4\kap}   
    # write_parabola(0,    +M,   0,     1,   0, 0)  # {(1 + \lam + \kap)^2 = 1}   

    # Boundaries between sub-phases:
    write_segment(-M,    +M,  +M,    -M,  1)  # {\lam = -\kap}
    write_segment( 1,     0,  1+M,   -M,  1)  # {\lam = 1-\kap}
    
    # write_parabola(0,  1,   4,     0,  +1, 1)  # {(1 + \lam + \kap)^2 = 4\kap}   
    write_parabola(4, +M,   4,     0,  +1, 1)  # {(1 + \lam + \kap)^2 = 4\kap}   

    write_funnybola(0.0001, +M,  0, 1)  # {\kap(\kap+\lam)^2 = \kap-\lam}   
    write_funnybola(-M, -0.0001,  0, 1)  # {\kap(\kap+\lam)^2 = \kap-\lam}   

    # Region labels:
    write_label( -2.57, +1.25, 2, NIL, NIL, 16, "a"); #  1 
    write_label( -3.47, +3.10, 2, NIL, NIL, 16, "b"); #  2 
    write_label( -2.60, +3.12, 2, NIL, NIL, 16, "c"); #  3 
    write_label( -0.60, +0.93, 2, NIL, NIL, 16, "d"); #  4 
    write_label( +0.27, +0.30, 2, NIL, NIL, 16, "e"); #  5 
    write_label( +0.60, +0.30, 2, NIL, NIL, 16, "f"); #  6 
    write_label( +1.03, +0.33, 2, NIL, NIL, 16, "g"); #  7 
    write_label( +3.10, +0.53, 2, NIL, NIL, 16, "h"); #  8 
    write_label( +7.00, -2.60, 2, NIL, NIL, 16, "i"); #  9 
    write_label( +1.32, -0.85, 2, NIL, NIL, 16, "k"); # 11 
    write_label( -0.70, -0.87, 2, NIL, NIL, 16, "m"); # 12 
    write_label( +2.00, -0.10, 2, NIL, NIL, 16, "n"); # 13 
    
    # Sample F plots:
  
    if (dense) {
      write_dense_F_plots(NP,M,RF,xlo,xhi,ylo,yhi)
    } else {
      if ((region <0) || (region == 1)) {
        write_F_plot("a1", -4, +2,   -4.0, +2.8, RF)
        write_F_plot("a2", -4, +1,   -4.0, +0.2, RF)
        write_F_plot("a3", -3, +1,   -2.1, +1.9, RF)
        write_F_plot("a4", -2, +1,   -0.2, +0.1, RF)
      }
      
      if ((region <0) || (region == 2)) {
        write_F_plot("b1", -2, +4,   -4.0, +3.8, RF)
        write_F_plot("b2", -1, +4,   -4.0, +3.2, RF)
        write_F_plot("b3", 00, +4,   -2.4, +2.3, RF)
      }

      if ((region <0) || (region == 3)) {
        write_F_plot("c1", -3, +5,   -4.0, +4.2, RF)
        write_F_plot("c2", -2, +5,   -4.0, +4.8, RF)
        write_F_plot("c3", -1, +5,   -2.0, +2.2, RF)
        write_F_plot("c4", 00, +5,   -0.4, +1.3, RF)
      }

      if ((region <0) || (region == 4)) {
        write_F_plot("d1", 00, +2,   -1.6, +1.7, RF)
        write_F_plot("d2", +1, +2,   -0.1, +0.9, RF)
        write_F_plot("d3", +2, +2,   -0.1, +0.2, RF)
      }

      if ((region <0) || (region == 5)) {
        if (region < 0) { 
          xq = 1.00; dxq = 1.00; yq = 3.00; RQ = RF  
        } else { 
          xq = 0.30; dxq = 0.50; yq = 0.75; RQ = 0.20
        }
        write_F_plot("e1", xq+0*dxq, yq, +0.05, +0.80, RQ)
        write_F_plot("e2", xq+1*dxq, yq, +0.10, +0.10, RQ)
        write_F_plot("e3", xq+2*dxq, yq, +0.40, +0.10, RQ)
      }

      if ((region <0) || (region == 6)) {
        write_F_plot("f1", +1, +4,   +0.4, +0.6, RF)
        write_F_plot("f2", +2, +4,   +0.6, +0.2, RF)
        write_F_plot("f3", +3, +4,   +0.8, +0.1, RF)
      }

      if ((region <0) || (region == 7)) {
        write_F_plot("g1", +1, +5,   +0.2, +0.7, RF)
        write_F_plot("g2", +2, +5,   +1.1, +0.2, RF)
        write_F_plot("g3", +3, +5,   +1.8, +0.1, RF)
      }

      if ((region <0) || (region == 8)) {
        write_F_plot("h1", +4, +2,   +0.3, +0.9, RF)
        write_F_plot("h2", +5, +2,   +2.1, +0.2, RF)
        write_F_plot("h3", +6, +2,   +5.0, +0.8, RF)
        write_F_plot("h4", +7, +2,   +5.0, +0.2, RF)
      }

      if ((region <0) || (region == 9)) {
        write_F_plot("i1", +7, +1,   +7.0, -2.6, RF)
      }

      if ((region <0) || (region == 10)) {
        write_F_plot("j1", +6, -4,   +1.3, -0.2, RF)
        write_F_plot("j2", +7, -4,   +4.0, -2.8, RF)
        write_F_plot("j3", +8, -4,   +4.0, -1.2, RF)
      }

      if ((region <0) || (region == 11)) {
        write_F_plot("k1", -3, -4,   +0.4, -0.2, RF)
        write_F_plot("k2", -2, -4,   +0.9, -0.2, RF)
        write_F_plot("k3", -1, -4,   +2.0, -1.2, RF)
        write_F_plot("k4", 00, -4,   +2.0, -1.8, RF)
      }

      if ((region <0) || (region == 12)) {
        write_F_plot("m1", -4, -2,   -4.0, -0.2, RF)
        write_F_plot("m2", -3, -2,   -0.1, -0.2, RF)
        write_F_plot("m3", -2, -2,   +2.0, -2.2, RF)
      }

      if ((region <0) || (region == 13)) {
        write_F_plot("n1", +5, -1,   +2.0, -0.1, RF)
      }
    }
}

function write_segment(x0,y0,x1,y1,code,  k,r,x,y) {
  # Writes points of segment from {(x0,y0)} to {(x1,y1)} with given {code}.
  for (k = 0; k <= NP; k++) {
    r = (k + 0.0)/NP
    # Compute point {r} of the line segment:
    x = (1-r)*x0 + r*x1
    y = (1-r)*y0 + r*y1
    printf "%12.5f %12.5f %d @\n", x, y, code
  }
  printf "\n" # Blank line to break the curve in gnuplot.
}

function write_ratbound(x0,x1,n,cmp,code,  a,k,r,x,y) {
  # Writes points of the boundary between 1(n-1) and 1(n) 
  # from {x0} to {x1} with given {code}.
  # If {cmp} is true, draws the boundary between 2(n-1) and 2(n) instead.
  for (k = 0; k <= NP; k++) {
    r = (k + 0.0)/NP
    # Compute the {\kappa} parameter:
    x = (1-r)*x0 + r*x1
    if (x > 0.0) {
      # Compute the corresponding {\lambda} parameter::
      a = 1/x
      y = (1 - x)/(1 - exp(log(a)*n))
      # Complement if requested:
      if (cmp == 1) { y = 1 - x - y; }
      printf "%12.5f %12.5f %d @\n", x, y, code
    }
  }
  printf "\n" # Blank line to break the curve in gnuplot.
}

function write_parabola(x0,x1,A,B,sgn,code,   smin,smax,s,k,r,x,Q,R,y) {
  # Writes data to plot the parabola {(1 + x + y)^2 = A*x + B} between {x0} and {x1}.
  # If {sgn} is zero, plots both branches in ccw order.
  # If {sgn} is {+1} or {-1}, plots only that branch
  # Assumes that {x0} is where the two branches come together.
  if (sgn == 0) { smin = -1; smax = +1} else { smin = sgn; smax = sgn }
  
  for (s = smin; s <= smax; s = s+2) {
    # Plots the branch with sign {s}
    for (k = 0; k <= NP; k++) {
      r = (k + 0.0)/NP
      if (s < 0) { 
        # Plot the negative branch in reverse order:
        r = 1-r;
      }
      # Compute point {r} of the {x} range:
      x = (1-r)*x0 + r*x1
      # Solve {(1 + x + y)^2 = A*x + B} for {y}:
      Q = A*x + B;
      if (Q >= 0) {
        R = s*sqrt(Q);
        y = R - 1 - x
        printf "%12.5f %12.5f %d @\n", x, y, code
      }
    }
  }
  printf "\n" # Blank line to break the curve in gnuplot.
}

function write_funnybola(x0,x1,sgn,code,   smin,smax,s,k,r,x,Q,R,y) {
  # Writes data to plot the curve {x(x + y)^2 = x - y} between {x0} and {x1}.
  # If {sgn} is zero, plots both branches in ccw order.
  # If {sgn} is {+1} or {-1}, plots only that branch
  # Assumes that the branches never meet.
  if (sgn == 0) { smin = -1; smax = +1} else { smin = sgn; smax = sgn }
  
  for (s = smin; s <= smax; s = s+2) {
    # Plots the branch with sign {s}
    for (k = 0; k <= NP; k++) {
      r = (k + 0.0)/NP
      if (s < 0) { 
        # Plot the negative branch in reverse order:
        r = 1-r;
      }
      # Compute point {r} of the {x} range:
      x = (1-r)*x0 + r*x1
      # Solve {x(x + y)^2 = x - y} for {y}:
      Q = 1.0/(2.0*x);
      R = sqrt(2 + Q)
      y = - (x + Q) + s*R
      printf "%12.5f %12.5f %d @\n", x, y, code
    }
    printf "\n" # Blank line to break the curve in gnuplot.
  }
}

function write_grid_lines(M,  m,k,code) { 
  # Writes data to plot grid lines from {-M} to {+M}:
  m = int(M)+1
  for (k = -m; k <= +m; k++) {
    code = (k == 0 ? 8 : 9)
    write_segment(-M,  k,   +M, k,  code)
    write_segment( k, -M,   k, +M,  code)
  }
}

function write_label(xt,yt,codet,xp,yp,codep,lab) {
  # Writes data to plot a label {lab}.
  # The label will be at {(xt,yt)} with given plot {codet}.
  # If {(xp,yp)} is not {NIL} and different from {(xt,yt)},
  # writes a line segment from {(xp,yp)} to {(xt,yt)} with plot code {codep}.
  if ((xp != NIL) && (yp != NIL) && ((xp != xt) || (yp != yt))) {
    printf "%12.5f %12.5f %d @\n", xt, yt, codep
    printf "%12.5f %12.5f %d @\n", xp, yp, codep
    printf "\n" # Blank line to break the curve in gnuplot.
  }
  printf "%12.5f %12.5f %d %s\n", xt, yt, codet, lab
}

function write_fill_corner(x,y,code) {
  # Writes data to plot a corner of a filled curve with given {code} at {(x,y)}
  printf "%12.5f %12.5f %d @\n", x, y, code
}

function write_dense_F_plots(NP,M,RF,xlo,xhi,ylo,yhi, np,kx,ky,x,y,skip) {
  np = 10 # Number of F-plots along each axis
  for (kx = 0; kx < np; kx++) {
    x = xlo + (kx + 0.5)/np*(xhi-xlo)
    for (ky = 0; ky < np; ky++) {
      y = ylo + (ky + 0.5)/np*(yhi-ylo)
      lab = sprintf("%d.%d", kx, ky)
      # Exclude selected regions:
      skip = 0
      # if ((y <= 0) && ((x + y <= 1) || (Delta(x,y) <= 0))) { skip = 1; }
      # if ((y >= 1) && (x + y >= 1)) { skip = 1; }
      # if ((y >= 0) && (y <= 1) && (x + y >= 1)) { skip = 1; }
      if (skip == 0) { write_F_plot(lab, NIL, NIL, x, y, RF) }
    }
  }
}

function write_F_plot(lab,xt,yt,xc,yc,R,   n,k,x,y) {
  # Writes a miniature plot of {F} with half-side {R} for parameters
  # {\kappa=xc}, {\lambda=yc}, with label {lab}.
  # The plot codes 3,4,5,7,20,21 are used for the 
  # plot frame, the diagonals, the graph of F, the graph of {F^{\circ 2}}, 
  # the plot background, and the label, respectively.  
  
  # If {xt,yt} are NIL or equal to {xc,yc}, the plot is centered at
  # {(xc,yc)}. Otherwise it is centered at {(xt,yt)} and a line is drawn
  # between the two points, with plot code 16.
  
  if ((xc != NIL) && (xt != NIL) && ((xt != xc) || (yt != yc))) {
    printf "%12.5f %12.5f %d @\n", xt, yt, 16
    printf "%12.5f %12.5f %d @\n", xc, yc, 16
    printf "\n" # Blank line to break the curve in gnuplot.
  } else {
    xt = xc; yt = yc
  }
  
  # Paint the plot background:
  printf "%12.5f %12.5f %d @\n", xt-R, yt-R, 20
  printf "%12.5f %12.5f %d @\n", xt+R, yt-R, 20
  printf "%12.5f %12.5f %d @\n", xt+R, yt+R, 20
  printf "%12.5f %12.5f %d @\n", xt-R, yt+R, 20
  printf "%12.5f %12.5f %d @\n", xt-R, yt-R, 20
  printf "\n"  # Blank line to break the curve in gnuplot.

  # Draw the frame:
  printf "%12.5f %12.5f %d @\n", xt-R, yt-R, 3
  printf "%12.5f %12.5f %d @\n", xt+R, yt-R, 3
  printf "%12.5f %12.5f %d @\n", xt+R, yt+R, 3
  printf "%12.5f %12.5f %d @\n", xt-R, yt+R, 3
  printf "%12.5f %12.5f %d @\n", xt-R, yt-R, 3
  printf "\n"  # Blank line to break the curve in gnuplot.
  
  # Draw the diagonals:
  printf "%12.5f %12.5f %d @\n", xt-R, yt-R, 4
  printf "%12.5f %12.5f %d @\n", xt+R, yt+R, 4
  printf "\n"  # Blank line to break the curve in gnuplot.
  
  printf "%12.5f %12.5f %d @\n", xt-R, yt+R, 4
  printf "%12.5f %12.5f %d @\n", xt+R, yt-R, 4
  printf "\n"  # Blank line to break the curve in gnuplot.
  
  # Plot the graph of {F^{\circ 2}}:
  for (k = 0; k <= NP; k++) {
    r = (k + 0.0)/NP
    s = F(r, xc, yc, dead, 2)
    x = xt + (2*r-1)*R
    y = yt + (2*s-1)*R
    printf "%12.5f %12.5f %d @\n", x, y, 7
  }
  printf "\n"  # Blank line to break the curve in gnuplot.
  
  # Plot the graph of {F^{\circ 8}}:
  for (k = 0; k <= NP; k++) {
    r = (k + 0.0)/NP
    s = F(r, xc, yc, dead, 8)
    x = xt + (2*r-1)*R
    y = yt + (2*s-1)*R
    printf "%12.5f %12.5f %d @\n", x, y, 22
  }
  printf "\n"  # Blank line to break the curve in gnuplot.

  # Plot the graph of {F}:
  for (k = 0; k <= NP; k++) {
    r = (k + 0.0)/NP
    s = F(r, xc, yc, dead, 1)
    x = xt + (2*r-1)*R
    y = yt + (2*s-1)*R
    printf "%12.5f %12.5f %d @\n", x, y, 5
  }
  
  # Print the label:
  x = xt;
  y = yt + 0.7*R;
  printf "%12.5f %12.5f %d %s\n", x, y, 21, lab
  
  printf "\n"  # Blank line to break the curve in gnuplot.
}

function Delta(x,y, B) {
  # Discriminant of equation {FQ(z) = z} for quadratic part.
  B = 1 - x + y
  return B*B + 4*x*y
}