text:
const GERG2008_names=[
        "methane",
        "nitrogen",
        "carbon dioxide",
        "ethane",
        "propane",
        "butane",
        "isobutane",
        "pentane",
        "isopentane",
        "hexane",
        "heptane",
        "octane",
        "nonane",
        "decane",
        "hydrogen",
        "oxygen",
        "carbon monoxide",
        "water",
        "hydrogen sulfide",
        "helium",
        "argon"
        ]

function getparams_gerg2008()

    components =[
            "methane",
            "nitrogen",
            "carbon dioxide",
            "ethane",
            "propane",
            "butane",
            "isobutane",
            "pentane",
            "isopentane",
            "hexane",
            "heptane",
            "octane",
            "nonane",
            "decane",
            "hydrogen",
            "oxygen",
            "carbon monoxide",
            "water",
            "hydrogen sulfide",
            "helium",
            "argon"
            ]

    lb_v = SingleParam("lower bound volume",components,[
        2.952789230664869e-5
        2.595917396536852e-5
        3.0350520997504374e-5
        4.181265338895899e-5
        5.7773265228561896e-5
        7.778071972263543e-5
        7.529511394064539e-5
        9.03473213517436e-5
        9.138272716465976e-5
        0.00010746846161549492
        0.00012307769607878738
        0.00014149360905504386
        0.00013613745434844812
        0.00015142994339404644
        1.6005013101722185e-5
        2.3229544237311946e-5
        2.9381748144372254e-5
        1.732307538792257e-5
        3.154736660208043e-5
        1.2820394828598075e-5
        2.511877404958246e-5
    ])
    acentricfactor = SingleParam("acentric factor",components,[
            0.0115478,
            0.0377215,
            0.223621,
            0.099493,
            0.152291,
            0.183521,
            0.200164,
            0.227875,
            0.251506,
            0.301261,
            0.349469,
            0.399552,
            0.44346,
            0.492328,
            -0.215993,
            0.0221798,
            0.0481621,
            0.344861,
            0.0941677,
            -0.390032,
            0.0,
        ])
        Mw = SingleParam("Molecular weight",components,[
            16.04246,
            28.0134,
            44.0095,
            30.06904,
            44.09562,
            58.1222,
            58.1222,
            72.14878,
            72.14878,
            86.17536,
            100.20194,
            114.22852,
            128.2551,
            142.28168,
            2.01588,
            31.9988,
            28.0101,
            18.01528,
            34.08088,
            4.002602,
            39.948,
        ])
        Tc = SingleParam("Critical Temperature (K)",components,[
            190.564,
            126.192,
            304.1282,
            305.322,
            369.825,
            425.125,
            407.817,
            469.7,
            460.35,
            507.82,
            540.13,
            569.32,
            594.55,
            617.7,
            33.19,
            154.595,
            132.86,
            647.096,
            373.1,
            5.1953,
            150.687,
        ])
        rhoc = SingleParam("Critical density (mol/dm3)",components,[
            10.139342719,
            11.1839,
            10.624978698,
            6.87085454,
            5.000043088,
            3.920016792,
            3.86014294,
            3.215577588,
            3.271,
            2.705877875,
            2.315324434,
            2.056404127,
            1.81,
            1.64,
            14.94,
            13.63,
            10.85,
            17.87371609,
            10.19,
            17.399,
            13.407429659,
        ])
        pc = SingleParam("Critical pressure (pa)",components,[
            4599000.,
            3400000.,
            7383000.,
            4872000.,
            4248000.,
            3640000.,
            3796000.,
            3380000.,
            3370000.,
            3025000.,
            2740000.,
            2490000.,
            2290000.,
            2110000.,
            1313000.,
            5043000.,
            3499000.,
            22064000.,
            8962910.,
            117000.,
            4898000.,
        ])

    n0 = SingleParam("n⁰",components,[
        (19.597508817, -83.959667892, 3.00088, 0.76315, 0.00460, 8.74432, -4.46921), #methane
        (11.083407489, -22.202102428, 2.50031, 0.13732, -0.14660, 0.90066, 0.0),     #nitrogen
        (11.925152758, -16.118762264, 2.50002, 2.04452, -1.06044, 2.03366, 0.01393),  #carbon dioxide
        (24.675437527, -77.425313760, 3.00263, 4.33939, 1.23722, 13.19740, -6.01989), #ethane
        (31.602908195, -84.463284382, 3.02939, 6.60569, 3.19700, 19.19210, -8.37267), #propane
        (20.884143364, -91.638478026, 3.33944, 9.44893, 6.89406, 24.46180, 14.78240), #n-butane
        (20.413726078, -94.467620036, 3.06714, 8.97575, 5.25156, 25.14230, 16.13880), #iso-butane
        (14.536611217, -89.919548319, 3.00000, 8.95043, 21.83600, 33.40320, 0.0),     #n-pentane
        (15.449907693, -101.298172792, 3.00000, 11.76180, 20.11010, 33.16880, 0.0),   #isopentane
        (14.345969349, -96.165722367, 3.00000, 11.69770, 26.81420, 38.61640, 0.0),    #n-hexane
        (15.063786601, -97.345252349, 3.00000, 13.72660, 30.47070, 43.55610, 0.0),    #n-heptane
        (15.864687161, -97.370667555, 3.00000, 15.68650, 33.80290, 48.17310, 0.0),   #n-octane
        (16.313913248, -102.160247463, 3.00000, 18.02410, 38.12350, 53.34150, 0.0),  #n-nonane
        (15.870791919, -108.858547525, 3.00000, 21.00690, 43.49310, 58.36570, 0.0),       #n-decane
        (13.796443393, -175.864487294, 1.47906, 0.95806, 0.45444, 1.56039, -1.37560), #hydrogen
        (10.001843586, -14.996095135, 2.50146, 1.07558, 1.01334, 0.0, 0.0),           #oxygen
        (10.813340744, -19.834733959, 2.50055, 1.02865, 0.00493, 0.0, 0.0),           #carbon monoxide
        (8.203520690, -11.996306443, 3.00392, 0.01059, 0.98763, 3.06904, 0.0),            #water
        (9.336197742, -16.266508995, 3.00000, 3.11942, 1.00243, 0.0, 0.0),            #hydrogen sulfide
        (13.628409737, -143.470759602, 1.50000, 0.0, 0.0, 0.0, 0.0),                  #helium
        (8.316631500, -4.946502600, 1.50000, 0.0, 0.0, 0.0, 0.0),                      #argon
    ])

    theta = SingleParam("ϑ",components,[
        (4.306474465, 0.936220902, 5.577233895, 5.722644361),                  #methane
        (5.251822620, -5.393067706, 13.788988208, 0.),                        #nitrogen
        (3.022758166, -2.844425476, 1.589964364, 1.121596090),               #carbon dioxide
        (1.831882406, 0.731306621, 3.378007481, 3.508721939),                #ethane
        (1.297521801, 0.543210978, 2.583146083, 2.777773271),                #propane
        (1.101487798, 0.431957660, 4.502440459, 2.124516319),                #n-butane
        (1.074673199, 0.485556021, 4.671261865, 2.191583480),                #isobutane
        (0.380391739, 1.789520971, 3.777411113, 0.),                          #n-pentane
        (0.635392636, 1.977271641, 4.169371131, 0.),                          #isopentane
        (0.359036667, 1.691951873, 3.596924107, 0.),                          #n-hexane
        (0.314348398, 1.548136560, 3.259326458, 0.),                          #n-heptane
        (0.279143540, 1.431644769, 2.973845992, 0.),                          #n-octane
        (0.263819696, 1.370586158, 2.848860483, 0.),                          #n-nonane
        (0.267034159, 1.353835195, 2.833479035, 0.),                          #n-decane
        (6.891654113, 9.847634830, 49.765290750, 50.367279301),              #hydrogen
        (14.461722565, 7.223325463, 0., 0.),                                   #oxygen
        (11.669802800, 5.302762306, 0., 0.),                                   #carbon monoxide
        (0.415386589, 1.763895929, 3.874803739, 0.),                          #water
        (4.914580541, 2.270653980, 0., 0.),                                    #hydrogen sulfide
        (0.0, 0., 0.0, 0.0),                                                  #helium
        (0.0, 0., 0.0, 0.0),                                                  #argon
    ])

    gamma_T = PairParam("γ_T",components,
    [1.0       0.979273  0.975665  1.04971   1.09866   1.17161   1.16112   1.24917   1.1889    1.33057   1.37985   1.44925   1.52853  1.56823   1.35264  0.95      0.987473  0.775811  0.961156  3.15978  0.989843
    0.979273  1.0       1.10765   1.09865   1.20126   1.30489   1.28446   1.41903   1.38177   1.47261   1.52098   1.73328   1.74912  1.82216   1.1757   0.993425  0.968808  0.995157  0.960174  1.47183  0.990035
    0.975665  1.10765   1.0       0.90095   0.908772  0.911498  0.929983  0.979217  0.961218  1.03868   1.04617   1.07446   1.14067  1.14551   1.78292  1.0       1.0       0.775454  0.926019  3.20746  1.05097
    1.04971   1.09865   0.90095   1.0       1.01473   1.03483   1.03328   1.06667   1.02115   1.09218   1.13453   1.16864   1.05603  1.23719   1.90201  1.0       1.06922   1.0       0.902737  1.0      1.0
    1.09866   1.20126   0.908772  1.01473   1.0       1.00739   1.00525   1.00834   1.00327   1.02566   1.05004   1.06369   1.10997  1.14091   2.30822  1.0       1.19756   0.600341  0.905829  1.0      1.0
    1.17161   1.30489   0.911498  1.03483   1.00739   1.0       1.00143   1.00214   1.00079   1.00916   1.02128   1.03318   1.0141   1.07647   2.50926  1.0       1.17392   0.615513  0.926157  1.0      1.24504
    1.16112   1.28446   0.929983  1.03328   1.00525   1.00143   1.0       1.0025    1.00184   1.00602   1.00989   1.01395   1.01782  1.02162   1.89531  1.0       1.16139   1.0       0.937131  1.0      1.0
    1.24917   1.41903   0.979217  1.06667   1.00834   1.00214   1.0025    1.0       1.00005   1.00076   1.00244   1.01642   1.10342  1.04904   2.01386  1.0       1.20604   0.447666  0.959066  1.0      1.0
    1.1889    1.38177   0.961218  1.02115   1.00327   1.00079   1.00184   1.00005   1.0       1.0012    1.00319   1.00565   1.00819  1.01083   1.99639  1.0       1.19933   1.0       0.982652  1.0      1.0
    1.33057   1.47261   1.03868   1.09218   1.02566   1.00916   1.00602   1.00076   1.0012    1.0       0.999763  1.00163   1.05537  1.02894   3.0212   1.0       1.23327   0.569681  0.956076  1.0      1.0
    1.37985   1.52098   1.04617   1.13453   1.05004   1.02128   1.00989   1.00244   1.00319   0.999763  1.0       0.998793  1.00115  1.00223   3.16914  1.0       1.25612   1.0       1.01345   1.0      1.0
    1.44925   1.73328   1.07446   1.16864   1.06369   1.03318   1.01395   1.01642   1.00565   1.00163   0.998793  1.0       1.00024  1.00719   2.19156  1.0       1.27657   0.662072  1.0       1.0      1.0
    1.52853   1.74912   1.14067   1.05603   1.10997   1.0141    1.01782   1.10342   1.00819   1.05537   1.00115   1.00024   1.0      1.00018   2.23435  1.0       1.29407   1.0       1.08656   1.0      1.0
    1.56823   1.82216   1.14551   1.23719   1.14091   1.07647   1.02162   1.04904   1.01083   1.02894   1.00223   1.00719   1.00018  1.0       3.786    1.0       1.80357   0.740416  1.10369   1.0      1.0
    1.35264   1.1757    1.78292   1.90201   2.30822   2.50926   1.89531   2.01386   1.99639   3.0212    3.16914   2.19156   2.23435  3.786     1.0      1.0       1.3775    1.0       1.0       1.0      1.0
    0.95      0.993425  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0      1.0       1.0       0.964768  1.0       1.0      0.99043
    0.987473  0.968808  1.0       1.06922   1.19756   1.17392   1.16139   1.20604   1.19933   1.23327   1.25612   1.27657   1.29407  1.80357   1.3775   1.0       1.0       1.0       1.02275   1.0      0.954216
    0.775811  0.995157  0.775454  1.0       0.600341  0.615513  1.0       0.447666  1.0       0.569681  1.0       0.662072  1.0      0.740416  1.0      0.964768  1.0       1.0       0.940587  1.0      1.07094
    0.961156  0.960174  0.926019  0.902737  0.905829  0.926157  0.937131  0.959066  0.982652  0.956076  1.01345   1.0       1.08656  1.10369   1.0      1.0       1.02275   0.940587  1.0       1.0      1.0
    3.15978   1.47183   3.20746   1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0      1.0       1.0       1.0       1.0       1.0      1.0
    0.989843  0.990035  1.05097   1.0       1.0       1.24504   1.0       1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0      0.99043   0.954216  1.07094   1.0       1.0      1.0])

    gamma_v = PairParam("γ_v",components,[
    1.0       1.01395   1.00281   1.00662  1.03847   1.04538  1.05432   1.12451   1.34369   1.05264   1.15666  1.11655   1.1419    1.14609   1.0187    1.0       1.0061    1.58502   1.04016   0.881406  1.01468
    1.01395   1.0       1.04758   1.04235  1.08103   1.14695  1.10058   1.07888   1.15414   1.19595   1.40455  1.18607   1.10041   1.0       0.970115  1.00869   1.09475   0.997082  1.25684   0.93263   1.00221
    1.00281   1.04758   1.0       1.03288  1.0476    1.22244  1.08191   1.06841   1.11679   0.851344  1.16459  1.10404   0.973386  1.18339   1.15279   1.0       1.0       1.54233   1.02409   0.864142  1.02921
    1.00662   1.04235   1.03288   1.0      1.00303   1.00618  1.00662   1.02609   1.04544   1.1697    1.05767  1.07192   1.14353   1.09836   1.10607   1.0       1.20142   1.0       1.03099   1.0       1.0
    1.03847   1.08103   1.0476    1.00303  1.0       1.00326  1.00116   1.01992   0.999432  1.05787   1.07965  1.10276   1.19977   1.05304   1.07401   1.0       1.10814   1.01176   1.01059   1.0       1.0
    1.04538   1.14695   1.22244   1.00618  1.00326   1.0      1.00041   1.01816   1.00273   1.035     1.01917  1.04691   1.04922   1.02785   1.23294   1.0       1.08474   1.22364   1.03337   1.0       1.21464
    1.05432   1.10058   1.08191   1.00662  1.00116   1.00041  1.0       1.00278   1.00228   1.01049   1.02167  1.03281   1.0473    1.06024   1.1476    1.0       1.08727   1.0       0.988591  1.0       1.0
    1.12451   1.07888   1.06841   1.02609  1.01992   1.01816  1.00278   1.0       1.00002   1.00248   1.00897  1.06922   1.03491   1.01637   1.18833   1.0       1.11995   0.956677  1.07654   1.0       1.0
    1.34369   1.15414   1.11679   1.04544  0.999432  1.00273  1.00228   1.00002   1.0       1.003     1.00993  1.01788   1.02899   1.03937   1.18434   1.0       1.11669   1.0       0.835763  1.0       1.0
    1.05264   1.19595   0.851344  1.1697   1.05787   1.035    1.01049   1.00248   1.003     1.0       1.00151  1.00627   1.02076   1.01351   1.24346   1.0       1.15515   1.17022   1.33928   1.0       1.0
    1.15666   1.40455   1.16459   1.05767  1.07965   1.01917  1.02167   1.00897   1.00993   1.00151   1.0      1.00677   1.00137   1.00297   1.15913   1.0       1.19035   1.0       1.08796   1.0       1.0
    1.11655   1.18607   1.10404   1.07192  1.10276   1.04691  1.03281   1.06922   1.01788   1.00627   1.00677  1.0       1.00136   1.00255   1.30525   1.0       1.21921   0.599484  1.0       1.0       1.0
    1.1419    1.10041   0.973386  1.14353  1.19977   1.04922  1.0473    1.03491   1.02899   1.02076   1.00137  1.00136   1.0       1.00081   1.34265   1.0       1.25215   1.0       1.08291   1.0       1.0
    1.14609   1.0       1.18339   1.09836  1.05304   1.02785  1.06024   1.01637   1.03937   1.01351   1.00297  1.00255   1.00081   1.0       1.12023   1.0       0.870185  0.551405  1.17171   1.0       1.0
    1.0187    0.970115  1.15279   1.10607  1.07401   1.23294  1.1476    1.18833   1.18434   1.24346   1.15913  1.30525   1.34265   1.12023   1.0       1.0       1.12142   1.0       1.0       1.0       1.0
    1.0       1.00869   1.0       1.0      1.0       1.0      1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.14317   1.0       1.0       0.993907
    1.0061    1.09475   1.0       1.20142  1.10814   1.08474  1.08727   1.11995   1.11669   1.15515   1.19035  1.21921   1.25215   0.870185  1.12142   1.0       1.0       1.0       1.10173   1.0       1.15972
    1.58502   0.997082  1.54233   1.0      1.01176   1.22364  1.0       0.956677  1.0       1.17022   1.0      0.599484  1.0       0.551405  1.0       1.14317   1.0       1.0       1.01483   1.0       1.03899
    1.04016   1.25684   1.02409   1.03099  1.01059   1.03337  0.988591  1.07654   0.835763  1.33928   1.08796  1.0       1.08291   1.17171   1.0       1.0       1.10173   1.01483   1.0       1.0       1.0
    0.881406  0.93263   0.864142  1.0      1.0       1.0      1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
    1.01468   1.00221   1.02921   1.0      1.0       1.21464  1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0       0.993907  1.15972   1.03899   1.0       1.0       1.0])

    beta_T = PairParam("β_T",components,[
    1.0       0.998099  1.02262   0.996337  0.98968  0.994175  0.980316  0.992128  1.0       0.981845  0.977432  0.957474  0.947717  0.937778  1.0       1.0       0.987412  1.06333   1.01109   1.0       0.990954
    1.0019    1.0       1.00589   1.00767   1.00268  0.994515  0.992868  1.0       1.0       1.0       1.0       1.0       0.956379  0.957934  0.946134  1.0       1.0       0.997191  1.00469   0.692869  0.99907
    0.977876  0.99414   1.0       1.01387   1.03362  1.01817   1.02334   1.027     1.01918   1.0       1.01181   1.02969   1.00769   1.02003   0.94232   1.0       1.0       0.997372  1.01603   0.768378  0.996513
    1.00368   0.992387  0.986319  1.0       0.9962   0.999131  1.0       0.998689  1.0       1.0       1.0       0.984068  1.0       0.970918  0.93297   1.0       1.0       1.0       0.990197  1.0       1.0
    1.01043   0.99733   0.967473  1.00381   1.0      1.00031   0.998012  0.996484  0.994364  1.0       1.0       1.0       1.0       0.985331  1.0       1.0       1.0       1.0       0.992574  1.0       1.0
    1.00586   1.00552   0.982153  1.00087   0.99969  1.0       1.00008   1.0       1.0       1.0       1.0       1.0       1.0       0.993688  1.0       1.0       1.0       1.0       0.985963  1.0       1.0
    1.02008   1.00718   0.977192  1.0       1.00199  0.999922  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       0.974551  1.0       1.0
    1.00793   1.0       0.973709  1.00131   1.00353  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       0.962007  1.0       1.0
    1.0       1.0       0.98118   1.0       1.00567  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
    1.01849   1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       0.997641  1.0       1.0       1.0       1.0       0.985891  1.0       1.0
    1.02309   1.0       0.988331  1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       0.988937  1.0       1.0
    1.04442   1.0       0.971165  1.01619   1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
    1.05517   1.04561   0.99237   1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
    1.06635   1.04391   0.980364  1.02995   1.01489  1.00635   1.0       1.0       1.0       1.00236   1.0       1.0       1.0       1.0       1.06482   1.0       1.04959   0.897162  0.973091  1.0       1.0
    1.0       1.05693   1.06121   1.07185   1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       0.939128  1.0       1.0       1.0       1.0       1.0       1.0       1.0
    1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.00002
    1.01275   1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       0.952749  1.0       1.0       1.0       1.0       1.02554   1.0       1.0
    0.940438  1.00282   1.00263   1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.11463   1.0       1.0       1.0       1.0       1.0       1.0       1.0
    0.989032  0.99533   0.984218  1.0099    1.00748  1.01424   1.02611   1.03949   1.0       1.01431   1.01119   1.0       1.0       1.02765   1.0       1.0       0.975099  1.0       1.0       1.0       1.0
    1.0       1.44327   1.30144   1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0
    1.00913   1.00093   1.0035    1.0       1.0      1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       0.999977  1.0       1.0       1.0       1.0       1.0])

    beta_v = PairParam("β_v",components,[
    1.0       0.998721  0.999518  0.997548  1.00483   0.979106  1.01124   0.94833   1.0      0.958015  0.962051  0.994741  1.00285  1.03309   1.0       1.0      0.997341  1.01278   1.0126    1.0       1.03463
    1.00128   1.0       0.977795  0.97888   0.974425  0.996083  0.986416  1.0       1.0      1.0       1.0       1.0       1.0      1.0       0.972532  1.0      1.0       0.999522  0.910394  0.969501  1.00417
    1.00048   1.02271   1.0       1.00253   0.996898  1.17476   1.07655   1.02431   1.06079  1.0       1.20547   1.02617   1.0      1.00015   0.904142  1.0      1.0       0.949056  0.906631  0.846648  1.00839
    1.00246   1.02158   0.997481  1.0       0.997607  0.999157  1.0       0.993851  1.0      1.0       1.0       1.00747   1.0      0.995676  0.925367  1.0      1.0       1.0       1.01082   1.0       1.0
    0.995196  1.02625   1.00311   1.0024    1.0       0.999796  0.999243  1.04492   1.04046  1.0       1.0       1.0       1.0      0.984104  1.0       1.0      1.0       1.0       0.936811  1.0       1.0
    1.02134   1.00393   0.851237  1.00084   1.0002    1.0       1.00088   1.0       1.0      1.0       1.0       1.0       1.0      0.976952  1.0       1.0      1.0       1.0       0.908113  1.0       1.0
    0.988885  1.01377   0.928892  1.0       1.00076   0.99912   1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       1.01299   1.0       1.0
    1.05449   1.0       0.976266  1.00619   0.957012  1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       0.984613  1.0       1.0
    1.0       1.0       0.942691  1.0       0.961114  1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0
    1.04382   1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.00152   1.0       1.0      1.0       1.0       0.754474  1.0       1.0
    1.03945   1.0       0.829552  1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       0.828967  1.0       1.0
    1.00529   1.0       0.974498  0.992586  1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0
    0.997156  1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0
    0.967973  1.0       0.999849  1.00434   1.01615   1.02359   1.0       1.0       1.0      0.998486  1.0       1.0       1.0      1.0       1.69536   1.0      1.0       1.0       0.975188  1.0       1.0
    1.0       1.02824   1.10602   1.08065   1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      0.589846  1.0       1.0      1.0       1.0       1.0       1.0       1.0
    1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       1.0       1.0       0.999747
    1.00267   1.0       1.0       1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       0.79566   1.0       1.0
    0.987378  1.00048   1.05368   1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0
    0.987558  1.09843   1.10299   0.989298  1.06745   1.10118   0.987172  1.01563   1.0      1.32543   1.20632   1.0       1.0      1.02544   1.0       1.0      1.25682   1.0       1.0       1.0       1.0
    1.0       1.03146   1.18113   1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.0      1.0       1.0       1.0       1.0       1.0
    0.966529  0.995851  0.991677  1.0       1.0       1.0       1.0       1.0       1.0      1.0       1.0       1.0       1.0      1.0       1.0       1.00025  1.0       1.0       1.0       1.0       1.0])

    n = SingleParam("nᵢ",components, [
    [0.57335704239162, -1.676068752373, 0.23405291834916, -0.21947376343441, 0.016369201404128, 0.01500440638928, 0.098990489492918, 0.58382770929055, -0.7478686756039, 0.30033302857974, 0.20985543806568, -0.018590151133061, -0.15782558339049, 0.12716735220791, -0.032019743894346, -0.068049729364536, 0.024291412853736, 0.0051440451639444, -0.019084949733532, 0.0055229677241291, -0.0044197392976085, 0.040061416708429, -0.033752085907575, -0.0025127658213357],
    [0.59889711801201, -1.6941557480731, 0.24579736191718, -0.23722456755175, 0.017954918715141, 0.014592875720215, 0.10008065936206, 0.73157115385532, -0.88372272336366, 0.31887660246708, 0.20766491728799, -0.019379315454158, -0.16936641554983, 0.13546846041701, -0.033066712095307, -0.060690817018557, 0.012797548292871, 0.0058743664107299, -0.018451951971969, 0.0047226622042472, -0.0052024079680599, 0.043563505956635, -0.036251690750939, -0.0028974026866543],
    [0.52646564804653, -1.4995725042592, 0.27329786733782, 0.12949500022786, 0.15404088341841, -0.58186950946814, -0.18022494838296, -0.095389904072812, -0.0080486819317679, -0.03554775127309, -0.28079014882405, -0.082435890081677, 0.010832427979006, -0.0067073993161097, -0.0046827907600524, -0.028359911832177, 0.019500174744098, -0.21609137507166, 0.43772794926972, -0.22130790113593, 0.015190189957331, -0.0153809489533],
    [0.63596780450714, -1.7377981785459, 0.28914060926272, -0.33714276845694, 0.022405964699561, 0.015715424886913, 0.11450634253745, 1.0612049379745, -1.2855224439423, 0.39414630777652, 0.31390924682041, -0.021592277117247, -0.21723666564905, -0.28999574439489, 0.42321173025732, 0.04643410025926, -0.13138398329741, 0.011492850364368, -0.033387688429909, 0.015183171583644, -0.0047610805647657, 0.046917166277885, -0.039401755804649, -0.0032569956247611],
    [1.0403973107358, -2.8318404081403, 0.84393809606294, -0.076559591850023, 0.09469737305728, 0.00024796475497006, 0.2774376042287, -0.043846000648377, -0.2699106478435, -0.06931341308986, -0.029632145981653, 0.01404012675138],
    [1.0626277411455, -2.862095182835, 0.88738233403777, -0.12570581155345, 0.10286308708106, 0.00025358040602654, 0.32325200233982, -0.037950761057432, -0.32534802014452, -0.079050969051011, -0.020636720547775, 0.005705380933475],
    [1.04293315891, -2.8184272548892, 0.8617623239785, -0.10613619452487, 0.098615749302134, 0.00023948208682322, 0.3033000485695, -0.041598156135099, -0.29991937470058, -0.080369342764109, -0.029761373251151, 0.01305963030314],
    [1.0968643098001, -2.9988888298061, 0.99516886799212, -0.16170708558539, 0.11334460072775, 0.00026760595150748, 0.40979881986931, -0.040876423083075, -0.38169482469447, -0.10931956843993, -0.03207322332799, 0.016877016216975],
    [1.0963, -3.0402, 1.0317, -0.1541, 0.11535, 0.00029809, 0.39571, -0.045881, -0.35804, -0.10107, -0.035484, 0.018156],
    [1.0553238013661, -2.6120615890629, 0.7661388296726, -0.29770320622459, 0.11879907733358, 0.00027922861062617, 0.46347589844105, 0.011433196980297, -0.48256968738131, -0.093750558924659, -0.0067273247155994, -0.0051141583585428],
    [1.0543747645262, -2.6500681506144, 0.81730047827543, -0.30451391253428, 0.122538687108, 0.00027266472743928, 0.4986582568167, -0.00071432815084176, -0.5423689552545, -0.13801821610756, -0.0061595287380011, 0.00048602510393022],
    [1.0722544875633, -2.4632951172003, 0.65386674054928, -0.36324974085628, 0.12713269626764, 0.0003071357277793, 0.5265685698754, 0.019362862857653, -0.58939426849155, -0.14069963991934, -0.0078966330500036, 0.0033036597968109],
    [1.1151, -2.702, 0.83416, -0.38828, 0.1376, 0.00028185, 0.62037, 0.015847, -0.61726, -0.15043, -0.012982, 0.0044325],
    [1.0461, -2.4807, 0.74372, -0.52579, 0.15315, 0.00032865, 0.84178, 0.055424, -0.73555, -0.18507, -0.020775, 0.012335],
    [5.3579928451252, -6.2050252530595, 0.13830241327086, -0.071397954896129, 0.015474053959733, -0.14976806405771, -0.026368723988451, 0.056681303156066, -0.060063958030436, -0.45043942027132, 0.424788402445, -0.021997640827139, -0.01049952137453, -0.0028955902866816],
    [0.88878286369701, -2.4879433312148, 0.59750190775886, 0.0096501817061881, 0.07197042871277, 0.00022337443000195, 0.18558686391474, -0.03812936803576, -0.15352245383006, -0.026726814910919, -0.025675298677127, 0.0095714302123668],
    [0.90554, -2.4515, 0.53149, 0.024173, 0.072156, 0.00018818, 0.19405, -0.043268, -0.12778, -0.027896, -0.034154, 0.016329],
    [0.82728408749586, -1.8602220416584, -1.1199009613744, 0.15635753976056, 0.87375844859025, -0.36674403715731, 0.053987893432436, 1.0957690214499, 0.053213037828563, 0.013050533930825, -0.41079520434476, 0.1463744334412, -0.055726838623719, -0.0112017741438, -0.0066062758068099, 0.0046918522004538],
    [0.87641, -2.0367, 0.21634, -0.050199, 0.066994, 0.00019076, 0.20227, -0.0045348, -0.2223, -0.034714, -0.014885, 0.0074154],
    [-0.45579024006737, 1.2516390754925, -1.5438231650621, 0.020467489707221, -0.34476212380781, -0.020858459512787, 0.016227414711778, -0.057471818200892, 0.019462416430715, -0.03329568012302, -0.010863577372367, -0.022173365245954],
    [0.85095714803969, -2.400322294348, 0.54127841476466, 0.016919770692538, 0.068825965019035, 0.00021428032815338, 0.17429895321992, -0.033654495604194, -0.13526799857691, -0.016387350791552, -0.024987666851475, 0.0088769204815709],
    ])
    t = SingleParam("tᵢ",components,[
    [0.125, 1.125, 0.375, 1.125, 0.625, 1.5, 0.625, 2.625, 2.75, 2.125, 2.0, 1.75, 4.5, 4.75, 5.0, 4.0, 4.5, 7.5, 14.0, 11.5, 26.0, 28.0, 30.0, 16.0],
    [0.125, 1.125, 0.375, 1.125, 0.625, 1.5, 0.625, 2.625, 2.75, 2.125, 2.0, 1.75, 4.5, 4.75, 5.0, 4.0, 4.5, 7.5, 14.0, 11.5, 26.0, 28.0, 30.0, 16.0],
    [0.0, 1.25, 1.625, 0.375, 0.375, 1.375, 1.125, 1.375, 0.125, 1.625, 3.75, 3.5, 7.5, 8.0, 6.0, 16.0, 11.0, 24.0, 26.0, 28.0, 24.0, 26.0],
    [0.125, 1.125, 0.375, 1.125, 0.625, 1.5, 0.625, 2.625, 2.75, 2.125, 2.0, 1.75, 4.5, 4.75, 5.0, 4.0, 4.5, 7.5, 14.0, 11.5, 26.0, 28.0, 30.0, 16.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.5, 0.625, 0.375, 0.625, 1.125, 2.625, 0.0, 0.25, 1.375, 4.0, 4.25, 5.0, 8.0, 8.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.5, 1.25, 1.875, 0.125, 1.5, 1.0, 0.75, 1.5, 0.625, 2.625, 5.0, 4.0, 4.5, 3.0, 4.0, 6.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    [0.0, 0.125, 0.75, 1.0, 0.75, 2.625, 0.125, 1.25, 2.0, 1.0, 4.5, 5.0],
    [0.25, 1.125, 1.5, 1.375, 0.25, 0.875, 0.625, 1.75, 3.625, 3.625, 14.5, 12.0],
    ])

    d = SingleParam("dᵢ",components,
    [
    [1, 1, 2, 2, 4, 4, 1, 1, 1, 2, 3, 6, 2, 3, 3, 4, 4, 2, 3, 4, 5, 6, 6, 7],
    [1, 1, 2, 2, 4, 4, 1, 1, 1, 2, 3, 6, 2, 3, 3, 4, 4, 2, 3, 4, 5, 6, 6, 7],
    [1, 1, 2, 3, 3, 3, 4, 5, 6, 6, 1, 4, 1, 1, 3, 3, 4, 5, 5, 5, 5, 5],
    [1, 1, 2, 2, 4, 4, 1, 1, 1, 2, 3, 6, 2, 3, 3, 4, 4, 2, 3, 4, 5, 6, 6, 7],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 2, 2, 4, 1, 5, 5, 5, 1, 1, 2, 5, 1],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 2, 2, 3, 4, 1, 5, 5, 1, 2, 4, 4, 1, 1],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    [1, 1, 1, 4, 1, 3, 5, 5, 5, 2, 1, 2],
    [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4],
    ])

    c = SingleParam("cᵢ",components,
    [[1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 6, 6, 6, 6],
    [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 6, 6, 6, 6],
    [1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 5, 5, 5, 6, 6],
    [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 6, 6, 6, 6],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 1, 1, 2, 2, 3, 3, 5],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 1, 2, 2, 2, 3, 5, 5],
    [1, 1, 2, 2, 3, 3],
    [1, 1, 1, 1, 1, 2, 3, 3],
    [1, 1, 2, 2, 3, 3]])

    idx1 = [1, 1, 1, 1, 1, 1, 1, 2, 2, 4, 4, 4, 5, 5, 6]
    idx2 = [2, 3, 4, 5, 6, 7, 15, 3, 4, 5, 6, 7, 6, 7, 7]

    Fij = PairParam("Fᵢⱼ",components,
    sparse(idx1,idx2,[
    1.0,
    1.0,
    1.0,
    1.0,
    1.0,
    0.771035405688,
    1.0,
    1.0,
    1.0,
    0.130424765150,
    0.281570073085,
    0.260632376098,
    0.312572600489e-1,
    -0.551609771024e-1,
    -0.551240293009e-1,
    ],21,21))

    dij = PairParam("dᵢⱼ",components,
    sparse(idx1,idx2,
    [
    [1, 4, 1, 2, 2, 2, 2, 2, 3],
    [1, 2, 3, 1, 2, 3],
    [3, 4, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3],
    [3, 3, 4, 4, 4, 1, 1, 1, 2],
    [2, 3, 1, 1, 1, 2],
    [2, 2, 3, 1, 2, 2],
    [1, 3, 3, 4],
    [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
    [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
    [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
    [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
    [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
    [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
    [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
    [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
    ],21,21
    ))

    tij = PairParam("tᵢⱼ",components,
    sparse(idx1,idx2,[
    [0.000, 1.850, 7.850, 5.400, 0.000, 0.750, 2.800, 4.450, 4.250],
    [2.600, 1.950, 0.000, 3.950, 7.950, 8.000],
    [0.650, 1.550, 3.100, 5.900, 7.050, 3.350, 1.200, 5.800, 2.700, 0.450, 0.550, 1.950],
    [1.850, 3.950, 0.000, 1.850, 3.850, 5.250, 3.850, 0.200, 6.500],
    [1.850, 1.400, 3.200, 2.500, 8.000, 3.750],
    [0.000, 0.050, 0.000, 3.650, 4.900, 4.450],
    [2.0, -1.0, 1.75, 1.4],
    [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
    [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
    [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
    [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
    [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
    [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
    [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
    [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
    ]
    ,21,21))

    nij = PairParam("nᵢⱼ",components,
    sparse(idx1,idx2,
    [
    [-0.0098038985517335, 0.00042487270143005, -0.034800214576142, -0.13333813013896, -0.011993694974627, 0.069243379775168, -0.31022508148249, 0.24495491753226, 0.22369816716981],
    [-0.10859387354942, 0.080228576727389, -0.0093303985115717, 0.040989274005848, -0.24338019772494, 0.23855347281124],
    [-0.00080926050298746, -0.00075381925080059, -0.041618768891219, -0.23452173681569, 0.14003840584586, 0.063281744807738, -0.034660425848809, -0.23918747334251, 0.0019855255066891, 6.1777746171555, -6.9575358271105, 1.0630185306388],
    [0.013746429958576, -0.0074425012129552, -0.0045516600213685, -0.0054546603350237, 0.0023682016824471, 0.18007763721438, -0.44773942932486, 0.0193273748882, -0.30632197804624],
    [0.28661625028399, -0.10919833861247, -1.137403208227, 0.76580544237358, 0.0042638000926819, 0.17673538204534],
    [-0.47376518126608, 0.48961193461001, -0.0057011062090535, -0.1996682004132, -0.69411103101723, 0.69226192739021],
    [-0.25157134971934, -0.0062203841111983, 0.088850315184396, -0.035592212573239],
    [2.5574776844118, -7.9846357136353, 4.7859131465806, -0.73265392369587, 1.3805471345312, 0.28349603476365, -0.49087385940425, -0.10291888921447, 0.11836314681968, 5.5527385721943e-5],
    [2.5574776844118, -7.9846357136353, 4.7859131465806, -0.73265392369587, 1.3805471345312, 0.28349603476365, -0.49087385940425, -0.10291888921447, 0.11836314681968, 5.5527385721943e-5],
    [2.5574776844118, -7.9846357136353, 4.7859131465806, -0.73265392369587, 1.3805471345312, 0.28349603476365, -0.49087385940425, -0.10291888921447, 0.11836314681968, 5.5527385721943e-5],
    [2.5574776844118, -7.9846357136353, 4.7859131465806, -0.73265392369587, 1.3805471345312, 0.28349603476365, -0.49087385940425, -0.10291888921447, 0.11836314681968, 5.5527385721943e-5],
    [2.5574776844118, -7.9846357136353, 4.7859131465806, -0.73265392369587, 1.3805471345312, 0.28349603476365, -0.49087385940425, -0.10291888921447, 0.11836314681968, 5.5527385721943e-5],
    [2.5574776844118, -7.9846357136353, 4.7859131465806, -0.73265392369587, 1.3805471345312, 0.28349603476365, -0.49087385940425, -0.10291888921447, 0.11836314681968, 5.5527385721943e-5],
    [2.5574776844118, -7.9846357136353, 4.7859131465806, -0.73265392369587, 1.3805471345312, 0.28349603476365, -0.49087385940425, -0.10291888921447, 0.11836314681968, 5.5527385721943e-5],
    [2.5574776844118, -7.9846357136353, 4.7859131465806, -0.73265392369587, 1.3805471345312, 0.28349603476365, -0.49087385940425, -0.10291888921447, 0.11836314681968, 5.5527385721943e-5],
    ],21,21))

    eta_ij = PairParam("ηᵢⱼ",components,
    sparse(idx1,idx2,
    [
    [1.000, 1.000, 0.250, 0.000, 0.000, 0.000, 0.000],
    [1.000, 0.500, 0.000],
    [1.000, 1.000, 1.000, 0.875, 0.750, 0.500, 0.000, 0.000, 0.000, 0.000],
    [0.250, 0.250, 0.000, 0.000],
    [0.250, 0.250, 0.000, 0.000],
    [1.000, 1.000, 0.875],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    ],
     21, 21)
    )

    epsilon_ij = PairParam("εᵢⱼ",components,
    sparse(idx1,idx2,
    [
    [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    ],
    21, 21)
    )

    beta_ij = PairParam("εᵢⱼ",components,
    sparse(idx1,idx2,
    [
    [1.000, 1.000, 2.500, 3.000, 3.000, 3.000, 3.000],
    [1.000, 2.000, 3.000],
    [1.000, 1.000, 1.000, 1.250, 1.500, 2.000, 3.000, 3.000, 3.000, 3.000],
    [0.750, 1.000, 2.000, 3.000],
    [0.750, 1.000, 2.000, 3.000],
    [1.000, 1.000, 1.250],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    ],
    21, 21)
    )
    gamma_ij = PairParam("γᵢⱼ",components,
    sparse(idx1, idx2,
    [
    [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5, 0.5],
    [0.5, 0.5, 0.5],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    Float64[],
    ],
    21, 21)
    )
    return Dict{Symbol,Any}(
        :Mw => Mw,
        :rhoc =>rhoc,
        :Tc =>Tc,
        :pc =>pc,
        :lb_v =>lb_v,
        :n0 =>n0,
        :theta =>theta,
        :gamma_T =>gamma_T,
        :gamma_v =>gamma_v,
        :beta_T =>beta_T,
        :beta_v =>beta_v,
        :n =>n,
        :t =>t,
        :d =>d,
        :c =>c,
        :nij =>nij,
        :tij =>tij,
        :dij =>dij,
        :Fij =>Fij,
        :eta_ij =>eta_ij,
        :epsilon_ij =>epsilon_ij,
        :beta_ij =>beta_ij,
        :gamma_ij =>gamma_ij,
    )
end

function getnames_gerg2008(components::Vector{String})
    comps = normalisestring.(components)
    names = normalisestring.(GERG2008_names)
    @assert allunique(comps)
    res = indexin(comps,names)
    if any(isnothing,res)
        idx = findfirst(isnothing,res)
        throw(error("component $(components[idx]) could not be found"))
    end
    return res
end

function f(s)
    state = zeros(Int, 2)
    state = state
    code = Meta.parse.(s)
    for c in code
        println(@macroexpand(@eval($c)))

    end
    return state
end


text:
__precompile__()

module MittLeffLite

import SpecialFunctions.gamma

export mittleff
export mittlefferr
export mittleffderiv

include("mittag_leffler.jl")

end # module

ext: jl


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2DAD.jl/blob/master/LICENSE

using ForwardDiff: value

mutable struct Contact2DAD <: BoundaryProblem
    master_elements :: Vector{Element}
    rotate_normals :: Bool
    dual_basis :: Bool
    iteration :: Int
    contact_state_in_first_iteration :: Symbol
    always_in_contact :: Set{Int64}
    use_scaling :: Bool
    alpha :: Dict{Int64, Real}
    beta :: Dict{Int64, Real}
    nadj_nodes :: Dict{Int64, Real}
    scaling_factors :: Dict{Int64, Real}
    print_summary :: Bool
end

function Contact2DAD()
    return Contact2DAD([], false, true, 0, :AUTO, Set(),
                        true, Dict(), Dict(), Dict(), Dict(), false)
end

function FEMBase.add_elements!(::Problem{Contact2DAD}, ::Any)
    error("use `add_slave_elements!` and `add_master_elements!` to add ",
          "elements to the Mortar2D problem.")
end

function FEMBase.add_slave_elements!(problem::Problem{Contact2DAD}, elements)
    for element in elements
        push!(problem.elements, element)
    end
end

function FEMBase.add_master_elements!(problem::Problem{Contact2DAD}, elements)
    for element in elements
        push!(problem.properties.master_elements, element)
    end
end

function FEMBase.get_slave_elements(problem::Problem{Contact2DAD})
    return problem.elements
end

function FEMBase.get_master_elements(problem::Problem{Contact2DAD})
    return problem.properties.master_elements
end

function get_slave_dofs(problem::Problem{Contact2DAD})
    dofs = Int64[]
    for element in get_slave_elements(problem)
        append!(dofs, get_gdofs(problem, element))
    end
    return sort(unique(dofs))
end

function get_master_dofs(problem::Problem{Contact2DAD})
    dofs = Int64[]
    for element in get_master_elements(problem)
        append!(dofs, get_gdofs(problem, element))
    end
    return sort(unique(dofs))
end

""" Find segment from slave element corresponding to master element nodes.

Parameters
----------
x1_, n1_
    slave element geometry and normal direction
x2
    master element node to project onto slave

Returns
-------
xi
    dimensionless coordinate on slave corresponding to
    projected master

"""
function project_from_master_to_slave_ad(
    slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector;
    tol=1.0e-10, max_iterations=20, debug=false) where E<:MortarElements2D

    """ Multiply basis / dbasis at `xi` with field. """
    function mul(func, xi, field)
        B = func(slave_element, [xi], time)
        return sum(B[i]*field[i] for i=1:length(B))
    end

    x1(xi1) = mul(get_basis, xi1, x1_)
    dx1(xi1) = mul(get_dbasis, xi1, x1_)
    n1(xi1) = mul(get_basis, xi1, n1_)
    dn1(xi1) = mul(get_dbasis, xi1, n1_)
    cross2(a, b) = cross([a; 0], [b; 0])[3]
    R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
    dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))

    xi1 = 0.0
    xi1_next = 0.0
    dxi1 = 0.0
    for i=1:max_iterations
        dxi1 = -R(xi1)/dR(xi1)
        dxi1 = clamp.(dxi1, -0.3, 0.3)
        xi1_next = clamp.(xi1 + dxi1, -1.0, 1.0)
        if norm(xi1_next - xi1) < tol
            return xi1_next
        end
        if debug
            @info("Debug data for iteration $i")
            @info("xi1 = $(value(xi1))")
            @info("R(xi1) = $(value(R(xi1)))")
            @info("dR(xi1) = $(value(dR(xi1)))")
            @info("dxi1 = $(value(dxi1))")
            @info("norm = $(value(norm(xi1_next - xi1)))")
            @info("xi1_next = $(value(xi1_next))")
        end
        xi1 = xi1_next
    end

    @info("Projection algorithm failed with the following data:")
    a, b = x1_.data
    @info("x11 = $(map(value, a))")
    @info("x12 = $(map(value, b))")
    a, b = n1_.data
    @info("n11 = $(map(value, a))")
    @info("n12 = $(map(value, b))")
    @info("x2 = $(map(value, x2))")
    @info("xi1 = $(value(xi1)), dxi1 = $(value(dxi1))")
    @info("-R(xi1) = $(value(-R(xi1)))")
    @info("dR(xi1) = $(value(dR(xi1)))")
    error("find projection from master to slave: did not converge")

end

function project_from_slave_to_master_ad(
    master_element::Element{E}, x1, n1, x2_;
    tol=1.0e-10, max_iterations=20) where E<:MortarElements2D

    x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_)
    dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_)
    cross2(a, b) = cross([a; 0], [b; 0])[3]
    R(xi2) = cross2(x2(xi2)-x1, n1)
    dR(xi2) = cross2(dx2(xi2), n1)

    xi2 = 0.0
    dxi2 = 0.0
    for i=1:max_iterations
        dxi2 = -R(xi2) / dR(xi2)
        xi2 += dxi2
        if norm(dxi2) < tol
            return xi2
        end
    end

    error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")

end

"""
Frictionless 2d finite sliding contact with forwarddiff.

true/false flags: finite_sliding, friction, use_forwarddiff
"""
function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Assembly,
                                    elements::Vector{Element{Seg2}}, time::Float64)

    props = problem.properties
    props.iteration += 1

    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    slave_elements = get_slave_elements(problem)

    function calculate_interface(x::Vector)

        ndofs = round(Int, length(x)/2)
        nnodes = round(Int, ndofs/field_dim)
        u = reshape(x[1:ndofs], field_dim, nnodes)
        la = reshape(x[ndofs+1:end], field_dim, nnodes)
        fc = zero(u)
        gap = zero(u)
        C = zero(la)
        S = Set{Int64}()

        # 1. update nodal normals for slave elements
        Q = [0.0 -1.0; 1.0 0.0]
        normals = zero(u)
        for element in slave_elements
            conn = get_connectivity(element)
            push!(S, conn...)
            gdofs = get_gdofs(problem, element)
            X_el = element("geometry", time)
            x_el = tuple( (X_el[i] + u[:,j] for (i,j) in enumerate(conn))... )
            #=
            for ip in get_integration_points(element, 3)
                dN = get_dbasis(element, ip, time)
                N = element(ip, time)
                t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
                normals[:, conn] += ip.weight*Q*t'*N
            end
            =#
            dN = get_dbasis(element, [0.0], time)
            t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
            n = Q*t'
            n /= norm(n)
            for c in conn
                normals[:,c] += n
            end
        end
        for i in 1:size(normals,2)
            normals[:,i] /= norm(normals[:,i])
        end
        # swap element normals in 2d if they point to inside of body
        if props.rotate_normals
            for i=1:size(normals,2)
                normals[:,i] = -normals[:,i]
            end
        end
        normals2 = Dict()
        for j in S
            normals2[j] = normals[:,j]
        end
        update!(slave_elements, "normal", time => normals2)

        alpha = empty!(problem.properties.alpha)
        beta = empty!(problem.properties.beta)
        nadj_nodes = empty!(problem.properties.nadj_nodes)

        # 2. loop all slave elements
        for slave_element in slave_elements

            slave_element_nodes = get_connectivity(slave_element)
            X1 = slave_element("geometry", time)
            u1 = ((u[:,i] for i in slave_element_nodes)...,)
            x1 = ((Xi+ui for (Xi,ui) in zip(X1,u1))...,)
            la1 = ((la[:,i] for i in slave_element_nodes)...,)
            n1 = ((normals[:,i] for i in slave_element_nodes)...,)
            nnodes = size(slave_element, 2)

            Ae = Matrix(1.0I, nnodes, nnodes)

            if problem.properties.dual_basis # construct dual basis

                De = zeros(nnodes, nnodes)
                Me = zeros(nnodes, nnodes)
                ae = zeros(nnodes)
                be = zeros(nnodes)

                for ip in get_integration_points(slave_element, 3)
                    dN = get_dbasis(slave_element, ip, time)
                    j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
                    detj = norm(j)
                    w = ip.weight*detj
                    N1 = slave_element(ip, time)
                    ae += w*vec(N1)/detj
                end

                nsegments = 0

                for master_element in get_master_elements(problem)

                    master_element_nodes = get_connectivity(master_element)
                    X2 = master_element("geometry", time)
                    u2 = ((u[:,i] for i in master_element_nodes)...,)
                    x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...,)

                    # calculate segmentation: we care only about endpoints
                    xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1])
                    xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2])
                    xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
                    l = 1/2*abs(xi1[2]-xi1[1])
                    isapprox(l, 0.0) && continue # no contribution in this master element

                    nsegments += 1
                    for ip in get_integration_points(slave_element, 3)
                        # jacobian of slave element in deformed state
                        dN = get_dbasis(slave_element, ip, time)
                        j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
                        detj = norm(j)
                        w = ip.weight*norm(j)*l
                        xi = ip.coords[1]
                        xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                        N1 = get_basis(slave_element, xi_s, time)
                        De += w*Matrix(Diagonal(vec(N1)))
                        Me += w*N1'*N1
                        be += w*vec(N1)/detj
                    end
                end

                for (i, j) in enumerate(get_connectivity(slave_element))
                    alpha[j] = get(alpha, j, 0.0) + ae[i]
                    beta[j] = get(beta, j, 0.0) + be[i]
                    nadj_nodes[j] = get(nadj_nodes, j, 0) + 1
                end

                if nsegments == 0
                    continue
                end

                Ae = De*inv(Me)

            end

            # 3. loop all master elements
            for master_element in get_master_elements(problem)

                master_element_nodes = get_connectivity(master_element)
                X2 = master_element("geometry", time)
                u2 = ((u[:,i] for i in master_element_nodes)...,)
                x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...,)

                #x1_midpoint = 1/2*(x1[1]+x1[2])
                #x2_midpoint = 1/2*(x2[1]+x2[2])
                #distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint))
                #distance > props.maximum_distance && continue

                # calculate segmentation: we care only about endpoints
                xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1])
                xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2])
                xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
                l = 1/2*abs(xi1[2]-xi1[1])
                isapprox(l, 0.0) && continue # no contribution in this master element

                slave_dofs = get_gdofs(problem, slave_element)
                master_dofs = get_gdofs(problem, master_element)

                # 4. loop integration points of segment
                for ip in get_integration_points(slave_element, 3)
                    # jacobian of slave element in deformed state
                    dN = get_dbasis(slave_element, ip, time)
                    j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
                    w = ip.weight*norm(j)*l

                    # project gauss point from slave element to master element
                    xi = ip.coords[1]
                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                    N1 = vec(get_basis(slave_element, xi_s, time))
                    x_s = interpolate(N1, x1) # coordinate in gauss point
                    n_s = interpolate(N1, n1) # normal direction in gauss point
                    t_s = Q'*n_s # tangent direction in gauss point
                    xi_m = project_from_slave_to_master_ad(master_element, x_s, n_s, x2)
                    N2 = vec(get_basis(master_element, xi_m, time))
                    x_m = interpolate(N2, x2)
                    Phi = Ae*N1

                    la_s = interpolate(Phi, la1) # traction force in gauss point
                    gn = -dot(n_s, x_s - x_m) # normal gap

                    fc[:,slave_element_nodes] += w*la_s*N1'
                    fc[:,master_element_nodes] -= w*la_s*N2'
                    gap[1,slave_element_nodes] += w*gn*Phi
                    #gap[1,slave_element_nodes] += w*gn*N1'

                end # done integrating segment

            end # master elements done

        end # slave elements done

        # at this point we have calculated contact force fc and gap for all slave elements.
        # next task is to find out are they in contact or not and remove inactive nodes

        state = problem.properties.contact_state_in_first_iteration
        if problem.properties.iteration == 1
            @info("First contact iteration, initial contact state = $state")
            if state == :AUTO
                avg_gap = ForwardDiff.value(mean([gap[1, j] for j in S]))
                std_gap = ForwardDiff.value(std([gap[1, j] for j in S]))
                if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12)
                    state = :ACTIVE
                else
                    state = :UNKNOWN
                end
                @info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
            end
        end

        is_active = Dict{Int, Bool}()
        condition = Dict()

        for j in S
            if j in props.always_in_contact
                is_active[j] = true
                continue
            end
            lan = dot(normals[:,j], la[:,j])
            condition[j] = ForwardDiff.value(lan - gap[1, j])
            is_active[j] = condition[j] > 0
        end

        if problem.properties.iteration == 1 && state == :ACTIVE
            for j in S
                is_active[j] = true
            end
        end

        if problem.properties.iteration == 1 && state == :INACTIVE
            for j in S
                is_active[j] = false
            end
        end

        if problem.properties.print_summary
            @info("Summary of nodes")
            for j in sort(collect(keys(is_active)))
                n = map(ForwardDiff.value, normals[:,j])
                @info("$j, c=$(condition[j]), s=$(is_active[j]), n=$n")
            end
        end

        for j in S

            if is_active[j]
                n = normals[:,j]
                t = Q'*n
                lan = dot(n, la[:,j])
                lat = dot(t, la[:,j])
                C[1,j] -= gap[1, j]
                C[2,j] += lat
            else
                C[:,j] = la[:,j]
            end

        end

        # Apply scaling to contact virtual work and contact constraints
        if problem.properties.use_scaling
            scaling = empty!(problem.properties.scaling_factors)
            for j in S
                isapprox(beta[j], 0.0) && continue
                #is_active[j] || continue
                #haskey(alpha, j) || continue
                #haskey(beta, j) || continue
                #haskey(nadj_nodes, j) || continue
                scaling[j] = alpha[j] / (nadj_nodes[j] * beta[j])
                C[1,j] *= scaling[j]
                fc[:,j] *= scaling[j]
            end
            update!(slave_elements, "contact scaling", time => scaling)
        end

        return vec([fc C])

    end

    # x doesn't mean deformed configuration here
    x = [problem.assembly.u; problem.assembly.la]
    if length(x) == 0
        error("2d autodiff contact problem: initialize problem.assembly.u & la before solution")
    end

    A = ForwardDiff.jacobian(calculate_interface, x)
    b = calculate_interface(x)
    A = sparse(A)
    b = -sparse(b)
    SparseArrays.droptol!(A, 1.0e-9)
    SparseArrays.droptol!(b, 1.0e-9)

    ndofs = round(Int, length(x)/2)
    K = A[1:ndofs,1:ndofs]
    C1 = A[1:ndofs,ndofs+1:end]
    C2 = A[ndofs+1:end,1:ndofs]
    D = A[ndofs+1:end,ndofs+1:end]
    f = b[1:ndofs]
    g = b[ndofs+1:end]

    f += C1*problem.assembly.la
    g += D*problem.assembly.la

    #=
    if !haskey(problem, "contact force")
        problem.fields["contact force"] = Field(time => f)
    else
        update!(problem.fields["contact force"], time => f)
    end

    fc = problem.fields["contact force"]

    if length(fc) > 1
        # kick in generalized alpha rule for time integration
        alpha = 0.5
        info("Applying Generalized alpha time integration")
        K = (1-alpha)*K
        C1 = (1-alpha)*C1
        f = alpha*fc[end-1].data
    end
    =#

    problem.assembly.K = K
    problem.assembly.C1 = copy(transpose(C1))
    problem.assembly.C2 = C2
    problem.assembly.D = D
    problem.assembly.f = sparse(f)
    problem.assembly.g = sparse(g)

end


text:
using Compat.Test
include("../ReflectBlock.jl")

@testset "Reflect" begin
    reg0 = rand_state(3)
    mirror = randn(1<<3)*im; mirror[:]/=norm(mirror)
    rf = Reflect(mirror)
    reg = copy(reg0)
    apply!(reg, rf)

    v0, v1 = vec(reg.state), vec(reg0.state)
    @test rf.state'*v0 ≈ rf.state'*v1
    @test v0-rf.state'*v0*rf.state ≈ -(v1-rf.state'*v1*rf.state)
end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
	using JellyFish

	@show jaro_winkler("jellyfish", "smellyfish")
	@show jaro_distance("jellyfish", "smellyfish")
	@show hamming_distance("apple", "pear")
	@show levenshtein_distance("jellyfish", "smellyfish")
	@show damerau_levenshtein_distance("jellyfish", "jellyfihs")
	@show soundex("jellyfish")
	@show metaphone("jellyfish")
	@show match_rating_codex("jellyfish")
	@show match_rating_comparison("jellyfish", "jellyfihs")
	@show nysiis("jellyfish")


text:
using Distributed
using SharedArrays

## fill table in parallel using shared arrays

function myrange_y(t::SharedArray,p)
    idx = indexpids(t)
    if idx == 0 # This worker is not assigned a piece
       return 1:0, 1:0
    end
    nchunks = length(procs(t))
    splits = [round(Int, s) for s in range(0, stop=size(t,1), length=nchunks+1)]
    years = [round(Int, y) for y in range(0, stop=size(p,1), length=nchunks+1)]
    splits[idx]+1:splits[idx+1], years[idx]+1:12:years[idx+1]
end

function fill_table!(t::SharedArray, pt, qt)
    split, years = myrange_y(t,pt)
    for (i,y) in zip(split, years)
        t[i,1] == 0.0 || continue # don't do previously done evals
        @show (myid(),i ,y)
        p = pt[y:y+11]
        q = qt[y:y+11,:]
        mchain, status = gchain(p,q)
        status || continue # don't do if infeasible point
        t[i,1], t[i,2], t[i,3] = RP_Cournot_test(mchain)
    end
    t
end

function shared_fill_table!(t, p, q)
    @sync begin
        for w in procs(t)
            @async remotecall_wait(fill_table!, w, t, p, q)
        end
    end
    t
end;

## parallel chain

# Semesters
function myrange_sem(t::SharedArray,p)
    idx = indexpids(t)
    if idx == 0 # This worker is not assigned a piece
       return 1:0, 1:0
    end
    nchunks = length(procs(t))
    splits = [round(Int, s) for s in range(0, stop=size(t,1), length=nchunks+1)]
    years = [round(Int, y) for y in range(0, stop=size(p,1), length=nchunks+1)]
    splits[idx]+1:splits[idx+1], years[idx]+1:6:years[idx+1]
end

function fill_chain_sem!(t::SharedArray, pt, qt)
    split, years = myrange_sem(t,pt)
    for (i,y) in zip(split, years)
        t[i,1,1] == 0.0 || continue # don't do previously done evals
        @show (myid(),i ,y)
        p = pt[y:y+5]
        q = qt[y:y+5,:]
        t[i,:,:], status = gchain(p,q)
        # status || continue # don't do if infeasible point
        # t[i,1], t[i,2] = RP_Cournot_test(mchain)
    end
    t
end

function fill_chain_sem_algo!(t::SharedArray, pt, qt)
    split, years = myrange_sem(t,pt)
    for (i,y) in zip(split, years)
        t[i,1,1] == 0.0 || continue # don't do previously done evals
        @show (myid(),i ,y)
        p = pt[y:y+5]
        q = qt[y:y+5,:]
        t[i,:,:], status = gchain_algo(p,q)
        # status || continue # don't do if infeasible point
        # t[i,1], t[i,2] = RP_Cournot_test(mchain)
    end
    t
end

function fill_chain_sem_alt!(t::SharedArray, pt, qt)
    split, years = myrange_sem(t,pt)
    for (i,y) in zip(split, years)
        t[i,1,1] == 0.0 || continue # don't do previously done evals
        @show (myid(),i ,y)
        p = pt[y:y+5]
        q = qt[y:y+5,:]
        t[i,:,:], status = gchain_alt(p,q)
        # status || continue # don't do if infeasible point
        # t[i,1], t[i,2] = RP_Cournot_test(mchain)
    end
    t
end

function shared_fill_chain!(t, p, q,f::Function)
    @sync begin
        for w in procs(t)
            @async remotecall_wait(f, w, t, p, q)
        end
    end
    t
end;

## parallel obj funs

function myrange(q::SharedArray)
    idx = indexpids(q)
    if idx == 0 # This worker is not assigned a piece
       return 1:0, 1:0
    end
    nchunks = length(procs(q))
    splits = [round(Int, s) for s in range(0, stop=size(q,1), length=nchunks+1)]
    splits[idx]+1:splits[idx+1]
end

function hM_chunk!(hM, mchain, gamma)
   # @show (irange, jrange, trange)  # display so we can see what's happening
   # GC.gc()
   # hM=zeros(size(mchain,2))
   T, N, nfast = size(mchain)
   geta=mchain[:,:,1]
   range = myrange(hM)
   for t in range, r=1:nfast
       val=0.0
       gtry=mchain[t,:,r]
       for j=1:N
           val=+gtry[j].*gamma[j]-geta[t,j].*gamma[j]
       end
       geta[t,:]=log.(rand()) < val[1] ? gtry : geta[t,:]
       hM[t,:]=+geta[t,:]./nfast
       # mod(r,1000)==0 && GC.gc()
       # return hM
   end
   # GC.gc()
   return hM
end

function hM_chunk_alt!(hM, mchain, gamma)
   # @show (irange, jrange, trange)  # display so we can see what's happening

   # hM=zeros(size(mchain,2))
   T, N, nfast = size(mchain)
   geta=mchain[:,:,1]
   range = myrange(hM)
   for t in range, r=1:nfast
       val=0.0
       gtry=mchain[t,:,r]
       for j=1:N
           val=+gtry[j].*gamma[j]-geta[t,j].*gamma[j]
       end
       geta[t,:]=rand() < 1.4^(val[1]) ? gtry : geta[t,:]
       hM[t,:]=+geta[t,:]./nfast
       # return hM
   end
   hM
end

# function hM_dist_chunk!(hM, mchain, gamma)
# Tried to use @distributed but it actually took more time
#    T, N, nfast = size(mchain)
#    geta=mchain[:,:,1]
#    range = myrange(hM)
#    for t in range
#        hM[t,:] = @distributed (+) for r=1:nfast
#            val=0.0
#            gtry=mchain[t,:,r]
#            for j=1:N
#                val=+gtry[j].*gamma[j]-geta[t,j].*gamma[j]
#            end
#            geta[t,:]=log.(rand()) < val[1] ? gtry : geta[t,:]
#            hM[t,:]=+geta[t,:]
#        end
#    end
#    hM./nfast
# end

function shared_hM_fill!(hM, m, g, fun::Function)
    @sync begin
        for w in procs(hM)
            @async remotecall_wait(fun, w, hM, m, g)
        end
    end
    hM
end;

##

function objfun_par(gamma0...; mchain=mchain_op)
    # @unpack nsim, nfast = Par
    gamma = gamma0
    T, N, nfast = size(mchain)

    hM = SharedArray{Float64,2}(zeros(T,N))
    shared_hM_fill!(hM, mchain, gamma, hM_chunk!)

    mhM=sum(hM, dims=1)./T
    omega=(hM.-mhM)'*(hM.-mhM)

    obj=mhM*inv(omega)*mhM'

    return obj[1]
end

function objfun_bbo_par(gamma0; mchain=mchain_op)

    gamma = gamma0
    T,N,nfast = size(mchain)

    hM = SharedArray{Float64,2}(zeros(T,N))
    shared_hM_fill!(hM, mchain, gamma, hM_chunk!)

    mhM=sum(hM, dims=1)./T
    omega=(hM.-mhM)'*(hM.-mhM)

    obj=mhM*inv(omega)*mhM'

    return obj[1]
end

function objfun_par_alt(gamma0...; mchain=mchain_op)
    # @unpack nsim, nfast = Par
    gamma = gamma0
    T, N, nfast = size(mchain)

    hM = SharedArray{Float64,2}(zeros(T,N))
    shared_hM_fill!(hM, mchain, gamma, hM_chunk_alt!)

    mhM=sum(hM, dims=1)./T
    omega=(hM.-mhM)'*(hM.-mhM)

    obj=mhM*inv(omega)*mhM'

    return obj[1]
end

function objfun_bbo_par_alt(gamma0; mchain=mchain_op)

    gamma = gamma0
    T,N,nfast = size(mchain)

    hM = SharedArray{Float64,2}(zeros(T,N))
    shared_hM_fill!(hM, mchain, gamma, hM_chunk_alt!)

    mhM=sum(hM, dims=1)./T
    omega=(hM.-mhM)'*(hM.-mhM)

    obj=mhM*inv(omega)*mhM'

    return obj[1]
end

## fill table in parallel using shared arrays
# Use mchain by semesters, test by decade

function myrange_dec(t::SharedArray)
    idx = indexpids(t)
    if idx == 0 # This worker is not assigned a piece
       return 1:0, 1:0
    end
    nchunks = length(procs(t))
    splits = collect(1:size(t,1))
    dec = [0, 16, 36, 56, 72]
    if idx ≤ 4
        splits[idx], dec[idx]+1:dec[idx+1]
    else
        1:0, 1:0
    end
end

function fill_table_dec!(t::SharedArray, m)
    i, split = myrange_dec(t)
    @show (myid(), i, split)
    mchain = m[split, :, :]
    t[i,1], t[i,2] = RP_Cournot_test_par(mchain)
    return t
end

function shared_fill_table_dec!(t, m)
    @sync begin
        for w in procs(t)
            @async remotecall_wait(fill_table_dec!, w, t, m)
        end
    end
    t
end;

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
module CausalInference
using LightGraphs
using LightGraphs.SimpleGraphs
using Combinatorics
using Base.Iterators

export dsep, skeleton, gausscitest, dseporacle, partialcor
export unshielded, pcalg, vskel
export cpdag
export digraph, vpairs, skel_oracle, pc_oracle, randdag
export cmitest, kl_entropy, kl_renyi, kl_mutual_information
export kl_cond_mi, kl_perm_mi_test, kl_perm_cond_mi_test
export n_ball
export fcialg, is_collider, is_triangle, is_parent
export is_discriminating_path, has_marks, set_marks!, is_uncovered_circle_path
export is_uncovered_PD_path, @arrow_str
export plot_pc_graph, plot_fci_graph
export orient_unshielded, orientable_unshielded, apply_pc_rules

include("graphs.jl")
include("combinations_without.jl")
include("klentropy.jl")
include("skeleton.jl")
include("dsep.jl")
include("pc.jl")
include("cpdag.jl")
include("fci.jl")
include("misc.jl")
include("recantingwitness.jl")
include("backdoor.jl")

end # module


text:
type FigureFigure{T<:Integer}
    r::Array{T,1}
    rnmax::T
    snmax::T
    snext::T
end

function grow!{T<:Integer}(ff::FigureFigure{T}, rnmax::T=100)
    ff.rnmax < rnmax || return nothing
    append!(ff.r, zeros(T, (rnmax-ff.rnmax)))
    snext = ff.snext
    for i in (ff.rnmax+1):rnmax
        ff.r[i] = ff.r[i-1] + snext
        snext += 1
        while snext in ff.r
            snext += 1
        end
    end
    ff.rnmax = rnmax
    ff.snmax = ff.r[end] - rnmax
    ff.snext = snext
    return nothing
end

function FigureFigure{T<:Integer}(rnmax::T=10)
    ff = FigureFigure([1], 1, 0, 2)
    grow!(ff, rnmax)
    return ff
end

function FigureFigure{T<:Integer}(rnmax::T, snmax::T)
    ff = FigureFigure(rnmax)
    while ff.snmax < snmax
        grow!(ff, 2ff.rnmax)
    end
    return ff
end

function make_ffr{T<:Integer}(nmax::T=10)
    ff = FigureFigure(nmax)
    function ffr{T<:Integer}(n::T)
        if n > ff.rnmax
            grow!(ff, 2n)
        end
        ff.r[n]
    end
end

function make_ffs{T<:Integer}(nmax::T=100)
    ff = FigureFigure(13, nmax)
    function ffs{T<:Integer}(n::T)
        while ff.snmax < n
            grow!(ff, 2ff.rnmax)
        end
        s = n
        for r in ff.r
            r <= s || return s
            s += 1
        end
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    QemuRunner

A `QemuRunner` represents an "execution context", an object that bundles all
necessary information to run commands within the container that contains
our crossbuild environment.  Use `run()` to actually run commands within the
`QemuRunner`, and `runshell()` as a quick way to get an interactive shell
within the crossbuild environment.
"""
mutable struct QemuRunner <: Runner
    qemu_cmd::Cmd
    append_line::String
    sandbox_cmd::Cmd
    env::Dict{String, String}
    platform::Platform

    comm_task
    exit_code::Cint
end

qemu_url = "https://github.com/Keno/QemuBuilder/releases/download/v2.12.50-0/Qemu.v2.12.50.x86_64-apple-darwin14.tar.gz"
qemu_hash = "98210f754e39008d982f447cad60043f75ed8831b3ff36d464856ac24d7c603c"
kernel_url = "https://github.com/Keno/LinuxBuilder/releases/download/v4.15-rc2/linux.x86_64-linux-gnu.tar.gz"
kernel_hash = "0d02413529e635d4af6d2122f6aba22d261f43616bb286da3855325988f9ac3b"

qemu_dependencies = [
    "https://github.com/staticfloat/PixmanBuilder/releases/download/v0.34.0-0/Pixman.x86_64-apple-darwin14.tar.gz" =>
        "7d10eba60ac7824093199d557c7fb9ec0750cd2d3f2314d5a2480677fc1f88aa",
    "https://github.com/staticfloat/GlibBuilder/releases/download/v2.54.2-2/Glib.x86_64-apple-darwin14.tar.gz" =>
        "29b52a6a95ef48a0da40e2c3a9385dd2fc93d472a965763f932a15fcd96a100b",
    "https://github.com/staticfloat/PcreBuilder/releases/download/v8.41-0/pcre.x86_64-apple-darwin14.tar.gz" =>
        "e9b99a17525ed9d694e312be7e556e0666f64515c99c7cb87257a411abf8fb9e",
    "https://github.com/staticfloat/GettextBuilder/releases/download/v0.19.8-0/libgettext.x86_64-apple-darwin14.tar.gz" =>
        "b7f64f9c34705c157f6c6a241454a2984cc80cefd11f46c63345341791cc012f",
    "https://github.com/staticfloat/LibffiBuilder/releases/download/v3.2.1-0/libffi.x86_64-apple-darwin14.tar.gz" =>
        "287f289b428cfae58084aa4bab21cb07cc405a7c5535b62a39dc8b68d163d250",
    "https://github.com/staticfloat/ZlibBuilder/releases/download/v1.2.11-3/Zlib.x86_64-apple-darwin14.tar.gz" =>
        "d921114efca229f53b4faaac3d6d153862559aa38bc6df1b702732036394c634",
]

qemu_path() = storage_dir("qemu", "bin", "qemu-system-x86_64")
kernel_path() = storage_dir("qemu", "bzImage")

"""
    update_qemu(;verbose::Bool = false)

Update our QEMU and Linux kernel installations, downloading and installing both
into `<storage_dir>/qemu`.
"""
function update_qemu(;verbose::Bool = false)
	should_delete = false
    for (url, hash) in [qemu_dependencies; qemu_url=>qemu_hash; kernel_url=>kernel_hash]
        should_delete |= !download_verify(
            url,
            hash,
            storage_dir("downloads", basename(url));
            verbose=verbose,
            force=true
        )
    end
    qemu_cache = storage_dir("qemu")
    if should_delete
        rm(qemu_cache; recursive=true, force=true)
    end
    if !isdir(qemu_cache)
        for (url, hash) in [qemu_dependencies; qemu_url=>qemu_hash]
            unpack(storage_dir("downloads", basename(url)), qemu_cache; verbose=verbose)
        end
    end
    if !isfile(kernel_path())
        unpack(storage_dir("downloads", basename(kernel_url)), qemu_cache; verbose=verbose)
    end
end

function QemuRunner(workspace_root::String;
                    cwd = nothing,
                    platform::Platform = platform_key_abi(),
                    workspaces::Vector = [],
                    extra_env=Dict{String, String}(),
                    verbose::Bool = false)
    global use_ccache

    # Choose and prepare our shards
    shards = choose_shards(platform)
    prepare_shard.(shards; mount_squashfs = false, verbose=verbose)

    # QEMU can use the .squashfs files directly, so we don't use `mount_path()`.
    qemu_mount_path(cs) = cs.archive_type == :squashfs ? download_path(cs) : mount_path(cs)

    # Download/unpack Qemu (TODO: Make this a Qemu.jll package and just rely on it.)
    update_qemu(;verbose=verbose)

    # Construct environment variables we'll use from here on out
    envs = merge(platform_envs(platform), extra_env)

    # Resource limits for QEMU.  For some reason, it seems to crash if we use too many CPUs....
    memory_limit = Int64(div(Sys.total_memory()*3, 4*1024^2))
    core_limit = min(Sys.CPU_THREADS, 3)

    qemu_cmd = ```
        $(qemu_path()) -kernel $(kernel_path())
        -m $(memory_limit)M
        -cpu host
        -smp 2
        -M accel=$(Sys.islinux() ? "kvm" : "hvf")
        -nographic
        -nodefaults
        -rtc base=utc,driftfix=slew
        -device virtio-serial
        -chardev stdio,id=charconsole0
        -device virtconsole,chardev=charconsole0,id=console0
        -device virtserialport,chardev=charcomm0,id=comm0
        -device virtio-net-pci,netdev=networking
        -netdev user,id=networking
    ```
    append_line = "quiet console=hvc0 root=/dev/vda rootflags=ro rootfstype=squashfs noinitrd init=/sandbox"

    # First shard is always the rootfs
    qemu_cmd = `$qemu_cmd -drive if=virtio,file=$(qemu_mount_path(shards[1])),format=raw`

    sandbox_cmd = ``
    if verbose
        # We put --verbose at the beginning here so that we can see argument parsing
        sandbox_cmd = `$sandbox_cmd --verbose`
    end

    if cwd != nothing
        sandbox_cmd = `$sandbox_cmd --cd $cwd`
    end

    # We always include the destdir and  workspace as one of our read-write workspace mountings, and always as the first
    insert!(workspaces, 1, workspace_root => "/workspace")

    # If we're enabling ccache, then mount in a read-writeable volume at /root/.ccache
    if use_ccache
        if !isdir(ccache_dir())
            mkpath(ccache_dir())
        end
        push!(workspaces, ccache_dir() => "/root/.ccache")
    end

    devchr = 'b'
    devnum = 0
    # Mount in read-writable workspaces over plan9
    for (outside, inside) in workspaces
        qemu_cmd = `$qemu_cmd -fsdev local,security_model=none,id=fsdev$(devnum),path=$(outside) -device virtio-9p-pci,id=fs$(devnum),fsdev=fsdev$(devnum),mount_tag=workspace$(devnum)`
        sandbox_cmd = `$sandbox_cmd --workspace 9p/workspace$(devnum):$(inside)`
        devnum += 1
    end

    # Mount in compiler shards (excluding the rootfs shard)
    for shard in shards[2:end]
        # Squashfs files get mounted in directly, actual directories have to be shared over plan9:
        if shard.archive_type == :squashfs
            qemu_cmd = `$qemu_cmd -drive if=virtio,file=$(qemu_mount_path(shard)),format=raw`
            sandbox_cmd = `$sandbox_cmd --map /dev/vd$(devchr):$(map_target(shard))`
            devchr = Char(Int(devchr) + 1)
        else
            qemu_cmd = `$qemu_cmd -fsdev local,security_model=none,id=fsdev$(devnum),path=$(qemu_mount_path(shard)),readonly -device virtio-9p-pci,id=fs$(devnum),fsdev=fsdev$(devnum),mount_tag=mapping$(devnum)`
            sandbox_cmd = `$sandbox_cmd --map 9p/mapping$(devnum):$(map_target(shard))`
            devnum += 1
        end
    end

    return QemuRunner(qemu_cmd, append_line, sandbox_cmd, envs, platform, nothing, 0)
end

function show(io::IO, x::QemuRunner)
    p = x.platform
    # Displays as, e.g., Linux x86_64 (glibc) QemuRunner
    write(io, "$(typeof(p).name.name)", " ", arch(p), " ",
          Sys.islinux(p) ? "($(p.libc)) " : "",
          "QemuRunner")
end

function qemu_gen_cmd(qr::QemuRunner, cmd::Cmd, comm_socket_path::String)
    qr.comm_task = @async begin
        # This is what we'll send to `sandbox`
        sandbox_args = String.(`$(qr.sandbox_cmd) -- $(cmd)`.exec)
        sandbox_env = ["$k=$v" for (k, v) in qr.env]

        # Wait until QEMU starts up the communications channel
        start_time = time()
        while !issocket(comm_socket_path)
            sleep(0.01)
            if time() - start_time > 2
                start_time = time()
                @warn("Unable to establish communication with QEMU on $(comm_socket_path); quitting QEMU comms task")
                qr.exit_code = -5040
                return
            end
        end

        # Once it has, open it up,
        commsock = connect(comm_socket_path)

        # We're going to spit out:
        #  [4 bytes] number of sandbox args
        #    for each sandbox arg:
        #      [4 bytes] length of arg
        #      [N bytes] arg literal string
        #  [4 bytes] number of environment mappings
        #    for each environment mapping
        #      [4 bytes] length of mapping
        #      [N bytes] mapping literal string
        write(commsock, UInt32(length(sandbox_args)))
        for idx in 1:length(sandbox_args)
            write(commsock, UInt32(length(sandbox_args[idx])))
            write(commsock, sandbox_args[idx])
        end

        write(commsock, UInt32(length(sandbox_env)))
        for idx in 1:length(sandbox_env)
            write(commsock, UInt32(length(sandbox_env[idx])))
            write(commsock, sandbox_env[idx])
        end

        # Wait for acknowledgement
        ack = read(commsock, UInt8)
        if ack != 0
            @warn("Acknowledgement from sandbox nonzero! $(ack)")
        end

        # Read exit status
        qr.exit_code = read(commsock, Cint)

        # Close communications socket
        close(commsock)
    end

    long_cmd = ```
        $(qr.qemu_cmd)
        -device virtio-serial -chardev socket,path=$(comm_socket_path),server,nowait,id=charcomm0
        -append $(qr.append_line)
    ```

    return `$long_cmd`
end

function qemu_wait(qr, oc)
    # Wait for QEMU to finish running
    did_succeed = wait(oc)

    if !did_succeed
        @warn("QEMU failed to run")
        return false
    end

    # Wait for our communications task to finish
    wait(qr.comm_task)

    # Check return code
    return qr.exit_code == 0
end

function Base.run(qr::QemuRunner, cmd, logpath::AbstractString; verbose::Bool = false, tee_stream=stdout)
    return temp_prefix() do prefix
        comm_socket_path = joinpath(prefix.path, "qemu_comm.socket")
        # Launch QEMU
        cmd = qemu_gen_cmd(qr, cmd, comm_socket_path)

        oc = OutputCollector(cmd; verbose=verbose, tee_stream=tee_stream)
        did_succeed = qemu_wait(qr, oc)

        if !isempty(logpath)
            # Write out the logfile, regardless of whether it was successful
            mkpath(dirname(logpath))
            open(logpath, "w") do f
                # First write out the actual command, then the command output
                println(f, cmd)
                print(f, merge(oc))
            end
        end

        # Return whether we succeeded or not
        return did_succeed
    end
end

function run_interactive(qr::QemuRunner, cmd::Cmd; stdin = nothing, stdout = nothing, stderr = nothing)
    temp_prefix() do prefix
        comm_socket_path = joinpath(prefix.path, "qemu_comm.socket")
        cmd = qemu_gen_cmd(qr, cmd, comm_socket_path)
        if stdin isa Base.AbstractCmd || stdin isa Base.TTY
            cmd = pipeline(cmd, stdin=stdin)
        end
        if stdout isa Base.AbstractCmd || stdout isa Base.TTY
            cmd = pipeline(cmd, stdout=stdout)
        end
        if stderr isa Base.AbstractCmd || stderr isa Base.TTY
            cmd = pipeline(cmd, stderr=stderr)
        end

        if stdout isa IOBuffer
            if !(stdin isa IOBuffer)
                stdin = devnull
            end

            process = open(cmd, "r", stdin)
            @async begin
                while !eof(process)
                    write(stdout, read(process))
                end
            end
            wait(process)
        else
            @static if Sys.isapple()
                println("Setting parent shell interrupt to ^] (use ^C as usual, ^] to interrupt Julia)")
                run(`stty intr ^\]`)
            end
            run(cmd)
            @static if Sys.isapple()
                run(`stty intr ^\C`)
            end
        end
    end

    # Qemu appears to mess with our terminal
    if stdin !== nothing
        Base.reseteof(stdin)
    end
end

function runshell(qr::QemuRunner, args...; kwargs...)
    run_interactive(qr, `/bin/bash`, args...; kwargs...)
end


text:
using LinearAlgebra
using NewRimuOperators
using Rimu
using Test

using NewRimuOperators:
    ParticleCountTerm, OnsiteInteractionTerm, NeighbourOneBodyTerm, FullOneBodyTerm,
    MomentumTwoBodyTerm, FullTwoBodyTerm, MomentumThreeBodyTerm

@testset "terms" begin
    function one_component_changed_or_zero(term, add)
        # Only one address changes at a time
        for (new_add, val) in offdiagonals(term, add)
            if new_add.components[1] ≠ add.components[1]
                new_add.components[2] == add.components[2] || return false
            elseif new_add.components[2] ≠ add.components[2]
                new_add.components[1] == add.components[1] || return false
            else
                val == 0 || return false
            end
        end
        return true
    end

    function good_structure(term, add)
        if LOStructure(term) == IsHermitian() || LOStructure(term) == IsDiagonal()
            return ishermitian(Matrix(term, add))
        elseif LOStructure(term) == AdjointKnown()
            return Matrix(term, add)' == Matrix(term', add)
        end
    end

    @testset "ParticleCountTerm" begin
        add = CompositeFS(
            BoseFS((0, 1, 2, 0)),
            FermiFS((1, 1, 0, 0)),
        )
        term = ParticleCountTerm((σ, p) -> σ * p)
        @test length(offdiagonals(term, add)) == 0
        @test diagonal_element(term, add) == 2 + 6 + 2 * (1 + 2)
        @test good_structure(term, add)
    end

    @testset "OnsiteInteractionTerm" begin
        add = CompositeFS(
            BoseFS((0, 1, 3, 0)),
            FermiFS((1, 1, 1, 0)),
        )
        term = OnsiteInteractionTerm((σ, τ) -> σ == τ ? 10 : 2)
        @test length(offdiagonals(term, add)) == 0
        @test diagonal_element(term, add) == 30 + 2 * (1 + 3)
        @test good_structure(term, add)
    end

    @testset "NeighbourOneBodyTerm" begin
        add = CompositeFS(
            BoseFS((0, 0, 3, 0)),
            FermiFS((1, 1, 1, 0)),
        )
        term = NeighbourOneBodyTerm(identity)
        @test length(offdiagonals(term, add)) == 2 + 6
        @test one_component_changed_or_zero(term, add)
        @test diagonal_element(term, add) == 0
        @test good_structure(term, add)
    end

    @testset "FullOneBodyTerm" begin
        add = CompositeFS(
            BoseFS((0, 0, 3, 0)),
            FermiFS((1, 1, 0, 0)),
        )
        term = FullOneBodyTerm((σ, q, p) -> σ == 1 ? (p - q * im) : (p * q + im))
        @test num_offdiagonals(term, add) == 3 * 3
        @test diagonal_element(term, add) == 14 - 7im
        @test one_component_changed_or_zero(term, add)
        @test good_structure(term, add)
    end

    @testset "MomentumTwoBodyTerm" begin
        add = CompositeFS(
            BoseFS((0, 0, 3, 0)),
            FermiFS((1, 1, 0, 0)),
        )
        term = MomentumTwoBodyTerm((σ,τ,p,q,k) -> σ * τ * mod1(k,num_modes(add)) + (p - q))
        # Broken with bosons due to how momentum_transfer_excitation is implemented.
        @test_broken good_structure(term, add)
    end

    @testset "FullTwoBodyTerm" begin
        add = CompositeFS(
            FermiFS((0, 0, 1, 0, 1, 0)),
            FermiFS((1, 1, 0, 0, 0, 0)),
        )
        term = FullTwoBodyTerm((σ,τ,p,q,r,s) -> σ * τ * (p + q + r + s))
        @test_broken good_structure(term, add)
    end

    @testset "MomentumThreeBodyTerm" begin
        add = CompositeFS(
            FermiFS((0, 0, 1, 0, 1, 0)),
            FermiFS((1, 1, 0, 0, 0, 0)),
        )
        term = MomentumThreeBodyTerm((σ,τ,k,l) -> σ * τ * (k + l))
    end

end


text:
struct Proposal
    W::Int
    S::Float64
end

rand_proposal() = Proposal(rand(5:15), rand(Uniform(0, .1)))

function crowdingDistance(arr::Array{Float64,2})
    rs, cs = size(arr)
    c_sum = zeros(rs)
    for k ∈ 1:cs
        a = arr[:, k]
        cd = Vector{Float64}(undef, rs)
        cd[1] = cd[end] = Inf
        sc = sort(a)
        ind = [ sort(collect(enumerate(a)), by = x -> x[2])[i][1] for i in 1:rs ]
        dist = euclidean(extrema(a)...)
        @inbounds for i in 2:rs - 1
            cd[i] = (sc[i + 1] - sc[i - 1])/dist
        end
        c_sum += cd[ind, :]
    end
    return c_sum
end

function dominates(x, y)
    strict_inequality_found = false
    for i ∈ eachindex(x)
        y[i] < x[i] && return false
        strict_inequality_found |= x[i] < y[i]
    end
    return strict_inequality_found
end

function fast_nds(arr)
    sz = size(arr, 2)
    fronts = Vector{Int64}[]
    ind = collect(axes(arr, 1))
    a = SVector{sz}.(eachrow(arr))
    N = length(a)
    ii = 1
    while ii <= N
        start = ii
        @inbounds for i in ii:N
            for j in start:N
                dominates(a[j], a[i]) && @goto l1
            end
            ind[ii], ind[i] = ind[i], ind[ii]
            a[ii], a[i] = a[i], a[ii]
            ii += 1
            @label l1
        end
        push!(fronts, ind[start:ii - 1])
    end
    return fronts
end

function mix(p::Proposal, q::Proposal)
    w = p.W < q.W ? rand(p.W:q.W) : rand(q.W:p.W)
    s = p.S < q.S ? rand(Uniform(p.S, q.S)) : rand(Uniform(q.S, p.S))
    return Proposal(w, s)
end

function crossover(v::Array{Proposal,1})
    parent_pop_size = length(v)
    child_arr = Array{Proposal,1}(undef, parent_pop_size)
    @inbounds for i in 1:parent_pop_size
        x, y = sample(1:parent_pop_size, 2, replace = false)
	child_arr[i] = mix(v[x], v[y])
    end
    return child_arr
end

function new_gen(v::Array{Proposal,1}; prob = 0.05)
    pop = vcat(v, crossover(v))
    s = sample(1:length(pop), floor(Int, 100*prob), replace = false)
    @inbounds for i ∈ s
	pop[i] = rand_proposal()
    end
    return pop
end

function intervene(ps::ParameterSetting, p::Proposal, numb_updates::Int; init_w = .5)
    @unpack max_tspan, pr_quick_rec, trans_rate, drop, pr_death, N, λ, hhs = ps
    ppl, grph = create_model(N, λ, hhs)
    window, scale = p.W, p.S
    daily_inf = Int[]
    weights = Float64[]
    cfs = Float64[]

    @inbounds for i in 1:window
        model_update!(ppl, grph, init_w, max_tspan, pr_quick_rec, trans_rate, drop, pr_death)
        push!(daily_inf, sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :I)))
        push!(weights, init_w)
    end

    pq = .0
    @inbounds for i in 1:numb_updates - window
        x = [ones(window) 1:window]
        y = daily_inf[end - (window - 1):end]
        _, cf = x \ y
        pq = pq + (cf * scale)
        w = rich_scaled(pq)
        model_update!(ppl, grph, w, max_tspan, pr_quick_rec, trans_rate, drop, pr_death)
        push!(daily_inf, sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :I)))
        push!(weights, w)
    end

    mw = mean(weights)
    undercap = sum(max.((daily_inf ./ 100) .- (length(ppl) / 1000), 0))
    tot_dec = sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :D))

    return tot_dec, mw, undercap
end

function sim_fnc_nsga(ps::ParameterSetting, p::Proposal, numb_updates::Int; n_test=5)
    dec = 0
    uc = 0.
    wgts = 0.

    i = 1
    while i <= n_test
        res = intervene(ps, p, numb_updates)
        if first(res) < 4 # epidemic didn't get off the ground
            i -= 1
        else
            dec += first(res)
            uc += last(res)
            wgts += res[2]
        end
        i += 1
    end
    return dec / n_test, uc / n_test, wgts / n_test
end

function select_best(ps::ParameterSetting, arr::Array{Proposal,1}, numb_updates::Int, parent_pop_size::Int)
    scores = Matrix(DataFrame(pmap(i -> sim_fnc_nsga(ps, arr[i], numb_updates), 1:length(arr))))
    cdist = crowdingDistance(scores)
    ranking = fast_nds(scores)
    rnkInd = Vector{Int64}(undef, 2*parent_pop_size)
    @inbounds for i in 1:length(ranking)
        rnkInd[ranking[i]] .= i
    end
    ns = hcat(rnkInd, cdist[:], 1:2*parent_pop_size)
    p2 = sortperm(view(ns, :, 2))
    s1 = ns[p2, :]
    p1 = sortperm(view(s1, :, 1))
    out = Int.(s1[p1, 3][1:parent_pop_size])
    return arr[out], scores[out, :]
end

extract_pams(p::Proposal) = (p.W, p.S)

function run_evo(ps::ParameterSetting, numb_updates::Int, numb_gen::Int, parent_pop_size::Int; full=false)
    pms_ar = Array{Tuple{Int64,Float64},1}[]
    start_pop = [ rand_proposal() for _ in 1:parent_pop_size ]
    push!(pms_ar, extract_pams.(start_pop))
    scrs = Array{Float64,3}(undef, parent_pop_size, 3, numb_gen)
    for i in 1:numb_gen
        new_pop = [ Proposal(pms_ar[i][j]...) for j in 1:parent_pop_size ]
        out_pop, out_scrs = select_best(ps, new_gen(new_pop), numb_updates, parent_pop_size)
        push!(pms_ar, extract_pams.(out_pop))
        scrs[:, :, i] = out_scrs
    end

    pf = fast_nds(scrs[:, :, end])
    pareto_front = pms_ar[end][pf[1]]

    if full
        return pms_ar, scrs, pareto_front
    else
        return pareto_front
    end
end


text:
__precompile__(false)

# ls -l ~/.julia/v0.6  | grep PCL | awk '{print $9}' | sort | grep -v "^PCL$"|\
#   awk '{print "-", "["$1"]""(https://github.com/JuliaPCL/"$1".jl)"}'
"""
The Julia wrapper for Point Cloud Library (PCL)

With the packages, e.g., you can visualize a point cloud in 5 lines:

```julia
using PCLCommon, PCLIO, PCLVisualization
cloud = PointCloud{PointXYZRGB}("your_pcd_file.pcd")
viewer = PCLVisualizer("pcl visualizer")
addPointCloud(viewer, cloud, id="input")
spin(viewer) # you will see the PCLVisualizer
```

![](assets/pcl_visualizer_5lines.png)

## Why?

PCL.jl was stated to give us **fast prototyping** using PCL in interactive
environment, without much loss of computational efficiency and flexibility.

You might find that there is already a python binding
([strawlab/python-pcl](https://github.com/strawlab/python-pcl)) for the purpose,
however, it lacks flexibility in particular; it only supports `PointXYZ` for
point clouds at the moment. I guess the reason why it doesn't support arbitary
point types for point clouds is that there is not elegant and simple way to
expose C++ templated classes in python
([actually I tried](https://github.com/r9y9/pypcl), but couldn't figure it out).
Since I wanted to use arbiraty point types for point clouds, I decided to create
a new one.

The reasons why I started writing the [Julia](http://julialang.org/) (not python)
binding are:

- **The Julia C++ interface [Keno/Cxx.jl](https://github.com/Keno/Cxx.jl) is quite powerful**: It enables us to call C++ functions/methods directly from Julia *without any glue code* (unlike cython and swig), embed C++ into Julia and vise versa.
- **Julia's types can be parameterized and exposing C++ templated classes is quite straightfood**: e.g. C++ `pcl::PointCloud<PointT>` can be represented as `PointCloud{PointT}` in Julia.

### Comparison to python-pcl

Pros:

- Support *arbitaty* point types for point clouds, whereas python-pcl only supports `PointXYZ`
- Support PCL 1.8 or later
- More consistent APIs
- Can write mixed C++/Julia
- PCLVisualizer: Jupyter integration using off-screen rendering

Cons:

- Only works on osx for now
- Hard to build entire dependencies
- Sometime segfaults in Cxx.jl and Julia when doing wrong

### Dependencies (in short)

- [pcl](https://github.com/PointCloudLibrary/pcl)
- [Julia](https://github.com/JuliaLang/julia)
- [Cxx.jl](https://github.com/Keno/Cxx.jl)

## API guidelines

- Function names should be exactly same between Julia and C++.
- C++ template classes are available in Julia as templated types
- C++ dereferences which sometimes needed in C++, are hidden in implementation in Julia

e.g.

In C++:

```cpp
pcl::PassThrough<pcl::PointXYZ> pass;
pass.setInputCloud (cloud);
pass.setFilterFieldName ("z");
pass.setFilterLimits (0.0, 1.0);
pass.filter (*cloud_filtered);
```

In Julia:

```julia
pass = PassThrough{PointXYZ}()
setInputCloud(pass, cloud)
setFilterFieldName(pass, "z")
setFilterLimits(pass, 0.0, 1.0)
filter(pass, cloud_filtered)
```

## Package structure

To simplify development and minimize dependencies, the Julia wrapper consists of
the packages below:

- [LibPCL.jl](@ref)
- [PCLCommon.jl](@ref)
- [PCLFeatures.jl](@ref)
- [PCLFilters.jl](@ref)
- [PCLIO.jl](@ref)
- [PCLKDTree.jl](@ref)
- [PCLKeyPoints.jl](@ref)
- [PCLOctree.jl](@ref)
- [PCLRecognition.jl](@ref)
- [PCLRegistration.jl](@ref)
- [PCLSampleConsensus.jl](@ref)
- [PCLSearch.jl](@ref)
- [PCLSegmentation.jl](@ref)
- [PCLSurface.jl](@ref)
- [PCLTracking.jl](@ref)
- [PCLVisualization.jl](@ref)

following the [PCL module structure](http://docs.pointclouds.org/trunk/modules.html)
except for the `LibPCL.jl` which manages binary dependencies (i.e. search
installed PCL shared libraries or build and install them if not found).
"""
module PCL

using Reexport
using DocStringExtensions

@reexport using PCLCommon
@reexport using PCLIO
@reexport using PCLFeatures
@reexport using PCLFilters
@reexport using PCLVisualization
@reexport using PCLKDTree
@reexport using PCLOctree
@reexport using PCLRecognition
@reexport using PCLSampleConsensus
@reexport using PCLSearch
@reexport using PCLSegmentation
@reexport using PCLTracking
@reexport using PCLRegistration
@reexport using PCLKeyPoints
@reexport using PCLSurface
@reexport using LibPCL

end # module PCL