text:
export B767FlutterCondition

"""
    B767FlutterCondition

"""
struct B767FlutterCondition <: ControlBenchmarkProblem end

@addbenchmark B767FlutterCondition "ifac06" "Boeing 767" "Boeing 767 at flutter condition" Set([:linear, :unstable, :continuous, :aerospace, :nonnormal, :nonminimum_phase, :ifac])

"""
    controlbenchmark(::B767FlutterCondition) -> sys
"""
function controlbenchmark( ::B767FlutterCondition )

    # The state transition matrix
    #           R    C        Value
    Aparts = [   1    1       0.1015;
                 2    1        19.77;
                 1    2       -19.77;
                 2    2       0.1015;
                 3    3    -0.023202;
                 4    3     0.092543;
                 3    4    -0.092543;
                 4    4    -0.023202;
                 5    5      -0.3165;
                 6    5        14.33;
                 5    6       -14.33;
                 6    6      -0.3165;
                 7    7      -0.3892;
                 8    7        22.29;
                 7    8       -22.29;
                 8    8      -0.3892;
                 9    9      -0.9883;
                10    9        36.16;
                 9   10       -36.16;
                10   10      -0.9883;
                11   11       -1.342;
                12   11        2.457;
                11   12       -2.547;
                12   12       -1.342;
                13   13       -2.312;
                14   13        21.51;
                13   14       -21.51;
                14   14       -2.312;
                15   15       -2.792;
                16   15        26.71;
                15   16       -26.71;
                16   16       -2.792;
                17   17       -2.796;
                18   17        63.79;
                17   18       -63.79;
                18   18       -2.796;
                19   19       -3.332;
                20   19        56.37;
                19   20       -56.37;
                20   20       -3.332;
                21   21       -3.418;
                22   21        88.67;
                21   22       -88.67;
                22   22       -3.418;
                23   23       -4.441;
                24   23        69.45;
                23   24       -69.45;
                24   24       -4.441;
                25   25       -5.108;
                26   25        53.39;
                25   26       -53.39;
                26   26       -5.108;
                27   27       -5.199;
                28   27        50.64;
                27   28       -50.64;
                28   28       -5.199;
                29   29       -5.301;
                30   30        -5.67;
                31   30        92.96;
                30   31       -92.96;
                31   31        -5.67;
                32   32       -6.197;
                33   32         40.0;
                32   33        -40.0;
                33   33       -6.197;
                34   34       -8.177;
                35   34        16.61;
                34   35       -16.61;
                35   35       -8.177;
                36   36       -8.221;
                37   36        139.1;
                36   37       -139.1;
                37   37       -8.221;
                38   38       -10.35;
                39   38        163.0;
                38   39       -163.0;
                39   39       -10.35;
                40   40       -11.72;
                41   40        109.3;
                40   41       -109.3;
                41   41       -11.72;
                42   42       -11.79;
                43   42        304.6;
                42   43       -304.6;
                43   43       -11.79;
                44   44       -33.27;
                45   45       -221.2;
                 1   46      -1970.0;
                 2   46      -2665.0;
                 3   46        586.9;
                 4   46       -501.1;
                 5   46        134.5;
                 6   46       -164.5;
                 7   46        608.5;
                 8   46       -77.54;
                 9   46      -1117.0;
                10   46       1790.0;
                11   46     -13858.0;
                12   46      15010.0;
                13   46     -21697.0;
                14   46      19966.0;
                15   46      -2861.0;
                16   46       6583.0;
                17   46      -2212.0;
                18   46       6321.0;
                19   46      13964.0;
                20   46       2382.0;
                21   46       -303.7;
                22   46        542.8;
                23   46       -527.1;
                24   46       2658.0;
                25   46        25341;
                26   46       6538.0;
                27   46      42048.0;
                28   46        210.6;
                30   46        994.2;
                31   46        769.2;
                32   46      42702.0;
                33   46     -76982.0;
                34   46       -288.8;
                35   46      -2237.0;
                36   46     -16942.0;
                37   46     -12513.0;
                38   46       5748.0;
                39   46     -14657.0;
                40   46      -4089.0;
                41   46       7940.0;
                42   46       9582.0;
                43   46       9474.0;
                48   46    -800000.0;
                 1   47        6.796;
                 2   47        7.196;
                 3   47    0.0017334;
                 4   47    0.0041488;
                 5   47       -0.269;
                 6   47       0.4627;
                 7   47       -1.835;
                 8   47       0.1595;
                 9   47        4.615;
                10   47       -3.752;
                11   47        0.555;
                12   47       -3.138;
                13   47        61.46;
                14   47       -53.56;
                15   47        13.14;
                16   47       -19.32;
                17   47       -12.28;
                18   47       -27.17;
                19   47       -17.24;
                20   47       -20.08;
                21   47        2.624;
                22   47       0.5789;
                23   47       -36.01;
                24   47        3.812;
                25   47       -36.45;
                26   47       -28.56;
                27   47       -59.43;
                28   47       -39.51;
                30   47        13.57;
                31   47         10.6;
                32   47       -120.7;
                33   47        140.6;
                34   47      -0.9428;
                35   47        8.455;
                36   47        94.43;
                37   47        68.02;
                38   47       -35.78;
                39   47        216.8;
                40   47       -258.1;
                41   47        309.0;
                42   47         37.3;
                43   47         41.4;
                46   47          1.0;
                48   47     -60800.0;
                 1   48      -0.1334;
                 2   48      -0.1568;
                 3   48  -4.3119e-05;
                 4   48  -8.0848e-05;
                 5   48    0.0049694;
                 6   48    -0.008432;
                 7   48     0.036978;
                 8   48   -0.0020525;
                 9   48    -0.084286;
                10   48     0.071871;
                11   48    -0.012727;
                12   48     0.065016;
                13   48       -1.286;
                14   48        1.121;
                15   48      -0.2426;
                16   48       0.3841;
                17   48       0.1229;
                18   48       0.2049;
                19   48     0.035851;
                20   48       0.3605;
                21   48      -0.1208;
                22   48       0.1194;
                23   48       0.3072;
                24   48    -0.037101;
                25   48       0.4514;
                26   48       0.3244;
                27   48       0.7558;
                28   48       0.5076;
                30   48    -0.021971;
                31   48    -0.092043;
                32   48         1.94;
                33   48       -2.796;
                34   48     0.023831;
                35   48      -0.1565;
                36   48       -1.144;
                37   48      -0.6821;
                38   48       0.2067;
                39   48      -0.6615;
                40   48       -1.059;
                41   48        1.175;
                42   48    -0.031126;
                43   48    -0.050827;
                47   48          1.0;
                48   48      -1060.0;
                 1   49      -1429.0;
                 2   49      -2222.0;
                 3   49        119.0;
                 4   49       -105.9;
                 5   49       -146.8;
                 6   49        196.0;
                 7   49        175.1;
                 8   49        9.616;
                 9   49       -34.14;
                10   49       -222.7;
                11   49      -2014.0;
                12   49       2668.0;
                13   49      -3652.0;
                14   49       2881.0;
                15   49      -1080.0;
                16   49      -3883.0;
                17   49      -1647.0;
                18   49      -1045.0;
                19   49        562.0;
                20   49        902.6;
                21   49       1129.0;
                22   49      -1258.0;
                23   49      -1267.0;
                24   49       -49.11;
                25   49       4360.0;
                26   49       3790.0;
                27   49       2501.0;
                28   49       -470.4;
                30   49        828.4;
                31   49        810.0;
                32   49      -3077.0;
                33   49       -348.3;
                34   49      -1499.0;
                35   49      -6578.0;
                36   49       -60.89;
                37   49       -29.37;
                38   49       -39.31;
                39   49        123.4;
                40   49        33.03;
                41   49        61.04;
                42   49        140.9;
                43   49        156.3;
                51   49    -800000.0;
                 1   50        18.03;
                 2   50       -2.684;
                 3   50   0.00048611;
                 4   50    0.0010274;
                 5   50      -0.7097;
                 6   50        1.345;
                 7   50       -1.082;
                 8   50        0.637;
                 9   50        6.425;
                10   50        1.899;
                11   50    -0.050439;
                12   50        -1.06;
                13   50        5.074;
                14   50       -26.77;
                15   50        48.49;
                16   50       -12.93;
                17   50        -5.13;
                18   50       -15.59;
                19   50       -25.48;
                20   50       -1.623;
                21   50        4.399;
                22   50       -6.993;
                23   50        7.583;
                24   50       -4.347;
                25   50       -131.3;
                26   50        90.71;
                27   50       -17.39;
                28   50       -18.75;
                30   50         12.9;
                31   50        19.07;
                32   50         35.9;
                33   50       -38.56;
                34   50       0.9758;
                35   50         49.5;
                36   50       0.2221;
                37   50       0.4639;
                38   50       -0.683;
                39   50        9.043;
                40   50        -2.04;
                41   50        2.311;
                42   50        2.627;
                43   50        3.098;
                49   50          1.0;
                51   50     -60800.0;
                 1   51      -0.1517;
                 2   51    -0.050292;
                 3   51  -7.4279e-06;
                 4   51  -1.1816e-05;
                 5   51    0.0012317;
                 6   51   -0.0039645;
                 7   51     0.011391;
                 8   51   -0.0028823;
                 9   51    -0.049071;
                10   51    -0.015192;
                11   51    0.0010263;
                12   51     0.011647;
                13   51      -0.1277;
                14   51        0.254;
                15   51      -0.3712;
                16   51     0.010825;
                17   51     0.020187;
                18   51     0.064588;
                19   51       0.1072;
                20   51     0.077575;
                21   51     -0.02724;
                22   51     0.029297;
                23   51   -0.0065626;
                24   51   -0.0084969;
                25   51       0.8394;
                26   51      -0.3495;
                27   51       0.1113;
                28   51       0.1395;
                30   51    -0.036967;
                31   51    -0.081012;
                32   51      -0.2467;
                33   51       0.1981;
                34   51    -0.041943;
                35   51       -0.498;
                36   51    -0.011955;
                37   51   -0.0099611;
                38   51    0.0078189;
                39   51    -0.087434;
                40   51     0.023536;
                41   51    -0.028729;
                42   51    -0.018542;
                43   51    -0.022171;
                50   51          1.0;
                51   51      -1060.0;
                53   52          1.0;
                 1   53       -3.158;
                 2   53        1.779;
                 3   53    0.0053132;
                 4   53    0.0080111;
                 5   53       0.1066;
                 6   53      -0.4078;
                 7   53       0.1001;
                 8   53      -0.1571;
                 9   53        -1.06;
                10   53       0.4854;
                11   53       -2.733;
                12   53        2.067;
                13   53        3.001;
                14   53        3.544;
                15   53       -7.199;
                16   53        7.207;
                17   53       0.6186;
                18   53        1.076;
                19   53        2.821;
                20   53      -0.7866;
                21   53      -0.1111;
                22   53        1.509;
                23   53      -0.3659;
                24   53       0.8635;
                25   53        7.433;
                26   53       -15.36;
                27   53         4.57;
                28   53      -0.2771;
                29   53       -1.251;
                30   53      -0.5031;
                31   53      -0.4524;
                32   53       -1.592;
                33   53        1.956;
                34   53       -5.981;
                35   53       -10.88;
                36   53       -1.121;
                37   53       0.1532;
                38   53      -0.4254;
                39   53      -0.7886;
                40   53      -0.1053;
                41   53       -2.007;
                42   53       0.1886;
                43   53    0.0092049;
                44   53       -17.07;
                45   53       -37.83;
                52   53      -0.2668;
                53   53       -1.033;
                48   54  -16000000.0;
                54   54        -20.0;
                51   55  -16000000.0;
                55   55        -20.0 ]

    A = sparse( Aparts[:, 1], Aparts[:, 2], Aparts[:, 3], 55, 55 )

    # The disturbance input mapping matrix
    #           R    C        Value
    B1parts = [ 52   1         0.3713;
                53   1          1.245;
                48   2       800000.0;
                54   2            1.0;
                51   3       800000.0;
                55   3            1.0 ]

    B1 = sparse( B1parts[:, 1], B1parts[:, 2], B1parts[:, 3], 55, 3 )

    # The control input mapping matrix
    #           R    C        Value
    B2parts = [ 48   1       800000.0;
                51   2       800000.0 ]

    B2 = sparse( B2parts[:, 1], B2parts[:, 2], B2parts[:, 3], 55, 2 )

    # The output measurement matrix
    #           R    C        Value
    C1parts = [  1   1     4.4247e-05;
                 2   1          35.92;
                 1   2     4.3403e-05;
                 2   2         -12.46;
                 1   3     4.9713e-05;
                 2   3        -0.5332;
                 1   4     4.5556e-05;
                 2   4        -0.5659;
                 1   5    -4.1311e-06;
                 2   5        -0.9664;
                 1   6     7.9201e-06;
                 2   6          12.75;
                 1   7     7.2496e-06;
                 2   7         -21.45;
                 1   8    -6.7385e-05;
                 2   8          6.822;
                 1   9    -0.00018236;
                 2   9          33.88;
                 1   10    2.0657e-05;
                 2   10            17;
                 1   11    0.00013194;
                 2   11        -0.138;
                 1   12   -0.00015849;
                 2   12       -0.8803;
                 1   13   -6.7655e-05;
                 2   13         5.465;
                 1   14   -8.7358e-06;
                 2   14        -1.614;
                 1   15   -9.8119e-05;
                 2   15         55.18;
                 1   16    3.1119e-05;
                 2   16          30.1;
                 1   17   -5.7479e-05;
                 2   17        -4.878;
                 1   18      0.000316;
                 2   18         14.36;
                 1   19    -0.0005775;
                 2   19         62.54;
                 1   20       0.00071;
                 2   20         55.18;
                 1   21   -9.6843e-05;
                 2   21        -138.8;
                 1   22    2.9934e-05;
                 2   22        -104.4;
                 1   23    0.00054014;
                 2   23         6.831;
                 1   24    -0.0004762;
                 2   24        -11.78;
                 1   25    -4.163e-06;
                 2   25        -55.04;
                 1   26   -0.00021609;
                 2   26        -30.66;
                 1   27   -0.00010541;
                 2   27         1.402;
                 1   28   -3.1883e-05;
                 2   28        -30.33;
                 1   29   -0.00015062;
                 2   29         4.869;
                 1   30   -0.00027714;
                 2   30        -74.28;
                 1   31   -0.00024361;
                 2   31         74.14;
                 1   32     2.799e-05;
                 2   32          16.3;
                 1   33   -0.00016592;
                 2   33         3.118;
                 1   34    1.2748e-05;
                 2   34         8.728;
                 1   35   -2.0762e-05;
                 2   35        -74.53;
                 1   36   -0.00012546;
                 2   36          1.46;
                 1   37    -4.436e-05;
                 2   37        -3.334;
                 1   38   -1.0147e-05;
                 2   38        -7.367;
                 1   39    0.00015343;
                 2   39        -1.032;
                 1   40    2.1956e-05;
                 2   40         2.867;
                 1   41   -2.9386e-05;
                 2   41         6.583;
                 1   42    6.5513e-05;
                 2   42        -1.361;
                 1   43    6.1813e-05;
                 2   43          1.35;
                 1   44    1.9995e-05;
                 2   44         10.71;
                 1   45    8.9674e-06;
                 2   45        -25.91;
                 2   46         48.52;
                 2   47        -7.932;
                 2   48      0.079915;
                 2   49        -25618;
                 2   50        -167.4;
                 2   51        0.2213;
                 2   53          1.88 ]

    C1 = sparse( C1parts[:, 1], C1parts[:, 2], C1parts[:, 3], 2, 55 )

    # The output performance matrix (regulated)
    #           R    C        Value
    C2parts = [  1   1      -0.00196156517750494;
                 2   1       -0.0129504190712116;
                 3   1      -0.00576793734796071;
                 1   2      -0.00178666088824936;
                 2   2      -0.00659743772384401;
                 3   2       0.00212717146699555;
                 1   3       0.00383083053052468;
                 2   3         0.000520641604177;
                 3   3      0.000855251280034119;
                 1   4       0.00417773050112139;
                 2   4      0.000964484966186617;
                 3   4      0.000948383511296986;
                 1   5       0.00271334336979307;
                 2   5       0.00603316894923389;
                 3   5      0.000112181294363187;
                 1   6       0.00343099389588498;
                 2   6       0.00288047491396818;
                 3   6      -0.00240891603216052;
                 1   7     -0.000578619067781213;
                 2   7        0.0125202768040487;
                 3   7       0.00399138133983713;
                 1   8     -0.000744991440487741;
                 2   8       0.00414241786038058;
                 3   8     -0.000938835501033062;
                 1   9      -0.00335052241025187;
                 2   9        0.0012720384899837;
                 3   9        0.0032244100235547;
                 1  10      0.000254794053305802;
                 2  10       0.00330629051960048;
                 3  10       0.00282214139440248;
                 1  11       0.00638400655782871;
                 2  11       0.00138777668592609;
                 3  11       0.00123153679928778;
                 1  12       0.00655735958952992;
                 2  12       0.00137912336286498;
                 3  12       0.00142672317304374;
                 1  13     -0.000209807584228979;
                 2  13      -0.00121471022470382;
                 3  13      -0.00140825325154959;
                 1  14      0.000596458529656505;
                 2  14      -0.00375806609441612;
                 3  14      -0.00084585979452862;
                 1  15       -0.0075373665495052;
                 2  15      -0.00171696351737595;
                 3  15      -0.00223360830272454;
                 1  16      -0.00707644914866206;
                 2  16      -0.00037245344675543;
                 3  16      -0.00256700603816976;
                 1  17      -0.00136471883345984;
                 2  17     -3.18478344161734e-07;
                 3  17      -0.00108487310064357;
                 1  18       0.00216827024662518;
                 2  18       0.00166540413413682;
                 3  18       0.00532966802437075;
                 1  19       0.00155164536876182;
                 2  19       0.00124319407977998;
                 3  19       0.00306005902720846;
                 1  20       0.00200693946097036;
                 2  20      0.000469803331192958;
                 3  20      0.000830833417719822;
                 1  21       -0.0111834035238652;
                 2  21      -0.00587272191747575;
                 3  21      -0.00193934175688557;
                 1  22       -0.0113474490102401;
                 2  22      0.000330773274011066;
                 3  22       -0.0026781386166515;
                 1  23      -0.00187663382640301;
                 2  23      -0.00235081943160252;
                 3  23      -0.00835372635714146;
                 1  24      -0.00340733721818079;
                 2  24      -0.00159365366375508;
                 3  24      -0.00927409193668038;
                 1  25     -0.000802000155610957;
                 2  25        0.0050308256946549;
                 3  25        0.0121838871958009;
                 1  26       -9.559624678825e-05;
                 2  26       0.00394267031971987;
                 3  26       0.00977027542088758;
                 1  27      2.96096285690989e-05;
                 2  27      -0.00467676055940434;
                 3  27        -0.013112079013261;
                 1  28     -0.000463632101563298;
                 2  28      0.000816657363892594;
                 3  28      -0.00333397735445219;
                 1  29       -0.0104723458423793;
                 2  29      -0.00241463768917823;
                 3  29       -0.0032118880428807;
                 1  30       0.00272730468778243;
                 2  30       -0.0239635754421163;
                 3  30        0.0131856456497208;
                 1  31      -0.00352988656497627;
                 2  31        0.0215622782926573;
                 3  31       -0.0125736338442791;
                 1  32      -0.00081809445273758;
                 2  32      -0.00171227630071785;
                 3  32      0.000249813514446276;
                 1  33     -0.000313547933177688;
                 2  33      -0.00127095682460106;
                 3  33       -0.0015855958028451;
                 1  34       0.00815476705614575;
                 2  34      -0.00221669292415526;
                 3  34      -0.00684639293350886;
                 1  35       0.00651722080031051;
                 2  35       0.00897421712462987;
                 3  35        0.0152032497858188;
                 1  36      1.85821264294483e-05;
                 2  36      4.96123855503845e-05;
                 3  36     -6.51769094081639e-05;
                 1  37     -0.000102828985135515;
                 2  37     -0.000575085428436506;
                 3  37     -1.22026701668119e-05;
                 1  38     -8.21972596623513e-05;
                 2  38     -0.000689381403868715;
                 3  38      3.14771289192645e-05;
                 1  39     -1.56037119250517e-05;
                 2  39     -7.34090239684468e-05;
                 3  39      3.93503242680413e-05;
                 1  40      0.000112563126289207;
                 2  40      0.000877230625320389;
                 3  40      8.89843251646041e-05;
                 1  41      0.000129413661473586;
                 2  41       0.00112997976973041;
                 3  41     -0.000317588734844295;
                 1  42     -3.99254913056809e-06;
                 2  42     -6.15107047594157e-05;
                 3  42      2.52474435339501e-06;
                 1  43      5.50696431802495e-06;
                 2  43      7.29763578153911e-05;
                 3  43      2.33691464328503e-06;
                 1  44      -0.00630528023694427;
                 2  44      -0.00523489989684617;
                 3  44       -0.0116673554929984;
                 1  45     -0.000260805176328999;
                 2  45     -0.000249215704160071;
                 3  45     -0.000737544661698531 ]

    C2 = sparse( C2parts[:, 1], C2parts[:, 2], C2parts[:, 3], 5, 55 )

    # Feedforward term for the disturbance to the regulated output
    D11 = spzeros( 5, 3 )

    # Feedforward term for the control inputs to the regulated output
    D12 = spzeros( 5, 2 )

    # Feedforward term for the disturbance to the measurments
    D21 = spzeros( 2, 3 )

    # Feedforward term for the control inputs to the measurments
    D22 = spzeros( 2, 2 )

    return (; :sys => HeteroStateSpace( A, B2, C1, D22, Continuous() ) )
end


text:
# export evalElliptic, init, testOperatorGradientInterpolation, dataPrep, interpolateFace

import PDESolver.evalResidual

function init(mesh::AbstractMesh{Tmsh},
              sbp::AbstractOperator,
              eqn::AbstractEllipticData{Tsol, Tres},
              opts,
              pmesh=mesh) where {Tmsh, Tsol, Tres}
  #
  # TODO: initMPIStructures(mesh, opts)
  #
  getBCFunctors(mesh, sbp, eqn, opts)
  getBCFunctors(pmesh, sbp, eqn, opts)
  # getFluxFunctors(mesh, sbp, eqn, opts)
  getSrcFuntors(mesh, sbp, eqn, opts)
  getDiffnFunc(mesh, sbp, eqn, opts)
  if haskey(opts, "Functional")
    getFunctional(mesh, sbp, eqn, opts)
  end

  calcDiffn(mesh, sbp, eqn, eqn.diffusion_func, eqn.lambda)
  dim = size(mesh.coords, 1)

  # interpolate diffusion from elements to faces
  for d2 = 1 : dim
    for d1 = 1 : dim
      lambda_elem = Base.view(eqn.lambda, d1, d2, :, :, :)
      lambda_face = Base.view(eqn.lambda_face, d1, d2, :, :, :, :)
      lambda_bndry = Base.view(eqn.lambda_bndry, d1, d2, :, :, :)
      interpolateFace(mesh, sbp, eqn, opts, lambda_elem, lambda_face)
      interpolateBoundary(mesh, sbp, eqn, opts, lambda_elem, lambda_bndry)

      # lambda_elem[:,:,:] = eqn.lambda[d1, d2, :,:,:]
      # interpolateFace(mesh, sbp, eqn, opts, lambda_elem, lambda_face)
      # interpolateBoundary(mesh, sbp, eqn, opts, lambda_elem, lambda_bndry)
      # eqn.lambda_face[d1, d2, :,:,:,:] = lambda_face[:,:,:,:]
      # eqn.lambda_bndry[d1, d2, :,:,:] = lambda_bndry[:,:,:]
    end
  end
  return nothing
end

function evalResidual(mesh::AbstractMesh,
                      sbp::AbstractOperator,
                      eqn::EllipticData,
                      opts::Dict,
                      t=0.0)
  dataPrep(mesh, sbp, eqn, opts)
  evalVolumeIntegrals(mesh, sbp, eqn, opts)
  evalFaceIntegrals(mesh, sbp, eqn, opts)
  evalBoundaryIntegrals(mesh, sbp, eqn)

  for e = 1:mesh.numEl
    for n = 1:mesh.numNodesPerElement
      for v = 1:mesh.numDofPerNode
        # eqn.res[v,n,e] = -1.0*eqn.res[v,n,e]
      end
    end
  end

  # if haskey(opts, "TimeAdvance") && opts["TimeAdvance"] == "BDF2"
  # BDF2(mesh, sbp, eqn, opts)
  # end
  #
  # TODO: parallel commmunication
  #
  # if mesh.commsize > 1
  #     evalSharedFaceIntegrals(mesh, sbp, eqn, opts)
  # end if
  return nothing
end

function dataPrep(mesh::AbstractMesh{Tmsh},
                         sbp::AbstractOperator,
                         eqn::AbstractEllipticData{Tsol, Tres},
                         opts) where {Tmsh, Tsol, Tres}
  #
  # TODO: Not sure if filtering is needed
  #

  fill!(eqn.res, 0.0)
  # fill!(eqn.res_edge, 0.0)

  # source term is multiplied 1/|J|
  calcSource(mesh, sbp, eqn, eqn.src_func, eqn.src)
  # calcDiffn(mesh, sbp, eqn, eqn.diffusion_func, eqn.lambda)
  calcGradient(mesh, sbp, eqn, eqn.q, eqn.q_grad)

  fill!(eqn.q_face, 0.0)
  fill!(eqn.q_grad_face, 0.0)
  fill!(eqn.q_bndry, 0.0)
  fill!(eqn.q_grad_bndry, 0.0)

  # interpolate diffusion from elements to faces
  # dim = size(mesh.coords, 1)
  # lambda_elem = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
  # lambda_face = Array{Tsol}(mesh.numDofPerNode, 2, mesh.numNodesPerFace, mesh.numInterfaces)
  # lambda_bndry = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numBoundaryFaces)

  # for d2 = 1 : dim
  # for d1 = 1 : dim
  # lambda_elem = sview(eqn.lambda, d1, d2, :, :, :)
  # lambda_face = sview(eqn.lambda_face, d1, d2, :, :, :, :)
  # lambda_bndry = sview(eqn.lambda_bndry, d1, d2, :, :, :)
  # interpolateFace(mesh, sbp, eqn, opts, lambda_elem, lambda_face)
  # interpolateBoundary(mesh, sbp, eqn, opts, lambda_elem, lambda_bndry)

  # # lambda_elem[:,:,:] = eqn.lambda[d1, d2, :,:,:]
  # # interpolateFace(mesh, sbp, eqn, opts, lambda_elem, lambda_face)
  # # interpolateBoundary(mesh, sbp, eqn, opts, lambda_elem, lambda_bndry)
  # # eqn.lambda_face[d1, d2, :,:,:,:] = lambda_face[:,:,:,:]
  # # eqn.lambda_bndry[d1, d2, :,:,:] = lambda_bndry[:,:,:]
  # end
  # end

  # interpolate q from volume to interior edge and boundary edge
  interpolateBoundary(mesh, sbp, eqn, opts, eqn.q, eqn.q_bndry)
  interpolateFace(mesh, sbp, eqn, opts, eqn.q, eqn.q_face)
  # interpolate grad(q) from volume to interior edge and boundary edge
  interpolateBoundary(mesh, sbp, eqn, opts, eqn.q_grad, eqn.q_grad_bndry)
  interpolateFace(mesh, sbp, eqn, opts, eqn.q_grad, eqn.q_grad_face)

  fill!(eqn.xflux_face, 0.0)
  fill!(eqn.yflux_face, 0.0)
  fill!(eqn.flux_face,  0.0)
  fill!(eqn.xflux_bndry, 0.0)
  fill!(eqn.yflux_bndry, 0.0)
  fill!(eqn.flux_bndry,  0.0)
  # eqn.flux_func = FluxType["SIPG"]
  calcFaceFlux(mesh, sbp, eqn, opts, mesh.interfaces,
               eqn.xflux_face, eqn.yflux_face, eqn.flux_face)
  getBCFluxes(mesh, sbp, eqn, opts)

  return nothing
end

function evalVolumeIntegrals(mesh::AbstractMesh{Tmsh},
                             sbp::AbstractOperator,
                             eqn::EllipticData{Tsol, Tres, Tdim},
                             opts) where {Tmsh, Tsol, Tres, Tdim}
  #
  # integral of term-1
  #
  weakdifferentiate2!(mesh, sbp, eqn, eqn.q_grad, eqn.res)
  # integral of term-6, i.e., source term
  strongintegrate(mesh, sbp, eqn.src, eqn.res)
end

function strongintegrate(mesh::AbstractMesh{Tmsh},
                         sbp::AbstractOperator{Tsbp},
                         f::AbstractArray{Tflx,3},
                         res::AbstractArray{Tres,3}) where {Tmsh, Tsbp, Tflx, Tres}
  for elem = 1:mesh.numEl
    jac = sview(mesh.jac, :, elem)
    for i = 1:mesh.numNodesPerElement
      res[:, i, elem] -= sbp.w[i]/jac[i]*f[:, i, elem]
    end
  end

end
#
# Volume integration of term looking like:
#     \int ∇v⋅(Λ∇f) dΩ
# where Λ is matrix coefficient
#
function weakdifferentiate2!(mesh::AbstractMesh{Tmsh},
                                    sbp::AbstractOperator{Tsbp},
                                    eqn::EllipticData{Tsol, Tres, Tdim},
                                    q_grad::AbstractArray{Tflx,4},
                                    res::AbstractArray{Tres,3}) where {Tmsh, Tsbp,Tflx,Tsol, Tres, Tdim}
  @assert (sbp.numnodes == size(q_grad, 2))
  @assert (sbp.numnodes == size(res,2))
  @assert (size(q_grad, 1) == size(res,1))    # numDofPerNode
  @assert (size(mesh.dxidx, 1) == size(q_grad, 4))
  dim           = size(q_grad, 4)
  nElems        = size(q_grad, 3)
  nNodesPerElem = size(q_grad, 2)
  nDofsPerNode  = size(q_grad, 1)
  Qx = Array{Tsbp}(nNodesPerElem, nNodesPerElem, dim)
  lambda_dqdx = Array{Tres}(dim, nDofsPerNode, nNodesPerElem)
  w = sview(sbp.w, :)
  for elem=1:nElems
    # compute Λ∇q
    lambda = sview(eqn.lambda, :,:,:,:, elem)
    lambda_dqdx[:,:,:] = 0.0
    for d1 = 1:dim
      for d2 = 1:dim
        for node = 1:nNodesPerElem
          for ivar = 1:nDofsPerNode
            lambda_dqdx[d1, ivar, node] += lambda[d1, d2, ivar, node]*q_grad[ivar, node, elem, d2]
          end
        end
      end
    end
    jac = sview(mesh.jac, :, elem)
    calcQx(mesh, sbp, Int(elem), Qx)
    for d = 1:dim
      for i = 1:sbp.numnodes
        for j = 1:sbp.numnodes
          for ivar = 1:nDofsPerNode
            # res[ivar, i, elem] += Qx[j,i,d]*q_grad[ivar, j, elem, d]/jac[j]
            res[ivar, i, elem] += Qx[j,i,d]*lambda_dqdx[d, ivar, j]/jac[j]
          end
        end
      end
    end
  end
end

#
# Integration over interior interfaces. The integration is categorized into 2
# parts.
#
# TODO: put the integral of ∫∇v ⋅ (fx, fy) dΓ into a function,
#       let's say interiorfaceintegrate2
#
function evalFaceIntegrals(
                                  mesh::AbstractDGMesh{Tmsh},
                                  sbp::AbstractOperator,
                                  eqn::EllipticData{Tsol, Tres, Tdim},
                                  opts) where {Tmsh, Tsol, Tres, Tdim}

  #
  # Integrate face_flux
  #
  interiorfaceintegrate!(mesh.sbpface, mesh.interfaces, eqn.flux_face, eqn.res)

  #
  # Integrate face_flux
  #
  sbpface = mesh.sbpface
  DxL = Array{Tmsh}(mesh.numNodesPerElement, mesh.numNodesPerElement, Tdim)
  DxR = Array{Tmsh}(mesh.numNodesPerElement, mesh.numNodesPerElement, Tdim)

  GxL = Array{Tmsh}(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numNodesPerElement, Tdim)
  GxR = Array{Tmsh}(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numNodesPerElement, Tdim)
  # interpolation matrix, keep in mind R is transfered for efficiency
  R = sview(sbpface.interp[:,:])
  w = sview(sbpface.wface, :)
  res = sview(eqn.res, :,:,:)

  nv = mesh.numNodesPerElement    # number of Nodes per elemet
  nf = mesh.numNodesPerFace       # number of nodes on interfaces
  ns = sbpface.stencilsize        # size of stencil for interpolation

  # loop over all the boundaries
  nfaces = length(mesh.interfaces)
  for f = 1:nfaces
    face = mesh.interfaces[f]

    xflux = sview(eqn.xflux_face, :,:,f)
    yflux = sview(eqn.yflux_face, :,:,f)
    elemL = face.elementL
    elemR = face.elementR
    faceL = face.faceL
    faceR = face.faceR
    pL = sview(sbpface.perm, :, faceL)
    pR = sview(sbpface.perm, :, faceR)

    calcDx(mesh, sbp, Int(elemL), DxL)
    calcDx(mesh, sbp, Int(elemR), DxR)
    #
    # Compute operator G = Λ∇
    #
    lambdaL = sview(eqn.lambda, :, :, :, :, elemL)
    lambdaR = sview(eqn.lambda, :, :, :, :, elemR)
    GxL[:,:,:,:] = 0.0
    GxR[:,:,:,:] = 0.0
    for d = 1:Tdim
      for i = 1:mesh.numNodesPerElement
        for j  = 1:mesh.numNodesPerElement
          GxL[:, i, j, d] = lambdaL[d, 1, :, i]*DxL[i,j,1] + lambdaL[d, 2, :, i]*DxL[i,j,2]
          GxR[:, i, j, d] = lambdaR[d, 1, :, i]*DxR[i,j,1] + lambdaR[d, 2, :, i]*DxR[i,j,2]
        end
      end
    end

    for i = 1 : nv
      for j = 1 : ns
        for k = 1 : nf
          kR = sbpface.nbrperm[k, face.orient]
          for ivar = 1:mesh.numDofPerNode
            res[ivar, i, elemL] += GxL[ivar, pL[j], i, 1]*R[j, k]*w[k]*xflux[ivar, k]
            res[ivar, i, elemL] += GxL[ivar, pL[j], i, 2]*R[j, k]*w[k]*yflux[ivar, k]
            res[ivar, i, elemR] += GxR[ivar, pR[j], i, 1]*R[j,kR]*w[k]*xflux[ivar, k]
            res[ivar, i, elemR] += GxR[ivar, pR[j], i, 2]*R[j,kR]*w[k]*yflux[ivar, k]
          end
        end
      end
    end

  end
  return nothing
end

function evalBoundaryIntegrals(mesh::AbstractMesh{Tmsh},
                               sbp::AbstractOperator,
                               eqn::EllipticData{Tsol, Tres, Tdim}) where {Tmsh, Tsol, Tres, Tdim}
  boundaryintegrate!(mesh.sbpface, mesh.bndryfaces, eqn.flux_bndry, eqn.res)

  sbpface = mesh.sbpface
  Dx = Array{Tmsh}((mesh.numNodesPerElement, mesh.numNodesPerElement, Tdim))
  # interpolation matrix, keep in mind R is transfered for efficiency
  R = sview(sbpface.interp[:,:])
  w = sview(sbpface.wface, :)
  res = sview(eqn.res, :,:,:)

  nv = mesh.numNodesPerElement
  nf = mesh.numNodesPerFace
  ns = sbpface.stencilsize
  Gx = Array{Tmsh}((mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numNodesPerElement, Tdim))
  # loop over all the boundaries
  for bc = 1:mesh.numBC
    indx0 = mesh.bndry_offsets[bc]
    indx1 = mesh.bndry_offsets[bc+1] - 1
    for f = indx0:indx1
      xflux = sview(eqn.xflux_bndry, :,:,f)
      yflux = sview(eqn.yflux_bndry, :,:,f)
      bndry = mesh.bndryfaces[f]
      elem = bndry.element
      face = bndry.face
      p = sview(sbpface.perm, :, face)

      calcDx(mesh, sbp, Int(elem), Dx)

      lambda = sview(eqn.lambda, :, :, :, :, elem)
      fill!(Gx, 0.0)
      for d = 1:Tdim
        for i = 1:mesh.numNodesPerElement
          for j  = 1:mesh.numNodesPerElement
            Gx[:, i, j, d] = lambda[d, 1, :, i]*Dx[i,j,1] + lambda[d, 2, :, i]*Dx[i,j,2]
          end
        end
      end

      for i = 1 : nv      # loop over volume nodes
        for j = 1 : ns      # loop over stencils
          for k = 1 : nf      # loop over face nodes
            for ivar = 1:mesh.numDofPerNode
              res[ivar, i, elem] += Gx[ivar, p[j], i, 1]*R[j ,k]*w[k]*xflux[ivar, k]
              res[ivar, i, elem] += Gx[ivar, p[j], i, 2]*R[j ,k]*w[k]*yflux[ivar, k]
            end
          end
        end
      end
    end
  end

  return nothing
end  # end evalBoundaryIntegrals
#
# Compute Dx
#
function calcDx(mesh::AbstractMesh{Tmsh},
                sbp::AbstractOperator,
                elem::Int,
                Dx::AbstractArray{Tmsh, 3}) where Tmsh
  @assert(size(Dx, 1) == mesh.numNodesPerElement)
  @assert(size(Dx, 2) == mesh.numNodesPerElement)
  @assert(size(Dx, 3) == size(mesh.dxidx, 1))
  dxidx = sview(mesh.dxidx, :,:,:,elem) # (dim, dim, numNodesPerElement)
  jac = mesh.jac[:, elem]
  Dx[:,:,:] = 0.0
  dim = size(Dx, 3)
  for d=1:dim
    for dd=1:dim
      for l=1:mesh.numNodesPerElement
        for m = 1:mesh.numNodesPerElement
          # Since dxidx is scaled by 1/|J|, we need to get it back,
          # that's why jac is here
          Dx[l,m,d] += dxidx[dd,d,l]*jac[l]*sbp.Q[l,m,dd]
        end
      end
    end
    # Until here Dx stores Qx, we need to left multiply H^(-1)
    for l=1:mesh.numNodesPerElement
      Dx[l,:,d] /= sbp.w[l]
    end
  end
end

function calcQx(mesh::AbstractMesh{Tmsh},
                sbp::AbstractOperator,
                elem::Int,
                Qx::AbstractArray{Tmsh, 3}) where Tmsh
  @assert(size(Qx, 1) == mesh.numNodesPerElement)
  @assert(size(Qx, 2) == mesh.numNodesPerElement)
  @assert(size(Qx, 3) == size(mesh.dxidx, 1))
  dxidx = sview(mesh.dxidx, :,:,:,elem) # (dim, dim, numNodesPerElement)
  jac = mesh.jac[:, elem]
  Qx[:,:,:] = 0.0
  dim = size(Qx, 3)
  for d=1:dim
    for dd=1:dim
      for l=1:mesh.numNodesPerElement
        # Since dxidx is scaled by 1/|J|, we need to get it back,
        # that's why jac is here
        for m = 1:mesh.numNodesPerElement
          Qx[l,m,d] += dxidx[dd,d,l]*jac[l]*sbp.Q[l,m,dd]
        end
      end
    end
  end
end

ext: jl


text:
using Base.Broadcast: Broadcasted

using CUDA: CuArrayStyle

# TODO: support more types (SparseVector, SparseMatrixCSC, COO, BSR)

## sparse broadcast style

# broadcast container type promotion for combinations of sparse arrays and other types
struct CuSparseVecStyle <: Broadcast.AbstractArrayStyle{1} end
struct CuSparseMatStyle <: Broadcast.AbstractArrayStyle{2} end
Broadcast.BroadcastStyle(::Type{<:CuSparseVector}) = CuSparseVecStyle()
Broadcast.BroadcastStyle(::Type{<:CuSparseMatrix}) = CuSparseMatStyle()
const SPVM = Union{CuSparseVecStyle,CuSparseMatStyle}

# CuSparseVecStyle handles 0-1 dimensions, CuSparseMatStyle 0-2 dimensions.
# CuSparseVecStyle promotes to CuSparseMatStyle for 2 dimensions.
# Fall back to DefaultArrayStyle for higher dimensionality.
CuSparseVecStyle(::Val{0}) = CuSparseVecStyle()
CuSparseVecStyle(::Val{1}) = CuSparseVecStyle()
CuSparseVecStyle(::Val{2}) = CuSparseMatStyle()
CuSparseVecStyle(::Val{N}) where N = Broadcast.DefaultArrayStyle{N}()
CuSparseMatStyle(::Val{0}) = CuSparseMatStyle()
CuSparseMatStyle(::Val{1}) = CuSparseMatStyle()
CuSparseMatStyle(::Val{2}) = CuSparseMatStyle()
CuSparseMatStyle(::Val{N}) where N = Broadcast.DefaultArrayStyle{N}()

Broadcast.BroadcastStyle(::CuSparseVecStyle, ::CuArrayStyle{1}) = CuSparseVecStyle()
Broadcast.BroadcastStyle(::CuSparseVecStyle, ::CuArrayStyle{2}) = CuSparseMatStyle()
Broadcast.BroadcastStyle(::CuSparseMatStyle, ::CuArrayStyle{2}) = CuSparseMatStyle()

# don't wrap sparse arrays with Extruded
Broadcast.extrude(x::CuSparseVecOrMat) = x

## detection of zero-preserving functions

# modified from SparseArrays.jl

# capturescalars takes a function (f) and a tuple of broadcast arguments, and returns a
# partially-evaluated function and a reduced argument tuple where all scalar operations have
# been applied already.
@inline function capturescalars(f, mixedargs)
    let (passedsrcargstup, makeargs) = _capturescalars(mixedargs...)
        parevalf = (passed...) -> f(makeargs(passed...)...)
        return (parevalf, passedsrcargstup)
    end
end
# Work around losing Type{T}s as DataTypes within the tuple that makeargs creates
@inline capturescalars(f, mixedargs::Tuple{Ref{Type{T}}, Vararg{Any}}) where {T} =
    capturescalars((args...)->f(T, args...), Base.tail(mixedargs))
@inline capturescalars(f, mixedargs::Tuple{Ref{Type{T}}, Ref{Type{S}}, Vararg{Any}}) where {T, S} =
    # This definition is identical to the one above and necessary only for
    # avoiding method ambiguity.
    capturescalars((args...)->f(T, args...), Base.tail(mixedargs))
@inline capturescalars(f, mixedargs::Tuple{CuSparseVecOrMat, Ref{Type{T}}, Vararg{Any}}) where {T} =
    capturescalars((a1, args...)->f(a1, T, args...), (mixedargs[1], Base.tail(Base.tail(mixedargs))...))
@inline capturescalars(f, mixedargs::Tuple{Union{Ref,AbstractArray{<:Any,0}}, Ref{Type{T}}, Vararg{Any}}) where {T} =
    capturescalars((args...)->f(mixedargs[1], T, args...), Base.tail(Base.tail(mixedargs)))

scalararg(::Number) = true
scalararg(::Any) = false
scalarwrappedarg(::Union{AbstractArray{<:Any,0},Ref}) = true
scalarwrappedarg(::Any) = false

@inline function _capturescalars()
    return (), () -> ()
end
@inline function _capturescalars(arg, mixedargs...)
    let (rest, f) = _capturescalars(mixedargs...)
        if scalararg(arg)
            return rest, @inline function(tail...)
                (arg, f(tail...)...)
            end # add back scalararg after (in makeargs)
        elseif scalarwrappedarg(arg)
            return rest, @inline function(tail...)
                (arg[], f(tail...)...) # TODO: This can put a Type{T} in a tuple
            end # unwrap and add back scalararg after (in makeargs)
        else
            return (arg, rest...), @inline function(head, tail...)
                (head, f(tail...)...)
            end # pass-through to broadcast
        end
    end
end
@inline function _capturescalars(arg) # this definition is just an optimization (to bottom out the recursion slightly sooner)
    if scalararg(arg)
        return (), () -> (arg,) # add scalararg
    elseif scalarwrappedarg(arg)
        return (), () -> (arg[],) # unwrap
    else
        return (arg,), (head,) -> (head,) # pass-through
    end
end

@inline _iszero(x) = x == 0
@inline _iszero(x::Number) = Base.iszero(x)
@inline _iszero(x::AbstractArray) = Base.iszero(x)
@inline _zeros_eltypes(A) = (zero(eltype(A)),)
@inline _zeros_eltypes(A, Bs...) = (zero(eltype(A)), _zeros_eltypes(Bs...)...)

## iteration helpers

"""
    CSRIterator{Ti}(row, args...)

A GPU-compatible iterator for accessing the elements of a single row `row` of several CSR
matrices `args` in one go. The row should be in-bounds for every sparse argument. Each
iteration returns a 2-element tuple: The current column, and each arguments' pointer index
(or 0 if that input didn't have an element at that column). The pointers can then be used to
access the elements themselves.

For convenience, this iterator can be passed non-sparse arguments as well, which will be
ignored (with the returned `col`/`ptr` values set to 0).
"""
struct CSRIterator{Ti,N,ATs}
    row::Ti
    col_ends::NTuple{N, Ti}
    args::ATs
end

function CSRIterator{Ti}(row, args::Vararg{<:Any, N}) where {Ti,N}
    # check that `row` is valid for all arguments
    @boundscheck begin
        ntuple(Val(N)) do i
            arg = @inbounds args[i]
            arg isa CuSparseDeviceMatrixCSR && checkbounds(arg, row, 1)
        end
    end

    col_ends = ntuple(Val(N)) do i
        arg = @inbounds args[i]
        if arg isa CuSparseDeviceMatrixCSR
            @inbounds(arg.rowPtr[row+1i32])
        else
            zero(Ti)
        end
    end

    CSRIterator{Ti, N, typeof(args)}(row, col_ends, args)
end

@inline function Base.iterate(iter::CSRIterator{Ti,N}, state=nothing) where {Ti,N}
    # helper function to get the column of a sparse array at a specific pointer
    @inline function get_col(i, ptr)
        arg = @inbounds iter.args[i]
        if arg isa CuSparseDeviceMatrixCSR
            col_end = @inbounds iter.col_ends[i]
            if ptr < col_end
                return @inbounds arg.colVal[ptr] % Ti
            end
        end
        typemax(Ti)
    end

    # initialize the state
    # - ptr: the current index into the colVal/nzVal arrays
    # - col: the current column index (cached so that we don't have to re-read each time)
    state = something(state,
        ntuple(Val(N)) do i
            arg = @inbounds iter.args[i]
            if arg isa CuSparseDeviceMatrixCSR
                ptr = @inbounds iter.args[i].rowPtr[iter.row] % Ti
                col = @inbounds get_col(i, ptr)
            else
                ptr = typemax(Ti)
                col = typemax(Ti)
            end
            (; ptr, col)
        end
    )

    # determine the column we're currently processing
    cols = ntuple(i -> @inbounds(state[i].col), Val(N))
    cur_col = min(cols...)
    cur_col == typemax(Ti) && return

    # fetch the pointers (we don't look up the values, as the caller might want to index
    # the sparse array directly, e.g., to mutate it). we don't return `ptrs` from the state
    # directly, but first convert the `typemax(Ti)` to a more convenient zero value.
    # NOTE: these values may end up unused by the caller (e.g. in the count_nnzs kernels),
    #       but LLVM appears smart enough to filter them away.
    ptrs = ntuple(Val(N)) do i
        ptr, col = @inbounds state[i]
        col == cur_col ? ptr : zero(Ti)
    end

    # advance the state
    new_state = ntuple(Val(N)) do i
        ptr, col = @inbounds state[i]
        if col == cur_col
            ptr += one(Ti)
            col = get_col(i, ptr)
        end
        (; ptr, col)
    end

    return (cur_col, ptrs), new_state
end

struct CSCIterator{Ti,N,ATs}
    col::Ti
    row_ends::NTuple{N, Ti}
    args::ATs
end

function CSCIterator{Ti}(col, args::Vararg{<:Any, N}) where {Ti,N}
    # check that `col` is valid for all arguments
    @boundscheck begin
        ntuple(Val(N)) do i
            arg = @inbounds args[i]
            arg isa CuSparseDeviceMatrixCSR && checkbounds(arg, 1, col)
        end
    end

    row_ends = ntuple(Val(N)) do i
        arg = @inbounds args[i]
        x = if arg isa CuSparseDeviceMatrixCSC
            @inbounds(arg.colPtr[col+1i32])
        else
            zero(Ti)
        end
        x
    end

    CSCIterator{Ti, N, typeof(args)}(col, row_ends, args)
end

@inline function Base.iterate(iter::CSCIterator{Ti,N}, state=nothing) where {Ti,N}
    # helper function to get the column of a sparse array at a specific pointer
    @inline function get_col(i, ptr)
        arg = @inbounds iter.args[i]
        if arg isa CuSparseDeviceMatrixCSC
            col_end = @inbounds iter.row_ends[i]
            if ptr < col_end
                return @inbounds arg.rowVal[ptr] % Ti
            end
        end
        typemax(Ti)
    end

    # initialize the state
    # - ptr: the current index into the rowVal/nzVal arrays
    # - row: the current row index (cached so that we don't have to re-read each time)
    state = something(state,
        ntuple(Val(N)) do i
            arg = @inbounds iter.args[i]
            if arg isa CuSparseDeviceMatrixCSC
                ptr = @inbounds iter.args[i].colPtr[iter.col] % Ti
                row = @inbounds get_col(i, ptr)
            else
                ptr = typemax(Ti)
                row = typemax(Ti)
            end
            (; ptr, row)
        end
    )

    # determine the row we're currently processing
    rows = ntuple(i -> @inbounds(state[i].row), Val(N))
    cur_row = min(rows...)
    cur_row == typemax(Ti) && return

    # fetch the pointers (we don't look up the values, as the caller might want to index
    # the sparse array directly, e.g., to mutate it). we don't return `ptrs` from the state
    # directly, but first convert the `typemax(Ti)` to a more convenient zero value.
    # NOTE: these values may end up unused by the caller (e.g. in the count_nnzs kernels),
    #       but LLVM appears smart enough to filter them away.
    ptrs = ntuple(Val(N)) do i
        ptr, row = @inbounds state[i]
        row == cur_row ? ptr : zero(Ti)
    end

    # advance the state
    new_state = ntuple(Val(N)) do i
        ptr, row = @inbounds state[i]
        if row == cur_row
            ptr += one(Ti)
            row = get_col(i, ptr)
        end
        (; ptr, row)
    end

    return (cur_row, ptrs), new_state
end

# helpers to index a sparse or dense array
function _getindex(arg::Union{CuSparseDeviceMatrixCSR,CuSparseDeviceMatrixCSC}, I, ptr)
    if ptr == 0
        zero(eltype(arg))
    else
        @inbounds arg.nzVal[ptr]
    end
end
_getindex(arg, I, ptr) = Broadcast._broadcast_getindex(arg, I)

## sparse broadcast implementation

# TODO: unify CSC/CSR kernels

# kernel to count the number of non-zeros in a row, to determine the row offsets
function compute_offsets_kernel(::Type{<:CuSparseMatrixCSR}, offsets::AbstractVector{Ti},
                                args...) where Ti
    # every thread processes an entire row
    row = threadIdx().x + (blockIdx().x - 1i32) * blockDim().x
    row > length(offsets)-1 && return
    iter = @inbounds CSRIterator{Ti}(row, args...)

    # count the nonzero columns of all inputs
    accum = zero(Ti)
    for (col, vals) in iter
        accum += one(Ti)
    end

    # the way we write the nnz counts is a bit strange, but done so that the result
    # after accumulation can be directly used as the rowPtr array of a CSR matrix.
    @inbounds begin
        if row == 1
            offsets[1] = 1
        end
        offsets[row+1] = accum
    end

    return
end
function compute_offsets_kernel(::Type{<:CuSparseMatrixCSC}, offsets::AbstractVector{Ti},
                                args...) where Ti
    # every thread processes an entire columm
    col = threadIdx().x + (blockIdx().x - 1i32) * blockDim().x
    col > length(offsets)-1 && return
    iter = @inbounds CSCIterator{Ti}(col, args...)

    # count the nonzero columns of all inputs
    accum = zero(Ti)
    for (col, vals) in iter
        accum += one(Ti)
    end

    # the way we write the nnz counts is a bit strange, but done so that the result
    # after accumulation can be directly used as the colPtr array of a CSC matrix.
    @inbounds begin
        if col == 1
            offsets[1] = 1
        end
        offsets[col+1] = accum
    end

    return
end

# broadcast kernels that iterate the elements of sparse arrays
function sparse_to_sparse_broadcast_kernel(f, output::CuSparseDeviceMatrixCSR{<:Any,Ti},
                                           offsets::Union{AbstractVector,Nothing},
                                           args...) where {Ti}
    # every thread processes an entire row
    row = threadIdx().x + (blockIdx().x - 1i32) * blockDim().x
    row > size(output, 1) && return
    iter = @inbounds CSRIterator{Ti}(row, args...)

    # fetch the row offset, and write it to the output
    @inbounds begin
        output_ptr = output.rowPtr[row] = offsets[row]
        if row == size(output, 1)
            output.rowPtr[row+1i32] = offsets[row+1i32]
        end
    end

    # set the values for this row
    for (col, ptrs) in iter
        I = CartesianIndex(row, col)
        vals = ntuple(Val(length(args))) do i
            arg = @inbounds args[i]
            ptr = @inbounds ptrs[i]
            _getindex(arg, I, ptr)
        end

        @inbounds output.colVal[output_ptr] = col
        @inbounds output.nzVal[output_ptr] = f(vals...)
        output_ptr += one(Ti)
    end

    return
end
function sparse_to_sparse_broadcast_kernel(f, output::CuSparseDeviceMatrixCSC{<:Any,Ti},
                                           offsets::Union{AbstractVector,Nothing},
                                           args...) where {Ti}
    # every thread processes an entire column
    col = threadIdx().x + (blockIdx().x - 1i32) * blockDim().x
    col > size(output, 2) && return
    iter = @inbounds CSCIterator{Ti}(col, args...)

    # fetch the column offset, and write it to the output
    @inbounds begin
        output_ptr = output.colPtr[col] = offsets[col]
        if col == size(output, 2)
            output.colPtr[col+1i32] = offsets[col+1i32]
        end
    end

    # set the values for this col
    for (row, ptrs) in iter
        I = CartesianIndex(col, row)
        vals = ntuple(Val(length(args))) do i
            arg = @inbounds args[i]
            ptr = @inbounds ptrs[i]
            _getindex(arg, I, ptr)
        end

        @inbounds output.rowVal[output_ptr] = row
        @inbounds output.nzVal[output_ptr] = f(vals...)
        output_ptr += one(Ti)
    end

    return
end
function sparse_to_dense_broadcast_kernel(::Type{<:CuSparseMatrixCSR}, f,
                                          output::CuDeviceArray, args...)
    # every thread processes an entire row
    row = threadIdx().x + (blockIdx().x - 1i32) * blockDim().x
    row > size(output, 1) && return
    iter = @inbounds CSRIterator{Int}(row, args...)

    # set the values for this row
    for (col, ptrs) in iter
        I = CartesianIndex(row, col)
        vals = ntuple(Val(length(args))) do i
            arg = @inbounds args[i]
            ptr = @inbounds ptrs[i]
            _getindex(arg, I, ptr)
        end

        @inbounds output[I] = f(vals...)
    end

    return
end
function sparse_to_dense_broadcast_kernel(::Type{<:CuSparseMatrixCSC}, f,
                                          output::CuDeviceArray, args...)
    # every thread processes an entire column
    col = threadIdx().x + (blockIdx().x - 1i32) * blockDim().x
    col > size(output, 2) && return
    iter = @inbounds CSCIterator{Int}(col, args...)

    # set the values for this col
    for (row, ptrs) in iter
        I = CartesianIndex(row, col)
        vals = ntuple(Val(length(args))) do i
            arg = @inbounds args[i]
            ptr = @inbounds ptrs[i]
            _getindex(arg, I, ptr)
        end

        @inbounds output[I] = f(vals...)
    end

    return
end

function Broadcast.copy(bc::Broadcasted{<:Union{CuSparseVecStyle,CuSparseMatStyle}})
    # on 1.6, ntuple closures often fail to infer
    VERSION < v"1.7" && @warn "Sparse broadcast is only supported on Julia 1.7 or higher" maxlog=1

    # find the sparse inputs
    bc = Broadcast.flatten(bc)
    sparse_args = findall(bc.args) do arg
        arg isa AbstractCuSparseArray
    end
    sparse_types = unique(map(i->nameof(typeof(bc.args[i])), sparse_args))
    if length(sparse_types) > 1
        error("broadcast with multiple types of sparse arrays ($(join(sparse_types, ", "))) is not supported")
    end
    sparse_typ = typeof(bc.args[first(sparse_args)])
    sparse_typ <: Union{CuSparseMatrixCSR,CuSparseMatrixCSC} ||
        error("broadcast with sparse arrays is currently only implemented for CSR and CSC matrices")
    Ti = if sparse_typ <: CuSparseMatrixCSR
        reduce(promote_type, map(i->eltype(bc.args[i].rowPtr), sparse_args))
    elseif sparse_typ <: CuSparseMatrixCSC
        reduce(promote_type, map(i->eltype(bc.args[i].colPtr), sparse_args))
    end

    # determine the output type
    Tv = Broadcast.combine_eltypes(bc.f, eltype.(bc.args))
    if !Base.isconcretetype(Tv)
        error("""GPU sparse broadcast resulted in non-concrete element type $Tv.
                 This probably means that the function you are broadcasting contains an error or type instability.""")
    end

    # partially-evaluate the function, removing scalars.
    parevalf, passedsrcargstup = capturescalars(bc.f, bc.args)
    # check if the partially-evaluated function preserves zeros. if so, we'll only need to
    # apply it to the sparse input arguments, preserving the sparse structure.
    if all(arg->isa(arg, AbstractSparseArray), passedsrcargstup)
        fofzeros = parevalf(_zeros_eltypes(passedsrcargstup...)...)
        fpreszeros = _iszero(fofzeros)
    else
        fpreszeros = false
    end

    # the kernels below parallelize across rows or cols, not elements, so it's unlikely
    # we'll launch many threads. to maximize utilization, parallelize across blocks first.
    rows, cols = size(bc)
    function compute_launch_config(kernel)
        config = launch_configuration(kernel.fun)
        if sparse_typ <: CuSparseMatrixCSR
            threads = min(rows, config.threads)
            blocks = max(cld(rows, threads), config.blocks)
            threads = cld(rows, blocks)
        elseif sparse_typ <: CuSparseMatrixCSC
            threads = min(cols, config.threads)
            blocks = max(cld(cols, threads), config.blocks)
            threads = cld(cols, blocks)
        end
        (; threads, blocks)
    end

    # allocate the output container
    if !fpreszeros
        # either we have dense inputs, or the function isn't preserving zeros,
        # so use a dense output to broadcast into.
        output = CuArray{Tv}(undef, size(bc))

        # since we'll be iterating the sparse inputs, we need to pre-fill the dense output
        # with appropriate values (while setting the sparse inputs to zero). we do this by
        # re-using the dense broadcast implementation.
        nonsparse_args = map(bc.args) do arg
            # NOTE: this assumes the broadcst is flattened, but not yet preprocessed
            if arg isa AbstractCuSparseArray
                zero(eltype(arg))
            else
                arg
            end
        end
        broadcast!(bc.f, output, nonsparse_args...)
    elseif length(sparse_args) == 1
        # we only have a single sparse input, so we can reuse its structure for the output.
        # this avoids a kernel launch and costly synchronization.
        sparse_arg = bc.args[first(sparse_args)]
        if sparse_typ <: CuSparseMatrixCSR
            offsets = rowPtr = sparse_arg.rowPtr
            colVal = similar(sparse_arg.colVal)
            nzVal = similar(sparse_arg.nzVal, Tv)
            output = CuSparseMatrixCSR(rowPtr, colVal, nzVal, size(bc))
        elseif sparse_typ <: CuSparseMatrixCSC
            offsets = colPtr = sparse_arg.colPtr
            rowVal = similar(sparse_arg.rowVal)
            nzVal = similar(sparse_arg.nzVal, Tv)
            output = CuSparseMatrixCSC(colPtr, rowVal, nzVal, size(bc))
        end
        # NOTE: we don't use CUSPARSE's similar, because that copies the structure arrays,
        #       while we do that in our kernel (for consistency with other code paths)
    else
        # determine the number of non-zero elements per row so that we can create an
        # appropriately-structured output container
        offsets = if sparse_typ <: CuSparseMatrixCSR
            CuArray{Ti}(undef, rows+1)
        elseif sparse_typ <: CuSparseMatrixCSC
            CuArray{Ti}(undef, cols+1)
        end
        let
            args = (sparse_typ, offsets, bc.args...)
            kernel = @cuda launch=false compute_offsets_kernel(args...)
            threads, blocks = compute_launch_config(kernel)
            kernel(args...; threads, blocks)
        end

        # accumulate these values so that we can use them directly as row pointer offsets,
        # as well as to get the total nnz count to allocate the sparse output array.
        # cusparseXcsrgeam2Nnz computes this in one go, but it doesn't seem worth the effort
        accumulate!(Base.add_sum, offsets, offsets)
        total_nnz = @allowscalar last(offsets[end]) - 1

        output = if sparse_typ <: CuSparseMatrixCSR
            colVal = CuArray{Ti}(undef, total_nnz)
            nzVal = CuArray{Tv}(undef, total_nnz)
            CuSparseMatrixCSR(offsets, colVal, nzVal, size(bc))
        elseif sparse_typ <: CuSparseMatrixCSC
            rowVal = CuArray{Ti}(undef, total_nnz)
            nzVal = CuArray{Tv}(undef, total_nnz)
            CuSparseMatrixCSC(offsets, rowVal, nzVal, size(bc))
        end
    end

    # perform the actual broadcast
    if output isa AbstractCuSparseArray
        args = (bc.f, output, offsets, bc.args...)
        kernel = @cuda launch=false sparse_to_sparse_broadcast_kernel(args...)
        threads, blocks = compute_launch_config(kernel)
        kernel(args...; threads, blocks)
    else
        args = (sparse_typ, bc.f, output, bc.args...)
        kernel = @cuda launch=false sparse_to_dense_broadcast_kernel(args...)
        threads, blocks = compute_launch_config(kernel)
        kernel(args...; threads, blocks)
    end

    return output
end


text:
export NLfunForward

"""
"""
function NLfunForward(A, x, g, params)

    e = params["numr"]*params["nr"]
    L = x[1:e]
    R = x[(e+1):end]
    L = reshape(L, params["numr"], params["nr"])
    R = reshape(R, params["numc"], params["nr"])

    if isempty(g)
        f1 = params["afun"](L*R', params)
        f2 = 0.
    else
        fp = params["afunT"](g)
        f1 = [vec(fp*R); vec(fp'*L)]
        f2 = vec(fp)
    end

    return f1,f2
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module Jentropy

using DataFrames
using StatsBase

export atleast2d,
  bayessupport,
  decimalise,
  calcentropy,
  checkinputs,
  condentropy,
  condprob,
  prob,
  requireents,
  requireprobs

include("atleast2d.jl")
include("bayessupport.jl")
include("decimalise.jl")
include("calcentropy.jl")
include("checkinputs.jl")
include("condentropy.jl")
include("condprob.jl")
include("prob.jl")
include("requireents.jl")
include("requireprobs.jl")
end


text:
# Buffer
# ======

const INIT_BUFFER_SIZE = Ref(256)

mutable struct Buffer
    # UTF-8 encoded sequence + 0-3 incomplete bytes
    data::Vector{UInt8}
    # the current reading position (the first byte of the UTF-8 sequence)
    p::Int
    # the end position of the UTF-8 sequence
    p_end::Int
    # the end position of the filled bytes
    p_fill::Int
    # the end position of the stream (-1 if yet found)
    p_eof::Int

    function Buffer()
        data = Vector{UInt8}(undef, INIT_BUFFER_SIZE[])
        return new(data, 1, 0, 0, -1)
    end
end

function fillbuffer!(input::IO, buffer::Buffer)
    if buffer.p_eof ≥ 0
        return 0
    end
    if buffer.p < buffer.p_fill
        # move data
        shift = buffer.p - 1
        copyto!(buffer.data, 1, buffer.data, buffer.p, buffer.p_fill - buffer.p + 1)
        buffer.p = 1
        buffer.p_end -= shift
        buffer.p_fill -= shift
    end
    if eof(input)
        # found EOF
        buffer.p_end = buffer.p_eof = buffer.p_fill
        return 0
    end
    n_avail = bytesavailable(input)
    n_free = length(buffer.data) - buffer.p_fill
    if n_free == 0  # no space to fill data
        resize!(buffer.data, length(buffer.data) * 2)
        n_free = length(buffer.data) - buffer.p_fill
    end
    # fill data into the buffer
    n = min(n_avail, n_free)
    @assert n > 0
    unsafe_read(input, pointer(buffer.data, buffer.p_fill + 1), n)
    buffer.p_fill += n
    # align p_end to a UTF-8 boundary
    p_new = buffer.p_end + 1
    p_end = buffer.p_fill
    if buffer.data[p_end] ≤ 0b01111111
        # ascii
        buffer.p_end = p_end
    else
        for i in 1:3
            if buffer.data[p_end] >> 6 == 0b10
                p_end -= 1
            else
                break
            end
            if p_end < buffer.p
                @goto utf8error
            end
        end
        if buffer.p_fill == p_end + leading_ones(buffer.data[p_end])
            buffer.p_end = buffer.p_fill
        else
            buffer.p_end = p_end - 1
        end
    end
    # validate UTF-8 encoding
    if buffer.p_end > p_new && !is_valid_utf8(buffer.data, p_new, buffer.p_end)
        @label utf8error
        throw(ErrorException("invalid UTF-8 sequence"))
    end
    return n
end

function is_valid_utf8(data::Vector{UInt8}, from::Int, to::Int)
    return ccall(:u8_isvalid, Cint, (Ptr{UInt8}, Csize_t), pointer(data, from), to - from + 1) > 0
end

function ensurebytes!(input::IO, buffer::Buffer, n::Int)
    if buffer.p_end - buffer.p + 1 ≥ n
        # ok
        return true
    end
    fillbuffer!(input, buffer)
    return buffer.p_end - buffer.p + 1 ≥ n
end

@inline function peekchar(input::IO, buffer::Buffer; offset::Int=0)
    ensurebytes!(input, buffer, offset+4)
    @inbounds if buffer.p + offset ≤ buffer.p_end
        # NOTE: UTF-8 encoding is validated in fillbuffer!
        b = buffer.data[buffer.p+offset]
        u = UInt32(b) << 24
        if b < 0b10000000
            return reinterpret(Char, u), 1
        elseif b ≤ 0b11011111
            u |= UInt32(buffer.data[buffer.p+offset+1]) << 16
            return reinterpret(Char, u), 2
        elseif b ≤ 0b11101111
            u |= UInt32(buffer.data[buffer.p+offset+1]) << 16
            u |= UInt32(buffer.data[buffer.p+offset+2]) <<  8
            return reinterpret(Char, u), 3
        else @assert b ≤ 0b11110111
            u |= UInt32(buffer.data[buffer.p+offset+1]) << 16
            u |= UInt32(buffer.data[buffer.p+offset+2]) <<  8
            u |= UInt32(buffer.data[buffer.p+offset+3])
            return reinterpret(Char, u), 4
        end
    else
        # no more character
        return '\0', 0
    end
end

function taketext!(buffer::Buffer, size::Int)
    @assert buffer.p + size - 1 ≤ buffer.p_end
    # This allocates less than `String(buffer.data[buffer.p:buffer.p+size-1])`.
    text = unsafe_string(pointer(buffer.data, buffer.p), size)
    buffer.p += size
    return text
end

function takeback!(buffer::Buffer, text::String)
    size = sizeof(text)
    @assert size ≤ buffer.p
    copyto!(buffer.data, buffer.p - size, text, 1, size)
    buffer.p -= size
    return buffer
end

function consume!(buffer::Buffer, size::Int)
    buffer.p += size
    return
end

function scanwhile(f::Function, input::IO, buffer::Buffer; offset::Int=0)
    o = offset
    while (char_n = peekchar(input, buffer, offset=o))[2] > 0
        char, n = char_n
        if !f(char)
            break
        end
        o += n
    end
    return o - offset
end

function scanbytes(f::Function, input::IO, buffer::Buffer)
    n = 0
    @label scan
    while buffer.p + n ≤ buffer.p_end
        @inbounds if !f(buffer.data[buffer.p+n])
            break
        end
        n += 1
    end
    if buffer.p + n > buffer.p_end && buffer.p_eof < 0 && fillbuffer!(input, buffer) > 0
        @goto scan
    end
    return n
end

scanwhitespace(input::IO, buffer::Buffer) =
    scanbytes(
        # space or tab
        b -> (b == 0x20) | (b == 0x09),
        input, buffer)
scanbarekey(input::IO, buffer::Buffer) =
    scanbytes(
        # - 0-9 A-Z _ a-z
        b -> (b == 0x2d) | (0x30 ≤ b ≤ 0x39) | (0x41 ≤ b ≤ 0x5a) | (b == 0x5f) | (0x61 ≤ b ≤ 0x7a),
        input, buffer)


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# ╔═╡ ea52d932-9304-11eb-0c41-4d06311308a6
begin
	using Pkg
	Pkg.activate(temp = true)

	Pkg.add("CairoMakie")

	using CairoMakie
	CairoMakie.activate!(type = "png")
end

# ╔═╡ 70c50868-5465-4b03-918b-a83c83d224c5
using PlutoUI; TableOfContents()

# ╔═╡ 528b77d3-2672-418a-b059-bdbbea41b54f
md"""
## Default attributes
"""

# ╔═╡ 38df7208-a6ab-4824-acc7-42c5f22f970a
default_orientation(legend_position) = legend_position in [:top, :bottom] ? :horizontal : :vertical

# ╔═╡ 0e38adf6-8853-43a7-b44e-1df517ed8bd2
default_nbanks(orientation, has_colorbar) = has_colorbar && (orientation == :horizontal) ? 2 : 1

# ╔═╡ 172b31a8-d853-4f64-a89c-9e61129b3455
default_titleposition(orientation) = orientation == :horizontal ? :left : :top

# ╔═╡ 59723752-2d54-4b8f-b55a-02f165fc2b6c
md"""
## Compose
"""

# ╔═╡ 0986d286-af46-4373-937f-3f73b28430f3
md"""
## Test plots
"""

# ╔═╡ 25912b01-7a74-44cc-9336-e50ccc84edb8
function positions(fig, has_legend, has_colorbar, vertical, legend_position)
	if has_legend || has_colorbar
		if legend_position == :bottom
			ax_pos 	 = fig[1,1]
			legs_pos = fig[2,1]
		elseif legend_position == :top
			ax_pos 	 = fig[2,1]
			legs_pos = fig[1,1]
		elseif legend_position == :right
			ax_pos   = fig[1,1]
			legs_pos = fig[1,2]
		elseif legend_position == :left
			ax_pos   = fig[1,2]
			legs_pos = fig[1,1]
		end
		i = 1
	else
		ax_pos = fig[1,1]
	end

	if has_legend
		leg_pos = vertical ? legs_pos[i,1] : legs_pos[1,i]
		i += 1
	else
		leg_pos = nothing
	end
	if has_colorbar
		cb_pos = vertical ? legs_pos[i,1] : legs_pos[1,i]
	else
		cb_pos = nothing
	end

	(; ax_pos, leg_pos, cb_pos)
end

# ╔═╡ 276a78cc-7a76-4901-a3ea-edcf60bb7d3b
function make_legend(fig, legend_title, n = 2)
	markersizes = [5, 10, 15, 20]

	group_size = [MarkerElement(marker = :circle, color = :black, strokecolor = :transparent,
    	markersize = ms) for ms in markersizes]

	legend = Legend(fig,
    	fill(group_size, n),
    	fill(string.(markersizes), n),
    	fill(legend_title, n)
	)
end

# ╔═╡ 8d63a603-6ed6-4d9c-8c50-19c7535324b7
function adjust_attributes!(legend; attr...)
	for (a,b) in pairs(attr)
		setproperty!(legend, a, b)
	end
end

# ╔═╡ 4dacaaee-05eb-488e-857b-9c7b95b6c181
md"""
## Legend
"""

# ╔═╡ 46084f3f-844b-4b62-9fe5-812ca917236e
function legend_attributes(orientation, has_colorbar, titleposition, nbanks)
	vertical = orientation == :vertical
	horizontal = !vertical

	stretch_height = has_colorbar && vertical
	stretch_width  = has_colorbar && horizontal

	(; nbanks, orientation, titleposition,
	   tellwidth = vertical    || stretch_width,
	   tellheight = horizontal || stretch_height,
	   framevisible = false)
end

# ╔═╡ b5318ac4-91f2-44ab-b6ea-f12a63fe68de
md"""
## Colorbar
"""

# ╔═╡ 72c5cdcf-d6ee-435b-b018-dbb2755be00c
function colorbar_titleposition(cbpos, titlepos, has_legend)
	if titlepos == :top
		cbtitlepos = cbpos[0,1]
	elseif titlepos == :left
		cbtitlepos = cbpos[1,0]
	end
end

# ╔═╡ 711d8f7a-16fd-42ab-9e7a-c3c91f76fe1e
function add_colorbar(colorbar, cbpos, titlepos, has_legend, cb_attributes)
	cb = Colorbar(cbpos[1,1]; colorbar.kwargs..., cb_attributes...)

	squeeze_label_height = titlepos == :top
	squeeze_label_width  = titlepos == :left

	if !isnothing(colorbar.title)
		cbtitlepos = colorbar_titleposition(cbpos, titlepos, has_legend)

		Label(cbtitlepos, colorbar.title, tellheight = squeeze_label_height, tellwidth = squeeze_label_width)
	end
end

# ╔═╡ f0fbd3f9-7e44-4e76-9c5e-f5347863cf8c
function colorbar_attributes(orientation)
	vertical = orientation == :vertical
	horizontal = !vertical

	if horizontal
		cb_attributes = (height = 18, width = Relative(1.0),
			tellwidth = true,
			vertical = vertical,
			valign = :center,
			ticksize = 5,
  			ticklabelpad = 1.5)
	else
		cb_attributes = (width = 18, height = Relative(1.0),
			vertical = vertical,
			halign = :left,
	)

	end
end

# ╔═╡ d9cc8536-9304-11eb-1a23-dd648f1d3ef7
function add_legend_and_colorbar(fig, legend_title, colorbar;
		legend_position = :top,
		orientation = default_orientation(legend_position),
		titleposition =  default_titleposition(orientation),
		nbanks = default_nbanks(orientation, !isnothing(colorbar))
		)

	has_colorbar = !isnothing(colorbar)
	has_legend = true

	vertical = orientation == :vertical

	ax_pos, legpos, cbpos = positions(fig, has_legend, has_colorbar, vertical, legend_position)

	ax = Axis(ax_pos)

	if has_legend
		legend = make_legend(legpos, legend_title)
		leg_attributes = legend_attributes(orientation, has_colorbar, titleposition, nbanks)
		adjust_attributes!(legend; leg_attributes...)
	end

	if has_colorbar
		cb_attributes = colorbar_attributes(orientation)
		add_colorbar(colorbar, cbpos, titleposition, has_legend, cb_attributes)
	end

	fig
end

# ╔═╡ edd7503c-ef60-4bd6-8b2a-1f0c517216a7
function test(legend_title, colorbar_title; kwargs...
		)

	my_colorbar = (args = (),
			  kwargs = (limits = (-1, 1), ),
			  title = "nice title")

	fig = Figure()

	add_legend_and_colorbar(fig, legend_title, my_colorbar; kwargs...)

end

# ╔═╡ 0303e9a8-9305-11eb-0737-c187ab85885b
test("bla bla", "bla", legend_position = :bottom)

# ╔═╡ 0b209f44-d7c1-424d-8fa1-0d5a46657058
test("bla bla", "bla", legend_position = :top, titleposition = :top)

# ╔═╡ 390a3394-63e9-479d-bd54-fa5088d16abb
test("bla bla", "bla", legend_position = :right, nbanks = 2)

# ╔═╡ d0ca4cf2-96f5-4d48-90d6-62909628a6ee
test("bla bla", "bla",  legend_position = :left, titleposition = :left)

# ╔═╡ Cell order:
# ╠═ea52d932-9304-11eb-0c41-4d06311308a6
# ╠═edd7503c-ef60-4bd6-8b2a-1f0c517216a7
# ╟─528b77d3-2672-418a-b059-bdbbea41b54f
# ╠═38df7208-a6ab-4824-acc7-42c5f22f970a
# ╠═0e38adf6-8853-43a7-b44e-1df517ed8bd2
# ╠═172b31a8-d853-4f64-a89c-9e61129b3455
# ╟─59723752-2d54-4b8f-b55a-02f165fc2b6c
# ╠═d9cc8536-9304-11eb-1a23-dd648f1d3ef7
# ╟─0986d286-af46-4373-937f-3f73b28430f3
# ╠═0303e9a8-9305-11eb-0737-c187ab85885b
# ╠═0b209f44-d7c1-424d-8fa1-0d5a46657058
# ╠═390a3394-63e9-479d-bd54-fa5088d16abb
# ╠═d0ca4cf2-96f5-4d48-90d6-62909628a6ee
# ╠═25912b01-7a74-44cc-9336-e50ccc84edb8
# ╠═276a78cc-7a76-4901-a3ea-edcf60bb7d3b
# ╠═8d63a603-6ed6-4d9c-8c50-19c7535324b7
# ╟─4dacaaee-05eb-488e-857b-9c7b95b6c181
# ╠═46084f3f-844b-4b62-9fe5-812ca917236e
# ╟─b5318ac4-91f2-44ab-b6ea-f12a63fe68de
# ╠═72c5cdcf-d6ee-435b-b018-dbb2755be00c
# ╠═711d8f7a-16fd-42ab-9e7a-c3c91f76fe1e
# ╠═f0fbd3f9-7e44-4e76-9c5e-f5347863cf8c
# ╠═70c50868-5465-4b03-918b-a83c83d224c5

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
using SHA

name = "libpg_query"
version = v"0.0.6"

# Collection of sources required to build libpg_query
sources = [
    GitSource(
        "https://github.com/relationalai-oss/libpg_query.git",
        "c18afd7395a19a35043e1af94ac01360b8730106",
    ),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libpg_query/
make build
mkdir -p $prefix/lib
cp libpg_query.* $prefix/lib/

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    MacOS(:x86_64)
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libpg_query", :libpg_query)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

version_string = string(version)

function compute_file_sha(file_name)
    return open(joinpath(".", "products", file_name)) do file
        bytes2hex(sha2_256(file))
    end
end

build_jl_file = """
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libpg_query"], :libpg_query),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/RelationalAI-oss/libpg_query_binary_builder/releases/download/$(version_string)"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:i686, libc=:glibc) => ("\$bin_prefix/libpg_query.v$(version_string).i686-linux-gnu.tar.gz", "$(compute_file_sha("libpg_query.v$(version_string).i686-linux-gnu.tar.gz"))"),
    Linux(:i686, libc=:musl) => ("\$bin_prefix/libpg_query.v$(version_string).i686-linux-musl.tar.gz", "$(compute_file_sha("libpg_query.v$(version_string).i686-linux-musl.tar.gz"))"),
    MacOS(:x86_64) => ("\$bin_prefix/libpg_query.v$(version_string).x86_64-apple-darwin14.tar.gz", "$(compute_file_sha("libpg_query.v$(version_string).x86_64-apple-darwin14.tar.gz"))"),
    Linux(:x86_64, libc=:glibc) => ("\$bin_prefix/libpg_query.v$(version_string).x86_64-linux-gnu.tar.gz", "$(compute_file_sha("libpg_query.v$(version_string).x86_64-linux-gnu.tar.gz"))"),
    Linux(:x86_64, libc=:musl) => ("\$bin_prefix/libpg_query.v$(version_string).x86_64-linux-musl.tar.gz", "$(compute_file_sha("libpg_query.v$(version_string).x86_64-linux-musl.tar.gz"))"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\\"\$(Sys.MACHINE)\\", parsed as \\"\$(triplet(platform_key_abi()))\\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)

"""

open(joinpath(".", "products", "build_libpg_query.v$(version_string).jl"), "w") do io
    write(io, build_jl_file)
end;


text:
"""

$(TYPEDEF)

Structure to contain the names of the output files.

$(TYPEDFIELDS)

"""
@with_kw mutable struct OutputFiles

    output::String
    solute_atoms::String
    solvent_atoms::String

end


text:
using ITensors
using Test
import Random: seed!

seed!(12345)

@testset "ITensor combine contract order $(2*N) tensors" for N in 1:5
  d = 1
  i = Index(QN(0) => d, QN(1) => d)
  is = IndexSet(ntuple(n -> settags(i, "i$n"), Val(N)))
  A = randomITensor(is'..., dag(is)...)
  B = randomITensor(is'..., dag(is)...)

  @test !ITensors.using_combine_contract()

  C_contract = A' * B

  ITensors.enable_combine_contract()
  @test ITensors.using_combine_contract()

  C_combine_contract = A' * B

  ITensors.disable_combine_contract()
  @test !ITensors.using_combine_contract()

  @test C_contract ≈ C_combine_contract
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 63952e48-fda6-11ea-29d5-c9e748d213cf
using Pkg, DrWatson

# ╔═╡ 639575b0-fda6-11ea-2629-8fdf4e78c434
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 3effb80c-fda4-11ea-1b76-33e901d4c0ad
md"## Clip-05-35-37s.jl"

# ╔═╡ 639604e4-fda6-11ea-0184-af7e26ad59fa
begin
	df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
	df = filter(row -> !(row[:neocortex_perc] == "NA"), df);
	df[!, :neocortex_perc] = parse.(Float64, df[:, :neocortex_perc])
	df[!, :lmass] = log.(df[:, :mass])
	scale!(df, [:kcal_per_g, :neocortex_perc, :lmass])
end;

# ╔═╡ 63a35040-fda6-11ea-0016-4bf41ecaa531
md"### snippet 5.35"

# ╔═╡ 63a85bda-fda6-11ea-1664-4d81fe6e0bc2
stan5_5 = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] K; // Outcome
 vector[N] NC; // Predictor
}

parameters {
 real a; // Intercept
 real bN; // Slope (regression coefficients)
 real < lower = 0 > sigma;    // Error SD
}

model {
  vector[N] mu;               // mu is a vector
  a ~ normal(0, 0.2);           //Priors
  bN ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bN * NC;
  K ~ normal(mu , sigma);     // Likelihood
}
";

# ╔═╡ 63b25188-fda6-11ea-2e79-47306de1cf27
md"## Define the SampleModel, etc."

# ╔═╡ 63b3736c-fda6-11ea-35fe-5b5f3d67aef6
begin
	m5_5s = SampleModel("m5.5", stan5_5);
	m5_5_data = Dict("N" => size(df, 1), "NC" => df[!, :neocortex_perc_s],
		"K" => df[!, :kcal_per_g_s]);
	rc5_5s = stan_sample(m5_5s, data=m5_5_data);
end;

# ╔═╡ 63bc7bba-fda6-11ea-1e4e-3deb1786b246
if success(rc5_5s)
  post5_5s_df = read_samples(m5_5s, :dataframe)
  title = "Kcal_per_g vs. neocortex_perc" * "\nshowing predicted and hpd range"
  plotbounds(
    df, :neocortex_perc, :kcal_per_g,
    post5_5s_df, [:a, :bN, :sigma];
    title=title
  )
end

# ╔═╡ 63bce55a-fda6-11ea-1485-dd20ccc99814
md"## End of clip-05-35-37s.jl"

# ╔═╡ Cell order:
# ╟─3effb80c-fda4-11ea-1b76-33e901d4c0ad
# ╠═63952e48-fda6-11ea-29d5-c9e748d213cf
# ╠═639575b0-fda6-11ea-2629-8fdf4e78c434
# ╠═639604e4-fda6-11ea-0184-af7e26ad59fa
# ╟─63a35040-fda6-11ea-0016-4bf41ecaa531
# ╠═63a85bda-fda6-11ea-1664-4d81fe6e0bc2
# ╟─63b25188-fda6-11ea-2e79-47306de1cf27
# ╠═63b3736c-fda6-11ea-35fe-5b5f3d67aef6
# ╠═63bc7bba-fda6-11ea-1e4e-3deb1786b246
# ╟─63bce55a-fda6-11ea-1485-dd20ccc99814


text:
using Test, Random, SparseArrays, LinearAlgebra
import LuxurySparse: IMatrix, PermMatrix

Random.seed!(2)

p1 = IMatrix{4}()
sp = sprand(ComplexF64, 4,4, 0.5)
ds = rand(ComplexF64, 4,4)
pm = PermMatrix([2,3,4,1], randn(4))
v = [0.5, 0.3im, 0.2, 1.0]
dv = Diagonal(v)

@testset "kron" begin
    for source in [p1, sp, ds, dv, pm]
        for target in [p1, sp, ds, dv, pm]
            lres = kron(source, target)
            rres = kron(target, source)
            flres = kron(Matrix(source), Matrix(target))
            frres = kron(Matrix(target), Matrix(source))
            @test lres == flres
            @test rres == frres
            @test eltype(lres) == eltype(flres)
            @test eltype(rres) == eltype(frres)
            if !(target === ds && source === ds)
                @test !(typeof(lres) <: StridedMatrix)
                @test !(typeof(rres) <: StridedMatrix)
            end
        end
    end
end


text:
module Testing
using DynamicalBilliards, Test
using DynamicalBilliards: isperiodic

export tag, omnibilliard, finitehorizon, testparticles, isperiodic
export acclevel

using DynamicalBilliards: distance, cellsize
export distance, cellsize

CONCAVE = Union{Semicircle}

tag(p::Particle) = "[STRAIGHT]"
tag(p::MagneticParticle) = "[MAGNETIC]"

function tag(bd::Billiard)
    s = DynamicalBilliards.isperiodic(bd) ? "[PERIODIC]" : ""
    if length(unique(map(x -> typeof(x), bd))) > 7
        s *= "[OMNI]"
    elseif any(o -> typeof(o) <: CONCAVE, bd)
        s *= "[CONCAVE]"
    end
    return s
end

tag(t::Union{<: RaySplitter, <: Tuple}) = "[RAYSPLIT]"

tag(args...) = join([tag(a) for a in args])

acclevel(::Particle{T}) where {T} = eps(T)^(3/4)
acclevel(::MagneticParticle{T}) where {T} = eps(T)^(3/4)
acclevel(bd::Billiard{T}) where {T} = isperiodic(bd) ? eps(T)^(3/5) : eps(T)^(3/4)
acclevel(t::Union{<: RaySplitter, <: Tuple}) = eps()^(3/4)
acclevel(args...) = maximum(acclevel(a) for a in args)

function omnibilliard()
    pref = 3 # how big is the frame with respect to the Disks
    h = pref*1.0; α = pref*0.8; r = 0.3; off = pref*0.25

    cos6 = cos(π/6)
    β = (h - cos6*α)/cos6
    t = α + 2β
    center_of_mass = [0.0, √3*t/6]
    startloc = [-α/2, 0.0]

    # create directions of the hexagonal 6:
    hexvert = [(cos(2π*i/6), sin(2π*i/6)) for i in 1:6]
    dirs = [SVector{2}(hexvert[i] .- hexvert[mod1(i+1, 6)]) for i in 1:6]

    frame = Obstacle{Float64}[]

    sp = startloc
    ep = startloc + α*dirs[1]
    normal = (w = ep .- sp; [-w[2], w[1]])
    push!(frame, Semicircle((sp+ep)/2, α/2, normal))

    for i in 2:6
        s = iseven(i) ? β : α
        x = mod(i, 4)
        T = if x == 1
            RandomWall
        elseif x == 2
            FiniteWall
        elseif x == 3
            SplitterWall
        elseif x == 0
            InfiniteWall
        end
        sp = i != 2 ? frame[i-1].ep : startloc + α*dirs[1]
        ep = sp + s*dirs[i]
        normal = (w = ep .- sp; [-w[2], w[1]])
        push!(frame, T(sp, ep, normal))
    end

    # Radii of circles that compose the Julia logo
    offset = [0.0, off]

    R = [cos(2π/3) -sin(2π/3);
         sin(2π/3)  cos(2π/3)]

    green = Disk(center_of_mass .+ offset, r, "green")
    red = Antidot(center_of_mass .+ R*offset, r, "red")
    purple = RandomDisk(center_of_mass .+ R*R*offset, r, "purple")

    bd = Billiard(green, red, purple, frame...)
end

finitehorizon(r = 0.5) = billiard_hexagonal_sinai(0.5; setting = "periodic")

function testparticles()
    p = Particle(0.31, 0.51, 2π*rand())
    mp = MagneticParticle(0.31, 0.51, 2π*rand(), 0.5)
    return p, mp
end

function all_tests(f, args...)
    omni_tests(f, args...)
    periodic_tests(f, args...)
end

function simple_tests(f, args...)
    ergodic_tests(f, args...)
    periodic_tests(f, args...)
end

function omni_tests(f, args...)
    bd = omnibilliard()
    p, mp = testparticles()
    e = Ellipse([0.68, 1.53], 0.3, 0.2)
    bd2 = Billiard(bd..., e)
    f(p, bd2, args...)
    f(mp, bd, args...)
end

function omni_tests_noellipse(f, args...)
    bd = omnibilliard()
    p, mp = testparticles()
    f(p, bd, args...)
    f(mp, bd, args...)
end
export omni_tests_noellipse

function periodic_tests(f, args...)
    bd = finitehorizon()
    p, mp = testparticles()
    f(p, bd, args...)
    f(mp, bd, args...)
end

function ergodic_tests(f, args...)
    for bd in (billiard_sinai(), billiard_stadium())
        p = Particle(0.1, 0.1, 2π*rand())
        mp = MagneticParticle(p, 1.0)
        f(p, bd, args...)
        f(mp, bd, args...)
    end
end

"""
    billiards_testset(description, f, args...; caller = all_tests)
Wrap a testset around `caller(f, args...)` which times the result
and separates the test from others (`println`).
"""
function billiards_testset(d, f, args...; caller = all_tests)
    println(">> TEST: $(d)")
    t = time()
    @testset "$(d)" begin
        caller(f, args...)
    end
    println("Required time: $(round(time()-t, digits=3)) sec.")
    separator()
end
separator() = println("\n", "- "^30, "\n")

export all_tests, simple_tests, omni_tests, periodic_tests, billiards_testset
export ergodic_tests

function basic_ray(ellipse)
    bd = billiard_rectangle()
    d = Antidot([0.25, 0.5], 0.15)
    if ellipse
        e = Ellipse([0.75, 0.5], 0.15, 0.25)
        bd = Billiard(bd..., d, e)
    else
        bd = Billiard(bd..., d)
    end
    sa = (φ, pflag, ω) -> pflag ? 0.5φ : 2.0φ
    newoantidot = ((x, bool) -> bool ? -0.5x : -2.0x)
    Tp = (p) -> (φ, pflag, ω) -> begin
        if pflag
            p
        else
            abs(φ) < π/4 ? (1-p) : 0.0
        end
    end
    ray = RaySplitter(ellipse ? [5, 6] : [5], Tp(0.5), sa, newoantidot)

    return bd, (ray,)
end

function extreme_ray(ellipse)
    x = 2.0; y = 1.0
    bdr =  billiard_rectangle(2.0, 1.0)
    sw = SplitterWall([x/2, 0.0], [x/2,y], [-1,0], true)
    a1 = Antidot([x/4, y/2], 0.25, "Left Antidot")
    if ellipse
        a2 = Ellipse([3x/4, y/2], 0.25, 0.2, true, "Ellipse")
    else
        a2 = Antidot([3x/4, y/2], 0.25, "Right Antidot")
    end

    sa(φ, pflag, ω) = pflag ? 2.0φ : 0.5φ
    function Tp(φ, pflag, ω)
        if pflag
            abs(φ) < π/4 ? 0.9 : 0.0
        else
            0.9
        end
    end

    newoantidot = ((x, bool) -> bool ? -0.5x : -2.0x)

    raywall = RaySplitter([3], Tp, sa, newoantidot)
    raya = RaySplitter([1, 2], Tp, sa, newoantidot)

    return Billiard(a1, a2, sw, bdr...), (raywall, raya)
end

function ray_tests(f, args...)
    bd, ray = basic_ray(true)
    p = randominside(bd)
    @testset "$(tag(p, bd, ray))" begin; f(p, bd, ray, args...); end
    bd, ray = basic_ray(false)
    p = randominside(bd, 1.0)
    @testset "$(tag(p, bd, ray))" begin; f(p, bd, ray, args...); end
end

export basic_ray, extreme_ray, ray_tests

end