text:
using SparseArrays
using LinearAlgebra

# DIAGONAL_SBP_D1 creates a diagonal norm SBP operator
# (D, Hinv, H, r) = diagonal_sbp_D1(p, N; xc = (-1,1))
#
# inputs:
#   p: sbp interior accuracy
#   N: finite difference grid size is N+1
#   xc: (keyword) grid span [default: (-1, 1)]
#
# outputs:
#   D:   difference operator Hinv*(-M+BS)
#   HI:  inverse of the SBP norm
#   H:   the SBP norm
#   r:   grid from xc[1] to xc[2]
#
#   References:
#   Operators for order 2, 4, 6 are from
#   @book{gustafsson2008high,
#     title={High order difference methods for time dependent PDE},
#     author={Gustafsson, Bertil},
#     year={2008},
#     publisher={Springer},
#     series={Springer Series in Computational Mathematics},
#     volume={38}
#   }
#
#   Order 8 operator is from
#   @article{MattssonSvardShoeybe2008JCP,
#     title={Stable and accurate schemes for the compressible Navier--Stokes equations},
#     author={Mattsson, Ken and Sv{\"a}rd, Magnus and Shoeybi, Mohammad},
#     journal={Journal of Computational Physics},
#     volume={227},
#     number={4},
#     pages={2293--2316},
#     year={2008},
#     doi={10.1016/j.jcp.2007.10.018}
#   }
#   with
#      x1 =  663/950
#      x2 = -229/1900
#      x3 =  415/547
#
#   Order 10 operator is from
#   @Article{MattssonAlmquist2013JCP,
#     author = {K. Mattsson and M. Almquist},
#     title = {A solution to the stability issues with block norm summation by parts operators},
#     journal = {Journal of Computational Physics},
#     volume = {253},
#     pages = {418--442},
#     year = {2013},
#     doi = {10.1016/j.jcp.2013.07.013}
#   }
#{{{
function diagonal_sbp_D1(p, N; xc = (-1, 1))

  if p == 2
    bhinv = [2]
    d  = [-1/2 0 1/2]
    bd  = [-1 1]
  elseif p == 4
    bhinv = [48/17 48/59 48/43 48/49]
    d  = [1/12 -2/3 0 2/3 -1/12]
    bd = [-24/17  59/34  -4/17  -3/34  0     0;
           -1/2    0      1/2    0     0     0;
            4/43 -59/86   0     59/86 -4/43  0;
            3/98   0    -59/98   0    32/49 -4/49]
  elseif p == 6
    bhinv = [43200/13649 8640/12013 4320/2711 4320/5359 8640/7877 43200/43801];

    d = [-1/60 3/20 -3/4 0 3/4 -3/20 1/60];

    bd = [ -21600/13649   104009/54596   30443/81894  -33311/27298    16863/27298  -15025/163788        0           0         0;
          -104009/240260       0          -311/72078   20229/24026   -24337/48052   36661/360390        0           0         0;
           -30443/162660     311/32532       0        -11155/16266    41287/32532  -21999/54220         0           0         0;
            33311/107180  -20229/21436     485/1398        0           4147/21436    25427/321540      72/5359      0         0;
           -16863/78770    24337/31508  -41287/47262   -4147/15754        0         342523/472620   -1296/7877    144/7877    0;
            15025/525612  -36661/262806  21999/87602  -25427/262806 -342523/525612       0          32400/43801 -6480/43801 720/43801];
  elseif p == 8
    bhinv = [5080320/1498139 725760/1107307 80640/20761 725760/1304999 725760/299527 80640/103097 725760/670091 5080320/5127739];
    d = [1/280 -4/105 1/5 -4/5 0 4/5 -1/5 4/105 -1/280];

    bd = [    -2540160/1498139        699846290793/311403172540  -10531586157/311403172540  -145951651079/186841903524    398124597/15570158627    39152001/113858564  -80631892229/934209517620  -6230212503/311403172540        0                 0               0               0;
          -24132630717/55557028660               0                 2113176981/23016483302      5686186719/11508241651   -3408473341/138098899812  -39291999/210388330     607046586/11508241651    3460467023/483346149342        0                 0               0               0;
            3510528719/90623010660      -704392327/1294614438               0                   503511235/1294614438      354781619/2589228876       407439/3944590        -2986842/16597621        169381493/3020767022          0                 0               0               0;
          145951651079/1139279786988   -5686186719/13562854607    -1510533705/27125709214               0                6763379967/54251418428    13948923/49589962    -1603900430/40688563821   -3742312557/189879964498        0                 0               0               0;
            -398124597/21790888777      3408473341/37355809332    -1064344857/12451936444     -6763379967/12451936444             0                  763665/1198108     -1282435899/12451936444    7822226819/261490665324    -2592/299527          0               0               0;
              -1864381/23506116           13097333/58765290           -407439/19588430           -4649641/11753058          -254555/1237164               0                 5346432/9794215           -923328/9794215          3072/103097       -288/103097        0               0;
           11518841747/417855345780     -607046586/6964255763         1192698/23768791         1603900430/20892767289    1282435899/27857023052   -48117888/63658645              0                 301190400/366539777     -145152/670091      27648/670091    -2592/670091        0;
            6230212503/1065851828540   -3460467023/319755548562   -1524433437/106585182854     3742312557/106585182854  -7822226819/639511097124   58169664/487135205   -2108332800/2804873233              0               4064256/5127739  -1016064/5127739  193536/5127739  -18144/5127739];
  elseif p == 10
    bhinv = [18289152000/5261271563 1828915200/2881040311 406425600/52175551 6096384/11662993 87091200/50124587 72576000/50124587 87091200/148333439 152409600/63867949 16257024/20608675 1828915200/1704508063 18289152000/18425967263];

    d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260];

    bd = [  -1.7380923775745425e+00   2.3557601935237220e+00  -1.5328406598563976e-01  -5.7266565770416333e-01  -1.8308103515008173e-01   1.8186748267946842e-01   2.0034232582598244e-01   2.2678007363666621e-02  -1.1782459320459637e-01  -3.0591175636402144e-02   3.4890895862586133e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
            -4.3020203737210871e-01   0.0000000000000000e+00   1.1837297346927406e-01   3.3928601158526644e-01   1.3241927733034406e-01  -8.7495003780608913e-02  -1.1750484124279399e-01  -1.6401912273575153e-02   6.2537843443041474e-02   1.7143274696828435e-02  -1.8155585855667674e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
             3.4348531361887280e-01  -1.4525207124434036e+00   0.0000000000000000e+00   2.9011513992277767e+00  -2.2419288742360557e+00  -5.4662873578741478e-01   1.2908050607446131e+00   6.1514504292452719e-02  -4.2442625460011202e-01   1.5579158905288801e-02   5.2969140277981920e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
             8.6111387816878188e-02  -2.7937273515056432e-01  -1.9467880944770807e-01   0.0000000000000000e+00   2.0170150914578375e-01   2.4269917331475005e-01  -7.7261988327590472e-02   5.0649247607525059e-02  -7.4775049946661561e-03  -4.0978487203372188e-02   1.8608207238964152e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
             9.1509035082611684e-02  -3.6243526359648576e-01   5.0007055839856984e-01  -6.7045605191055857e-01   0.0000000000000000e+00  -1.7807807859119628e-02   7.5000761407401195e-01  -2.2979723229714316e-01  -1.2521154324370892e-01   6.8278284106004450e-02  -4.1575927541817690e-03   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
            -7.5752056274147259e-02   1.9956355926115746e-01   1.0160630736447970e-01  -6.7227694623145351e-01   1.4839839882599690e-02   0.0000000000000000e+00   5.4091068834671807e-01  -1.2712520372174399e-01  -8.9292453564020990e-02   1.6181541970619609e-01  -5.4289154769785249e-02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
            -3.3838029883391296e-02   1.0867927550524317e-01  -9.7293058702223670e-02   8.6783825404790446e-02  -2.5344131542932297e-01  -2.1934035945002228e-01   0.0000000000000000e+00   2.7184438867288430e-01   1.9102691945078512e-01  -4.8646826827046824e-02  -6.2407959378425991e-03   4.6597719614658163e-04   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
            -1.5567948806367624e-02   6.1656604470023607e-02  -1.8844858059892756e-02  -2.3122780265804038e-01   3.1560994521078772e-01   2.0951677187991255e-01  -1.1048784865195491e+00   0.0000000000000000e+00   1.1823059621092409e+00  -5.3610400867086083e-01   1.5931375952374752e-01  -2.3673846172827626e-02   1.8939076938262100e-03   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00;
             2.6737701764454301e-02  -7.7712278574126673e-02   4.2981266272823705e-02   1.1284579710276557e-02   5.6847566375570611e-02   4.8647834370398067e-02  -2.5665536068472994e-01  -3.9083324869946684e-01   0.0000000000000000e+00   6.5716944195909766e-01  -1.5822272208022428e-01   4.6954983762905661e-02  -7.8258306271509429e-03   6.2606645017207550e-04   0.0000000000000000e+00   0.0000000000000000e+00;
             9.4425181052687698e-03  -2.8976375375532045e-02  -2.1459742428921558e-03   8.4117843695442701e-02  -4.2165149106440383e-02  -1.1991463562335723e-01   8.8902467992349743e-02   2.4105392677971343e-01  -8.9388344421253152e-01   0.0000000000000000e+00   8.6496680152924643e-01  -2.5547312415382800e-01   6.3868281038457000e-02  -1.0644713506409501e-02   8.5157708051276015e-04   0.0000000000000000e+00;
            -9.9625965676187218e-03   2.8387641187789508e-02  -6.7495090936003027e-03  -3.5335033597892078e-02   2.3750992019053968e-03   3.7216380474824604e-02   1.0550378667904333e-02  -6.6265458456725809e-02   1.9908619649258188e-01  -8.0014409359906680e-01   0.0000000000000000e+00   8.2714572225493910e-01  -2.3632734921569687e-01   5.9081837303924217e-02  -9.8469728839873684e-03   7.8775783071898962e-04];
  else
    error(string("Operators for order ", p, " are not implemented"))
  end

  (bm, bn) = size(bd);
  Np = N+1

  if(Np < 2*bm || Np < bn)
    error("Grid not big enough to support the operator. Grid must have N >= ", max(bn,2*bm))
  end

  h = (xc[2] - xc[1]) / N
  @assert h > 0
  H_I = 1:Np
  H_V = ones(Np)
  H_V[1:bm] = 1 ./ bhinv[:]
  H_V[Np-bm+1:Np] = 1 ./ bhinv[end:-1:1]
  H = sparse(H_I, H_I, h * H_V)
  HI = sparse(H_I, H_I, 1 ./ (h * H_V))

  n = floor(Int64, length(d)/2);
  B_I1 = (bm+1:Np-bm) * ones(1,p+1)
  B_J1 = ones(Np-2bm,1) * (-div(p,2):div(p,2))' + B_I1
  B_V1 = ones(Np-2bm)
  B_V1 = kron(B_V1, d)

  B_I2 = (1:bm) * ones(1, bn)
  B_J2 = ones(bm) * (1:bn)'
  B_V2 = bd
  B_I3 = (Np+1) .- B_I2
  B_J3 = (Np+1) .- B_J2
  B_V3 = -B_V2
  D = sparse([B_I1[:];B_I2[:];B_I3[:]],
             [B_J1[:];B_J2[:];B_J3[:]],
             [B_V1[:];B_V2[:];B_V3[:]]/h, N+1, N+1)

  r = range(xc[1], stop=xc[2], length=N+1)

  (D, HI, H, r)
end
#}}}

#DIAGONAL_SBP_D2 creates a diagonal norm SBP operator for the 2nd drerivative
# (D, S0, SN, Hinv, H, r) = diagonal_sbp_D2(p, N; xc = (-1,1))
#
# inputs:
#   p: sbp interior accuracy
#   N: finite difference grid size is N+1
#   xc: (keyword) grid span [default: (-1, 1)]
#
# outputs:
#   D:   difference operator Hinv*(-M-S0+SN)
#   S0:  boundary derivative operator at first grid point
#   SN:  boundary derivative operator at last grid point
#   HI:  inverse of the SBP norm
#   H:   the SBP norm
#   r:   grid from xc[1] to xc[2]
#
#   References:
#   Operators for order 2, 4, 6, 8 are from
#   @book{gustafsson2008high,
#     title={High order difference methods for time dependent PDE},
#     author={Gustafsson, Bertil},
#     year={2008},
#     publisher={Springer},
#     series={Springer Series in Computational Mathematics},
#     volume={38}
#   }
#
#   Order 10 operator is from
#   @Article{MattssonAlmquist2013JCP,
#     author = {K. Mattsson and M. Almquist},
#     title = {A solution to the stability issues with block norm summation by
#     parts operators},
#     journal = {Journal of Computational Physics},
#     volume = {253},
#     pages = {418--442},
#     year = {2013},
#     doi = {10.1016/j.jcp.2013.07.013}
#   }
#{{{
function diagonal_sbp_D2(p, N; xc = (-1, 1))

  if p == 2
    bhinv = [2];

    d  = [1 -2 1];
    bd  = d;

    BS = [3/2 -2 1/2];
  elseif p == 4
    bhinv = [48/17 48/59 48/43 48/49];

    d  = [-1/12 4/3 -5/2 4/3 -1/12];

    bd = [ 2    -5       4     -1       0      0;
           1    -2       1      0       0      0;
          -4/43 59/43 -110/43  59/43   -4/43   0;
          -1/49  0      59/49 -118/49  64/49  -4/49];

    BS = [11/6 -3 3/2 -1/3];

  elseif p == 6
    bhinv = [43200/13649 8640/12013 4320/2711 4320/5359 8640/7877 43200/43801];

    d = [1/90 -3/20 3/2 -49/18 3/2 -3/20 1/90];

    bd = [  114170/40947   -438107/54596   336409/40947  -276997/81894     3747/13649     21035/163788      0           0         0
              6173/5860      -2066/879       3283/1758      -303/293       2111/3516       -601/4395        0           0         0
            -52391/81330    134603/32532   -21982/2711    112915/16266   -46969/16266     30409/54220       0           0         0
             68603/321540   -12423/10718   112915/32154   -75934/16077    53369/21436    -54899/160770     48/5359      0         0
             -7053/39385     86551/94524   -46969/23631    53369/15754   -87904/23631    820271/472620  -1296/7877     96/7877    0
             21035/525612   -24641/131403   30409/87602   -54899/131403  820271/525612  -117600/43801   64800/43801 -6480/43801 480/43801];

    BS = [25/12 -4 3 -4/3 1/4];

  elseif p == 8
    bhinv = [5080320/1498139 725760/1107307 80640/20761 725760/1304999 725760/299527 80640/103097 725760/670091 5080320/5127739];
    d = [-1/560 8/315 -1/5 8/5 -205/72 8/5 -1/5 8/315 -1/560]

    bd = zeros(8,12);
    bd[1, 1] =  4870382994799/1358976868290;
    bd[1, 2] =  -893640087518/75498714905  ;
    bd[1, 3] =   926594825119/60398971924  ;
    bd[1, 4] = -1315109406200/135897686829 ;
    bd[1, 5] =    39126983272/15099742981  ;
    bd[1, 6] =    12344491342/75498714905  ;
    bd[1, 7] =  -451560522577/2717953736580;

    bd[2, 1] =  333806012194/390619153855;
    bd[2, 2] = -154646272029/111605472530;
    bd[2, 3] =    1168338040/33481641759 ;
    bd[2, 4] =   82699112501/133926567036;
    bd[2, 5] =    -171562838/11160547253 ;
    bd[2, 6] =  -28244698346/167408208795;
    bd[2, 7] =   11904122576/167408208795;
    bd[2, 8] =   -2598164715/312495323084;

    bd[3, 1] =   7838984095/52731029988;
    bd[3, 2] =   1168338040/5649753213 ;
    bd[3, 3] =    -88747895/144865467  ;
    bd[3, 4] =    423587231/627750357  ;
    bd[3, 5] = -43205598281/22599012852;
    bd[3, 6] =   4876378562/1883251071 ;
    bd[3, 7] =  -5124426509/3766502142 ;
    bd[3, 8] =  10496900965/39548272491;

    bd[4, 1] =  -94978241528/828644350023;
    bd[4, 2] =   82699112501/157837019052;
    bd[4, 3] =    1270761693/13153084921 ;
    bd[4, 4] = -167389605005/118377764289;
    bd[4, 5] =   48242560214/39459254763 ;
    bd[4, 6] =  -31673996013/52612339684 ;
    bd[4, 7] =   43556319241/118377764289;
    bd[4, 8] =  -44430275135/552429566682;

    bd[5, 1] =   1455067816/21132528431;
    bd[5, 2] =   -171562838/3018932633 ;
    bd[5, 3] = -43205598281/36227191596;
    bd[5, 4] =  48242560214/9056797899 ;
    bd[5, 5] = -52276055645/6037865266 ;
    bd[5, 6] =  57521587238/9056797899 ;
    bd[5, 7] = -80321706377/36227191596;
    bd[5, 8] =   8078087158/21132528431;
    bd[5, 9] =        -1296/299527     ;

    bd[6, 1] =   10881504334/327321118845;
    bd[6, 2] =  -28244698346/140280479505;
    bd[6, 3] =    4876378562/9352031967  ;
    bd[6, 4] =  -10557998671/12469375956 ;
    bd[6, 5] =   57521587238/28056095901 ;
    bd[6, 6] = -278531401019/93520319670 ;
    bd[6, 7] =   73790130002/46760159835 ;
    bd[6, 8] = -137529995233/785570685228;
    bd[6, 9] =          2048/103097      ;
    bd[6,10] =          -144/103097      ;

    bd[7, 1] = -135555328849/8509847458140;
    bd[7, 2] =   11904122576/101307707835 ;
    bd[7, 3] =   -5124426509/13507694378  ;
    bd[7, 4] =   43556319241/60784624701  ;
    bd[7, 5] =  -80321706377/81046166268  ;
    bd[7, 6] =   73790130002/33769235945  ;
    bd[7, 7] = -950494905688/303923123505 ;
    bd[7, 8] =  239073018673/141830790969 ;
    bd[7, 9] =       -145152/670091       ;
    bd[7,10] =         18432/670091       ;
    bd[7,11] =        -1296/670091        ;

    bd[8, 1] =             0             ;
    bd[8, 2] =   -2598164715/206729925524;
    bd[8, 3] =   10496900965/155047444143;
    bd[8, 4] =  -44430275135/310094888286;
    bd[8, 5] =     425162482/2720130599  ;
    bd[8, 6] = -137529995233/620189776572;
    bd[8, 7] =  239073018673/155047444143;
    bd[8, 8] = -144648000000/51682481381 ;
    bd[8, 9] =       8128512/5127739     ;
    bd[8,10] =      -1016064/5127739     ;
    bd[8,11] =        129024/5127739     ;
    bd[8,12] =         -9072/5127739     ;

    BS = [4723/2100 -839/175 157/35 -278/105 103/140 1/175 -6/175];

  elseif p == 10
    bhinv = [18289152000/5261271563 1828915200/2881040311 406425600/52175551 6096384/11662993 87091200/50124587 72576000/50124587 87091200/148333439 152409600/63867949 16257024/20608675 1828915200/1704508063 18289152000/18425967263];

    M = zeros(11,16);

    M[ 1, 1] =  1.2056593789671863908;
    M[ 1, 2] = -1.3378814169347239658;
    M[ 1, 3] =  0.0036847309286546532061;
    M[ 1, 4] =  0.15698288365600946515;
    M[ 1, 5] = -0.0037472461482539197952;
    M[ 1, 6] = -0.0062491712449361657064;
    M[ 1, 7] = -0.029164045872729581661;
    M[ 1, 8] =  0.00054848184117832929161;
    M[ 1, 9] =  0.013613461413384884448;
    M[ 1,10] = -0.0025059220258337808220;
    M[ 1,11] = -0.00094113457993630916498;

    M[ 2, 2] =  2.1749807117105597139;
    M[ 2, 3] = -0.12369059547124894597;
    M[ 2, 4] = -0.83712574037924152603;
    M[ 2, 5] =  0.050065127254670973258;
    M[ 2, 6] =  0.0081045853127317536361;
    M[ 2, 7] =  0.097405846039248226536;
    M[ 2, 8] = -0.00068942461520402214720;
    M[ 2, 9] = -0.041326971493379188475;
    M[ 2,10] =  0.0075778529605774119402;
    M[ 2,11] =  0.0025800256160095691057;

    M[ 3, 3] =  0.18361596652499065332;
    M[ 3, 4] =  0.048289690013342693109;
    M[ 3, 5] = -0.19719621435164680412;
    M[ 3, 6] =  0.11406859029505842791;
    M[ 3, 7] = -0.029646295985488126964;
    M[ 3, 8] = -0.0016038463172861201306;
    M[ 3, 9] =  0.0032879841528337653050;
    M[ 3,10] = -0.00093242311589807387463;
    M[ 3,11] =  0.00012241332668787820533;

    M[ 4, 4] =  1.2886524606662484673;
    M[ 4, 5] = -0.14403037739488789185;
    M[ 4, 6] = -0.44846291607489015475;
    M[ 4, 7] = -0.10598334599408054277;
    M[ 4, 8] = -0.015873275740355918053;
    M[ 4, 9] =  0.073988493386459608166;
    M[ 4,10] = -0.012508848749152899785;
    M[ 4,11] = -0.0039290233894513005339;

    M[ 5, 5] =  0.51482665719685186210;
    M[ 5, 6] =  0.051199577887125103015;
    M[ 5, 7] = -0.36233561810883077365;
    M[ 5, 8] =  0.091356850268746392169;
    M[ 5, 9] =  0.0024195916108052419451;
    M[ 5,10] = -0.0018564214413731389338;
    M[ 5,11] = -0.00070192677320704413827;

    M[ 6, 6] =  0.68636003380365860083;
    M[ 6, 7] = -0.28358848290867614908;
    M[ 6, 8] = -0.13836006478253396528;
    M[ 6, 9] =  0.0076158070663111995297;
    M[ 6,10] =  0.011447010307180005164;
    M[ 6,11] = -0.0021349696610286552676;

    M[ 7, 7] =  1.5216081480839085990;
    M[ 7, 8] = -0.42653865162216293237;
    M[ 7, 9] = -0.42047484981879143123;
    M[ 7,10] =  0.019813359263872926304;
    M[ 7,11] =  0.019221397241190103344;

    M[ 8, 8] =  1.0656733504627815335;
    M[ 8, 9] = -0.66921872668484232217;
    M[ 8,10] =  0.12022033144141336599;
    M[ 8,11] = -0.030157881394591483631;

    M[ 9, 9] =  2.4064247712949611684;
    M[ 9,10] = -1.5150200315922263367;
    M[ 9,11] =  0.17373015320416595052;

    M[10,10] =  2.7682502485427255096;
    M[10,11] = -1.5975407111468405444;

    M[11,11] =  2.9033627686681129471;
    d     = [1/3150  -5/1008  5/126  -5/21  5/3 -5269/1800 5/3 -5/21 5/126 -5/1008 1/3150]
    for k = 1:5
      M[11-5+k,11 .+ (1:k)] = -d[k:-1:1];
    end

    M[1:11,1:11] = M[1:11,1:11]'+M[1:11,1:11]-diagm(0 => diag(M[1:11,1:11]));

    BS = zeros(1,16);

    BS[1:7] = -[-49/20 6 -15/2 20/3 -15/4 6/5 -1/6];
    e0 = zeros(11,1);
    e0[1] = 1;

    bd = diagm(0 => bhinv[:]) * (-M+e0*BS);

  else
    error(string("Operators for order ", p, " are not implemented"))
  end

  (bm, bn) = size(bd);

  M = N+1;
  if(M < 2*bm || M < bn)
    error("Grid not big enough to support the operator. Grid must have N >= ", max(bn,2*bm))
  end

  h = (xc[2] - xc[1]) / N
  @assert h > 0
  H_I = 1:M
  H_V = ones(M)
  H_V[1:bm] = bhinv[:]
  H_V[M-bm+1:M] = bhinv[end:-1:1]
  HI = sparse(H_I, H_I, H_V / h)
  H  = sparse(H_I, H_I, h ./ H_V)

  n = floor(Int64, length(d)/2);
  B_I1 = (bm+1:M-bm) * ones(1,p+1)
  B_J1 = ones(M-2bm,1) * (-div(p,2):div(p,2))' + B_I1
  B_V1 = ones(M-2bm)
  B_V1 = kron(B_V1, d)

  B_I2 = (1:bm) * ones(1, bn)
  B_J2 = ones(bm) * (1:bn)'
  B_V2 = bd
  B_I3 = (M+1) .- B_I2
  B_J3 = (M+1) .- B_J2
  B_V3 = B_V2
  D = sparse([B_I1[:];B_I2[:];B_I3[:]],
             [B_J1[:];B_J2[:];B_J3[:]],
             [B_V1[:];B_V2[:];B_V3[:]]/h^2, N+1, N+1)

  S0 = sparse(ones(length(BS)), 1:length(BS), -BS[:]/h, N+1, N+1);
  SN = sparse((N+1) * ones(length(BS)), (N+1):-1:(N+2-length(BS)),
              BS[:]/h, N+1, N+1);

  r = range(xc[1], stop=xc[2], length=N+1)

  (D, S0, SN, HI, H, r)

end
#}}}

#VARIABLE_DIAGONAL_SBP_D2 creates a diagonal norm SBP operator for the 2nd
# derivative with variable coefficients, e.g., approximates ∂/∂r( b(r) ∂u/∂r)
#
# (D, BS, Hinv, H, r) = variable_diagonal_sbp_D2(p, N, B; xc = (-1,1))
#{{{
function variable_diagonal_sbp_D2(p, N, B; xc = (-1, 1))
  r = range(xc[1], stop=xc[2], length=N+1)
  variable_diagonal_sbp_D2(p, N, B(collect(r));xc=xc)
end
function variable_diagonal_sbp_D2(p, N, B::T; xc = (-1, 1)) where T <: Number
  variable_diagonal_sbp_D2(p, N, B*ones(N+1);xc=xc)
end

function variable_diagonal_sbp_D2(p, N, B::Array{Float64,1}; xc = (-1, 1))
  @assert length(B) == N+1

  if p == 2
    bhinv = [2];

    I_M0 = [1]
    J_M0 = [1]
    V_M0 = [(B[1]+B[2])/2]

    I_MN = [N+1]
    J_MN = [N+1]
    V_MN = [(B[N]+B[N+1])/2]

    I_M = [2:N+1;
           2:N  ;
           1:N  ]
    J_M = [1:N  ;
           2:N  ;
           2:N+1]
    V_M = [-(B[1:N  ]+B[2:N+1])/2;
            (B[1:N-1]+2B[2:N]+B[3:N+1])/2;
           -(B[1:N  ]+B[2:N+1])/2];
    M = sparse([I_M0;I_M;I_MN], [J_M0;J_M;J_MN], [V_M0;V_M;V_MN])

    BS = [3/2 -2 1/2];
  elseif p == 4
    bhinv = [48/17 48/59 48/43 48/49];

    BS = [11/6 -3 3/2 -1/3];

    V_M0 = zeros(6,6)

    (b1,b2,b3,b4,b5,b6,b7,b8) = B[1:8]
    V_M0[1, 1] =                (12/17)b1 + (59/192)b2 + (27010400129/345067064608)b3 + (69462376031/2070402387648)b4
    V_M0[1, 2] = V_M0[2, 1] = - (59/68)b1 - (6025413881/21126554976)b3 - (537416663/7042184992)b4
    V_M0[1, 3] = V_M0[3, 1] =   (2/17)b1 - (59/192)b2 + (213318005/16049630912)b4 + (2083938599/8024815456)b3
    V_M0[1, 4] = V_M0[4, 1] =   (3/68)b1 - (1244724001/21126554976)b3 + (752806667/21126554976)b4
    V_M0[1, 5] = V_M0[5, 1] =   (49579087/10149031312)b3 - (49579087/10149031312)b4
    V_M0[1, 6] = V_M0[6, 1] = - (1/784)b4 + (1/784)b3
    V_M0[2, 2] =                (3481/3264)b1 + (9258282831623875/7669235228057664)b3 + (236024329996203/1278205871342944)b4
    V_M0[2, 3] = V_M0[3, 2] = - (59/408)b1 - (29294615794607/29725717938208)b3 - (2944673881023/29725717938208)b4
    V_M0[2, 4] = V_M0[4, 2] = - (59/1088)b1 + (260297319232891/2556411742685888)b3 - (60834186813841/1278205871342944)b4
    V_M0[2, 5] = V_M0[5, 2] = - (1328188692663/37594290333616)b3 + (1328188692663/37594290333616)b4
    V_M0[2, 6] = V_M0[6, 2] = - (8673/2904112)b3 + (8673/2904112)b4
    V_M0[3, 3] =                (1/51)b1 + (59/192)b2 + (13777050223300597/26218083221499456)b4 + (564461/13384296)b5 + (378288882302546512209/270764341349677687456)b3
    V_M0[3, 4] = V_M0[4, 3] =   (1/136)b1 - (125059/743572)b5 - (4836340090442187227/5525802884687299744)b3 - (17220493277981/89177153814624)b4
    V_M0[3, 5] = V_M0[5, 3] = - (10532412077335/42840005263888)b4 + (1613976761032884305/7963657098519931984)b3 + (564461/4461432)b5
    V_M0[3, 6] = V_M0[6, 3] = - (960119/1280713392)b4 - (3391/6692148)b5 + (33235054191/26452850508784)b3
    V_M0[4, 4] =                (3/1088)b1 + (507284006600757858213/475219048083107777984)b3 + (1869103/2230716)b5 + (1/24)b6 + (1950062198436997/3834617614028832)b4
    V_M0[4, 5] = V_M0[5, 4] = - (4959271814984644613/20965546238960637264)b3 - (1/6)b6 - (15998714909649/37594290333616)b4 - (375177/743572)b5
    V_M0[4, 6] = V_M0[6, 4] = - (368395/2230716)b5 + (752806667/539854092016)b3 + (1063649/8712336)b4 + (1/8)b6
    V_M0[5, 5] =                (8386761355510099813/128413970713633903242)b3 + (2224717261773437/2763180339520776)b4 + (5/6)b6 + (1/24)b7 + (280535/371786)b5
    V_M0[5, 6] = V_M0[6, 5] = - (35039615/213452232)b4 - (1/6)b7 - (13091810925/13226425254392)b3 - (1118749/2230716)b5 - (1/2)b6
    V_M0[6, 6] =                (3290636/80044587)b4 + (5580181/6692148)b5 + (5/6)b7 + (1/24)b8 + (660204843/13226425254392)b3 + (3/4)b6
    I_M0 = (1:6) * ones(1,6)
    J_M0 = I_M0'

    V_MN = zeros(6,6)
    (b1,b2,b3,b4,b5,b6,b7,b8) = B[N+1:-1:N-6]
    V_MN[6, 6] =                (12/17)b1 + (59/192)b2 + (27010400129/345067064608)b3 + (69462376031/2070402387648)b4
    V_MN[6, 5] = V_MN[5, 6] = - (59/68)b1 - (6025413881/21126554976)b3 - (537416663/7042184992)b4
    V_MN[6, 4] = V_MN[4, 6] =   (2/17)b1 - (59/192)b2 + (213318005/16049630912)b4 + (2083938599/8024815456)b3
    V_MN[6, 3] = V_MN[3, 6] =   (3/68)b1 - (1244724001/21126554976)b3 + (752806667/21126554976)b4
    V_MN[6, 2] = V_MN[2, 6] =   (49579087/10149031312)b3 - (49579087/10149031312)b4
    V_MN[6, 1] = V_MN[1, 6] = - (1/784)b4 + (1/784)b3
    V_MN[5, 5] =                (3481/3264)b1 + (9258282831623875/7669235228057664)b3 + (236024329996203/1278205871342944)b4
    V_MN[5, 4] = V_MN[4, 5] = - (59/408)b1 - (29294615794607/29725717938208)b3 - (2944673881023/29725717938208)b4
    V_MN[5, 3] = V_MN[3, 5] = - (59/1088)b1 + (260297319232891/2556411742685888)b3 - (60834186813841/1278205871342944)b4
    V_MN[5, 2] = V_MN[2, 5] = - (1328188692663/37594290333616)b3 + (1328188692663/37594290333616)b4
    V_MN[5, 1] = V_MN[1, 5] = - (8673/2904112)b3 + (8673/2904112)b4
    V_MN[4, 4] =                (1/51)b1 + (59/192)b2 + (13777050223300597/26218083221499456)b4 + (564461/13384296)b5 + (378288882302546512209/270764341349677687456)b3
    V_MN[4, 3] = V_MN[3, 4] =   (1/136)b1 - (125059/743572)b5 - (4836340090442187227/5525802884687299744)b3 - (17220493277981/89177153814624)b4
    V_MN[4, 2] = V_MN[2, 4] = - (10532412077335/42840005263888)b4 + (1613976761032884305/7963657098519931984)b3 + (564461/4461432)b5
    V_MN[4, 1] = V_MN[1, 4] = - (960119/1280713392)b4 - (3391/6692148)b5 + (33235054191/26452850508784)b3
    V_MN[3, 3] =                (3/1088)b1 + (507284006600757858213/475219048083107777984)b3 + (1869103/2230716)b5 + (1/24)b6 + (1950062198436997/3834617614028832)b4
    V_MN[3, 2] = V_MN[2, 3] = - (4959271814984644613/20965546238960637264)b3 - (1/6)b6 - (15998714909649/37594290333616)b4 - (375177/743572)b5
    V_MN[3, 1] = V_MN[1, 3] = - (368395/2230716)b5 + (752806667/539854092016)b3 + (1063649/8712336)b4 + (1/8)b6
    V_MN[2, 2] =                (8386761355510099813/128413970713633903242)b3 + (2224717261773437/2763180339520776)b4 + (5/6)b6 + (1/24)b7 + (280535/371786)b5
    V_MN[2, 1] = V_MN[1, 2] = - (35039615/213452232)b4 - (1/6)b7 - (13091810925/13226425254392)b3 - (1118749/2230716)b5 - (1/2)b6
    V_MN[1, 1] =                (3290636/80044587)b4 + (5580181/6692148)b5 + (5/6)b7 + (1/24)b8 + (660204843/13226425254392)b3 + (3/4)b6

    I_MN = (N-4:N+1) * ones(1,6)
    J_MN = I_MN'

    I_M = [7:N-3;
           7:N-4;
           7:N-5;
           6:N-5;
           5:N-5]
    J_M = [5:N-5;
           6:N-5;
           7:N-5;
           7:N-4;
           7:N-3]
    V_M = [(+(1/8)B[(5:N-5).+2] - (1/6)B[(5:N-5).+1] + (1/8)B[(5:N-5)]);
           (-(1/6)B[(6:N-5).+2] - (1/2)B[(6:N-5).+1] - (1/2)B[(6:N-5)] - (1/6)B[(6:N-5).-1]);
           (+(1/24)B[(7:N-5).+2] + (5/6)B[(7:N-5).+1] + (3/4)B[(7:N-5)] + (5/6)B[(7:N-5).-1] + (1/24)B[(7:N-5).-2]);
           (-(1/6)B[(7:N-4).+1] - (1/2)B[(7:N-4)] - (1/2)B[(7:N-4).-1] - (1/6)B[(7:N-4).-2]);
           (+(1/8)B[(7:N-3)] -(1/6)B[(7:N-3).-1] +(1/8)B[(7:N-3).-2]);
          ]

    M = sparse([I_M0[:];I_M;I_MN[:]], [J_M0[:];J_M;J_MN[:]],
               [V_M0[:];V_M;V_MN[:]])

  elseif p == 6

    bhinv = [43200/13649 8640/12013 4320/2711 4320/5359 8640/7877 43200/43801];

    BS = [25/12 -4 3 -4/3 1/4];

    V_M0 = zeros(9,9)

    (b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12) = B[1:12]
    V_M0[1, 1] = 0.79126675946955820939b1 + 0.29684720906380007429b2 + 0.0031855190887964290152b3 + 0.016324040425909519534b4 + 0.031603022440944150877b5 + 0.031679647480161052996b6 + 0.031485777339472539205b7
    V_M0[1, 2] = V_M0[2, 1] = -1.0166893393503381444b1 - 0.028456273704916113690b3 - 0.041280298383492988198b4 - 0.13922814516201405075b5 -0.11957773256112017666b6 -0.11942677565293334109b7
    V_M0[1, 3] = V_M0[3, 1] = 0.070756429372437150463b1 - 0.18454761060241510503b2 - 0.043641631471118923470b4 + 0.24323679072077324609b5 + 0.15821270735372154440b6 + 0.16023485783647863076b7
    V_M0[1, 4] = V_M0[4, 1] = 0.22519915328913532127b1 - 0.16627487110970548953b2 + 0.027105309616486712977b3 - 0.19166461859684399091b5 - 0.076841171601990145944b6 - 0.082195869498316975759b7
    V_M0[1, 5] = V_M0[5, 1] = -0.052244034642020563167b1 + 0.044400639485098762210b2 - 0.0010239765473093878745b3 + 0.074034846453161740905b4 + 0.012416255689984968954b6  +0.071886528478926012827b5 + 0.013793629971047355034b7
    V_M0[1, 6] = V_M0[6, 1] = -0.018288968138771973527b1 +0.0095746331632217580607b2 - 0.00081057845305764042779b3 - 0.0073488455877755196984b4 + 0.010636019497239069970b5 - 0.013159670383826183824b6 - 0.021179364788387535246b7
    V_M0[1, 7] = V_M0[7, 1] = 0.0019118885633161709274b4 - 0.040681303555291499361b5 + 0.013196749810737491670b6 + 0.025572665181237836763b7
    V_M0[1, 8] = V_M0[8, 1] = 0.015596528711367857640b5 -0.0064861841573315378995b6 -0.0091103445540363197401b7
    V_M0[1, 9] = V_M0[9, 1] = 0.00055939836966298630593b6 - 0.0013848225351007963723b5 + 0.00082542416543781006633b7
    V_M0[2, 2] =              1.3063321571116676286b1 +0.25420017604573457435b3 + 0.10438978280925626095b4 + 0.66723280210321129509b5 + 0.46818193597227494411b6 + 0.46764154101958369201b7
    V_M0[2, 3] = V_M0[3, 2] = -0.090914102699924646049b1 +0.11036113131714764253b4 -1.2903975449975188870b5 -0.66396052487350447871b6 -0.66159744640052061842b7
    V_M0[2, 4] = V_M0[4, 2] = -0.28935573956534316666b1 -0.24213200040645927216b3 +1.1876702550280310277b5 +0.39565981499041363328b6 +0.38600489217558000007b7
    V_M0[2, 5] = V_M0[5, 2] = 0.067127744758037639890b1 +0.0091471926820756301800b3 -0.18721961430038080217b4 -0.13193585588531745301b6 -0.48715757368119118874b5 -0.10475163122754481381b7
    V_M0[2, 6] = V_M0[6, 2] = 0.023499279745900688694b1 +0.0072409053835651813164b3 +0.018583789963916794487b4-0.092896161339386761743b5 +0.12235132704188076670b6 +0.11135203204362950339b7
    V_M0[2, 7] = V_M0[7, 2] = -0.0048347914064469075906b4 +0.23106838326878204031b5 -0.10807741421960079917b6-0.11815617764273433354b7
    V_M0[2, 8] = V_M0[8, 2] = -0.083681414344034553537b5 +0.040934994667670546616b6 +0.042746419676364006921b7
    V_M0[2, 9] = V_M0[9, 2] = -0.0035765451326969831434b6 +0.0073893991241210786821b5 -0.0038128539914240955387b7
    V_M0[3, 3] =              0.0063271611471368738078b1 +0.11473182007158685275b2 +0.11667405542796800075b4+2.7666108082854440372b5 +1.0709206899608171042b6 +1.0131613910329730572b7
    V_M0[3, 4] = V_M0[4, 3] = 0.020137694138847972466b1 +0.10337179946308864017b2 -2.9132216211517427243b5-0.87558073434822622598b6 -0.69099571834888124265b7
    V_M0[3, 5] = V_M0[5, 3] = -0.0046717510915754628683b1 -0.027603533656377128278b2 -0.19792902986208699745b4 +0.54029853383734330523b6 +1.2391775930319110779b5 +0.26280380502473582273b7
    V_M0[3, 6] = V_M0[6, 3] = -0.0016354308669218878195b1 -0.0059524752758832596197b2 +0.019646827777442752194b4+0.32366400126390466006b5 -0.46595166932288709739b6 -0.22172727209417368594b7
    V_M0[3, 7] = V_M0[7, 3] = -0.0051113531893524745496b4 -0.53558781637747543460b5 +0.33283351044897389336b6+0.20786565911785401579b7
    V_M0[3, 8] = V_M0[8, 3] = 0.18243281741342895622b5 -0.10598160301968184459b6 -0.076451214393747111630b7
    V_M0[3, 9] = V_M0[9, 3] = 0.0092090899634437994856b6 -0.015915028188724931671b5 +0.0067059382252811321853b7
    V_M0[4, 4] =              0.064092997759871869867b1 +0.093136576388046999489b2 +0.23063676246347492291b3+3.6894403082837166203b5 +1.1905503386876088738b6 +0.59124795468888565194b7
    V_M0[4, 5] = V_M0[5, 4] = -0.014868958192656041286b1 -0.024870405993901607642b2 -0.0087129289077117541871b3 -1.2635078373718242057b6 -0.30583173978439973269b7 -1.4706919260458029548b5
    V_M0[4, 6] = V_M0[6, 4] = -0.0052051474298559556576b1 -0.0053630987475285424890b2 -0.0068971427657906095463b3-0.78575245216674501017b5 +0.22911480054237346001b7 +0.99770643562927505292b6
    V_M0[4, 7] = V_M0[7, 4] = 0.66972974880676622652b5 -0.50132473560721279390b6 -0.17951612431066454373b7
    V_M0[4, 8] = V_M0[8, 4] = -0.20229090601117515652b5 +0.14534218580636584986b6 +0.056948720204809306656b7
    V_M0[4, 9] = V_M0[9, 4] = -0.012004296184410038337b6 -0.0047769156693859238415b7 +0.016781211853795962179b5
    V_M0[5, 5] =              0.0034494550959102336252b1 +0.0066411834994278261016b2 +0.00032915450832718628585b3+0.33577217075764772000b4 +2.0964133295790264390b6 +0.23173232041831268550b7 +0.0061078257643682645765b8 +0.71091258506833766956b5
    V_M0[5, 6] = V_M0[6, 5] = 0.0012075440723041938061b1 +0.0014321166657521476075b2 +0.00026055826461832559573b3-0.033329411132516353908b4 -0.28082416973855326836b7 -0.027209080835250836084b8+0.10458654356829219874b5 -1.3484369866671155432b6
    V_M0[5, 7] = V_M0[7, 5] = 0.0086710380841746926251b4 +0.17360734113554285637b6 +0.053313621252876254126b8 -0.24249352624045263018b5 +0.15690152576785882706b7
    V_M0[5, 8] = V_M0[8, 5] = -0.086316839802171222760b6 +0.026988423604709992435b7 +0.080981941477156510853b5 -0.032764636390806391639b8
    V_M0[5, 9] = V_M0[9, 5] = 0.0074620594845308550733b6 -0.00081216403616686789496b7 +0.00055227020881270902093b8-0.0072021656571766961993b5
    V_M0[6, 6] =              0.00042272261734493450425b1 +0.00030882419443789644048b2 +0.00020625757066474306202b3+0.0033083434042009682567b4 +0.58280470164050018158b5 +0.80541742203662154736b7+0.13383632334100334433b8 +0.0055555555555555555556b9 +1.1903620718618930511b6
    V_M0[6, 7] = V_M0[7, 6] = -0.00086070442526864133026b4 -0.17480747086739049893b5 -0.31325448501150501650b8 -0.025000000000000000000b9 -0.31691663053104292713b7 -0.66916070916479291611b6
    V_M0[6, 8] = V_M0[8, 6] = 0.033546617916933521087b5 -0.33436200223869714050b7 +0.050000000000000000000b9+0.21697906098076027508b6 +0.18383632334100334433b8
    V_M0[6, 9] = V_M0[9, 6] = 0.029125184768230046430b7 +0.022790919164749163916b8 -0.030689859975187405305b6 -0.0017817995133473605962b5 -0.030555555555555555556b9
    V_M0[7, 7] =              0.00022392237357715991790b4 +0.12754377854309566738b5 +1.0116994839296081646b6 +0.96988172751725752475b8 +0.12500000000000000000b9 +0.0055555555555555555556b10 +0.48231775430312815001b7
    V_M0[7, 8] = V_M0[8, 7] = -0.037841139730330129499b5 -0.29975568851348273616b6 -0.30000000000000000000b9 - 0.025000000000000000000b10 -0.39914868674468211784b7 -0.43825448501150501650b8
    V_M0[7, 9] = V_M0[9, 7] = 0.046981462180226839339b6 -0.29668637874712374587b8 +0.050000000000000000000b10 +0.17163557041460064817b7 +0.0030693461522962583624b5 +0.17500000000000000000b9
    V_M0[8, 8] =              0.012303289427168044554b5 +0.11836475296458983325b6 +0.94105118982279433342b7 +0.95000000000000000000b9 +0.12500000000000000000b10 +0.0055555555555555555556b11 +0.56994743445211445545b8
    V_M0[8, 9] = V_M0[9, 8] = -0.023080678926719163396b6 -0.29866250537751494972b7 -0.30000000000000000000b10 -0.025000000000000000000b11 -0.0010477348605150508026b5 -0.42720908083525083608b8 - 0.42500000000000000000b9
    V_M0[9, 9] =              0.0051393702211491099770b6 +0.12477232150094220014b7 +0.95055227020881270902b8 +0.95000000000000000000b10 +0.12500000000000000000b11 +0.0055555555555555555556b12 +0.000091593624651536418269b5 +0.56111111111111111111b9
    I_M0 = (1:9) * ones(1,9)
    J_M0 = I_M0'

    V_MN = zeros(9,9)
    (b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12) = B[N+1:-1:N-10]
    V_MN[9, 9] =              0.79126675946955820939b1 + 0.29684720906380007429b2 + 0.0031855190887964290152b3 + 0.016324040425909519534b4 + 0.031603022440944150877b5 + 0.031679647480161052996b6 + 0.031485777339472539205b7
    V_MN[9, 8] = V_MN[8, 9] = -1.0166893393503381444b1 - 0.028456273704916113690b3 - 0.041280298383492988198b4 - 0.13922814516201405075b5 -0.11957773256112017666b6 -0.11942677565293334109b7
    V_MN[9, 7] = V_MN[7, 9] = 0.070756429372437150463b1 - 0.18454761060241510503b2 - 0.043641631471118923470b4 + 0.24323679072077324609b5 + 0.15821270735372154440b6 + 0.16023485783647863076b7
    V_MN[9, 6] = V_MN[6, 9] = 0.22519915328913532127b1 - 0.16627487110970548953b2 + 0.027105309616486712977b3 - 0.19166461859684399091b5 - 0.076841171601990145944b6 - 0.082195869498316975759b7
    V_MN[9, 5] = V_MN[5, 9] = -0.052244034642020563167b1 + 0.044400639485098762210b2 - 0.0010239765473093878745b3 + 0.074034846453161740905b4 + 0.012416255689984968954b6  +0.071886528478926012827b5 + 0.013793629971047355034b7
    V_MN[9, 4] = V_MN[4, 9] = -0.018288968138771973527b1 +0.0095746331632217580607b2 -0.00081057845305764042779b3 -0.0073488455877755196984b4 +0.010636019497239069970b5 -0.013159670383826183824b6-0.021179364788387535246b7
    V_MN[9, 3] = V_MN[3, 9] = 0.0019118885633161709274b4 -0.040681303555291499361b5 +0.013196749810737491670b6 +0.025572665181237836763b7
    V_MN[9, 2] = V_MN[2, 9] = 0.015596528711367857640b5 -0.0064861841573315378995b6 -0.0091103445540363197401b7
    V_MN[9, 1] = V_MN[1, 9] = 0.00055939836966298630593b6 -0.0013848225351007963723b5 +0.00082542416543781006633b7
    V_MN[8, 8] =              1.3063321571116676286b1 +0.25420017604573457435b3 +0.10438978280925626095b4 +0.66723280210321129509b5 +0.46818193597227494411b6 +0.46764154101958369201b7
    V_MN[8, 7] = V_MN[7, 8] = -0.090914102699924646049b1 +0.11036113131714764253b4 -1.2903975449975188870b5 -0.66396052487350447871b6 -0.66159744640052061842b7
    V_MN[8, 6] = V_MN[6, 8] = -0.28935573956534316666b1 -0.24213200040645927216b3 +1.1876702550280310277b5 +0.39565981499041363328b6 +0.38600489217558000007b7
    V_MN[8, 5] = V_MN[5, 8] = 0.067127744758037639890b1 +0.0091471926820756301800b3 -0.18721961430038080217b4 -0.13193585588531745301b6 -0.48715757368119118874b5 -0.10475163122754481381b7
    V_MN[8, 4] = V_MN[4, 8] = 0.023499279745900688694b1 +0.0072409053835651813164b3 +0.018583789963916794487b4-0.092896161339386761743b5 +0.12235132704188076670b6 +0.11135203204362950339b7
    V_MN[8, 3] = V_MN[3, 8] = -0.0048347914064469075906b4 +0.23106838326878204031b5 -0.10807741421960079917b6-0.11815617764273433354b7
    V_MN[8, 2] = V_MN[2, 8] = -0.083681414344034553537b5 +0.040934994667670546616b6 +0.042746419676364006921b7
    V_MN[8, 1] = V_MN[1, 8] = -0.0035765451326969831434b6 +0.0073893991241210786821b5 -0.0038128539914240955387b7
    V_MN[7, 7] =               0.0063271611471368738078b1 +0.11473182007158685275b2 +0.11667405542796800075b4+2.7666108082854440372b5 +1.0709206899608171042b6 +1.0131613910329730572b7
    V_MN[7, 6] = V_MN[6, 7] = 0.020137694138847972466b1 +0.10337179946308864017b2 -2.9132216211517427243b5-0.87558073434822622598b6 -0.69099571834888124265b7
    V_MN[7, 5] = V_MN[5, 7] = -0.0046717510915754628683b1 -0.027603533656377128278b2 -0.19792902986208699745b4 +0.54029853383734330523b6 +1.2391775930319110779b5 +0.26280380502473582273b7
    V_MN[7, 4] = V_MN[4, 7] = -0.0016354308669218878195b1 -0.0059524752758832596197b2 +0.019646827777442752194b4+0.32366400126390466006b5 -0.46595166932288709739b6 -0.22172727209417368594b7
    V_MN[7, 3] = V_MN[3, 7] = -0.0051113531893524745496b4 -0.53558781637747543460b5 +0.33283351044897389336b6+0.20786565911785401579b7
    V_MN[7, 2] = V_MN[2, 7] = 0.18243281741342895622b5 -0.10598160301968184459b6 -0.076451214393747111630b7
    V_MN[7, 1] = V_MN[1, 7] = 0.0092090899634437994856b6 -0.015915028188724931671b5 +0.0067059382252811321853b7
    V_MN[6, 6] =              0.064092997759871869867b1 +0.093136576388046999489b2 +0.23063676246347492291b3+3.6894403082837166203b5 +1.1905503386876088738b6 +0.59124795468888565194b7
    V_MN[6, 5] = V_MN[5, 6] = -0.014868958192656041286b1 -0.024870405993901607642b2 -0.0087129289077117541871b3 -1.2635078373718242057b6 -0.30583173978439973269b7 -1.4706919260458029548b5
    V_MN[6, 4] = V_MN[4, 6] = -0.0052051474298559556576b1 -0.0053630987475285424890b2 -0.0068971427657906095463b3-0.78575245216674501017b5 +0.22911480054237346001b7 +0.99770643562927505292b6
    V_MN[6, 3] = V_MN[3, 6] = 0.66972974880676622652b5 -0.50132473560721279390b6 -0.17951612431066454373b7
    V_MN[6, 2] = V_MN[2, 6] = -0.20229090601117515652b5 +0.14534218580636584986b6 +0.056948720204809306656b7
    V_MN[6, 1] = V_MN[1, 6] = -0.012004296184410038337b6 -0.0047769156693859238415b7 +0.016781211853795962179b5
    V_MN[5, 5] =              0.0034494550959102336252b1 +0.0066411834994278261016b2 +0.00032915450832718628585b3+0.33577217075764772000b4 +2.0964133295790264390b6 +0.23173232041831268550b7 +0.0061078257643682645765b8 +0.71091258506833766956b5
    V_MN[5, 4] = V_MN[4, 5] = 0.0012075440723041938061b1 +0.0014321166657521476075b2 +0.00026055826461832559573b3-0.033329411132516353908b4 -0.28082416973855326836b7 -0.027209080835250836084b8+0.10458654356829219874b5 -1.3484369866671155432b6
    V_MN[5, 3] = V_MN[3, 5] = 0.0086710380841746926251b4 +0.17360734113554285637b6 +0.053313621252876254126b8 -0.24249352624045263018b5 +0.15690152576785882706b7
    V_MN[5, 2] = V_MN[2, 5] = -0.086316839802171222760b6 +0.026988423604709992435b7 +0.080981941477156510853b5 -0.032764636390806391639b8
    V_MN[5, 1] = V_MN[1, 5] = 0.0074620594845308550733b6 -0.00081216403616686789496b7 +0.00055227020881270902093b8-0.0072021656571766961993b5
    V_MN[4, 4] =              0.00042272261734493450425b1 +0.00030882419443789644048b2 +0.00020625757066474306202b3+0.0033083434042009682567b4 +0.58280470164050018158b5 +0.80541742203662154736b7+0.13383632334100334433b8 +0.0055555555555555555556b9 +1.1903620718618930511b6
    V_MN[4, 3] = V_MN[3, 4] = -0.00086070442526864133026b4 -0.17480747086739049893b5 -0.31325448501150501650b8 -0.025000000000000000000b9 -0.31691663053104292713b7 -0.66916070916479291611b6
    V_MN[4, 2] = V_MN[2, 4] = 0.033546617916933521087b5 -0.33436200223869714050b7 +0.050000000000000000000b9+0.21697906098076027508b6 +0.18383632334100334433b8
    V_MN[4, 1] = V_MN[1, 4] = 0.029125184768230046430b7 +0.022790919164749163916b8 -0.030689859975187405305b6 -0.0017817995133473605962b5 -0.030555555555555555556b9
    V_MN[3, 3] =              0.00022392237357715991790b4 +0.12754377854309566738b5 +1.0116994839296081646b6 +0.96988172751725752475b8 +0.12500000000000000000b9 +0.0055555555555555555556b10 +0.48231775430312815001b7
    V_MN[3, 2] = V_MN[2, 3] = -0.037841139730330129499b5 -0.29975568851348273616b6 -0.30000000000000000000b9 - 0.025000000000000000000b10 -0.39914868674468211784b7 -0.43825448501150501650b8
    V_MN[3, 1] = V_MN[1, 3] = 0.046981462180226839339b6 -0.29668637874712374587b8 +0.050000000000000000000b10 +0.17163557041460064817b7 +0.0030693461522962583624b5 +0.17500000000000000000b9
    V_MN[2, 2] =               0.012303289427168044554b5 +0.11836475296458983325b6 +0.94105118982279433342b7 +0.95000000000000000000b9 +0.12500000000000000000b10 +0.0055555555555555555556b11 +0.56994743445211445545b8
    V_MN[2, 1] = V_MN[1, 2] = -0.023080678926719163396b6 -0.29866250537751494972b7 -0.30000000000000000000b10 -0.025000000000000000000b11 -0.0010477348605150508026b5 -0.42720908083525083608b8 - 0.42500000000000000000b9
    V_MN[1, 1] =              0.0051393702211491099770b6 +0.12477232150094220014b7 +0.95055227020881270902b8 +0.95000000000000000000b10 +0.12500000000000000000b11 +0.0055555555555555555556b12 +0.000091593624651536418269b5 +0.56111111111111111111b9

    I_MN = (N-7:N+1) * ones(1,9)
    J_MN = I_MN'

    I_M = [10:N-5;
        10:N-6;
	10:N-7;
	10:N-8;
	9:N-8;
	8:N-8;
	7:N-8]
    J_M = [7:N-8;
        8:N-8;
	9:N-8;
	10:N-8;
	10:N-7;
	10:N-6;
	10:N-5]
    V_M =  [(-(11/360)B[(10:N-5).-3] + (1/40)B[(10:N-5).-2] + (1/40)B[(10:N-5).-1] - (11/360)B[(10:N-5)]);
    	    (+(1/20)B[(10:N-6).-3] + (7/40)B[(10:N-6).-2] - (3/10)B[(10:N-6).-1] + (7/40)B[(10:N-6)] + (1/20)B[(10:N-6).+1]);
            (-(1/40)B[(10:N-7).-3] -(3/10)B[(10:N-7).-2] - (17/40)B[(10:N-7).-1] - (17/40)B[(10:N-7)] - (3/10)B[(10:N-7).+1] - (1/40)B[(10:N-7).+2]);
	    (+(1/180)B[(10:N-8).-3] +(1/8)B[(10:N-8).-2] + (19/20)B[(10:N-8).-1] + (101/180)B[(10:N-8)] + (19/20)B[(10:N-8).+1] + (1/8)B[(10:N-8).+2]+ (1/180)B[(10:N-8).+3]);
	    (-(1/40)B[(9:N-8).-2] - (3/10)B[(9:N-8).-1] - (17/40)B[(9:N-8)] - (17/40)B[(9:N-8).+1] - (3/10)B[(9:N-8).+2] - (1/40)B[(9:N-8).+3]);
	    (+(1/20)B[(8:N-8).-1] + (7/40)B[(8:N-8)] - (3/10)B[(8:N-8).+1] + (7/40)B[(8:N-8).+2]+ (1/20)B[(8:N-8).+3] );
	    (-(11/360)B[(7:N-8)] + (1/40)B[(7:N-8).+1] + (1/40)B[(7:N-8).+2] - (11/360)B[(7:N-8).+3]);
          ]

    M = sparse([I_M0[:];I_M[:];I_MN[:]],
               [J_M0[:];J_M[:];J_MN[:]],
               [V_M0[:];V_M[:];V_MN[:]])

  else
    error(string("Operators for order ", p, " are not implemented"))
  end

  bm = length(bhinv)
  Np = N+1

  if(Np < 2*bm)
    error("Grid not big enough to support the operator. Grid must have N >= ", max(bn,2*bm))
  end

  h = (xc[2] - xc[1]) / N
  M = M/h
  @assert h > 0
  H_I = 1:Np
  H_V = ones(Np)
  H_V[1:bm] = bhinv[:]
  H_V[Np-bm+1:Np] = bhinv[end:-1:1]
  HI = sparse(H_I, H_I, H_V / h)
  H  = sparse(H_I, H_I, h ./ H_V)

  S0 = sparse(ones(length(BS)), 1:length(BS), -B[1]*BS[:]/h, N+1, N+1);
  SN = sparse((N+1) * ones(length(BS)), (N+1):-1:(N+2-length(BS)),
              B[N+1]*BS[:]/h, N+1, N+1);
  D = HI * (-M + SN - S0)

  r = range(xc[1], stop=xc[2], length=N+1)

  BS = SN - S0
  (D, BS, HI, H, r)
end
#}}}

nothing


text:
#!/usr/bin/env julia
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Adds the julia license line `new_license` to configured file extensions in `rootdirs`.
#
# Option `old_license` to remove an existing license first in case one wants to change the
# license text in the future.
#
# Checks also if somewhere else in the file the license text is found (`copy/past error`)
# and if possible deletes such lines - if other text is on the same line it raises an error.

### CONFIG HERE

const print_result = true  # prints files which where not processed.

const rootdirs = [
    "../base",
    "../cli",
    "../contrib",
    "../src",
    "../stdlib",
]

# to exclude whole sub directories
const excludedirs = [
    # see: https://github.com/JuliaLang/julia/pull/11073#issuecomment-98090053
    "../base/ryu",
    "../src/flisp",
    "../stdlib/TOML/test/testfiles",
    "../test/testhelpers/allocation_file.jl",
]

const skipfiles = [
    "../contrib/add_license_to_files.jl",
    "../contrib/asan/check.jl",
    # files to check - already copyright
    # see: https://github.com/JuliaLang/julia/pull/11073#issuecomment-98099389
    "../base/special/trig.jl",
    "../base/special/exp.jl",
    "../base/special/rem_pio2.jl",
    #
    "../src/abi_llvm.cpp",
    "../src/abi_ppc64le.cpp",
    "../src/abi_win32.cpp",
    "../src/abi_win64.cpp",
    "../src/abi_x86.cpp",
    "../src/abi_x86_64.cpp",
    "../src/disasm.cpp",
    "../src/getopt.c",
    "../src/getopt.h",
    "../src/support/END.h",
    "../src/support/ENTRY.amd64.h",
    "../src/support/ENTRY.i387.h",
    "../src/support/_setjmp.win32.S",
    "../src/support/MurmurHash3.c",
    "../src/support/MurmurHash3.h",
    "../src/support/asprintf.c",
    "../src/support/dirname.c",
    "../src/support/strptime.c",
    "../src/support/strtod.c",
    "../src/support/tzfile.h",
    "../src/support/utf8.c",
    "../src/crc32c.c",
]

const ext_prefix = Dict([
    (".jl", "# "),
    (".sh", "# "),
    (".h", "// "),
    (".c", "// "),
    (".cpp", "// "),
    (".S", "// "),
])

const new_license = "This file is a part of Julia. License is MIT: https://julialang.org/license"

# Old License text if such should be first removed - or empty string
const old_license = ""

### END CONFIG HERE

function check_lines!(
    path::AbstractString, lines::Vector, checktxt::AbstractString,
    prefix::AbstractString, oldcheck::Bool)
    remove = []
    for i in 1:length(lines)
        line = lines[i]
        if occursin(checktxt, line)
            if strip(line) == strip(prefix * checktxt) || strip(line) == strip(checktxt)
                push!(remove, i)
            else
                error(string("`path` contains an additional line with ",
                         oldcheck ? "old" : "new",
                         " license.\nlinenum: $(i): $(line) \n`path:` $(path)\n",
                         "Fix this first or add the file to `skipfiles`.\n\n"))
            end
        end
    end
    deleteat!(lines, remove)
end

license_linenum(line) = startswith(strip(line), "#!") ? 2 : 1

# Collects all absolute file paths in rootdir inclusive subdirs
function getfilespaths!(filepaths::Vector, rootdir::AbstractString)
    isdir(rootdir) || error(string("`rootdir` must be an directory. "))
    abs_rootdir = abspath(rootdir)
    for name in readdir(abs_rootdir)
        path = joinpath(abs_rootdir, name)
        islink(path) && continue
        if isdir(path)
            getfilespaths!(filepaths, path)
        else
            push!(filepaths, joinpath(abs_rootdir, name))
        end
    end
end

function add_license_line!(unprocessed::Vector, src::AbstractString, new_license::AbstractString,
                          old_license::AbstractString, ext_prefix::Dict, abs_excludedirs::Vector,
                                                                               skipfiles::Vector)

    for name in readdir(src)
        path = normpath(joinpath(src, name))
        islink(path) && continue
        if isdir(path)
            if path in abs_excludedirs
                getfilespaths!(unprocessed, path)
                continue
            else
                add_license_line!(unprocessed, path, new_license, old_license,
                                  ext_prefix, abs_excludedirs, skipfiles)
            end
        elseif path in skipfiles
            push!(unprocessed, path)
            continue
        else
            ext = splitext(path)[2]
            if ext in keys(ext_prefix)
                prefix = ext_prefix[ext]
                f = open(path, "r")
                lines = readlines(f, keep=true)
                close(f)
                isempty(lines) && (push!(unprocessed, path); continue)
                isempty(old_license) || check_lines!(path, lines, old_license, prefix, true)
                check_lines!(path, lines, new_license, prefix, false)
                isempty(lines) && continue  # file consisting of just license header
                # check shebang file
                linenum = license_linenum(lines[1])
                if !isempty(strip(lines[linenum]))
                    insert!(lines, linenum, string(prefix, new_license, "\n\n"))
                else
                    insert!(lines, linenum, string(prefix, new_license, "\n"))
                end
                open(path, "w") do f
                    for line in lines
                        write(f, line)
                    end
                end
            else
                push!(unprocessed, path)
            end
        end
    end
    return unprocessed
end

# Returns a new Vector with all absolute path: raises an error if path does not exist
function abspaths(A::Vector)
    abs_A = []
    for p in A
        abs_p = isabspath(p) ? normpath(p) : normpath(joinpath(@__DIR__, p))
        ispath(abs_p) || error(string("`abs_p` seems not to be an existing path. ",
                                      "Adjust your configuration: <", p, "> : ", abs_p, "\n"))
        push!(abs_A, abs_p)
    end
    return abs_A
end

function add_license(rootdirs::Vector, new_license::AbstractString, old_license::AbstractString,
                     ext_prefix::Dict, excludedirs::Vector, skipfiles::Vector, print_result::Bool)
    isempty(strip(new_license)) && error("`new_license` may not only contain white space.")
    abs_skipfiles = abspaths(skipfiles)
    abs_rootdirs = abspaths(rootdirs)
    abs_excludedirs = abspaths(excludedirs)
    for p in abs_rootdirs
        p in abs_excludedirs && error(string("Seems one of the `rootdirs` is also included in ",
                                             "the `excludedirs`. `rootdirs`: ", p))
        unprocessed = []
        add_license_line!(unprocessed, p, new_license, old_license, ext_prefix,
                                                abs_excludedirs, abs_skipfiles)
        if print_result
            println("\nUnprocessed files in rootdir: <$(p)>\n")
            for file in unprocessed
                println("  <", basename(file), "> : ", file)
            end
        end
    end
end

## ---------------

add_license(rootdirs, new_license, old_license, ext_prefix, excludedirs, skipfiles, print_result)

ext: jl

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T}
    m, n = length(vy), length(vx)
    gx = reshape(repeat(vx, inner = m, outer = 1), m, n)
    gy = reshape(repeat(vy, inner = 1, outer = n), m, n)
    return gx, gy
end

function plotedges(edges)::Tuple{Array{Float64,1}, Array{Float64,1}}
    x = Array{Float64,1}()
    y = Array{Float64,1}()
    for edge in edges
        push!(x, getx(geta(edge)))
        push!(x, getx(getb(edge)))
        push!(x, NaN)
        push!(y, gety(geta(edge)))
        push!(y, gety(getb(edge)))
        push!(y, NaN)
    end
    return x, y
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
#!/usr/bin/env julia
# del2z <delta.z@aliyun.com>

#module Word2Vec

using Flux
using Flux: glorot_normal, ADAM, OneHotMatrix, OneHotVector, onehot, onehotbatch, onecold, logitbinarycrossentropy, throttle
using JSON
using DataStructures
using StatsBase

const config = JSON.parsefile(joinpath(@__DIR__, "word2vec.json"))

include("prepro.jl")

corpus = segment!(loaddata("corpus.txt"))
wordfreq = genvocab(corpus, 1)

index2word = ["<UNK>"; "</s>"; collect(keys(wordfreq))]
vocabsize = length(index2word)
word2index = OrderedDict(zip(index2word, 1:vocabsize))

## model graph
embedsize = config["model"]["embed_size"]
W̃_in = param(glorot_normal(embedsize, vocabsize))
W̃_out = param(glorot_normal(vocabsize, embedsize))
b̃_out = param(zeros(vocabsize))

function dataset end
function model end
function loss end

window = config["model"]["window"]
@assert isodd(window)
negscale = config["model"]["neg_scale"]
if config["method"] == "cbow"
    dataset(corpus::Vector{Vector{String}}) = begin
        Xs = Vector{OneHotMatrix}()
        Ys = Vector{OneHotVector}()
        for seq in corpus
            size(seq, 1) >= window || continue
            for k in 1:(size(seq, 1) - window + 1)
                chunk = seq[k:(k - 1 + window)]
                in(chunk[(window + 1) ÷ 2], index2word) || continue
                push!(Xs, onehotbatch([chunk[i] for i in filter(x -> x != (window + 1) ÷ 2, 1:window)],
                                      index2word, "<UNK>"))
                push!(Ys, onehot(chunk[(window + 1) ÷ 2], index2word, "<UNK>"))
            end
        end
        Xs, Ys
    end

    model(x::OneHotMatrix) = begin
        W̃_out * sum(W̃_in * x, dims = 2) + b̃_out
    end

    loss(x::OneHotMatrix, y::OneHotVector) = begin
        posid = onecold(y)
        negids = sample(setdiff(1:size(y, 1), [posid]), negscale, replace = false, ordered = true)
        ŷ = model(x)
        ℒ = logitbinarycrossentropy(ŷ[posid], 1) + sum(logitbinarycrossentropy.(ŷ[negids], zeros(length(negids))))
        ℒ / (1 + length(negids))
    end
else
    dataset(corpus::Vector{Vector{String}}) = begin
        Xs = Vector{OneHotVector}()
        Ys = Vector{OneHotMatrix}()
        for seq in corpus
            size(seq, 1) >= window || continue
            for k in 1:(size(seq, 1) - window + 1)
                chunk = seq[k:(k - 1 + window)]
                in(chunk[(window + 1) ÷ 2], index2word) || continue
                push!(Xs, onehot(chunk[(window + 1) ÷ 2], index2word, "<UNK>"))
                push!(Ys, onehotbatch([chunk[i] for i in filter(x -> x != (window + 1) ÷ 2, 1:window)],
                                      index2word, "<UNK>"))
            end
        end
        Xs, Ys
    end

    model(x::OneHotVector) = begin
        W̃_out * W̃_in * x + b̃_out
    end

    loss(x::OneHotVector, y::OneHotMatrix) = begin
        posids = onecold(y)
        negids = sample(setdiff(1:size(y, 1), posids), negscale * length(posids), replace = false, ordered = true)
        ŷ = model(x)
        ℒ = sum(logitbinarycrossentropy.(ŷ[posids], ones(length(posids)))) + sum(logitbinarycrossentropy.(ŷ[negids], zeros(length(negids))))
        ℒ / (length(posids) + length(negids))
    end
end

## training
Xs, Ys = dataset(corpus)
@info size(Xs), size(Ys)
opt = ADAM(0.002)

for ep in 1:config["model"]["epochs"]
    Flux.train!(loss, params(W̃_in, W̃_out, b̃_out), zip(Xs, Ys), opt)
    @info sum(map(p -> loss(p...), zip(Xs, Ys)))
end

## testing

#end


text:
using quantumoptics

# System parameters for a single decaying spin
ω0 = 1.
γ = 1.
R = 0.5

# System operators
spinbasis = SpinBasis(1//2)
sz = sigmaz(spinbasis)
sm = sigmam(spinbasis)
sp = sigmap(spinbasis)

H = ω0 * sz
J = [γ*sm, R*sp]

# Initial state
Ψ₀ = spinup(spinbasis)
ρ₀ = Ψ₀ ⊗ dagger(Ψ₀)

# Integration options
Ntrajectories = 1000
T = [0:0.01:10;]

# Master time evolution
tout_master, ρt_master = timeevolution.master(T, ρ₀, H, J)

# Averaged MCWF time evolution
ρt_average = Operator[ρ₀*0. for i=1:length(T)]
for i=1:Ntrajectories
    tout, Ψt = timeevolution.mcwf(T, Ψ₀, H, J)
    for j=1:length(T)
        ρt_average[j] += (Ψt[j] ⊗ dagger(Ψt[j]))/Ntrajectories
    end
end

# Example trajectories
T_example = [0.,100.]
tout_example, Ψt_example = timeevolution.mcwf(T_example, Ψ₀, H, J; seed=UInt(1),
                                        display_beforeevent=true,
                                        display_afterevent=true)

using PyCall
@pyimport matplotlib.pyplot as plt

plt.subplot(2, 1, 1)
plt.plot(T, real(expect(sz, ρt_master)))
plt.plot(T, real(expect(sz, ρt_average)))
plt.xlabel("\$\\mathrm\{Time\} \$")
plt.ylabel("\$\\sigma^z\$")

plt.subplot(2, 1, 2)
plt.plot(tout_example, real(expect(sz, Ψt_example)))
plt.xlabel("\$\\mathrm\{Time\} \$")
plt.ylabel("\$\\sigma^z\$")

plt.show()


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__
cd(ProjDir) #do

d_string = "dag {A -> {E Z}; B -> {D Z}; Z -> {D E}; E -> D}"
dag = DAG("test_open_paths_02", d_string);

@test to_dagitty(dag.d) == "dag { D <- E; D <- Z; E <- Z; D <- B; Z <- B; E <- A; Z <- A }"

bs = basis_set(dag)
@test bs[1][1] == :B
@test bs[1][2] == :A
@test bs[1][3:end] == Symbol[]
@test bs[2][1] == :D
@test bs[2][2] == :A
@test length(bs[2]) == 5
@test bs[3][1] == :E
@test bs[3][2] == :B
@test bs[3][3:end] == [:A, :Z]

#=
fname = joinpath(ProjDir, "test_open_paths_02.dot")
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)
=#

ap  = all_paths(dag, :D, :E)
bp = backdoor_paths(dag, ap, :D)
adjs = adjustment_sets(dag, :E, :D)

@testset "Open_path_02" begin

  @test bp == [[:D, :Z, :E], [:D, :Z, :A, :E], [:D, :B, :Z, :E], [:D, :B, :Z, :A, :E]]
  @test adjs == [[:Z, :B], [:Z, :A]]

end


text:
using NumericalDataManipulation
using Test


text:
# load required tensor flow functionality

ENV["PYTHON"] = "/usr/bin/python3.6"

using Pkg
Pkg.build("PyCall")

using PyCall

py"""
import sys
sys.path.insert(0, "/home/singhk/learning-odometry/")
"""
# pushfirst!(PyVector(pyimport("sys")."path"), "/home/singhk/learning-odometry/" )

tf = pyimport("tensorflow")

tf.keras.models.load_model

Base.cd(ENV["USER"]*"/learning-odometry/")
mmodel = tf.keras.models.load_model("original_network0.h5")

mmodel.predict

A = [rand(4) for i in 1:25]

mmodel.predict([[A]])
# 1×2 Array{Float32,2}:
#  0.0169379  1.58055

models = []
# 0-element Array{Any,1}

using ProgressMeter

@showprogress "loading models" for i in 0:9
  push!(models, tf.keras.models.load_model("original_network$i.h5"))
end

models[1].predict([[A]])

preds = map(x->x.predict([[A]]), models)

## Using aggregate python functions

py"""
import sys
sys.path.insert(0, "/home/singhk/learning-odometry/")
"""

predsingle = pyimport("queryAllModels")["predictSinglePose"]

predsingle(A)
# returns  10×3 Array{Real,2}:


text:
"""
    ancilla_extraction(zxg::ZXGraph) -> ZXDiagram

Extract a quantum circuit from a general `ZXGraph` even without a gflow.
It will introduce post-selection operators.
"""
function ancilla_extraction(zxg::ZXGraph)
    nzxg = copy(zxg)
    simplify!(Rule(:scalar), nzxg)
    ins = copy(get_inputs(nzxg))
    outs = copy(get_outputs(nzxg))
    nbits = length(outs)
    gads = Dict{Int, Int}()
    for v in spiders(nzxg)
        if spider_type(nzxg, v) == SpiderType.Z && degree(nzxg, v) == 1
            v1 = neighbors(nzxg, v)[1]
            if phase(nzxg, v1) in (0, 1)
                gads[v1] = v
            end
        end
    end

    for i in eachindex(ins)
        v = ins[i]
        @inbounds u = neighbors(nzxg, v)[1]
        if !is_hadamard(nzxg, u, v)
            insert_spider!(nzxg, u, v)
        end
    end

    frontiers = copy(outs)
    circ = ZXDiagram(nbits)
    unextracts = Set(spiders(nzxg))
    qubit_map = Dict{Int, Int}()
    for i in eachindex(ins)
        v = ins[i]
        qubit_map[v] = i
        delete!(unextracts, v)
    end
    for i in eachindex(frontiers)
        @inbounds v = frontiers[i]
        qubit_map[v] = i
        delete!(unextracts, v)
        @inbounds v1 = neighbors(nzxg, v)[1]
        if !(v1 in ins)
            is_hadamard(nzxg, v, v1) && pushfirst_gate!(circ, Val(:H), i)
            frontiers[i] = v1
            qubit_map[v1] = i
            rem_edge!(nzxg, v, v1)
        end
    end

    while !isempty(unextracts)
        update_frontier_ancilla!(frontiers, nzxg, gads, qubit_map, unextracts, circ)
    end

    for v in frontiers
        if degree(nzxg, v) == 1
            @inbounds v1 = neighbors(nzxg, v)[1]
            if is_hadamard(nzxg, v1, v)
                pushfirst_gate!(circ, Val(:H), qubit_map[v])
            end
        end
    end

    M = biadjacency(nzxg, frontiers, ins)
    M, steps = normalize_perm(M)
    for step in steps
        @assert step.op === :swap
        # pushfirst_gate!(circ, Val(:SWAP), [step.r1, step.r2])
        pushfirst_gate!(circ, Val(:CNOT), step.r1, step.r2)
        pushfirst_gate!(circ, Val(:CNOT), step.r2, step.r1)
        pushfirst_gate!(circ, Val(:CNOT), step.r1, step.r2)
    end

    simplify!(Rule(:i1), circ)
    simplify!(Rule(:i2), circ)
    return circ
end

function update_frontier_ancilla!(frontiers, nzxg, gads, qubit_map, unextracts, circ)
    nbs = Int[]
    for i in 1:length(frontiers)
        v = frontiers[i]
        if phase(nzxg, v) != 0
            pushfirst_gate!(circ, Val(:Z), i, phase(nzxg, v))
            set_phase!(nzxg, v, zero(phase(nzxg, v)))
        end
        for j in (i+1):length(frontiers)
            u = frontiers[j]
            if has_edge(nzxg, u, v)
                pushfirst_gate!(circ, Val(:CZ), i, j)
                rem_edge!(nzxg, u, v)
            end
        end
    end
    for i in eachindex(frontiers)
        v = frontiers[i]
        nb_v = neighbors(nzxg, v)
        if length(nb_v) == 1
            delete!(unextracts, v)
            @inbounds u = nb_v[1]
            if spider_type(nzxg, u) == SpiderType.Z
                pushfirst_gate!(circ, Val(:H), i)
                frontiers[i] = u
                qubit_map[u] = i
                rem_edge!(nzxg, u, v)
                return frontiers
            end
        elseif length(nb_v) == 0
            delete!(unextracts, v)
        end
        for u in nb_v
            if haskey(gads, u)
                rewrite!(Rule(:pivot), nzxg, [u, gads[u], v])
                pushfirst_gate!(circ, Val(:H), i)
                delete!(unextracts, gads[u])
                delete!(gads, u)
                frontiers[i] = u
                qubit_map[u] = i
                return frontiers
            else
                spider_type(nzxg, u) == SpiderType.Z &&
                !(u in nbs) && push!(nbs, u)
            end
        end
    end

    length(nbs) == 0 && return frontiers

    M = biadjacency(nzxg, frontiers, nbs)
    M0, steps = gaussian_elimination(M)
    ws = Int[]
    @inbounds for i = 1:length(frontiers)
        if sum(M0[i,:]) == 1
            push!(ws, nbs[findfirst(isone, M0[i,:])])
        end
    end
    if length(ws) > 0
        @inbounds for e in findall(M .== 1)
            if has_edge(nzxg, frontiers[e[1]], nbs[e[2]])
                rem_edge!(nzxg, frontiers[e[1]], nbs[e[2]])
            end
        end
        @inbounds for e in findall(M0 .== 1)
            add_edge!(nzxg, frontiers[e[1]], nbs[e[2]])
        end

        @inbounds for step in steps
            if step.op == :addto
                ctrl = qubit_map[frontiers[step.r2]]
                loc = qubit_map[frontiers[step.r1]]
                pushfirst_gate!(circ, Val{:CNOT}(), loc, ctrl)
            elseif step.op == :swap
                q1 = qubit_map[frontiers[step.r1]]
                q2 = qubit_map[frontiers[step.r2]]

                # pushfirst_gate!(circ, Val{:SWAP}(), [q1, q2])
                pushfirst_gate!(circ, Val{:CNOT}(), q1, q2)
                pushfirst_gate!(circ, Val{:CNOT}(), q2, q1)
                pushfirst_gate!(circ, Val{:CNOT}(), q1, q2)
            end
        end
        return frontiers
    else
        w = nbs[1]
        push!(frontiers, w)
        add_ancilla!(circ, SpiderType.Z, SpiderType.Z)
        qubit_map[w] = length(frontiers)
        return frontiers
    end
end

"""
    circuit_extraction(zxg::ZXGraph)

Extract circuit from a graph-like ZX-diagram.
"""
function circuit_extraction(zxg::ZXGraph{T, P}) where {T, P}
    nzxg = copy(zxg)
    nbits = nqubits(nzxg)
    gads = Set{T}()
    for v in spiders(nzxg)
        if spider_type(nzxg, v) == SpiderType.Z && degree(nzxg, v) == 1
            v1 = neighbors(nzxg, v)[1]
            if phase(nzxg, v1) in (0, 1)
                push!(gads, v, v1)
            end
        end
    end

    Outs = get_outputs(nzxg)
    Ins = get_inputs(nzxg)
    if nbits == 0
        nbits = length(Outs)
    end
    cir = Chain()
    if length(Outs) != length(Ins)
        return cir
    end
    for v1 in Ins
        @inbounds v2 = neighbors(nzxg, v1)[1]
        if !is_hadamard(nzxg, v1, v2)
            insert_spider!(nzxg, v1, v2)
        end
    end
    @inbounds frontier = [neighbors(nzxg, v)[1] for v in Outs]
    qubit_map = Dict(zip(frontier, 1:nbits))

    for i = 1:nbits
        @inbounds w = neighbors(nzxg, Outs[i])[1]
        @inbounds if is_hadamard(nzxg, w, Outs[i])
            pushfirst_gate!(cir, Val{:H}(), i)
        end
        if phase(nzxg, w) != 0
            pushfirst_gate!(cir, Val{:Rz}(), i, phase(nzxg, w))
            set_phase!(nzxg, w, zero(P))
        end
        @inbounds rem_edge!(nzxg, w, Outs[i])
    end
    for i = 1:nbits
        for j = i+1:nbits
            @inbounds if has_edge(nzxg, frontier[i], frontier[j])
                if is_hadamard(nzxg, frontier[i], frontier[j])
                    pushfirst_gate!(cir, Val{:CZ}(), i, j)
                    rem_edge!(nzxg, frontier[i], frontier[j])
                end
            end
        end
    end

    old_frontier = copy(frontier)
    max_iter = 1000
    current_iter = 1
    while !isempty(frontier)
        update_frontier!(nzxg, gads, frontier, qubit_map, cir)
        if frontier != old_frontier
            old_frontier = copy(frontier)
            current_iter = 1
        else
            current_iter += 1
            if current_iter > max_iter
                error("Circuit extraction failed!")
            end
        end
    end

    frontier = T[]
    for v in Ins
        nb = neighbors(nzxg, v)
        if length(nb) > 0
            push!(frontier, nb[])
        end
    end
    sort!(frontier, by = (v->qubit_map[v]))
    M = biadjacency(nzxg, frontier, Ins)
    M, steps = gaussian_elimination(M)
    for step in steps
        if step.op == :addto
            ctrl = step.r2
            loc = step.r1
            pushfirst_gate!(cir, Val{:CNOT}(), loc, ctrl)
        elseif step.op == :swap
            q1 = step.r1
            q2 = step.r2
            pushfirst_gate!(cir, Val{:SWAP}(), q1, q2)
        end
    end

    simplify_swap!(cir)
    return cir
end

"""
    update_frontier!(zxg, frontier, qubit_map, cir)

Update frontier. This is an important step in the circuit extraction algorithm.
For more detail, please check the paper [arXiv:1902.03178](https://arxiv.org/abs/1902.03178).
"""
function update_frontier!(zxg::ZXGraph{T, P}, gads::Set{T}, frontier::Vector{T}, qubit_map::Dict{T, Int}, cir) where {T, P}
    # TODO: use inplace methods
    deleteat!(frontier, [spider_type(zxg, f) != SpiderType.Z || (degree(zxg, f)) == 0 for f in frontier])

    for i = 1:length(frontier)
        v = frontier[i]
        nb_v = neighbors(zxg, v)
        u = findfirst([u in gads for u in nb_v])
        if u !== nothing
            u = nb_v[u]
            gad_u = zero(T)
            for w in neighbors(zxg, u)
                if w in gads
                    gad_u = w
                    break
                end
            end
            rewrite!(Rule{:pivot}(), zxg, [u, gad_u, v])
            pop!(gads, u)
            pop!(gads, gad_u)
            frontier[i] = u
            qubit_map[u] = qubit_map[v]
            pushfirst_gate!(cir, Val(:H), qubit_map[u])
            delete!(qubit_map, v)
            for j = 1:length(frontier)
                for k = j+1:length(frontier)
                    if is_hadamard(zxg, frontier[j], frontier[k])
                        pushfirst_gate!(cir, Val(:CZ), qubit_map[frontier[j]], qubit_map[frontier[k]])
                        rem_edge!(zxg, frontier[j], frontier[k])
                    end
                end
            end

            return frontier
        end
    end

    SetN = Set{T}()
    for f in frontier
        union!(SetN, neighbors(zxg, f))
    end
    N = collect(SetN)

    M = biadjacency(zxg, frontier, N)
    M0, steps = gaussian_elimination(M)
    ws = T[]
    @inbounds for i = 1:length(frontier)
        if sum(M0[i,:]) == 1
            push!(ws, N[findfirst(isone, M0[i,:])])
        end
    end
    # M1 = biadjacency(zxg, frontier, ws)
    @inbounds for e in findall(M .== 1)
        if has_edge(zxg, frontier[e[1]], N[e[2]])
            rem_edge!(zxg, frontier[e[1]], N[e[2]])
        end
    end
    @inbounds for e in findall(M0 .== 1)
        add_edge!(zxg, frontier[e[1]], N[e[2]])
    end

    @inbounds for step in steps
        if step.op == :addto
            ctrl = qubit_map[frontier[step.r2]]
            loc = qubit_map[frontier[step.r1]]
            pushfirst_gate!(cir, Val{:CNOT}(), loc, ctrl)
        elseif step.op == :swap
            q1 = qubit_map[frontier[step.r1]]
            q2 = qubit_map[frontier[step.r2]]

            pushfirst_gate!(cir, Val{:SWAP}(), q1, q2)
        end
    end

    for i in 1:length(frontier)
        v = frontier[i]
        if degree(zxg, v) > 1
            continue
        end
        w = neighbors(zxg, v)[1]
        if is_hadamard(zxg, v, w)
            pushfirst_gate!(cir, Val(:H), qubit_map[v])
        end
        if spider_type(zxg, w) == SpiderType.Z
            qubit_map[w] = qubit_map[v]
            if phase(zxg, w) != 0
                pushfirst_gate!(cir, Val{:Rz}(), qubit_map[w], phase(zxg, w))
                set_phase!(zxg, w, zero(P))
            end
            rem_edge!(zxg, v, w)
        else
            rem_edge!(zxg, v, w)
            add_edge!(zxg, v, w, EdgeType.SIM)
        end
        frontier[i] = w
    end

    @inbounds for i1 = 1:length(frontier)
        for i2 = i1+1:length(frontier)
            if has_edge(zxg, frontier[i1], frontier[i2])
                pushfirst_gate!(cir, Val{:CZ}(), qubit_map[frontier[i1]],
                    qubit_map[frontier[i2]])
                rem_edge!(zxg, frontier[i1], frontier[i2])
            end
        end
    end
    return frontier
end

"""
    biadjacency(zxg, F, N)

Return the biadjacency matrix of `zxg` from vertices in `F` to vertices in `N`.
"""
function biadjacency(zxg::ZXGraph{T, P}, F::Vector{T}, N::Vector{T}) where {T, P}
    M = zeros(Int, length(F), length(N))

    for i = 1:length(F)
        for v2 in neighbors(zxg, F[i])
            if v2 in N
                M[i, findfirst(isequal(v2), N)] = 1
            end
        end
    end
    return M
end

"""
    GEStep

A struct for representing steps in the Gaussian elimination.
"""
struct GEStep
    op::Symbol
    r1::Int
    r2::Int
end

"""
    gaussian_elimination(M[, steps])

Return result and steps of Gaussian elimination of matrix `M`. Here we assume
that the elements of `M` is in binary field F_2 = {0,1}.
"""
function gaussian_elimination(M::Matrix{T}, steps::Vector{GEStep} = Vector{GEStep}(); rev = false) where {T<:Integer}
    M = copy(M)
    nr, nc = size(M)
    current_col = 1
    for i = 1:nr
        if sum(M[i,:]) == 0
            continue
        end
        while current_col <= nc
            rs = findall(!iszero, M[i:nr, current_col])
            if length(rs) > 0
                sort!(rs, by = k -> sum(M[k,:]), rev = rev)
                r0 = rs[1]
                r0 += i - 1
                r0 == i && break
                M_r0 = M[r0,:]
                M[r0,:] = M[i,:]
                M[i,:] = M_r0
                step = GEStep(:swap, r0, i)
                push!(steps, step)
                break
            else
                current_col += 1
            end
        end
        current_col > nc && break
        for j = 1:nr
            j == i && continue
            if M[j, current_col] == M[i, current_col]
                M[j,:] = M[j,:] .⊻ M[i,:]
                step = GEStep(:addto, i, j)
                push!(steps, step)
            end
        end
        current_col += 1
    end
    return M, steps
end

function normalize_perm(M::Matrix{T}, steps::Vector{GEStep} = Vector{GEStep}()) where {T<:Integer}
    nr, nc = size(M)
    @assert nc <= nr
    @assert all(sum(M; dims = 1) .<= 1) && all(sum(M; dims = 2) .<= 1)
    @assert sum(M) == nc

    cur_r = 1
    while cur_r <= nc
        cur_c = cur_r
        if M[cur_r, cur_c] != 1
            tar_r = findfirst(isone, M[:, cur_r])
            tar_c = findfirst(isone, M[cur_r, :])
            M[cur_r, cur_c] = 1
            M[tar_r, cur_c] = 0
            if tar_c !== nothing
                M[tar_r, tar_c] = 1
                M[cur_r, tar_c] = 0
            end
            push!(steps, GEStep(:swap, cur_r, tar_r))
        end
        cur_r += 1
    end

    return M, steps
end

# """
#     reverse_gaussian_elimination(M, steps)

# Apply back the operations in `steps` to `M`.
# """
# function reverse_gaussian_elimination(M, steps)
#     for i = length(steps):-1:1
#         s = steps[i]
#         op = s.op
#         r1 = s.r1
#         r2 = s.r2
#         if op == :addto
#             M[r2,:] = M[r2,:] .⊻ M[r1,:]
#         else
#             r_temp = M[r1,:]
#             M[r1,:] = M[r2,:]
#             M[r2,:] = r_temp
#         end
#     end
#     return M
# end


text:
export ess, actime

# Effective sample size (ESS)
actypes = (:bm, :imse, :ipse)

function ess(x::Vector{Float64}; vtype::Symbol=:imse, args...)
  @assert in(vtype, actypes) "Unknown ESS type $vtype"

  return length(x)*mcvar(x; vtype=:iid)/mcvar(x; vtype=vtype, args...)
end

ess(x::Matrix{Float64}, pars::Ranges=1:size(x, 2); vtype::Symbol=:imse, args...) =
  Float64[ess(x[:, pars[i]]; vtype=vtype, args...) for i = 1:length(pars)]

ess(x::Matrix{Float64}, par::Real; vtype::Symbol=:imse, args...) = ess(x, par:par; vtype=vtype, args...)

function ess(c::MCMCChain, pars::Ranges=1:size(c.samples, 2); vtype::Symbol=:imse, args...)
  if isa(c.task.sampler, ARS)
    indx = find(c.diagnostics["accept"])
  else
    indx = 1:size(c.samples, 1)
  end
  ess(c.samples[indx, :], pars; vtype=vtype, args...)
end

ess(c::MCMCChain, par::Real; vtype::Symbol=:imse, args...) = ess(c, par:par; vtype=vtype, args...)

# Integrated autocorrelation time
function actime(x::Vector{Float64}; vtype::Symbol=:imse, args...)
  @assert in(vtype, actypes) "Unknown actime type $vtype"

  return mcvar(x; vtype=vtype, args...)/mcvar(x; vtype=:iid)
end

actime(x::Matrix{Float64}, pars::Ranges=1:size(x, 2); vtype::Symbol=:imse, args...) =
  Float64[actime(x[:, pars[i]]; vtype=vtype, args...) for i = 1:length(pars)]

actime(x::Matrix{Float64}, par::Real; vtype::Symbol=:imse, args...) = actime(x, par:par; vtype=vtype, args...)

function actime(c::MCMCChain, pars::Ranges=1:size(c.samples, 2); vtype::Symbol=:imse, args...)
  if isa(c.task.sampler, ARS)
    indx = find(c.diagnostics["accept"])
  else
    indx = 1:size(c.samples, 1)
  end
  actime(c.samples[indx, :], pars; vtype=vtype, args...)
end

actime(c::MCMCChain, par::Real; vtype::Symbol=:imse, args...) = actime(c, par:par; vtype=vtype, args...)


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the ideal MHD equations
equations = IdealGlmMhdMulticomponentEquations1D(gammas           = (2.0, 2.0),
                                                 gas_constants    = (2.0, 2.0))

"""
    initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdMulticomponentEquations1D)

Compound shock tube test case for one dimensional ideal MHD equations. It is bascially an
MHD extension of the Sod shock tube. Taken from Section V of the article
- Brio and Wu (1988)
  An Upwind Differencing Scheme for the Equations of Ideal Magnetohydrodynamics
  [DOI: 10.1016/0021-9991(88)90120-9](https://doi.org/10.1016/0021-9991(88)90120-9)
"""
function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdMulticomponentEquations1D)
  # domain must be set to [0, 1], γ = 2, final time = 0.12
  if x[1] < 0.5
    rho       = 1.0
    prim_rho  = SVector{ncomponents(equations), real(equations)}(2^(i-1) * (1-2)/(1-2^ncomponents(equations)) * rho for i in eachcomponent(equations))
  else
    rho       = 0.125
    prim_rho  = SVector{ncomponents(equations), real(equations)}(2^(i-1) * (1-2)/(1-2^ncomponents(equations)) * rho for i in eachcomponent(equations))
  end
  v1 = 0.0
  v2 = 0.0
  v3 = 0.0
  p = x[1] < 0.5 ? 1.0 : 0.1
  B1 = 0.75
  B2 = x[1] < 0.5 ? 1.0 : -1.0
  B3 = 0.0

  prim_other = SVector(v1, v2, v3, p, B1, B2, B3)
  return prim2cons(vcat(prim_other, prim_rho), equations)
end
initial_condition = initial_condition_briowu_shock_tube

boundary_conditions = BoundaryConditionDirichlet(initial_condition)

surface_flux = flux_lax_friedrichs
volume_flux  = flux_hindenlang_gassner
basis = LobattoLegendreBasis(3)

indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.8,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = 0
coordinates_max = 1
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                n_cells_max=10_000,
                periodicity=false)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.12)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_errors=(:l2_error_primitive,
                                                            :linf_error_primitive))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_indicator = IndicatorHennemannGassner(semi,
                                          alpha_max=0.5,
                                          alpha_min=0.001,
                                          alpha_smooth=true,
                                          variable=density_pressure)
amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                      base_level=4,
                                      max_level=6, max_threshold=0.01)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=5,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=0.3)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        amr_callback, stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
module TerminalClock

using Dates
using TOML

using Configurations
using Preferences

export clock, stopwatch, countdown

include("types.jl")

default_tomlfile() = joinpath(dirname(@__FILE__), "dials", "ASCII.toml")

function set_dials(tomlfile::AbstractString)
    @set_preferences!("tomlfile" => tomlfile)
    @info("New dials set; restart your Julia session for this change to take effect!")
end

function clear_dials()
    @set_preferences!("tomlfile" => nothing)
    @info("Dials cleared; restart your Julia session for this change to take effect!")
end

function load_dials()
    tomlfile = @load_preference("tomlfile", default_tomlfile())
    toml = tomlfile |> TOML.parsefile
    return toml
end

const toml = load_dials()
const DIALS_LARGE = from_dict(Dials, toml["Large"])
const DIALS_SMALL = from_dict(Dials, get(toml, "Small", toml["Large"]))

colon(dial::Dials) = Dial(getfield(dial, :colon))

const COLON_LARGE = colon(DIALS_LARGE)
const COLON_SMALL = colon(DIALS_SMALL)

function n2d(dials::Dials, n::Int)
    n == 0 && return Dial(getfield(dials, :zero))
    n == 1 && return Dial(getfield(dials, :one))
    n == 2 && return Dial(getfield(dials, :two))
    n == 3 && return Dial(getfield(dials, :three))
    n == 4 && return Dial(getfield(dials, :four))
    n == 5 && return Dial(getfield(dials, :five))
    n == 6 && return Dial(getfield(dials, :six))
    n == 7 && return Dial(getfield(dials, :seven))
    n == 8 && return Dial(getfield(dials, :eight))
    n == 9 && return Dial(getfield(dials, :nine))
    DomainError("n should satisfy 0≤ n ≤9, actual $n")
end

n2d(n::Int, s::Symbol) = n2d(n, Val(s)) # symbol dispatcher
n2d(n::Int, ::Val{:Large}) = n2d(DIALS_LARGE, n)
n2d(n::Int, ::Val{:Small}) = n2d(DIALS_SMALL, n)
n2d(n::Int) = n2d(n, :Large)

function clearline(; move_up::Bool=false)
    buf = IOBuffer()
    print(buf, "\x1b[2K") # clear line
    print(buf, "\x1b[999D") # rollback the cursor
    move_up && print(buf, "\x1b[1A") # move up
    print(buf |> take! |> String)
end

function clearlines(H::Integer)
    for _ in 1:H
        clearline(move_up=true)
    end
end

clock(dt::DateTime) = clock(Time(dt))

function clock(t::Time)
    buf = IOBuffer()
    h = hour(t)
    m = minute(t)
    s = second(t)
    h1, h2 = n2d.(divrem(h, 10))
    m1, m2 = n2d.(divrem(m, 10))
    s1, s2 = n2d.(divrem(s, 10))
    str = hcat(h1, h2, COLON_LARGE, m1, m2, COLON_LARGE, s1, s2).str
    print(buf, str)
    return buf |> take! |> String
end

function stopwatch(t::Time)
    buf = IOBuffer()
    h = hour(t)
    m = minute(t)
    s = second(t)
    ms = millisecond(t)
    h1, h2 = n2d.(divrem(h, 10))
    m1, m2 = n2d.(divrem(m, 10))
    s1, s2 = n2d.(divrem(s, 10))
    ms12, ms3_ = divrem(ms, 10)
    ms3 = n2d.(ms3_, :Small)
    ms1, ms2 = n2d.(divrem(ms12, 10), :Small)

    str = hcat(m1, m2, COLON_LARGE, s1, s2, COLON_SMALL, ms1, ms2, ms3).str
    print(buf, str)
    return buf |> take! |> String
end

const Optional = Union{Nothing,Int}

function setup_timer(;
    hour=nothing::Optional,
    minute=nothing::Optional,
    second=nothing::Optional,
)
    if all(isnothing.([hour, second, minute]))
        hour = 0
        minute = 3
        second = 0
    else
        hour = isnothing(hour) ? 0 : hour
        minute = isnothing(minute) ? 0 : minute
        second = isnothing(second) ? 0 : second
    end
    return Time(hour, minute, second)
end

function countdown(;
    hour=nothing::Optional,
    minute=nothing::Optional,
    second=nothing::Optional,
)
    countdown(setup_timer(; hour, minute, second))
end

function countdown(t::Time)
    while true
        try
            str = clock(t)
            println(str)
            sleep(1)
            if Dates.hour(t) == Dates.minute(t) == Dates.second(t) == 0
                break
            end
            t -= Second(1)
            H = length(split(str, "\n"))
            clearlines(H)
        catch e
            isa(e, InterruptException) || error(e)
            break
        end
    end
end

function stopwatch(duration=0.1::AbstractFloat)
    start = Dates.now()
    sleep(0.001) # warmup
    while true
        try
            Δ = Dates.now() - start # millisecond
            periods =
                Dates.canonicalize(Dates.CompoundPeriod(Dates.Millisecond(Δ))).periods
            str = stopwatch(Time(periods...))
            println(str)
            sleep(duration)
            H = length(split(str, "\n"))
            clearlines(H)
        catch e
            isa(e, InterruptException) || error(e)
            break
        end
    end
end

function clock()
    while true
        try
            str = clock(Dates.now())
            println(str)
            sleep(0.5)
            H = length(split(str, "\n"))
            clearlines(H)
        catch e
            isa(e, InterruptException) || error(e)
            break
        end
    end
end

end # module


text:
CURRENTDIR, CURRENTFILE = splitdir(@__FILE__)
ROOTDIR = splitdir(splitdir(CURRENTDIR)[1])[1]
SRCDIR = joinpath(ROOTDIR, "src")
DATADIR = joinpath(ROOTDIR, "data")

include(joinpath(SRCDIR, "rv_model.jl"))
include(joinpath(SRCDIR, "utils_ex.jl"))

using RvModelKeplerian

srand(314159)

param_true = make_param_true_ex2()

# Make true values and one set of simulated data with noise
num_obs = 50
observation_timespan = 2*365.25
times = observation_timespan*sort(rand(num_obs));
model_true = map(t->calc_model_rv(param_true, t),times);
sigma_obs_scalar = 2.0
sigma_obs = sigma_obs_scalar*ones(num_obs);
jitter = RvModelKeplerian.num_jitters >=1 ? extract_jitter(param_true) : 0.
sigma_eff = sqrt(sigma_obs_scalar^2+jitter^2)

set_times(times);
set_obs( model_true .+ sigma_eff .* randn(length(times)) );
set_sigma_obs(sigma_obs);

writedlm(
  joinpath(DATADIR, "two_planets.csv"),
  zip(RvModelKeplerian.times,RvModelKeplerian.obs,RvModelKeplerian.sigma_obs),
  ','
)

# # Test
#
# using ForwardDiff
#
# function test_rv_example2()
#   plogtarget(param_true)
#   ForwardDiff.gradient(plogtarget,param_true)
#   ForwardDiff.hessian(plogtarget,param_true)
#   result = HessianResult(param_true)
#   ForwardDiff.hessian!(result, plogtarget, param_true);
# end


text:
s = "Hello"

println("$s, World!")

s = """
1
2
3
"""

println(s)

println("1 + 2 = 3" == "1 + 2 = $(1 + 2)")

println(repeat("abc ", 3))

arr = ["I'm", "Happy", "Today"]

println(join(arr, " "))

# regex
regex_pattern = r"^(a|b|c)d$"
println(typeof(regex_pattern))
# occursin is returns true if matched
println(occursin(regex_pattern, "ad"))
println(occursin(regex_pattern, "aad"))

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
import Distributions
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")

@halfwidth_plot_defaults()

plotlyjs()
function p()
ks_data = CSV.read((@__DIR__)*"/../data/errors/ks_error.csv",DataFrame)
ks_data_lwr = CSV.read((@__DIR__)*"/../data/errors/ks_error_lwr.csv",DataFrame)
ks_data_upr = CSV.read((@__DIR__)*"/../data/errors/ks_error_upr.csv",DataFrame)

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
markerstyles_vec = [:cross,:diamond,:circle,:dot]
colours = [1;3;2;4]
for (c,col) in enumerate(names(ks_data))
    if c!=2
        plot!(log10.(1:2:21),log10.(ks_data[:,col]),
            ribbon=(log10.(ks_data[:,col])-log10.(ks_data_lwr[:,col]),log10.(ks_data_upr[:,col])-log10.(ks_data[:,col])),
            label=col,linestyle=linestyles_vec[c],
            # marker=markerstyles_vec[c],
            linewidth=2,color=colours[c])
    end
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error")
# plot!(legend=(0.7,0.92))
plot!()
 @add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/ks_error_formatted.pdf")
end
p()

function p()
l1_cdf_data = CSV.read((@__DIR__)*"/../data/errors/l1_error.csv",DataFrame)
l1_cdf_data_lwr = CSV.read((@__DIR__)*"/../data/errors/l1_error_lwr.csv",DataFrame)
l1_cdf_data_upr = CSV.read((@__DIR__)*"/../data/errors/l1_error_upr.csv",DataFrame)
colours = [1;3;2;4]
plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
markerstyles_vec = [:cross,:diamond,:circle,:cross]
for (c,col) in enumerate(names(l1_cdf_data))
    if c!=2
        plot!(log10.(1:2:21),log10.(l1_cdf_data[:,col]),
            ribbon=(log10.(l1_cdf_data[:,col])-log10.(l1_cdf_data_lwr[:,col]),log10.(l1_cdf_data_upr[:,col])-log10.(l1_cdf_data[:,col])),
            label=col,linestyle=linestyles_vec[c],
            # marker=markerstyles_vec[c],
            linewidth=2,color=colours[c])
    end
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - CDF")
# plot!(legend=(0.225,0.5))
# plot!(ylims=(-2.5,-0.55))
plot!()
 @add_lines!(l1_cdf_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/l1_cdf_error_formatted.pdf")
end
p()

function p()
l1_cell_probs = CSV.read((@__DIR__)*"/../data/errors/l1_cell_probs_error.csv",DataFrame)
l1_cell_probs_lwr = CSV.read((@__DIR__)*"/../data/errors/l1_cell_probs_error_lwr.csv",DataFrame)
l1_cell_probs_upr = CSV.read((@__DIR__)*"/../data/errors/l1_cell_probs_error_upr.csv",DataFrame)
@halfwidth_plot_defaults()
colours = [1;3;2;4]
plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
markerstyles_vec = [:cross,:diamond,:circle,:cross]
for (c,col) in enumerate(names(l1_cell_probs))
    if c!=2
        plot!(log10.(1:2:21),log10.(l1_cell_probs[:,col]),
            ribbon=(log10.(l1_cell_probs[:,col])-log10.(l1_cell_probs_lwr[:,col]),log10.(l1_cell_probs_upr[:,col])-log10.(l1_cell_probs[:,col])),
            label=col,linestyle=linestyles_vec[c],
            # marker=markerstyles_vec[c],
            linewidth=2,color=colours[c])
    end
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - cell probabilities")
# plot!(legend=(0.225,0.5))
plot!()
 @add_lines!(l1_cell_probs,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/l1_cell_probs_error_formatted.pdf")
end
p()

function p()
@fullwidth_plot_defaults()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
dg_data = CSV.read((@__DIR__)*"/../data/first_return_cdf_approximations/order_21_model_dg1.csv",DataFrame)
qbdrap_data = CSV.read((@__DIR__)*"/../data/first_return_cdf_approximations/order_21_model_qbdrap4.csv",DataFrame)
sim_cdf = CSV.read("fluidfluid_discontinuous/data/sims_cdf_evaluated.csv",DataFrame)
plot(dg_data.x,dg_data.phase_4, label="DG",color=1;linewidth=2,linestyle=linestyles_vec[1])
plot!(dg_data.x,qbdrap_data.phase_4, label="QBD-RAP",color=4;linewidth=2,linestyle=linestyles_vec[4])
plot!(dg_data.x,sim_cdf.phase_4, label="Sim",color=5,linewidth=2,linestyle=:dot)
plot!(title="CDF of X at first return in phase 00")
plot!(xlabel="x")
plot!(ylabel="probability")
plot!(legend=(0.15,0.9))
savefig((@__DIR__)*"/phase_4_cdf.pdf")
end
p()


text:
module InlineExports

import Base: @__doc__

export @export

quote
    """
        @export

    Return the expression with all bindings exported.

    ```
    julia> module M
               using InlineExports
               @export begin
                   const a = 2
                   abstract type S <: Number end
                   struct T <: S
                       val
                   end
               end
               @export f(x::TT) where {TT<:S} = x.val^2
           end
    M

    julia> using .M

    julia> f(T(a))
    4
    ```
    """
    macro $(Symbol("export"))(expr::Expr)
        r = handle(expr)
        if r isa Symbol
            return quote
                export $(esc(r))
                @__doc__ $(esc(expr))
            end
        else
            return quote
                export $(map(esc, r)...)
                @__doc__ $(esc(expr))
            end
        end
    end
end |> eval

handle(::Any) = nothing
handle(x::Symbol) = x
handle(expr::Expr) = handle(Val(expr.head), expr)

handle(::Val{:block}, expr) = filter(x -> x !== nothing, map(handle, expr.args))
handle(::Val{:const}, expr) = handle(expr.args[1])
handle(::Val{:(::)}, expr) = handle(expr.args[1])
handle(::Val{:(=)}, expr) = handle(expr.args[1])
handle(::Val{:function}, expr) = handle(expr.args[1])
handle(::Val{:where}, expr) = handle(expr.args[1])
handle(::Val{:macro}, expr) = Symbol("@", handle(expr.args[1]))
handle(::Val{:struct}, expr) = handle(expr.args[2])
handle(::Union{Val{:abstract}, Val{:primitive}}, expr) = handle(expr.args[1])

handle(::Val{:<:}, expr) = handle(expr.args[1])
handle(::Val{:curly}, expr) = handle(expr.args[1])
handle(::Val{:call}, expr) = handle(expr.args[1])
handle(::Val{:macrocall}, expr) = filter(x -> x !== nothing, map(handle, expr.args[3:end]))

include("NoExport.jl")

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# helper function to extract raw data
# from uniform scaling objects
raw(a::UniformScaling) = a.λ
raw(a) = a

"""
    NestedVariogram(cs, γs)

A nested variogram model `γ = c₁γ₁ + c₂γ₂ + ⋯ + cₙγₙ` with
coefficients `cs = (c₁, c₂, ..., cₙ)` and variogram models
`γs = (γ₁, γ₂, ..., γₙ)`.
"""
struct NestedVariogram{CS,GS} <: Variogram
  cs::CS
  γs::GS

  function NestedVariogram{CS,GS}(cs, γs) where {CS,GS}
    @assert all(issymmetric, cs) "coefficients must be symmetric"
    new(cs, γs)
  end
end

NestedVariogram(cs, γs) = NestedVariogram{typeof(cs),typeof(γs)}(cs, γs)

(g::NestedVariogram)(h)                  = raw(sum(g.cs .* map(γ -> γ(h), g.γs)))
(g::NestedVariogram)(u::Point, v::Point) = raw(sum(g.cs .* map(γ -> γ(u, v), g.γs)))

sill(g::NestedVariogram)         = raw(sum(g.cs .* map(sill, g.γs)))
nugget(g::NestedVariogram)       = raw(sum(g.cs .* map(nugget, g.γs)))
Base.range(g::NestedVariogram)   = maximum(range(γ) for γ in g.γs)
isstationary(g::NestedVariogram) = all(isstationary, g.γs)
isisotropic(g::NestedVariogram) = all(isisotropic, g.γs)

# algebraic structure
*(c, γ::Variogram)                          = NestedVariogram((c,), (γ,))
*(c, γ::NestedVariogram)                    = NestedVariogram(map(x->c.*x, γ.cs), γ.γs)
+(γ₁::Variogram, γ₂::Variogram)             = NestedVariogram((1, 1), (γ₁, γ₂))
+(γ₁::NestedVariogram, γ₂::Variogram)       = NestedVariogram((γ₁.cs..., 1), (γ₁.γs..., γ₂))
+(γ₁::Variogram, γ₂::NestedVariogram)       = NestedVariogram((1, γ₂.cs...), (γ₁, γ₂.γs...))
+(γ₁::NestedVariogram, γ₂::NestedVariogram) = NestedVariogram((γ₁.cs..., γ₂.cs...), (γ₁.γs..., γ₂.γs...))

"""
    structures(γ)

Return the individual structures of a (possibly nested)
variogram as a tuple. The structures are the total nugget
`cₒ`, and the coefficients (or contributions) `cs` for the
remaining non-trivial structures `γs` after normalization
(i.e. sill=1, nugget=0).
"""
function structures(γ::Variogram)
  cₒ = nugget(γ)
  c  = sill(γ) - nugget(γ)
  T  = typeof(c)
  γ  = @set γ.sill = one(T)
  γ  = @set γ.nugget = zero(T)
  cₒ, (c,), (γ,)
end

function structures(γ::NestedVariogram)
  ks, gs = γ.cs, γ.γs

  # total nugget and contributions
  cₒ = raw(sum(@. ks * nugget(gs)))
  cs = @. raw(ks * (sill(gs) - nugget(gs)))

  # discard nugget effect terms
  inds = findall(g->!(g isa NuggetEffect), gs)
  cs, γs = cs[inds], gs[inds]

  # adjust sill and nugget
  γs = map(γs) do γ
    T = typeof(sill(γ))
    γ = @set γ.sill = one(T)
    γ = @set γ.nugget = zero(T)
  end

  cₒ, cs, γs
end

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, g::NestedVariogram)
  O = IOContext(io, :compact => true)
  coeffs = 1 .* raw.(g.cs)
  models = nameof.(typeof.(g.γs))
  lines  = ["$c × $γ" for (c, γ) in zip(coeffs, models)]
  print(O, join(lines, " + "))
end

function Base.show(io::IO, ::MIME"text/plain", g::NestedVariogram)
  O = IOContext(io, :compact => true)
  coeffs = 1 .* raw.(g.cs)
  models = g.γs
  nstruc = length(models)
  name   = "NestedVariogram{$nstruc}"
  header = isisotropic(g) ? name : name*" (anisotropic)"
  println(O, header)
  println(O, "  structures")
  lines = ["    └─$γ" for γ in models]
  print(O, join(lines, "\n"))
  println(O)
  println(O, "  coefficients")
  lines = ["    └─$c" for c in coeffs]
  print(O, join(lines, "\n"))
end