text:
#=
input:
- Julia version:
- Author: quangio
- Date: 2019-05-09
=#

module TestSet
    using DelimitedFiles
    dist_mat_48 = [
             0      4727      1205      6363      3657      3130      2414       563       463      5654      1713      1604      2368      2201      1290      1004      3833      2258      3419      2267      2957       720      1700      5279      2578      6076      3465      2654      3625      3115      1574      3951      1748      2142      6755      2383      3306      1029      3530       825      2188      4820      3489      1947      6835      1542      2379      3744
      4727         0      3588      2012      1842      6977      6501      5187      5028      2327      4148      4723      3635      3125      4907      3930      7463      6338      7243      5105      4043      4022      3677      2863      3106      1850      7173      6630      1204      6814      6001      3447      5253      2656      3123      6274      7183      5622      3085      4564      2756      1591      7027      6186      3472      5461      4390      2088
      1205      3588         0      5163      2458      3678      3071      1742      1444      4462      1184      1520      1498      1103      1501       951      4298      2903      3967      2169      2209       652       828      4136      1518      4873      3954      3254      2446      3581      2441      2960      1966       950      5564      2916      3878      2035      2482      1027      1395      3617      3891      2686      5661      2023      1867      2560
      6363      2012      5163         0      2799      8064      7727      6878      6581      1402      5366      5946      4679      4378      6225      5709      8417      7578      8296      6135      4802      5707      4982      2322      4178       320      8186      7800      2778      7859      7408      3763      6461      4223      1427      7451      8263      7131      3669      6011      4638      1681      7987      7502      1877      6758      5360      2844
      3657      1842      2458      2799         0      5330      4946      4200      3824      2012      2573      3157      1924      1580      3427      3179      5749      4793      5577      3409      2223      3066      2185      1860      1401      2491      5486      5035       894      5141      4611      1669      3677      1590      3113      4682      5533      4352      1252      3227      2426      1169      5313      4706      3241      3962      2651       304
      3130      6977      3678      8064      5330         0       743      3209      2670      6929      2831      2266      3407      3854      2178      4076       727       881       293      1930      3310      3672      3315      6199      3932      7745       365       482      5774       261      1659      4513      1746      4431      7910       769       207      2225      4435      2681      5053      6384       550      1224      7805      1670      2704      5230
      2414      6501      3071      7727      4946       743         0      2468      1952      6673      2380      1795      3051      3405      1604      3382      1469       168      1020      1681      3110      2993      2827      6009      3552      7412      1104       267      5300       821       916      4348      1270      3890      7698       332       900      1484      4185      2049      4415      6051      1219       482      7635      1054      2432      4884
       563      5187      1742      6878      4200      3209      2468         0       718      6203      2241      2051      2920      2762      1687      1304      3932      2331      3487      2669      3487      1175      2260      5840      3141      6596      3563      2728      4120      3240      1559      4507      2082      2658      7304      2512      3364       985      4091      1319      2544      5358      3632      1987      7391      1785      2879      4296
       463      5028      1444      6581      3824      2670      1952       718         0      5789      1602      1343      2330      2291       970      1451      3376      1796      2959      1951      2835      1112      1725      5346      2628      6285      3007      2193      3889      2661      1122      3920      1372      2391      6883      1927      2845       611      3543       676      2590      4993      3039      1486      6934      1112      2196      3876
      5654      2327      4462      1402      2012      6929      6673      6203      5789         0      4392      4947      3648      3501      5274      5183      7216      6535      7140      5022      3621      5077      4090       922      3207      1131      7014      6714      2437      6707      6477      2476      5432      3599      1102      6376      7121      6284      2497      5160      4318       937      6795      6507      1268      5773      4249      1914
      1713      4148      1184      5366      2573      2831      2380      2241      1602      4392         0       586       766      1029       883      2040      3353      2224      3100      1049      1246      1625       503      3841      1196      5054      3042      2488      2945      2676      2087      2331      1114      1650      5459      2132      3037      1958      1997       931      2513      3701      2923      2137      5459      1394       711      2534
      1604      4723      1520      5946      3157      2266      1795      2051      1343      4947       586         0      1299      1612       406      2208      2824      1639      2542       694      1586      1767      1050      4357      1770      5633      2498      1907      3520      2128      1558      2778       531      2171      6003      1552      2472      1538      2506       791      2912      4277      2403      1564      5983       827       892      3109
      2368      3635      1498      4679      1924      3407      3051      2920      2330      3648       766      1299         0       646      1642      2446      3840      2905      3655      1488       730      2096       697      3076       533      4363      3567      3122      2453      3219      2842      1592      1791      1480      4706      2772      3610      2721      1232      1656      2550      3001      3403      2860      4697      2126       756      1836
      2201      3125      1103      4378      1580      3854      3405      2762      2291      3501      1029      1612       646         0      1853      2026      4349      3247      4119      1997      1341      1753       606      3078       419      4070      4052      3517      1923      3690      3032      1866      2142       838      4593      3161      4060      2788      1380      1663      1932      2736      3915      3138      4647      2395      1351      1592
      1290      4907      1501      6225      3427      2178      1604      1687       970      5274       883       406      1642      1853         0      2029      2803      1438      2466       986      1987      1593      1253      4716      2072      5915      2454      1764      3710      2082      1204      3164       497      2287      6342      1419      2379      1134      2867       554      2885      4569      2405      1289      6338       555      1297      3406
      1004      3930       951      5709      3179      4076      3382      1304      1451      5183      2040      2208      2446      2026      2029         0      4759      3220      4368      2900      3151       442      1765      4960      2444      5443      4396      3610      2932      4034      2572      3891      2525      1590      6278      3313      4261      2033      3398      1476      1241      4287      4390      2928      6419      2428      2749      3337
      3833      7463      4298      8417      5749       727      1469      3932      3376      7216      3353      2824      3840      4349      2803      4759         0      1601       477      2359      3617      4345      3851      6433      4372      8098       370      1206      6267       726      2384      4754      2335      4991      8148      1452       609      2949      4752      3331      5687      6746       437      1948      8005      2334      3098      5618
      2258      6338      2903      7578      4793       881       168      2331      1796      6535      2224      1639      2905      3247      1438      3220      1601         0      1165      1563      2988      2829      2666      5882      3401      7263      1233       399      5138       923       794      4227      1117      3724      7565       286      1049      1348      4051      1881      4248      5903      1322       355      7508       887      2302      4736
      3419      7243      3967      8296      5577       293      1020      3487      2959      7140      3100      2542      3655      4119      2466      4368       477      1165         0      2170      3520      3965      3588      6393      4183      7977       202       767      6041       438      1932      4706      2027      4711      8107      1061       132      2503      4652      2972      5344      6617       486      1501      7989      1962      2939      5469
      2267      5105      2169      6135      3409      1930      1681      2669      1951      5022      1049       694      1488      1997       986      2900      2359      1563      2170         0      1430      2460      1547      4333      2019      5817      2079      1694      3910      1733      1813      2668       654      2694      6029      1366      2130      1991      2525      1474      3542      4455      1923      1641      5957      1071       777      3302
      2957      4043      2209      4802      2223      3310      3110      3487      2835      3621      1246      1586       730      1341      1987      3151      3617      2988      3520      1430         0      2779      1387      2905      1062      4482      3398      3119      2922      3087      3115      1240      1953      2175      4607      2796      3501      3119      1136      2173      3268      3136      3189      3029      4527      2355       711      2042
       720      4022       652      5707      3066      3672      2993      1175      1112      5077      1625      1767      2096      1753      1593       442      4345      2829      3965      2460      2779         0      1401      4781      2166      5427      3984      3212      2946      3620      2224      3603      2089      1496      6178      2906      3861      1719      3132      1040      1479      4211      3969      2553      6290      2012      2336      3189
      1700      3677       828      4982      2185      3315      2827      2260      1725      4090       503      1050       697       606      1253      1765      3851      2666      3588      1547      1387      1401         0      3621       903      4675      3537      2954      2475      3169      2427      2254      1578      1148      5177      2598      3521      2194      1833      1074      2054      3340      3423      2541      5213      1801      1077      2190
      5279      2863      4136      2322      1860      6199      6009      5840      5346       922      3841      4357      3076      3078      4716      4960      6433      5882      6393      4333      2905      4781      3621         0      2718      2042      6254      6024      2569      5966      5913      1687      4807      3384      1716      5699      6384      5787      1852      4687      4285      1272      6022      5892      1629      5178      3581      1639
      2578      3106      1518      4178      1401      3932      3552      3141      2628      3207      1196      1770       533       419      2072      2444      4372      3401      4183      2019      1062      2166       903      2718         0      3864      4097      3635      1932      3748      3274      1448      2284      1164      4286      3283      4136      3086       967      1973      2285      2507      3935      3331      4312      2589      1284      1340
      6076      1850      4873       320      2491      7745      7412      6596      6285      1131      5054      5633      4363      4070      5915      5443      8098      7263      7977      5817      4482      5427      4675      2042      3864         0      7866      7483      2515      7539      7101      3449      6146      3938      1375      7134      7944      6831      3349      5709      4397      1363      7667      7190      1798      6446      5041      2528
      3465      7173      3954      8186      5486       365      1104      3563      3007      7014      3042      2498      3567      4052      2454      4396       370      1233       202      2079      3398      3984      3537      6254      4097      7866         0       839      5973       374      2019      4569      1996      4669      7970      1085       305      2581      4532      2976      5339      6509       287      1581      7844      1974      2838      5369
      2654      6630      3254      7800      5035       482       267      2728      2193      6714      2488      1907      3122      3517      1764      3610      1206       399       767      1694      3119      3212      2954      6024      3635      7483       839         0      5427       558      1181      4349      1377      4044      7723       356       653      1744      4218      2241      4614      6121       955       743      7644      1231      2465      4957
      3625      1204      2446      2778       894      5774      5300      4120      3889      2437      2945      3520      2453      1923      3710      2932      6267      5138      6041      3910      2922      2946      2475      2569      1932      2515      5973      5427         0      5612      4824      2550      4050      1498      3476      5071      5980      4470      2096      3388      1911      1501      5831      4994      3704      4264      3209      1196
      3115      6814      3581      7859      5141       261       821      3240      2661      6707      2676      2128      3219      3690      2082      4034       726       923       438      1733      3087      3620      3169      5966      3748      7539       374       558      5612         0      1716      4280      1624      4298      7679       735       420      2263      4216      2606      4967      6179       400      1277      7567      1609      2501      5032
      1574      6001      2441      7408      4611      1659       916      1559      1122      6477      2087      1558      2842      3032      1204      2572      2384       794      1932      1813      3115      2224      2427      5913      3274      7101      2019      1181      4824      1716         0      4330      1180      3346      7545      1023      1808       578      4062      1438      3693      5763      2115       440      7537       763      2404      4603
      3951      3447      2960      3763      1669      4513      4348      4507      3920      2476      2331      2778      1592      1866      3164      3891      4754      4227      4706      2668      1240      3603      2254      1687      1448      3449      4569      4349      2550      4280      4330         0      3184      2510      3402      4031      4698      4281       533      3245      3612      2187      4339      4265      3296      3576      1941      1381
      1748      5253      1966      6461      3677      1746      1270      2082      1372      5432      1114       531      1791      2142       497      2525      2335      1117      2027       654      1953      2089      1578      4807      2284      6146      1996      1377      4050      1624      1180      3184         0      2685      6475      1022      1952      1341      2963      1050      3358      4787      1926      1086      6436       422      1244      3619
      2142      2656       950      4223      1590      4431      3890      2658      2391      3599      1650      2171      1480       838      2287      1590      4991      3724      4711      2694      2175      1496      1148      3384      1164      3938      4669      4044      1498      4298      3346      2510      2685         0      4697      3693      4636      2975      1981      1909      1124      2718      4565      3548      4830      2839      2140      1751
      6755      3123      5564      1427      3113      7910      7698      7304      6883      1102      5459      6003      4706      4593      6342      6278      8148      7565      8107      6029      4607      6178      5177      1716      4286      1375      7970      7723      3476      7679      7545      3402      6475      4697         0      7393      8097      7370      3515      6249      5379      2001      7738      7556       461      6829      5267      3013
      2383      6274      2916      7451      4682       769       332      2512      1927      6376      2132      1552      2772      3161      1419      3313      1452       286      1061      1366      2796      2906      2598      5699      3283      7134      1085       356      5071       735      1023      4031      1022      3693      7393         0       965      1542      3883      1913      4286      5772      1121       600      7322       902      2128      4608
      3306      7183      3878      8263      5533       207       900      3364      2845      7121      3037      2472      3610      4060      2379      4261       609      1049       132      2130      3501      3861      3521      6384      4136      7944       305       653      5980       420      1808      4698      1952      4636      8097       965         0      2380      4629      2877      5250      6583       570      1380      7986      1866      2904      5432
      1029      5622      2035      7131      4352      2225      1484       985       611      6284      1958      1538      2721      2788      1134      2033      2949      1348      2503      1991      3119      1719      2194      5787      3086      6831      2581      1744      4470      2263       578      4281      1341      2975      7370      1542      2380         0      3952      1127      3197      5518      2658      1002      7395       951      2429      4380
      3530      3085      2482      3669      1252      4435      4185      4091      3543      2497      1997      2506      1232      1380      2867      3398      4752      4051      4652      2525      1136      3132      1833      1852       967      3349      4532      4218      2096      4216      4062       533      2963      1981      3515      3883      4629      3952         0      2873      3080      2012      4324      4046      3478      3328      1755      1000
       825      4564      1027      6011      3227      2681      2049      1319       676      5160       931       791      1656      1663       554      1476      3331      1881      2972      1474      2173      1040      1074      4687      1973      5709      2976      2241      3388      2606      1438      3245      1050      1909      6249      1913      2877      1127      2873         0      2374      4392      2943      1659      6285      1012      1563      3254
      2188      2756      1395      4638      2426      5053      4415      2544      2590      4318      2513      2912      2550      1932      2885      1241      5687      4248      5344      3542      3268      1479      2054      4285      2285      4397      5339      4614      1911      4967      3693      3612      3358      1124      5379      4286      5250      3197      3080      2374         0      3386      5284      3997      5585      3386      3125      2664
      4820      1591      3617      1681      1169      6384      6051      5358      4993       937      3701      4277      3001      2736      4569      4287      6746      5903      6617      4455      3136      4211      3340      1272      2507      1363      6509      6121      1501      6179      5763      2187      4787      2718      2001      5772      6583      5518      2012      4392      3386         0      6314      5837      2205      5095      3680      1169
      3489      7027      3891      7987      5313       550      1219      3632      3039      6795      2923      2403      3403      3915      2405      4390       437      1322       486      1923      3189      3969      3423      6022      3935      7667       287       955      5831       400      2115      4339      1926      4565      7738      1121       570      2658      4324      2943      5284      6314         0      1676      7603      1964      2662      5184
      1947      6186      2686      7502      4706      1224       482      1987      1486      6507      2137      1564      2860      3138      1289      2928      1948       355      1501      1641      3029      2553      2541      5892      3331      7190      1581       743      4994      1277       440      4265      1086      3548      7556       600      1380      1002      4046      1659      3997      5837      1676         0      7521       744      2325      4670
      6835      3472      5661      1877      3241      7805      7635      7391      6934      1268      5459      5983      4697      4647      6338      6419      8005      7508      7989      5957      4527      6290      5213      1629      4312      1798      7844      7644      3704      7567      7537      3296      6436      4830       461      7322      7986      7395      3478      6285      5585      2205      7603      7521         0      6805      5208      3102
      1542      5461      2023      6758      3962      1670      1054      1785      1112      5773      1394       827      2126      2395       555      2428      2334       887      1962      1071      2355      2012      1801      5178      2589      6446      1974      1231      4264      1609       763      3576       422      2839      6829       902      1866       951      3328      1012      3386      5095      1964       744      6805         0      1644      3928
      2379      4390      1867      5360      2651      2704      2432      2879      2196      4249       711       892       756      1351      1297      2749      3098      2302      2939       777       711      2336      1077      3581      1284      5041      2838      2465      3209      2501      2404      1941      1244      2140      5267      2128      2904      2429      1755      1563      3125      3680      2662      2325      5208      1644         0      2532
      3744      2088      2560      2844       304      5230      4884      4296      3876      1914      2534      3109      1836      1592      3406      3337      5618      4736      5469      3302      2042      3189      2190      1639      1340      2528      5369      4957      1196      5032      4603      1381      3619      1751      3013      4608      5432      4380      1000      3254      2664      1169      5184      4670      3102      3928      2532         0
    ]

    loadVrpEil22() = readdlm("/home/quangio/academic/opt/data/abc")
    loadVrpEil13() = readdlm("/home/quangio/academic/opt/data/eil13")
    eil22_capacities = [
    0
 1100
  700
  800
 1400
 2100
  400
  800
  100
  500
  600
 1200
 1300
 1300
  300
  900
 2100
 1000
  900
 2500
 1800
  700
  ]
end

max_stars_repo_name: quangIO/mTSP

max_issues_repo_name: quangIO/mTSP

max_forks_repo_name: quangIO/mTSP


text:
lv(x) = GlobalRef(LoopVectorization, x)

"""
    Instruction

`Instruction` represents a function via its module and symbol. It is
similar to a `GlobalRef` and may someday be replaced by `GlobalRef`.
"""
struct Instruction
    mod::Symbol
    instr::Symbol
end
# lower(instr::Instruction) = Expr(:(.), instr.mod, QuoteNode(instr.instr))
# Base.convert(::Type{Expr}, instr::Instruction) = Expr(:(.), instr.mod, QuoteNode(instr.instr))
function callexpr(instr::Instruction)
    if instr.mod === :LoopVectorization
        Expr(:call, lv(instr.instr))
    else#if instr.mod === :Main
        Expr(:call, instr.instr)
    end
end
function callexpr(instr::Instruction, arg)
    ce = callexpr(instr)
    append!(ce.args, arg)
    ce
end
Base.hash(instr::Instruction, h::UInt64) = hash(instr.instr, hash(instr.mod, h))
# function Base.isless(instr1::Instruction, instr2::Instruction)
#     if instr1.mod === instr2.mod
#         isless(instr1.instr, instr2.instr)
#     else
#         isless(instr1.mod, instr2.mod)
#     end
# end
Base.isequal(ins1::Instruction, ins2::Instruction) = (ins1.instr === ins2.instr) && (ins1.mod === ins2.mod)

"""
    InstructionCost

Store parameters related to performance for individual CPU instructions.

$(TYPEDFIELDS)
"""
struct InstructionCost
    "A flag indicating how instruction cost scales with vector width (128, 256, or 512 bits)"
    scaling::Float64 # sentinel values: -3 == no scaling; -2 == offset_scaling, -1 == linear scaling, >0 ->  == latency == reciprocal throughput
    """The number of clock cycles per operation when many of the same operation are repeated in sequence.
    Think of it as the inverse of the flow rate at steady-state. It is typically ≤ the `scalar_latency`."""
    scalar_reciprocal_throughput::Float64
    """The minimum delay, in clock cycles, associated with the instruction.
    Think of it as the delay from turning on a faucet to when water starts coming out the end of the pipe.
    See also `scalar_reciprocal_throughput`."""
    scalar_latency::Int
    "Number of floating-point registered used"
    register_pressure::Int
end
InstructionCost(sl::Int, srt::Float64, scaling::Float64 = -3.0) = InstructionCost(scaling, srt, sl, 0)

nocost(c::InstructionCost) = c.scalar_reciprocal_throughput == 0.0
flatcost(c::InstructionCost) = c.scaling == -3.0
offsetscaling(c::InstructionCost) = c.scaling == -2.0
linearscaling(c::InstructionCost) = c.scaling == -1.0

function scalar_cost(ic::InstructionCost)#, ::Type{T} = Float64) where {T}
    @unpack scalar_reciprocal_throughput, scalar_latency, register_pressure = ic
    scalar_reciprocal_throughput, scalar_latency, register_pressure
end
function vector_cost(ic::InstructionCost, Wshift, sizeof_T)
    srt, sl, srp = scalar_cost(ic)
    if flatcost(ic) || Wshift == 0 || nocost(ic) # No scaling
        return srt, sl, srp
    elseif offsetscaling(ic) # offset scaling
        srt *= 1 << (Wshift + VectorizationBase.intlog2(sizeof_T) - 4)
        if (sizeof_T << Wshift) == 64 # VectorizationBase.register_size() # These instructions experience double latency with zmm
            sl += sl
        end
    elseif linearscaling(ic) # linear scaling
        W = 1 << Wshift
        extra_latency = sl - srt
        srt *= W
        sl = round(Int, srt + extra_latency)
    # else # we assume custom cost, and that latency == recip_throughput
    #     scaling = ic.scaling
    #     sl, srt = round(Int,scaling), scaling
    end
    srt, sl, srp
end

const OPAQUE_INSTRUCTION = InstructionCost(-1.0, 20, 20.0, 16)

instruction_cost(instruction::Instruction) = instruction.mod === :LoopVectorization ? COST[instruction.instr] : OPAQUE_INSTRUCTION
instruction_cost(instruction::Symbol) = get(COST, instruction, OPAQUE_INSTRUCTION)
scalar_cost(instr::Instruction) = scalar_cost(instruction_cost(instr))
vector_cost(instr::Instruction, Wshift, sizeof_T) = vector_cost(instruction_cost(instr), Wshift, sizeof_T)
# function cost(instruction::InstructionCost, Wshift, sizeof_T)
#     Wshift == 0 ? scalar_cost(instruction) : vector_cost(instruction, Wshift, sizeof_T)
# end

# function cost(instruction::Instruction, Wshift, sizeof_T)
#     cost( instruction_cost(instruction), Wshift, sizeof_T )
# end

# Just a semi-reasonable assumption; should not be that sensitive to anything other than loads

# Comments on setindex!
# 1. Not a part of dependency chains, so not really twice as expensive as getindex?
# 2. getindex loads a register, not setindex!, but we place cost on setindex!
#    as a heuristic means of approximating register pressure, since many loads can be
#    consolidated into a single register. The number of LICM-ed setindex!, on the other
#    hand, should indicate how many registers we're keeping live for the sake of eventually storing.
const COST = Dict{Symbol,InstructionCost}(
    :getindex => InstructionCost(-3.0,0.5,3,0),
    :conditionalload => InstructionCost(-3.0,0.5,3,0),
    :setindex! => InstructionCost(-3.0,1.0,3,0),
    :conditionalstore! => InstructionCost(-3.0,1.0,3,0),
    :zero => InstructionCost(1,0.5),
    :one => InstructionCost(3,0.5),
    :(+) => InstructionCost(4,0.5),
    :(-) => InstructionCost(4,0.5),
    :(*) => InstructionCost(4,0.5),
    :(/) => InstructionCost(13,4.0,-2.0),
    :vadd => InstructionCost(4,0.5),
    :vadd1 => InstructionCost(4,0.5),
    :add_fast => InstructionCost(4,0.5),
    :vsub => InstructionCost(4,0.5),
    :sub_fast => InstructionCost(4,0.5),
    # :vadd! => InstructionCost(4,0.5),
    # :vsub! => InstructionCost(4,0.5),
    # :vmul! => InstructionCost(4,0.5),
    :vmul => InstructionCost(4,0.5),
    :vmul_nsw => InstructionCost(4,0.5),
    :vadd_nsw => InstructionCost(4,0.5),
    :vsub_nsw => InstructionCost(4,0.5),
    :mul_fast => InstructionCost(4,0.5),
    # :vfdiv => InstructionCost(13,4.0,-2.0),
    # :vfdiv! => InstructionCost(13,4.0,-2.0),
    :rem_fast => InstructionCost(13,4.0,-2.0), # FIXME
    :div_fast => InstructionCost(13,4.0,-2.0),
    :vdiv_fast => InstructionCost(20,4.0,-2.0), # FIXME
    :÷ => InstructionCost(13,4.0,-2.0),
    # :evadd => InstructionCost(4,0.5),
    # :evsub => InstructionCost(4,0.5),
    # :evmul => InstructionCost(4,0.5),
    # :evfdiv => InstructionCost(13,4.0,-2.0),
    :vsum => InstructionCost(6,2.0),
    :vprod => InstructionCost(6,2.0),
    :reduced_add => InstructionCost(4,0.5),# ignoring reduction part of cost, might be nop
    :reduced_prod => InstructionCost(4,0.5),# ignoring reduction part of cost, might be nop
    :reduced_max => InstructionCost(4,0.5),# ignoring reduction part of cost, might be nop
    :reduced_min => InstructionCost(4,0.5),# ignoring reduction part of cost, might be nop
    :reduce_to_add => InstructionCost(0,0.0,0.0,0),
    :reduce_to_prod => InstructionCost(0,0.0,0.0,0),
    :abs => InstructionCost(1, 0.5),
    :abs2 => InstructionCost(4,0.5),
    :abs2_fast => InstructionCost(4,0.5),
    :(==) => InstructionCost(1, 0.5),
    :(!=) => InstructionCost(1, 0.5),
    :(≠) => InstructionCost(1, 0.5),
    :(===) => InstructionCost(1, 0.5),
    :(!==) => InstructionCost(1, 0.5),
    :(isnan) => InstructionCost(1, 0.5),
    :(isfinite) => InstructionCost(2, 1.0),
    :(isinf) => InstructionCost(2, 1.0),
    :isequal => InstructionCost(1, 0.5),
    :(!) => InstructionCost(1, 0.5),
    :(~) => InstructionCost(1, 0.5),
    :(&) => InstructionCost(1, 0.5),
    :(|) => InstructionCost(1, 0.5),
    :(⊻) => InstructionCost(1, 0.5),
    :(%) => InstructionCost(13, 4.0, -2.0),
    :(rem) => InstructionCost(13, 4.0, -2.0),
    :(>) => InstructionCost(1, 0.5),
    :(<) => InstructionCost(1, 0.5),
    :(>=) => InstructionCost(1, 0.5),
    :(<=) => InstructionCost(1, 0.5),
    :(≥) => InstructionCost(1, 0.5),
    :(≤) => InstructionCost(1, 0.5),
    :(>>) => InstructionCost(1, 0.5),
    :(>>>) => InstructionCost(1, 0.5),
    :(<<) => InstructionCost(1, 0.5),
    :isodd => InstructionCost(1, 0.5),
    :iseven => InstructionCost(1, 0.5),
    :max => InstructionCost(4,0.5),
    :min => InstructionCost(4,0.5),
    :max_fast => InstructionCost(4,0.5),
    :min_fast => InstructionCost(4,0.5),
    :relu => InstructionCost(4,0.5),
    # Instruction(:ifelse) => InstructionCost(1, 0.5),
    :ifelse => InstructionCost(1, 0.5),
    :inv => InstructionCost(13,4.0,-2.0,1),
    :inv_fast => InstructionCost(10,4.0,-2.0,1), # FIXME
    # :vinv => InstructionCost(13,4.0,-2.0,1),
    :muladd => InstructionCost(4,0.5), # + and * will fuse into this, so much of the time they're not twice as expensive
    :fma => InstructionCost(4,0.5), # + and * will fuse into this, so much of the time they're not twice as expensive
    :vmuladd_fast => InstructionCost(4,0.5), # + and * will fuse into this, so much of the time they're not twice as expensive
    :vfma_fast => InstructionCost(4,0.5), # + and * will fuse into this, so much of the time they're not twice as expensive
    :vfmadd => InstructionCost(4,0.5), # + and * will fuse into this, so much of the time they're not twice as expensive
    :vfmsub => InstructionCost(4,0.5), # - and * will fuse into this, so much of the time they're not twice as expensive
    :vfnmadd => InstructionCost(4,0.5), # + and -* will fuse into this, so much of the time they're not twice as expensive
    :vfnmsub => InstructionCost(4,0.5), # - and -* will fuse into this, so much of the time they're not twice as expensive
    :vfmadd_fast => InstructionCost(4,0.5), # + and * will fuse into this, so much of the time they're not twice as expensive
    :vfmsub_fast => InstructionCost(4,0.5), # - and * will fuse into this, so much of the time they're not twice as expensive
    :vfnmadd_fast => InstructionCost(4,0.5), # + and -* will fuse into this, so much of the time they're not twice as expensive
    :vfnmsub_fast => InstructionCost(4,0.5), # - and -* will fuse into this, so much of the time they're not twice as expensive
    :vfmadd231 => InstructionCost(4,0.5), # + and * will fuse into this, so much of the time they're not twice as expensive
    :vfmsub231 => InstructionCost(4,0.5), # - and * will fuse into this, so much of the time they're not twice as expensive
    :vfnmadd231 => InstructionCost(4,0.5), # + and -* will fuse into this, so much of the time they're not twice as expensive
    :vfnmsub231 => InstructionCost(4,0.5), # - and -* will fuse into this, so much of the time they're not twice as expensive
    :vfmaddsub => InstructionCost(4,0.5),
    :vfmsubadd => InstructionCost(4,0.5),
    :sqrt => InstructionCost(15,4.0,-2.0),
    :sqrt_fast => InstructionCost(15,4.0,-2.0),
    :log => InstructionCost(-3.0, 15, 30, 11),
    :log2 => InstructionCost(-3.0, 15, 30, 11),
    :log10 => InstructionCost(-3.0, 15, 30, 11),
    :log_fast => InstructionCost(-3.0, 15, 30, 11),
    :log2_fast => InstructionCost(-3.0, 15, 30, 11),
    :log10_fast => InstructionCost(-3.0, 15, 30, 11),
    :log1p => InstructionCost(-3.0, 15, 30, 11),
    :exp => InstructionCost(-3.0,13.0,26.0,14),
    :exp2 => InstructionCost(-3.0,10.0,40.0,14),
    :exp10 => InstructionCost(-3.0,13.0,26.0,14),
    :expm1 => InstructionCost(-3.0,30.0,60.0,19),
    :(^) => InstructionCost(-3.0,200.0,400.0,26), # FIXME
    :pow_fast => InstructionCost(-3.0,200.0,400.0,26), # FIXME
    :sin => InstructionCost(-3,30.0,60.0,23),
    :cos => InstructionCost(-3,27.0,60.0,26),
    :sincos => InstructionCost(-3,37.0,85.0,26),
    :sincos_fast => InstructionCost(-3,37.0,85.0,26),
    :sinpi => InstructionCost(18,15.0,68.0,23),
    :cospi => InstructionCost(18,15.0,68.0,26),
    :sincospi => InstructionCost(25,37.0,70.0,26),
    :log_fast => InstructionCost(20,20.0,40.0,20),
    :exp_fast => InstructionCost(20,20.0,20.0,18),
    :sin_fast => InstructionCost(18,15.0,68.0,23),
    :cos_fast => InstructionCost(18,15.0,68.0,26),
    :sinpi_fast => InstructionCost(18,15.0,68.0,23),
    :cospi_fast => InstructionCost(18,15.0,68.0,26),
    :sincospi_fast => InstructionCost(25,22.0,70.0,26),
    :tanh => InstructionCost(-3.0,80.0,160.0,26), # FIXME
    :tanh_fast => InstructionCost(-3.0,30.0,60.0,20), # FIXME
    :sigmoid_fast => InstructionCost(-3.0,16.0,66.0,15), # FIXME
    :identity => InstructionCost(0,0.0,0.0,0),
    :adjoint => InstructionCost(0,0.0,0.0,0),
    :conj => InstructionCost(0,0.0,0.0,0),
    :transpose => InstructionCost(0,0.0,0.0,0),
    :first => InstructionCost(0,0.0,0.0,0),
    :second => InstructionCost(0,0.0,0.0,0),
    :third => InstructionCost(0,0.0,0.0,0),
    :fourth => InstructionCost(0,0.0,0.0,0),
    :fifth => InstructionCost(0,0.0,0.0,0),
    :sixth => InstructionCost(0,0.0,0.0,0),
    :seventh => InstructionCost(0,0.0,0.0,0),
    :eighth => InstructionCost(0,0.0,0.0,0),
    :ninth => InstructionCost(0,0.0,0.0,0),
    :tenth => InstructionCost(0,0.0,0.0,0),
    :eleventh => InstructionCost(0,0.0,0.0,0),
    :twelfth => InstructionCost(0,0.0,0.0,0),
    :thirteenth => InstructionCost(0,0.0,0.0,0),
    :last => InstructionCost(0,0.0,0.0,0),
    :prefetch => InstructionCost(0,0.0,0.0,0),
    :prefetch0 => InstructionCost(0,0.0,0.0,0),
    :prefetch1 => InstructionCost(0,0.0,0.0,0),
    :prefetch2 => InstructionCost(0,0.0,0.0,0),
  :convert => InstructionCost(4,0.5),
  :vpermilps177 => InstructionCost(1, 1.0),
  :vmovsldup => InstructionCost(1, 1.0),
  :vmovshdup => InstructionCost(1, 1.0),
  :exponent => InstructionCost(8, 1.0),
  :significand => InstructionCost(8, 1.0)
)

# # @inline prefetch0(x::Ptr, i) = VectorizationBase.prefetch(x, Val{3}(), Val{0}())
# @inline prefetch0(x, i) = VectorizationBase.prefetch(gep(stridedpointer(x), (data(i),)), Val{3}(), Val{0}())
# @inline prefetch0(x, I::Tuple) = VectorizationBase.prefetch(gep(x, I), Val{3}(), Val{0}())
# @inline prefetch0(x, i, j) = VectorizationBase.prefetch(gep(stridedpointer(x), (data(i), data(j))), Val{3}(), Val{0}())
# # @inline prefetch0(x, i, j, oi, oj) = VectorizationBase.prefetch(gep(stridedpointer(x), (data(i) + data(oi) - 1, data(j) + data(oj) - 1)), Val{3}(), Val{0}())
# @inline prefetch1(x, i) = VectorizationBase.prefetch(gep(stridedpointer(x), (data(i),)), Val{2}(), Val{0}())
# @inline prefetch1(x, i, j) = VectorizationBase.prefetch(gep(stridedpointer(x), (data(i), data(j))), Val{2}(), Val{0}())
# # @inline prefetch1(x, i, j, oi, oj) = VectorizationBase.prefetch(gep(stridedpointer(x), (data(i) + data(oi) - 1, data(j) + data(oj) - 1)), Val{2}(), Val{0}())
# @inline prefetch2(x, i) = VectorizationBase.prefetch(gep(stridedpointer(x), (data(i),)), Val{1}(), Val{0}())
# @inline prefetch2(x, i, j) = VectorizationBase.prefetch(gep(stridedpointer(x), (data(i), data(j))), Val{1}(), Val{0}())
# # @inline prefetch2(x, i, j, oi, oj) = VectorizationBase.prefetch(gep(stridedpointer(x), (data(i) + data(oi) - 1, data(j) + data(oj) - 1)), Val{1}(), Val{0}())

Base.convert(::Type{Instruction}, instr::Symbol) = Instruction(instr)
# instruction(f::Symbol, m) = f ∈ keys(COST) ? Instruction(f) : Instruction(m, f)
# instruction(f::Symbol) = f ∈ keys(COST) ? Instruction(:LoopVectorization, f) : Instruction(Symbol(""), f)
function instruction(f::Symbol)
    # f === :ifelse && return Instruction(:LoopVectorization, :ifelse)
    # @assert f ∈ keys(COST)
    f ∈ keys(COST) ? Instruction(:LoopVectorization, f) : Instruction(Symbol(""), f)
end
# instruction(f::Symbol, m::Symbol) = f ∈ keys(COST) ? Instruction(:LoopVectorization, f) : Instruction(m, f)
Instruction(instr::Symbol) = instruction(instr)

struct IfElseOp{F}; f::F; end
struct IfElseReducer{F}; f::F; end
struct IfElseReduced{F}; f::F; end
struct IfElseReduceTo{F}; f::F; end
struct IfElseCollapser{F}; f::F; end
@inline (ieo::IfElseOp)(a, b) = ifelse(ieo.f(a, b), a, b)

@inline (ier::IfElseReducer)(a) = VectorizationBase.ifelse_reduce(ier.f, a)
@inline function (ier::IfElseReducer)(a, b)
  f = ier.f
  r = VectorizationBase.ifelse_reduce(f, b)
  ifelse(f(a,r), a, r)
end
@inline (ier::IfElseReducer)(a::VecUnroll) = VecUnroll(VectorizationBase.fmap(ier, VectorizationBase.data(a)))
@inline (ier::IfElseReducer)(a::VecUnroll, b::VecUnroll) = VecUnroll(VectorizationBase.fmap(ier, VectorizationBase.data(a), VectorizationBase.data(b)))

@inline (ier::IfElseReduced)(x::NativeTypes, y::NativeTypes) = ifelse(ier.f(x,y), x, y)
@inline (ier::IfElseReduced)(x::AbstractSIMD{W}, y::AbstractSIMD{W}) where {W} = ifelse(ier.f(x,y), x, y)
@inline function (ier::IfElseReduced)(x::AbstractSIMD, y::AbstractSIMD)
  r = IfElseReduceTo(ier.f)(x, y)
  ifelse(ier.f(r,y), r, y)
end
@inline (ier::IfElseReduced)(x::VecUnroll, y::VecUnroll) = VecUnroll(fmap(ier, getfield(x, :data), getfield(y, :data)))
@inline function (ier::IfElseReduced)(x::AbstractSIMD, y::NativeTypes)
  f = ier.f
  r = VectorizationBase.ifelse_reduce(f, x)
  ifelse(f(r, y), r, y)
end

@inline (ier::IfElseReduceTo)(a::NativeTypes, ::NativeTypes) = a
@inline (ier::IfElseReduceTo)(a::AbstractSIMD, ::NativeTypes) = VectorizationBase.ifelse_reduce(ier.f, a)
@inline (ier::IfElseReduceTo)(a::AbstractSIMD{W}, ::AbstractSIMD{W}) where {W} = a
@inline function (ier::IfElseReduceTo)(a::AbstractSIMD, b::AbstractSIMD)
  x, y = VectorizationBase.splitvector(a) # halve recursively
  ier(ifelse(ier.f(x,y), x, y), b)
end
@inline (ier::IfElseReduceTo)(a::VecUnroll, b::VecUnroll) = VecUnroll(VectorizationBase.fmap(ier, VectorizationBase.data(a), VectorizationBase.data(b)))

@inline (iec::IfElseCollapser)(a) = VectorizationBase.collapse(IfElseOp(iec.f), a)
@inline (iec::IfElseCollapser)(a, ::StaticInt{C}) where {C} = VectorizationBase.contract(IfElseOp(iec.f), a, StaticInt{C}())

struct IfElseOpMirror{F,A,B}; f::F; a::A; b::B; end
struct IfElseReducerMirror{F,A,B}; f::F; a::A; b::B; end
struct IfElseReducedMirror{F,A,B}; f::F; a::A; b::B; end
struct IfElseReduceToMirror{F,A}; f::F; a::A; end
struct IfElseCollapserMirror{F,A}; f::F; a::A; end
@inline (ieo::IfElseOpMirror)(a, b) = ifelse(ieo.f(ieo.a, ieo.b), a, b)

@inline _first_ifelse_reduce_mirror(f::F, a, b) where {F} = getfield(VectorizationBase.ifelse_reduce_mirror(f, a, b), 1, false)
@inline (ier::IfElseReducerMirror)(a) = _ifelse_reduce_mirror(ier.f, a, ier.a)
@inline function _ifelse_reduce_mirror(f::F, a, b, c, d) where {F}
  r, rm = VectorizationBase.ifelse_reduce_mirror(f, b, d)
  ifelse(f(c, rm), a, r)
end
@inline (ier::IfElseReducerMirror)(a, b) = _ifelse_reduce_mirror(ier.f, a, b, ier.a, ier.b)
@inline (ier::IfElseReducerMirror)(a::VecUnroll) = VecUnroll(VectorizationBase.fmap(_first_ifelse_reduce_mirror, ier.f, VectorizationBase.data(a), VectorizationBase.data(ier.a)))
@inline function (ier::IfElseReducerMirror)(a::VecUnroll, b::VecUnroll)
  VecUnroll(VectorizationBase.fmap(_ifelse_reduce_mirror, ier.f, VectorizationBase.data(a), VectorizationBase.data(b), VectorizationBase.data(ier.a), VectorizationBase.data(ier.b)))
end

@inline IfElseReducedMirror(f::F, a::A) where {F,A} = IfElseReducedMirror{F,A,Nothing}(f, a, nothing)
@inline (ier::IfElseReducedMirror)(x::NativeTypes, y::NativeTypes) = ifelse(ier.f(ier.a, ier.b), x, y)
@inline (ier::IfElseReducedMirror)(x::AbstractSIMD{W}, y::AbstractSIMD{W}) where {W} = ifelse(ier.f(ier.a, ier.b), x, y)
@inline function _reduce_mirror(f::F, x, y, a, b) where {F}
  r, rm = IfElseReduceToMirror(f, a, b)(x, y)
  ifelse(f(r,y), r, y)
end
@inline (ier::IfElseReducedMirror)(x::AbstractSIMD, y::AbstractSIMD) = _reduce_mirror(ier.f, x, y, ier.a, ier.b)
@inline (ier::IfElseReducedMirror)(x::VecUnroll, y::VecUnroll) = VecUnroll(fmap(_reduce_mirror, ier.f, getfield(x, :data), getfield(y, :data), getfield(ier.a, :data), getfield(ier.b, :data)))
@inline function (ier::IfElseReducedMirror)(x::AbstractSIMD, y::NativeTypes)
  f = ier.f
  r, rm = VectorizationBase.ifelse_reduce_mirror(f, x, ier.a)
  # @show rm, ier.b, r, y
  ifelse(f(rm, ier.b), r, y)
end

@inline (ier::IfElseReduceToMirror)(a::NativeTypes, ::NativeTypes) = a
@inline (ier::IfElseReduceToMirror)(a::AbstractSIMD, ::NativeTypes) = VectorizationBase.ifelse_reduce_mirror(ier.f, a, ier.a)
@inline (ier::IfElseReduceToMirror)(a::AbstractSIMD{W}, ::AbstractSIMD{W}) where {W} = a
@inline function (ier::IfElseReduceToMirror)(a::AbstractSIMD, b::AbstractSIMD)
  x, y = VectorizationBase.splitvector(a) # halve recursively
  w, z = VectorizationBase.splitvector(ier.a) # halve recursively
  f = ier.f
  fwz = f(w,z)
  IfElseReduceToMirror(f, ifelse(fwz, w, z))(ifelse(fwz, x, y), b)
end
@inline (ier::IfElseReduceToMirror)(a::VecUnroll, b::VecUnroll) = VecUnroll(VectorizationBase.fmap(ier, VectorizationBase.data(a), VectorizationBase.data(b)))

@inline (iec::IfElseCollapserMirror)(a) = getfield(VectorizationBase.ifelse_collapse_mirror(iec.f, a, iec.a), 1, false)
@inline (iec::IfElseCollapserMirror)(a, ::StaticInt{1}) = getfield(VectorizationBase.ifelse_collapse_mirror(iec.f, a, iec.a), 1, false)
# @inline function (iec::IfElseCollapserMirror)(a, ::StaticInt{C}) where {C}
#   VectorizationBase.contract(IfElseOp(iec.f), a, StaticInt{C}())
# end

const ADDITIVE_IN_REDUCTIONS = 1.0
const MULTIPLICATIVE_IN_REDUCTIONS = 2.0
const ANY = 3.0
const ALL = 4.0
const MAX = 5.0
const MIN = 6.0
const IFELSE = 7.0

const REDUCTION_CLASS = Dict{Symbol,Float64}(
    :+ => ADDITIVE_IN_REDUCTIONS,
    :- => ADDITIVE_IN_REDUCTIONS,
    :* => MULTIPLICATIVE_IN_REDUCTIONS,
    :vadd => ADDITIVE_IN_REDUCTIONS,
    :vsub => ADDITIVE_IN_REDUCTIONS,
    :add_fast => ADDITIVE_IN_REDUCTIONS,
    :sub_fast => ADDITIVE_IN_REDUCTIONS,
    :vadd_fast => ADDITIVE_IN_REDUCTIONS,
    :vsub_fast => ADDITIVE_IN_REDUCTIONS,
    # :vadd! => ADDITIVE_IN_REDUCTIONS,
    # :vsub! => ADDITIVE_IN_REDUCTIONS,
    :vmul => MULTIPLICATIVE_IN_REDUCTIONS,
    :mul_fast => MULTIPLICATIVE_IN_REDUCTIONS,
    :vmul_fast => MULTIPLICATIVE_IN_REDUCTIONS,
    # :vmul! => MULTIPLICATIVE_IN_REDUCTIONS,
    # :evadd => ADDITIVE_IN_REDUCTIONS,
    # :evsub => ADDITIVE_IN_REDUCTIONS,
    # :evmul => MULTIPLICATIVE_IN_REDUCTIONS,
    :& => ALL,
    :| => ANY,
    :muladd => ADDITIVE_IN_REDUCTIONS,
    :fma => ADDITIVE_IN_REDUCTIONS,
    :vmuladd_fast => ADDITIVE_IN_REDUCTIONS,
    :vfma_fast => ADDITIVE_IN_REDUCTIONS,
    :vfmadd => ADDITIVE_IN_REDUCTIONS,
    :vfmsub => ADDITIVE_IN_REDUCTIONS,
    :vfnmadd => ADDITIVE_IN_REDUCTIONS,
    :vfnmsub => ADDITIVE_IN_REDUCTIONS,
    :vfmadd_fast => ADDITIVE_IN_REDUCTIONS,
    :vfmsub_fast => ADDITIVE_IN_REDUCTIONS,
    :vfnmadd_fast => ADDITIVE_IN_REDUCTIONS,
    :vfnmsub_fast => ADDITIVE_IN_REDUCTIONS,
    :vfmadd231 => ADDITIVE_IN_REDUCTIONS,
    :vfmsub231 => ADDITIVE_IN_REDUCTIONS,
    :vfnmadd231 => ADDITIVE_IN_REDUCTIONS,
    :vfnmsub231 => ADDITIVE_IN_REDUCTIONS,
    # :vfmadd! => ADDITIVE_IN_REDUCTIONS,
    # :vfnmadd! => ADDITIVE_IN_REDUCTIONS,
    # :vfmsub! => ADDITIVE_IN_REDUCTIONS,
    # :vfnmsub! => ADDITIVE_IN_REDUCTIONS,
    :reduced_add => ADDITIVE_IN_REDUCTIONS,
    :reduced_prod => MULTIPLICATIVE_IN_REDUCTIONS,
    :reduced_all => ALL,
    :reduced_any => ANY,
    :max => MAX,
    :min => MIN,
    :max_fast => MAX,
  :min_fast => MIN,
  :vfmaddsub => ADDITIVE_IN_REDUCTIONS,
  :vfmsubadd => ADDITIVE_IN_REDUCTIONS
)
reduction_instruction_class(instr::Symbol) = get(REDUCTION_CLASS, instr, NaN)
reduction_instruction_class(instr::Instruction) = reduction_instruction_class(instr.instr)
function reduction_to_single_vector(x::Float64)
    if x == ADDITIVE_IN_REDUCTIONS
        :collapse_add
    elseif x == MULTIPLICATIVE_IN_REDUCTIONS
        :collapse_mul
    elseif x == MAX
        :collapse_max
    elseif x == MIN
        :collapse_min
    elseif x == ALL
        :collapse_and
    elseif x == ANY
        :collapse_or
    else
        throw("Reduction not found.")
    end
end
function reduce_to_onevecunroll(x::Float64)
    if x == ADDITIVE_IN_REDUCTIONS
        :+
    elseif x == MULTIPLICATIVE_IN_REDUCTIONS
        :*
    elseif x == MAX
        :max
    elseif x == MIN
        :min
    elseif x == ALL
        :&
    elseif x == ANY
        :|
    else
        throw("Reduction not found.")
    end
end
function reduce_number_of_vectors(x::Float64)
    if x == ADDITIVE_IN_REDUCTIONS
        :contract_add
    elseif x == MULTIPLICATIVE_IN_REDUCTIONS
        :contract_mul
    elseif x == MAX
        :contract_max
    elseif x == MIN
        :contract_min
    elseif x == ALL
        :contract_and
    elseif x == ANY
        :contract_or
    else
        throw("Reduction not found.")
    end
end
function reduction_to_scalar(x::Float64)
    if x == ADDITIVE_IN_REDUCTIONS
        :vsum
    elseif x == MULTIPLICATIVE_IN_REDUCTIONS
        :vprod
    elseif x == MAX
        :vmaximum
    elseif x == MIN
        :vminimum
    elseif x == ALL
        :vall
    elseif x == ANY
        :vany
    else
        throw("Reduction not found.")
    end
end
function reduction_scalar_combine(x::Float64)
    # x == 1.0 ? :reduced_add : x == 2.0 ? :reduced_prod : x == 3.0 ? :reduced_any : x == 4.0 ? :reduced_all : x == 5.0 ? :reduced_max : x == 6.0 ? :reduced_min : throw("Reduction not found.")
  if x == ADDITIVE_IN_REDUCTIONS
        :reduced_add
    elseif x == MULTIPLICATIVE_IN_REDUCTIONS
        :reduced_prod
    elseif x == MAX
        :reduced_max
    elseif x == MIN
        :reduced_min
    elseif x == ALL
        :reduced_all
    elseif x == ANY
        :reduced_any
    else
        throw("Reduction not found.")
    end
end

function reduction_zero(x::Float64)
    # x == 1.0 ? :zero : x == 2.0 ? :one : x == 3.0 ? :false : x == 4.0 ? :true : x == 5.0 ? :typemin : x == 6.0 ? :typemax : throw("Reduction not found.")
    if x == ADDITIVE_IN_REDUCTIONS
        :zero
    elseif x == MULTIPLICATIVE_IN_REDUCTIONS
        :one
    elseif x == MAX
        :typemin
    elseif x == MIN
        :typemax
    elseif x == ALL
        :max_mask
    elseif x == ANY
        :zero_mask
    else
        throw("Reduction not found.")
    end
end
function reduction_zero_class(x::Symbol)::Float64
  if x === :one
    MULTIPLICATIVE_IN_REDUCTIONS
  elseif x === :typemin
    MAX
  elseif x === :typemax
    MIN
  elseif x === :max_mask
    ALL
  elseif x === :zero_mask
    ANY
  elseif x === :zero#sorted last, as should go into preamble_zeros
    ADDITIVE_IN_REDUCTIONS
  else
    throw("Reduction not found.")
  end
end
reduction_zero(x) = reduction_zero(reduction_instruction_class(x))

function isreductcombineinstr(instr::Symbol)
    instr ∈ (:reduced_add, :reduced_prod, :reduce_to_add, :reduce_to_prod, :reduced_max, :reduced_min, :reduce_to_max, :reduce_to_min)
end
isreductcombineinstr(instr::Instruction) = isreductcombineinstr(instr.instr)

const FUNCTIONSYMBOLS = IdDict{Type{<:Function},Instruction}(
    typeof(+) => :add_fast,
    typeof(VectorizationBase.vadd) => :add_fast,
    typeof(add_fast) => :add_fast,
    typeof(-) => :sub_fast,
    typeof(VectorizationBase.vsub) => :sub_fast,
    typeof(sub_fast) => :sub_fast,
    typeof(*) => :mul_fast,
    typeof(VectorizationBase.vmul) => :mul_fast,
    typeof(mul_fast) => :mul_fast,
    typeof(/) => :div_fast,
    typeof(^) => :pow_fast,
    typeof(VectorizationBase.vdiv) => :div_fast,
    typeof(÷) => :(÷),
    typeof(div_fast) => :div_fast,
    typeof(rem_fast) => :rem_fast,
    typeof(==) => :(==),
    typeof(!=) => :(!=),
    typeof(isequal) => :isequal,
    typeof(isnan) => :isnan,
    typeof(isinf) => :isinf,
    typeof(isfinite) => :isfinite,
    typeof(abs) => :abs,
    typeof(abs2) => :abs2_fast,
    typeof(abs2_fast) => :abs2_fast,
    typeof(~) => :(~),
    typeof(!) => :(!),
    typeof(&) => :(&),
    typeof(|) => :(|),
    typeof(⊻) => :(⊻),
    typeof(>) => :(>),
    typeof(<) => :(<),
    typeof(>=) => :(>=),
    typeof(<=) => :(<=),
    typeof(inv) => :inv_fast,
    typeof(inv_fast) => :inv_fast,
    typeof(muladd) => :vmuladd_fast,
    typeof(fma) => :vfma_fast,
    typeof(VectorizationBase.vfma) => :vfma_fast,
    typeof(VectorizationBase.vmuladd) => :vmuladd_fast,
    typeof(VectorizationBase.vfmsub) => :vfmsub_fast,
    typeof(VectorizationBase.vfnmadd) => :vfnmadd_fast,
    typeof(VectorizationBase.vfnmsub) => :vfnmsub_fast,
    typeof(VectorizationBase.vfma_fast) => :vfma_fast,
    typeof(VectorizationBase.vmuladd_fast) => :vmuladd_fast,
    typeof(VectorizationBase.vfmsub_fast) => :vfmsub_fast,
    typeof(VectorizationBase.vfnmadd_fast) => :vfnmadd_fast,
    typeof(VectorizationBase.vfnmsub_fast) => :vfnmsub_fast,
    typeof(VectorizationBase.vfmadd231) => :vfmadd231,
    typeof(VectorizationBase.vfmsub231) => :vfmsub231,
    typeof(VectorizationBase.vfnmadd231) => :vfnmadd231,
    typeof(VectorizationBase.vfnmsub231) => :vfnmsub231,
    typeof(sqrt) => :sqrt_fast,
    typeof(sqrt_fast) => :sqrt_fast,
    typeof(log) => :log_fast,
    typeof(log2) => :log2_fast,
    typeof(log10) => :log10_fast,
    typeof(log_fast) => :log_fast,
    typeof(log1p) => :log1p,
    # typeof(VectorizationBase.vlog) => :log,
    typeof(SLEEFPirates.log) => :log,
    typeof(exp) => :exp,
    typeof(exp2) => :exp2,
    typeof(exp10) => :exp10,
    typeof(Base.FastMath.exp_fast) => :exp,
    typeof(expm1) => :expm1,
    # typeof(VectorizationBase.vexp) => :exp,
    typeof(SLEEFPirates.exp) => :exp,
    typeof(sin) => :sin_fast,
    typeof(Base.FastMath.sin_fast) => :sin_fast,
    typeof(SLEEFPirates.sin) => :sin_fast,
    typeof(cos) => :cos_fast,
    typeof(Base.FastMath.cos_fast) => :cos_fast,
    typeof(SLEEFPirates.cos) => :cos_fast,
    typeof(sincos) => :sincos_fast,
    typeof(Base.FastMath.sincos_fast) => :sincos_fast,
    typeof(SLEEFPirates.sincos) => :sincos_fast,
    typeof(tan) => :tan_fast,
    typeof(Base.FastMath.sincos_fast) => :sincos_fast,
    typeof(SLEEFPirates.sincos) => :sincos_fast,
    typeof(Base.tanh) => :tanh,
    typeof(tanh_fast) => :tanh_fast,
    typeof(sigmoid_fast) => :sigmoid_fast,
    typeof(max) => :max_fast,
    typeof(min) => :min_fast,
    typeof(max_fast) => :max_fast,
    typeof(min_fast) => :min_fast,
    typeof(relu) => :relu,
    typeof(<<) => :<<,
    typeof(>>) => :>>,
    typeof(>>>) => :>>>,
    typeof(%) => :(%),
    typeof(÷) => :(÷),
    typeof(Base.ifelse) => :ifelse,
    typeof(ifelse) => :ifelse,
    typeof(identity) => :identity,
    typeof(conj) => :identity,#conj,
    typeof(÷) => :vdiv_fast
    # typeof(zero) => :zero,
    # typeof(one) => :one,
    # typeof(axes) => :axes,
    # typeof(eltype) => :eltype
)

# implement whitelist for avx_support that package authors may use to conservatively guard `@turbo` application
for f ∈ keys(FUNCTIONSYMBOLS)
    @eval ArrayInterface.can_avx(::$f) = true
end

ext: jl


text:
using MD5: MD5
using Nettle: Nettle
using BenchmarkTools

data = randstring(10^4)

@show Nettle.hexdigest("md5", data)
display(@btime Nettle.hexdigest("md5", data))
println()

@show bytes2hex(MD5.md5(data))
display(@btime MD5.md5(data))
println()


text:
using Documenter
using HelFEM
using Plots

makedocs(
    sitename="HelFEM.jl",
    pages = [
        "Overview" => "index.md",
        "polynomial.md",
        "radialbasis.md",
        "FEM basis" => "fembasis.md",
        "utilities.md",
    ],
)

deploydocs(repo = "github.com/mortenpi/HelFEM.jl")


text:
function molecule(::Molecule"acetic_acid")
  return """
C	0.9843630	-0.9796550	0.0000000
C	0.0000000	0.1414350	0.0000000
O	0.2954430	1.3096060	0.0000000
H	1.9693630	-0.5399610	0.0000000
H	0.8492070	-1.6037890	0.8786330
H	0.8492070	-1.6037890	-0.8786330
O	-1.2695100	-0.2791920	0.0000000
H	-1.7814160	0.5335460	0.0000000
  """
end


text:
import UUIDs

"""Internal storage common to InfrastructureSystems types."""
mutable struct InfrastructureSystemsInternal
    uuid::Base.UUID
    ext::Union{Nothing, Dict{String, Any}}
end

"""
Creates PowerSystemInternal with a new UUID.
"""
InfrastructureSystemsInternal() = InfrastructureSystemsInternal(UUIDs.uuid4(), nothing)

"""
Creates PowerSystemInternal with an existing UUID.
"""
InfrastructureSystemsInternal(u::UUIDs.UUID) = InfrastructureSystemsInternal(u, nothing)

"""
Return a user-modifiable dictionary to store extra information.
"""
function get_ext(obj::InfrastructureSystemsInternal)
    if isnothing(obj.ext)
        obj.ext = Dict{String, Any}()
    end

    return obj.ext
end

"""
Clear any value stored in ext.
"""
function clear_ext(obj::InfrastructureSystemsInternal)
    obj.ext = nothing
end

"""
Gets the UUID for any PowerSystemType.
"""
function get_uuid(obj::InfrastructureSystemsType)::Base.UUID
    return obj.internal.uuid
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions related with the IERS (International Earth Orientation and
#   Reference Systems Service) EOP (Earth Orientation Parameters) data.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export EOPData_IAU1980, EOPData_IAU2000A
export get_iers_eop, get_iers_eop_iau_1980, get_iers_eop_iau_2000A
export read_iers_eop

################################################################################
#                       Private Structures and Variables
################################################################################

"""
    get_iers_eop(data_type::Symbol = :IAU1980; force_download = false)

Download and parse the IERS EOP C04 data. The data type is specified by
`data_type` symbol. Supported values are:

* `IAU1980`: Get IERS EOP C04 IAU1980 data.
* `IAU2000A`: Get IERS EOP C04 IAU2000A data.

If `data_type` is omitted, then it defaults to `IAU1980`.

The files are downloaded using the `RemoteFile` package with daily updates.
Hence, if one desires to force a download before the scheduled time, then set
the keyword `force_download` to `true`.

# Returns

A structure (`EOPData_IAU1980` or `EOPData_IAU2000A`, depending on `data_type`)
with the interpolations of the EOP parameters. Notice that the interpolation
indexing is set to the Julian Day.

"""
function get_iers_eop(data_type::Symbol = :IAU1980; force_download = false)
    if data_type == :IAU1980
        return get_iers_eop_iau_1980(force_download = force_download)
    elseif data_type == :IAU2000A
        return get_iers_eop_iau_2000A(force_download = force_download)
    else
        throw(ArgumentError("Unknow EOP type. It must be :IAU1980 or :IAU2000A."))
    end
end

"""
    get_iers_eop_iau_1980(url::String = "https://datacenter.iers.org/data/latestVersion/223_EOP_C04_14.62-NOW.IAU1980223.txt")

Get the IERS EOP C04 IAU1980 data from the URL `url`. If `url` is omitted, then
it defaults to https://datacenter.iers.org/data/latestVersion/223_EOP_C04_14.62-NOW.IAU1980223.txt

The file is downloaded using the `RemoteFile` package with daily updates. Hence,
if one desires to force a download before the scheduled time, then set the
keyword `force_download` to `true`.

# Returns

The structure `EOPData_IAU1980` with the interpolations of the EOP parameters.
Notice that the interpolation indexing is set to the Julian Day.

# Remarks

For every field in `EOPData_IAU1980` to interpolation between two points in the
grid is linear. If extrapolation is needed, then if will use the nearest value
(flat extrapolation).

"""
function get_iers_eop_iau_1980(url::String = "https://datacenter.iers.org/data/latestVersion/223_EOP_C04_14.62-NOW.IAU1980223.txt";
                               force_download = false)

    _eop_iau1980 = @RemoteFile(@eval($url),
                               file="EOP_IAU1980.TXT",
                               updates=:daily)
    download(_eop_iau1980; force = force_download)

    # Parse the data removing the header.
    eop = readdlm(path(_eop_iau1980); skipblanks=true, skipstart=14)

    # Return the parsed EOP data.
    parse_iers_eop_iau_1980(eop)
end

"""
    get_iers_eop_iau_2000A(url::String = "https://datacenter.iers.org/data/latestVersion/224_EOP_C04_14.62-NOW.IAU2000A224.txt"; force_download = false)

Get the IERS EOP C04 IAU2000A data from the URL `url`. If `url` is omitted, then
it defaults to https://datacenter.iers.org/data/latestVersion/224_EOP_C04_14.62-NOW.IAU2000A224.txt

The file is downloaded using the `RemoteFile` package with daily updates. Hence,
if one desires to force a download before the scheduled time, then set the
keyword `force_download` to `true`.

# Returns

The structure `EOPData_IAU2000A` with the interpolations of the EOP parameters.
Notice that the interpolation indexing is set to the Julian Day.

# Remarks

For every field in `EOPData_IAU2000A` to interpolation between two points in the
grid is linear. If extrapolation is needed, then if will use the nearest value
(flat extrapolation).

"""
function get_iers_eop_iau_2000A(url::String = "https://datacenter.iers.org/data/latestVersion/224_EOP_C04_14.62-NOW.IAU2000A224.txt";
                                force_download = false)

    _eop_iau2000A = @RemoteFile(@eval($url),
                                file="EOP_IAU2000A.TXT",
                                updates=:daily)
    download(_eop_iau2000A; force = force_download)

    # Parse the data removing the header.
    eop = readdlm(path(_eop_iau2000A); skipblanks=true, skipstart=14)

    # Return the parsed EOP data.
    parse_iers_eop_iau_2000A(eop)
end

"""
    read_iers_eop(filename::String, data_type::Symbol = :IAU1980)

Read IERS EOP Data from the file `filename`. The user must specify if the data
is related to the model IAU 1980 (`data_type = :IAU1980`), which is the default,
or to the model IAU 2000A (`data_type = :IAU2000A`).

# Returns

A structure (`EOPData_IAU1980` or `EOPData_IAU2000A`, depending on `data_type`)
with the interpolations of the EOP parameters. Notice that the interpolation
indexing is set to the Julian Day.

# Remarks

The input file **must be exactly the same** as provided by IERS. One can
download it using the following commands:

* IAU 1980

        curl -O https://datacenter.iers.org/data/latestVersion/223_EOP_C04_14.62-NOW.IAU1980223.txt
        wget https://datacenter.iers.org/data/latestVersion/223_EOP_C04_14.62-NOW.IAU1980223.txt

* IAU 2000A

        curl -O https://datacenter.iers.org/data/latestVersion/224_EOP_C04_14.62-NOW.IAU2000A224.txt
        wget https://datacenter.iers.org/data/latestVersion/224_EOP_C04_14.62-NOW.IAU2000A224.txt

"""
function read_iers_eop(filename::String, data_type::Symbol = :IAU1980)
    # Parse the data removing the header.
    eop = readdlm(filename; skipblanks=true, skipstart=14)

    if data_type == :IAU1980
        parse_iers_eop_iau_1980(eop)
    elseif data_type == :IAU2000A
        parse_iers_eop_iau_2000A(eop)
    else
        throw(ArgumentError("Unknow EOP type. It must be :IAU1980 or :IAU2000A."))
    end
end

################################################################################
#                              Private Functions
################################################################################

function parse_iers_eop_iau_1980(eop::Matrix)
    # Check if the dimension seems correct.
    if size(eop,2) != 16
        error("The input data does not have the correct dimension.")
    end

    # Create the EOP Data structure by creating the interpolations.
    #
    # The interpolation will be linear between two points in the grid. The
    # extrapolation will be flat, considering the nearest point.
	knots = (eop[:,4] .+ 2400000.5,)

    EOPData_IAU1980(
        extrapolate(interpolate(knots, eop[:, 5], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:, 6], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:, 7], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:, 8], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:, 9], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,10], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,11], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,12], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,13], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,14], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,15], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,16], Gridded(Linear())), Flat())
    )
end

function parse_iers_eop_iau_2000A(eop::Matrix)
    # Check if the dimension seems correct.
    if size(eop,2) != 16
        error("The input data does not have the correct dimension.")
    end

    # Create the EOP Data structure by creating the interpolations.
    #
    # The interpolation will be linear between two points in the grid. The
    # extrapolation will be flat, considering the nearest point.
	knots = (eop[:,4] .+ 2400000.5,)

    EOPData_IAU2000A(
        extrapolate(interpolate(knots, eop[:, 5], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:, 6], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:, 7], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:, 8], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:, 9], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,10], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,11], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,12], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,13], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,14], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,15], Gridded(Linear())), Flat()),
        extrapolate(interpolate(knots, eop[:,16], Gridded(Linear())), Flat())
    )
end


text:
using Images, FileIO

function floodfill!(img::Matrix{<:Color}, initnode::CartesianIndex{2}, target::Color, replace::Color)
    img[initnode] != target && return img
    # constants
    north = CartesianIndex(-1,  0)
    south = CartesianIndex( 1,  0)
    east  = CartesianIndex( 0,  1)
    west  = CartesianIndex( 0, -1)

    queue = [initnode]
    while !isempty(queue)
        node = pop!(queue)
        if img[node] == target
            wnode = node
            enode = node + east
        end
        # Move west until color of node does not match target color
        while checkbounds(Bool, img, wnode) && img[wnode] == target
            img[wnode] = replace
            if checkbounds(Bool, img, wnode + north) && img[wnode + north] == target
                push!(queue, wnode + north)
            end
            if checkbounds(Bool, img, wnode + south) && img[wnode + south] == target
                push!(queue, wnode + south)
            end
            wnode += west
        end
        # Move east until color of node does not match target color
        while checkbounds(Bool, img, enode) && img[enode] == target
            img[enode] = replace
            if checkbounds(Bool, img, enode + north) && img[enode + north] == target
                push!(queue, enode + north)
            end
            if checkbounds(Bool, img, enode + south) && img[enode + south] == target
                push!(queue, enode + south)
            end
            enode += east
        end
    end
    return img
end

img = Gray{Bool}.(load("data/unfilledcircle.png"))
floodfill!(img, CartesianIndex(100, 100), Gray(false), Gray(true))
save("data/filledcircle.png", img)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
```
measurement{T<:AbstractFloat}(m::Model990{T}, TTT::Matrix{T}, RRR::Matrix{T},
                              CCC::Matrix{T}; shocks::Bool = true)
```

Assign measurement equation
```
X_t = ZZ*S_t + DD + u_t
```
where
```
u_t = eta_t + MM*eps_t
var(eta_t) = EE
var(u_t) = HH = EE + MM*QQ*MM'
cov(eps_t,u_t) = VV = QQ*MM'
```
"""
function measurement{T<:AbstractFloat}(m::Model990{T},
                                       TTT::Matrix{T},
                                       RRR::Matrix{T},
                                       CCC::Matrix{T};
                                       shocks::Bool = true)
    endo = m.endogenous_states
    exo  = m.exogenous_shocks
    obs  = m.observables

    # If shocks = true, then return measurement equation matrices with rows and columns for
    # anticipated policy shocks
    if shocks
        _n_observables = n_observables(m)
        _n_states = n_states_augmented(m)
        _n_shocks_exogenous = n_shocks_exogenous(m)
        endo_new = m.endogenous_states_augmented
    else
        _n_observables = n_observables(m) - n_anticipated_shocks(m)
        _n_states = n_states_augmented(m) - n_anticipated_shocks(m)
        _n_shocks_exogenous = n_shocks_exogenous(m) - n_anticipated_shocks(m)
        endo_new = Dict(
            [(key,m.endogenous_states_augmented[key] - n_anticipated_shocks(m)) for key in keys(m.endogenous_states_augmented)])
    end

    ZZ = zeros(_n_observables, _n_states)
    DD = zeros(_n_observables, 1)
    MM = zeros(_n_observables, _n_shocks_exogenous)
    EE = zeros(_n_observables, _n_observables)
    QQ = zeros(_n_shocks_exogenous, _n_shocks_exogenous)

    ## Output growth - Quarterly!
    ZZ[obs[:obs_gdp], endo[:y_t]]      = 1.0
    ZZ[obs[:obs_gdp], endo_new[:y_t1]] = -1.0
    ZZ[obs[:obs_gdp], endo[:z_t]]      = 1.0
    DD[obs[:obs_gdp]]                  = 100*(exp(m[:z_star])-1)

    ## Hours growth
    ZZ[obs[:obs_hours], endo[:L_t]] = 1.0
    DD[obs[:obs_hours]]             = m[:Lmean]

    ## Labor Share/real wage growth
    ZZ[obs[:obs_wages], endo[:w_t]]      = 1.0
    ZZ[obs[:obs_wages], endo_new[:w_t1]] = -1.0
    ZZ[obs[:obs_wages], endo[:z_t]]      = 1.0
    DD[obs[:obs_wages]]                  = 100*(exp(m[:z_star])-1)

    ## Inflation (GDP Deflator)
    ZZ[obs[:obs_gdpdeflator], endo[:π_t]]            = m[:Γ_gdpdef]
    ZZ[obs[:obs_gdpdeflator], endo_new[:e_gdpdef_t]] = 1.0
    DD[obs[:obs_gdpdeflator]]                        = 100*(m[:π_star]-1) + m[:δ_gdpdef]

    ## Inflation (Core PCE)
    ZZ[obs[:obs_corepce], endo[:π_t]]             = 1.0
    ZZ[obs[:obs_corepce], endo_new[:e_corepce_t]] = 1.0
    DD[obs[:obs_corepce]]                         = 100*(m[:π_star]-1)

    ## Nominal interest rate
    ZZ[obs[:obs_nominalrate], endo[:R_t]] = 1.0
    DD[obs[:obs_nominalrate]]             = m[:Rstarn]

    ## Consumption Growth
    ZZ[obs[:obs_consumption], endo[:c_t]]      = 1.0
    ZZ[obs[:obs_consumption], endo_new[:c_t1]] = -1.0
    ZZ[obs[:obs_consumption], endo[:z_t]]      = 1.0
    DD[obs[:obs_consumption]]                  = 100*(exp(m[:z_star])-1)

    ## Investment Growth
    ZZ[obs[:obs_investment], endo[:i_t]]      = 1.0
    ZZ[obs[:obs_investment], endo_new[:i_t1]] = -1.0
    ZZ[obs[:obs_investment], endo[:z_t]]      = 1.0
    DD[obs[:obs_investment]]                  = 100*(exp(m[:z_star])-1)

    ## Spreads
    ZZ[obs[:obs_spread], endo[:ERtil_k_t]] = 1.0
    ZZ[obs[:obs_spread], endo[:R_t]]       = -1.0
    DD[obs[:obs_spread]]                   = 100*log(m[:spr])

    ## 10 yrs infl exp
    TTT10                          = (1/40)*((eye(size(TTT, 1)) - TTT)\(TTT - TTT^41))
    ZZ[obs[:obs_longinflation], :] = TTT10[endo[:π_t], :]
    DD[obs[:obs_longinflation]]    = 100*(m[:π_star]-1)

    ## Long Rate
    ZZ[obs[:obs_longrate], :]               = ZZ[6, :]*TTT10
    ZZ[obs[:obs_longrate], endo_new[:lr_t]] = 1.0
    DD[obs[:obs_longrate]]                  = m[:Rstarn]

    ## TFP
    ZZ[obs[:obs_tfp], endo[:z_t]]       = (1-m[:α])*m[:Iendoα] + 1*(1-m[:Iendoα])
    ZZ[obs[:obs_tfp], endo_new[:tfp_t]] = 1.0
    ZZ[obs[:obs_tfp], endo[:u_t]]       = m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα] )
    ZZ[obs[:obs_tfp], endo_new[:u_t1]]  = -(m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα]) )

    QQ[exo[:g_sh], exo[:g_sh]]           = m[:σ_g]^2
    QQ[exo[:b_sh], exo[:b_sh]]           = m[:σ_b]^2
    QQ[exo[:μ_sh], exo[:μ_sh]]           = m[:σ_μ]^2
    QQ[exo[:z_sh], exo[:z_sh]]           = m[:σ_z]^2
    QQ[exo[:λ_f_sh], exo[:λ_f_sh]]       = m[:σ_λ_f]^2
    QQ[exo[:λ_w_sh], exo[:λ_w_sh]]       = m[:σ_λ_w]^2
    QQ[exo[:rm_sh], exo[:rm_sh]]         = m[:σ_r_m]^2
    QQ[exo[:σ_ω_sh], exo[:σ_ω_sh]]       = m[:σ_σ_ω]^2
    QQ[exo[:μ_e_sh], exo[:μ_e_sh]]       = m[:σ_μ_e]^2
    QQ[exo[:γ_sh], exo[:γ_sh]]           = m[:σ_γ]^2
    QQ[exo[:π_star_sh], exo[:π_star_sh]] = m[:σ_π_star]^2
    QQ[exo[:lr_sh], exo[:lr_sh]]         = m[:σ_lr]^2
    QQ[exo[:zp_sh], exo[:zp_sh]]         = m[:σ_z_p]^2
    QQ[exo[:tfp_sh], exo[:tfp_sh]]       = m[:σ_tfp]^2
    QQ[exo[:gdpdef_sh], exo[:gdpdef_sh]] = m[:σ_gdpdef]^2
    QQ[exo[:corepce_sh], exo[:corepce_sh]]       = m[:σ_corepce]^2

    # These lines set the standard deviations for the anticipated shocks. They
    # are here no longer calibrated to the std dev of contemporaneous shocks,
    # as we had in 904
    if shocks
        for i = 1:n_anticipated_shocks(m)
            ZZ[obs[symbol("obs_nominalrate$i")], :]              = ZZ[obs[:obs_nominalrate], :]*(TTT^i)
            DD[obs[symbol("obs_nominalrate$i")]]                 = m[:Rstarn]
            QQ[exo[symbol("rm_shl$i")], exo[symbol("rm_shl$i")]] = m[symbol("σ_r_m$i")]^2
        end
    end

    # Adjustment to DD because measurement equation assumes CCC is the zero vector
    if any(CCC != 0)
        DD += ZZ*((UniformScaling(1) - TTT)\CCC)
    end

    HH    = EE + MM*QQ*MM'
    VV    = QQ*MM'
    VVall = [[RRR*QQ*RRR' RRR*VV];
             [VV'*RRR'    HH]]

    return Measurement(ZZ, DD, QQ, EE, MM, VVall)
end

type Measurement{T<:AbstractFloat}
    ZZ::Matrix{T}
    DD::Matrix{T}
    QQ::Matrix{T}
    EE::Matrix{T}
    MM::Matrix{T}
    VVall::Matrix{T}
end
function Base.getindex(M::Measurement, d::Symbol)
    if d in (:ZZ, :DD, :QQ, :EE, :MM, :VVall)
        return getfield(M, d)
    else
        throw(KeyError(d))
    end
end


text:
__precompile__()
@doc raw"""
    AFD2D_mod

# Finite-diference modeling of acoustic wave equation

```math
\frac{1}{c^2(\mathbf{x}} \frac{\partial P}{\partial t} = \nabla P + \dot{q}(t)
\delta(\mathbf{x} - \mathbf{x}_{s})
```
"""
module AFD2d_mod

#using Yeppp
#using Plots
using VegaLite
using DataFrames
using BenchmarkTools
using DelimitedFiles
using Printf
using Documenter
using StaticArrays

include("SU_mod.jl")
include("KaiserWindow_mod.jl")
include("LanczosInterpolation_mod.jl")
include("InvSmooth2d_mod.jl")
push!(LOAD_PATH, @__DIR__)
#push!(LOAD_PATH, pwd())

#import SU_mod
#import KaiserWindow_mod
#import LanczosInterpolation_mod
#import InvSmooth2d_mod

using .SU_mod
using .KaiserWindow_mod
using .LanczosInterpolation_mod
using .InvSmooth2d_mod
#using Profile

export AFD2d, Model2d, afd_mod

const Rabc=1.e-3
const nabc=49
const CFL=0.5
const ndrv1=2
const ndrv2=7
const dtsampling=0.004

const drv1=[ 9.0/8.0,
            -1.0/24.0]
#
const drv2=[-3.171097979234166,
      1.883703563645937,
     -0.392013173326410,
      0.126419594523541,
     -0.043436978346425,
      0.013633723622292,
     -0.002757740501852]
#
const KAISERWIND=4
#
@doc raw"""
    Grid2d(n1::Int64,
           n2::Int64,
           o1::Float64,
           o2::Float64,
           d1::Float64
           d2::Float64)

    Type to define a regular grid for Finite-Difference

##    Constructor Signature:

* n1: number of grid point along dimension 1
* n2: number of grid point along dimension 2
* o1: x1-coordinate of grid origin
* o2: x2-coordinate of grid origin
* d1: sampling interval along x1-coordinate
* d2: sampling interval along x2-coordinate
"""
mutable struct Grid2d
n1::Int64
n2::Int64
o1::Float64
o2::Float64
d1::Float64
d2::Float64
end
#

mutable struct Model2d
grd::Grid2d
n1fd::Int64
n2fd::Int64
vel::Array{Float64,2}
@doc raw"""
    Model2d(file::String)

    Type Model2d stores the 2D velocity model and
    extends its grid boundaries for FD computations

## Constructor signature:

* file: velocity model header file
"""
function Model2d(file::String)
        hdr=readdlm(file,' ')
        n1::Int64=hdr[1,1]
        n2::Int64=hdr[1,2]
        o1::Float64=hdr[1,3]
        o2::Float64=hdr[1,4]
        d1::Float64=hdr[1,5]
        d2::Float64=hdr[1,6]
        binfile::String=hdr[2,1]

        grd::Grid2d=Grid2d(n1,n2,o1,o2,d1,d2)

        n1fd::Int64=n1+2*nabc
        n2fd::Int64=n2+2*nabc
        vel=zeros(Float64,n1fd,n2fd)
        value::Float32=Float32(0.0)
        iob::IOStream=open(binfile,"r")
        @assert isopen(iob)
        for i2=1:n2
            for i1=1:n1
                    vel[nabc+i1,nabc+i2]=Float64(read(iob,Float32))
            end
        end
        close(iob)

# left boundary
        for i2=1:nabc
        for i1=nabc+1:nabc+n1
            vel[i1,i2]=vel[i1,nabc+1]
        end
        end
# right boundary
        for i2 in nabc+n2+1:n2fd
        for i1 in nabc+1:nabc+n1
            vel[i1,i2]=vel[i1,nabc+n2]
        end
        end

# top boundary
        for i2 in 1:n2fd
        for i1 in 1:nabc
            vel[i1,i2]=vel[nabc+1,i2]
        end
# bottom boundary
        for i1 in nabc+n1+1:n1fd
            vel[i1,i2]=vel[nabc+n1,i2]
        end
        end

        new(grd,n1fd,n2fd,vel)
end

function Model2d(n1::Int64,
                 n2::Int64,
                 o1::Float64,
                 o2::Float64,
                 d1::Float64,
                 d2::Float64
                )
grd::Grid2d=Grid2d(n1,n2,o1,o2,d1,d2)
n1fd::Int64=n1+2*nabc
n2fd::Int64=n2+2*nabc
vel=zeros(Float64,n1fd,n2fd)
new(grd,n1fd,n2fd,vel)
end
end

mutable struct AFD2dGeom
grd::Grid2d
n::Int64
i1::Array{Int64,1}
i2::Array{Int64,1}
wind1::Array{Float64,2}
wind2::Array{Float64,2}

function AFD2dGeom(grd::Grid2d,
                  x1stat::Array{Float64,1},x2stat::Array{Float64,1})
@assert length(x1stat) == length(x2stat)
# Kaiser window parameters
bKaiser::Float64=4.14
lKaiser::Float64=Float64(KAISERWIND)

nstat::Int64 = length(x1stat)
i1::Array{Int64,1}=zeros(Int64,nstat)
i2::Array{Int64,1}=zeros(Int64,nstat)
wind1::Array{Float64,2}=zeros(Float64,2*KAISERWIND+1,nstat)
wind2::Array{Float64,2}=zeros(Float64,2*KAISERWIND+1,nstat)
for is=1:nstat
    s1::Float64=(x1stat[is]-grd.o1)/grd.d1+Float64(nabc)
    s2::Float64=(x2stat[is]-grd.o2)/grd.d2+Float64(nabc)
    @assert s1 > 0.0 && s2 > 0.0
    i1[is]=floor(Int64,s1)
    i2[is]=floor(Int64,s2)

#=
Kaiser window coefficients
for each station
=#
    dx1::Float64=s1-Float64(i1[is])
    @simd for i=-KAISERWIND:KAISERWIND
          @inbounds  wind1[i+KAISERWIND+1,is]=KaiserWindow(Float64(i)+dx1,lKaiser,bKaiser)
    end
    dx2::Float64=s2-Float64(i2[is])
    @simd for i=-KAISERWIND:KAISERWIND
           @inbounds wind2[i+KAISERWIND+1,is]=KaiserWindow(Float64(i)+dx2,lKaiser,bKaiser)
    end
#= normalization of Kaiser window =#
    wsum::Float64=0.0
    for j=-KAISERWIND:KAISERWIND
    for i=-KAISERWIND:KAISERWIND
       @inbounds wsum += wind1[i+KAISERWIND+1,is]*wind2[j+KAISERWIND+1,is]
    end
    end
    wsum=1.0/sqrt(wsum)
    @inbounds wind1[:,is] *= wsum
    @inbounds wind2[:,is] *= wsum
end
new(grd,nstat,i1,i2,wind1,wind2)
end
end

# data structure for acoustic finite difference
mutable struct AFD2d
mdl::Model2d
freq::Float64
dt::Float64
h1::Float64
h2::Float64
abc1::Array{Float64,2}
abc2::Array{Float64,2}
abch1::Array{Float64,2}
abch2::Array{Float64,2}

function AFD2d(m::Model2d,freq::Float64)
vmin::Float64=m.vel[1,1]
vmax::Float64=m.vel[1,1]
@inbounds for i2 in 1:m.n2fd
@inbounds for i1 in 1:m.n1fd
    vmin=min(vmin,m.vel[i1,i2])
    vmax=max(vmax,m.vel[i1,i2])
end
end
# check dispersion
d::Float64=vmin/(6.0*freq)

if  min(m.grd.d1,m.grd.d2) > d
    println("vmin = ",vmin)
    println("d    = ",d)
    println("d1   = ",m.grd.d1)
    exit()
end
# check stability
arg::Float64=pi/sqrt(2.0)
cfl::Float64=CFL
sum::Float64=drv2[1]
for ic in 2:ndrv2
   @inbounds sum += 2.0*drv2[ic]*cos(Float64(ic)*arg)
end
mu::Float64 = cfl*sqrt(2.0/abs(sum))
dt::Float64 = Float64(floor(Int64,10000.0*mu*min(m.grd.d1,m.grd.d2)/vmax)/10000)

println("dt = ",dt)
ndt = floor(Int64,dtsampling/dt)
if  ndt > 1
    dt= dtsampling/Float64(ndt)
    ndt+=1
else
    dt=dtsampling
    ndt=1
end

println("d1 = ",m.grd.d1)
println("d2 = ",m.grd.d2)
println("dt = ",dt)
println("vmin = ",vmin)
println("vmax = ",vmax)
# fill up structure

mu = (dt*dt) #/(m.grd.d1*m.grd.d2)
@simd  for i2 in 1:m.n2fd
@simd for i1 in 1:m.n1fd
        @inbounds m.vel[i1,i2] *= mu*m.vel[i1,i2]
end
end

h1::Float64=1.0/m.grd.d1
h2::Float64=1.0/m.grd.d2

# CPML absorbing boundary conditions

abc1::Array{Float64,2}=zeros(Float64,2,m.n1fd)
abc2::Array{Float64,2}=zeros(Float64,2,m.n2fd)
abch1::Array{Float64,2}=zeros(Float64,2,m.n1fd)
abch2::Array{Float64,2}=zeros(Float64,2,m.n2fd)

d01::Float64 = -1.5*log(Rabc)*vmax/(Float64(nabc)*m.grd.d1)
for i in 1:nabc
    uabc::Float64=Float64(i)/Float64(nabc)
    alpha::Float64=pi*freq*(1.0-uabc)
    d1::Float64 = d01*(uabc-sin(2.0*pi*uabc)/(2.0*pi))
    aux1::Float64 = alpha+d1
    aux2::Float64 = exp(-aux1*dt)
    abc1[1,nabc+1-i] = aux2
    abc1[2,nabc+1-i] = -d1*(1.0-aux2)/aux1
    abc1[1,m.n1fd-nabc+i] = aux2
    abc1[2,m.n1fd-nabc+i] = -d1*(1.0-aux2)/aux1
end

for i in 1:nabc
    uabc::Float64=(Float64(i-1)+0.5)/Float64(nabc)
    alpha::Float64=pi*freq*(1.0-uabc)
    d1::Float64 = d01*(uabc-sin(2.0*pi*uabc)/(2.0*pi))
    aux1::Float64 = alpha+d1
    aux2::Float64 = exp(-aux1*dt)
    abch1[1,nabc+1-i] = aux2
    abch1[2,nabc+1-i] = -d1*(1.0-aux2)/aux1

    uabc =(Float64(i)+0.5)/Float64(nabc)
    alpha=pi*freq*(1.0-uabc)
    d1   = d01*(uabc-sin(2.0*pi*uabc)/(2.0*pi))
    aux1 = alpha+d1
    aux2 = exp(-aux1*dt)
    abch1[1,m.n1fd-nabc+i] = aux2
    abch1[2,m.n1fd-nabc+i] = -d1*(1.0-aux2)/aux1
end

d02::Float64 = -1.5*log(Rabc)*vmax/(Float64(nabc)*m.grd.d2)
for i in 1:nabc
    uabc::Float64=Float64(i)/Float64(nabc)
    alpha::Float64=pi*freq*(1.0-uabc)
    d2::Float64 = d02*(uabc-sin(2.0*pi*uabc)/(2.0*pi))
    aux1::Float64 = alpha+d2
    aux2::Float64 = exp(-aux1*dt)
    abc2[1,nabc+1-i] = aux2
    abc2[2,nabc+1-i] = -d2*(1.0-aux2)/aux1
    abc2[1,m.n2fd-nabc+i] = aux2
    abc2[2,m.n2fd-nabc+i] = -d2*(1.0-aux2)/aux1
end

for i in 1:nabc
    uabc::Float64=(Float64(i-1)+0.5)/Float64(nabc)
    alpha::Float64=pi*freq*(1.0-uabc)
    d2::Float64 = d02*(uabc-sin(2.0*pi*uabc)/(2.0*pi))
    aux1::Float64 = alpha+d2
    aux2::Float64 = exp(-aux1*dt)
    abch2[1,nabc+1-i] = aux2
    abch2[2,nabc+1-i] = -d2*(1.0-aux2)/aux1

    uabc=(Float64(i)+0.5)/Float64(nabc)
    alpha=pi*freq*(1.0-uabc)
    d2   = d02*(uabc-sin(2.0*pi*uabc)/(2.0*pi))
    aux1 = alpha+d2
    aux2 = exp(-aux1*dt)
    abch2[1,m.n2fd-nabc+i] = aux2
    abch2[2,m.n2fd-nabc+i] = -d2*(1.0-aux2)/aux1
end

println("vmin = ",vmin)
println("vmax = ",vmax)
new(m,freq,dt,h1,h2,abc1,abc2,abch1,abch2)
end

end

function afd_ricker(FD::AFD2d,t::Float64)
arg::Float64=(pi*FD.freq*(t-1.0/FD.freq))^2
pulse::Float64=exp(-arg)*(1.0-2.0*arg)
pulse
end

function afd_update(FD::AFD2d,p1::Array{Float64,2},
                              p2::Array{Float64,2},
                              phi1::Array{Float64,2},
                              phi2::Array{Float64,2},
                              psi1::Array{Float64,2},
                              psi2::Array{Float64,2},
        )
#=

 Finite difference
 time update

=#
    Threads.@threads  for i2 in ndrv2:FD.mdl.n2fd-ndrv2
        @inbounds for i1 in ndrv2:FD.mdl.n1fd-ndrv2
        dp1::Float64 =0.0
        dp2::Float64 =0.0
        for ic in 0:ndrv1-1
                dp1 += drv1[ic+1]*(p2[i1+ic+1,i2]-p2[i1-ic,i2])
        end
        dp1 *= FD.h1
        phi1[i1,i2]=FD.abch1[1,i1]*phi1[i1,i2]+FD.abch1[2,i1]*dp1

        for ic in 0:ndrv1-1
                dp2 += drv1[ic+1]*(p2[i1,i2+ic+1]-p2[i1,i2-ic])
        end
        dp2 *= FD.h2
        phi2[i1,i2]=FD.abch2[1,i2]* phi2[i1,i2]+FD.abch2[2,i2]*dp2
    end
    end

    Threads.@threads  for i2 in ndrv2:FD.mdl.n2fd-ndrv2
        @inbounds for i1 in ndrv2:FD.mdl.n1fd-ndrv2

        dphi1::Float64 =0.0
        dphi2::Float64 =0.0
        for ic in 0:ndrv1-1
                dphi1 += drv1[ic+1]*(phi1[i1+ic,i2]-phi1[i1-ic-1,i2])
        end
        dphi1 *= FD.h1

        for ic in 0:ndrv1-1
                dphi2 += drv1[ic+1]*(phi2[i1,i2+ic]-phi2[i1,i2-ic-1])
        end
        dphi2 *= FD.h2

        laplacian1::Float64=drv2[1]*p2[i1,i2]
        laplacian2::Float64=drv2[1]*p2[i1,i2]
        for ic in 1:ndrv2-1
             @inbounds laplacian1 += drv2[ic+1]*(p2[i1-ic,i2]+p2[i1+ic,i2])
        end
        laplacian1 *= FD.h1*FD.h1
        for ic in 1:ndrv2-1
             @inbounds laplacian2 += drv2[ic+1]*(p2[i1,i2-ic]+p2[i1,i2+ic])
        end
        laplacian2 *= FD.h2*FD.h2
        psi1[i1,i2]=FD.abc1[1,i1]* psi1[i1,i2]+FD.abc1[2,i1]*(laplacian1+dphi1)
        psi2[i1,i2]=FD.abc2[1,i2]* psi2[i1,i2]+FD.abc2[2,i2]*(laplacian2+dphi2)

        laplacian::Float64=laplacian1+laplacian2+dphi1+dphi2+psi1[i1,i2]+psi2[i1,i2]
        p1[i1,i2]=2.0*p2[i1,i2]-p1[i1,i2]+FD.mdl.vel[i1,i2]*laplacian
    end
    end

    Threads.@threads for i2 in 1:FD.mdl.n2fd
    @inbounds for i1 in 1:FD.mdl.n1fd
        swap::Float64=p1[i1,i2]
        p1[i1,i2]=p2[i1,i2]
        p2[i1,i2]=swap
    end
    end

    # Dirichlet boundary to avoid stability
    Threads.@threads for i1 in 1:FD.mdl.n1fd
            p2[i1,ndrv2]=0.0
            p2[i1,FD.mdl.n2fd-ndrv2]=0.0
    end

    Threads.@threads for i2 in 1:FD.mdl.n2fd
            p2[ndrv2,i2]=0.0
            p2[FD.mdl.n1fd-ndrv2,i2]=0.0
    end

end

function afd_src_injection(FD::AFD2d,Src::AFD2dGeom,t::Float64,p1::Array{Float64,2})
@inbounds  Threads.@threads  for is=1:Src.n
@inbounds        for iw2=-KAISERWIND:KAISERWIND
                     i2s::Int64=min(max(Src.i2[is]+iw2,1),FD.mdl.n2fd)
                     wsrc2::Float64=Src.wind2[iw2+KAISERWIND+1]
@inbounds        for iw1=-KAISERWIND:KAISERWIND
                     i1s::Int64=min(max(Src.i1[is]+iw1,1),FD.mdl.n1fd)
                     wsrc::Float64=Src.wind1[iw1+KAISERWIND+1]*wsrc2

                    p1[i1s,i2s] -= wsrc*FD.mdl.vel[i1s,i2s]*(afd_ricker(FD,t+FD.dt)+
                         afd_ricker(FD,t-FD.dt)-
                         2.0*afd_ricker(FD,t))/(FD.dt*FD.dt)
        end
        end
    end
    nothing
end

function afd_rcv_field(FD::AFD2d,Rcv::AFD2dGeom,is::Int64,p2::Array{Float64,2})
                 field::Float64 = 0.0
@inbounds        for iw2=-KAISERWIND:KAISERWIND
                     i2s::Int64=min(max(Rcv.i2[is]+iw2,1),FD.mdl.n2fd)
                     wrcv2::Float64=Rcv.wind2[iw2+KAISERWIND+1]
@inbounds        for iw1=-KAISERWIND:KAISERWIND
                     i1s::Int64=min(max(Rcv.i1[is]+iw1,1),FD.mdl.n1fd)
                     wrcv::Float64=Rcv.wind1[iw1+KAISERWIND+1]*wrcv2
                     field += wrcv*p2[i1s,i2s]
        end
        end
        return field
end

@doc raw"""
    afd_mod(FD::AFD2d,
            T::Float64,
            x1s::Array{Float64,1},
            x2s::Array{Float64,1};
            x1r::Array{Float64,1},
            x2r::Array{Float64,1},
            section::Array{SU_trc,1}
            )

    Acoustic wavefield forward modeling

## Signature:

* FD : finite-Difference type
* T  : simulation time
* x1s: source stations x1-coordinates
* x2s: source stations x2-coordinates
* x1r: receiver stations x1-coordinates
* x2r: receiver stations x2-coordinates
"""
function afd_mod!(FD::AFD2d,adj::Bool,
                 T::Float64,
                 x1s::Array{Float64,1},
                 x2s::Array{Float64,1};
                 x1r::Array{Float64,1},
                 x2r::Array{Float64,1},
                 section::Array{SU_trc,1}
                )

nt::Int64  = floor(Int64,T/FD.dt+1.0)
ndt::Int64 = floor(Int64,dtsampling/FD.dt)+1

#= set source stations =#
Src::AFD2dGeom=AFD2dGeom(FD.mdl.grd,x1s,x2s)
#= set receiver stations =#
Rcv::AFD2dGeom=AFD2dGeom(FD.mdl.grd,x1r,x2r)

ns::Int16=floor(Int16,T/dtsampling+1.0)

scale::Float64=10.0
scalco::Int16=Int16(scale)
scalel::Int16=Int16(scale)
@inbounds @simd for ir=1:length(x1r)
    section[ir].hdr.scalco=-floor(Int16,scalco)
    section[ir].hdr.scalel=-floor(Int16,scalel)
    section[ir].hdr.tracl=Int32(ir)
    section[ir].hdr.ns=ns
    section[ir].hdr.dt=Int16(1000000.0*dtsampling)
    section[ir].hdr.sx=Int32(x2s[1]*scale)
    section[ir].hdr.sdepth=Int32(x1s[1]*scale)
    section[ir].hdr.gx=Int32(x2r[ir]*scale)
    section[ir].hdr.gelev=Int32(-x1r[ir]*scale)
    section[ir].hdr.ntr=Int32(length(x1r))
end

irec::Int64 = 0
trec::Float64=0.0
dtrec::Float64=Float64(dtsampling)

# state field allocation and initialization
p1::Array{Float64,2}=zeros(Float64,FD.mdl.n1fd,FD.mdl.n2fd)
p2::Array{Float64,2}=zeros(Float64,FD.mdl.n1fd,FD.mdl.n2fd)

# memory field for ABC allocation and initialization
phi1::Array{Float64,2}=zeros(Float64,FD.mdl.n1fd,FD.mdl.n2fd)
phi2::Array{Float64,2}=zeros(Float64,FD.mdl.n1fd,FD.mdl.n2fd)
psi1::Array{Float64,2}=zeros(Float64,FD.mdl.n1fd,FD.mdl.n2fd)
psi2::Array{Float64,2}=zeros(Float64,FD.mdl.n1fd,FD.mdl.n2fd)

zm::Array{Float64,1}=zeros(Float64,FD.mdl.n1fd)
xm::Array{Float64,1}=zeros(Float64,FD.mdl.n2fd)

@inbounds for i1=1:FD.mdl.n1fd
        zm[i1]=FD.mdl.grd.o1 + Float64(i1-1)*FD.mdl.grd.d1
end
@inbounds for i2=1:FD.mdl.n2fd
        xm[i2]=FD.mdl.grd.o2 + Float64(i2-1)*FD.mdl.grd.d2
end

ndtsampling::Int64 = floor(Int64,dtsampling/FD.dt)+1

isnap::Int64=0
tsnap::Float64 = 0.0
dtsnap::Float64 = 50*dtsampling

        iob=open("fdmovie.bin","w")
        @assert isopen(iob)

if  adj == 0
@inbounds for it in 1:nt

    t::Float64=Float64(it-1)*FD.dt

    # source injection:

    afd_src_injection(FD,Src,t,p1)

    afd_update(FD,p1,p2,phi1,phi2,psi1,psi2)

    if abs(t-tsnap) < 0.5*FD.dt
       isnap += 1
       tsnap += dtsnap
       time::String=@sprintf("%.3f",Float64(isnap-1)*dtsnap)
       #=
      img=heatmap(xm,zm,FD.mdl.vel*maximum(p2)/maximum(FD.mdl.vel),
                  reuse=false,xlabel="x (m)",ylabel="y (m)",
                  yflip=true,colorbar=false)

      heatmap!(xm,zm,p2,xlabel="x (m)",ylabel="y
                       (m)",alpha=0.25,title="Snapshot at "*time*" s",yflip=true)
              #display(plot(img))
        =#

        for i2=1:size(p2,2)
            for i1=1:size(p2,1)
                write(iob,Float32(p2[i1,i2]))
            end
        end
        #=
        data=DataFrame(p2,:auto)
        select(data,[:x1,:x2]) |> @vlplot({:rect,color=data},x=:x1,y=:x2) |> display
        =#
    end

    if  abs(t-trec) < FD.dt
        irec += 1
        trec += dtrec
        @inbounds Threads.@threads for ir=1:length(x1r)
                section[ir].dat[irec] = Float32(afd_rcv_field(FD,Rcv,ir,p2))
        end
    end

end
        close(iob)

else
@inbounds for it in 1:nt

    t::Float64=Float64(it-1)*FD.dt

    # source injection:

    afd_src_injection(FD,Src,t,p1)

    afd_update(FD,p1,p2,phi1,phi2,psi1,psi2)

    if abs(t-tsnap) < 0.5*FD.dt
       isnap += 1
       tsnap += dtsnap
       time::String=@sprintf("%.3f",Float64(isnap-1)*dtsnap)
      img=heatmap(xm,zm,FD.mdl.vel*maximum(p2)/maximum(FD.mdl.vel),
                  reuse=false,xlabel="x (m)",ylabel="y (m)",
                  yflip=true,colorbar=false)

      heatmap!(xm,zm,p2,xlabel="x (m)",ylabel="y
                       (m)",alpha=0.25,title="Snapshot at "*time*" s",yflip=true)
           display(plot(img))

              for i2=1:size(p2,2)
                for i1=1:size(p2,1)
                write(iob,Float32(p2[i1,i2]))
            end
        end
    end

    if  abs(t-trec) < FD.dt
        irec += 1
        trec += dtrec
        @inbounds Threads.@threads for ir=1:length(x1r)
                section[ir].dat[irec] = Float32(afd_rcv_field(FD,Rcv,ir,p2))
        end
    end

end
        close(iob)

end
end

#function __init__()
function Main()

freq::Float64=20.0
T::Float64=3.0
adj::Bool=0

#mdl::Model2d=Model2d("marm_ii_vp_10m.hdr")
#mdl::Model2d=Model2d("velstep.hdr")
mdl::Model2d=Model2d("velbd.hdr")
println("n1, n2 ",mdl.n1fd," ",mdl.n2fd)
#=
vels::Array{Float64,2}=zeros(mdl.n1fd,mdl.n2fd)
println("n1, n2 ",size(vels,1)," ",size(vels,2))

Smooth2d!(mdl.vel,vels,0.1)
        binfile="vels.bin"
        iob=open(binfile,"w")
        @assert isopen(iob)
        for i2=1:size(vels,2)
                for i1=1:size(vels,1)
                write(iob,Float32(vels[i1,i2]))
            end
        end
        close(iob)
=#

AFD::AFD2d=AFD2d(mdl,freq)
x1s::Array{Float64,1}=zeros(Float64,1)
x2s::Array{Float64,1}=zeros(Float64,1)
x1r::Array{Float64,1}=zeros(Float64,AFD.mdl.grd.n2)
x2r::Array{Float64,1}=zeros(Float64,AFD.mdl.grd.n2)

x1s[1]=200*AFD.mdl.grd.d1 + AFD.mdl.grd.o1
x2s[1]=Float64(AFD.mdl.grd.n2-1)*AFD.mdl.grd.d2/2.0+AFD.mdl.grd.o2

@simd for ir=1:length(x1r)
     @inbounds x1r[ir] = mdl.grd.o1 + 120*mdl.grd.d1
end

@simd for ir=1:length(x2r)
    @inbounds x2r[ir] = mdl.grd.o2 + Float64(ir-1)*mdl.grd.d2
end

ns::Int64=floor(Int64,T/dtsampling)+1
seis::Array{SU_trc,1}=Array{SU_trc,1}(undef,mdl.grd.n2)
Threads.@threads for itrc=1:length(x1r)
    @inbounds seis[itrc]=SU_trc(ns)
end

println("x1s = ",x1s[1])
println("x2s = ",x2s[1])

in::IOStream=open("survey.su","r")
try
    isopen(in)
catch
   println("Error opening file survey.su file")
end
sutrace::SU_trc=SU_trc(ns)
SUGetTrace(in,sutrace)
close(in)

#=
#@assert isopen(io)
println("ntr = ",length(x1r))
println("ns  = ",ns)
@inbounds for itrc=1:length(x1r)
        seis[itrc].hdr.tracl=Int32(itrc)
        seis[itrc].hdr.ns=Int16(ns)
        seis[itrc].hdr.dt=Int16(2000)
        SUPutTrace(io,seis[itrc])
        #=
        @inbounds for is=1:ns
                write(io,Float32(seis[itrc].dat[is]))
        end
        =#
end

=#

@time afd_mod!(AFD,adj,T,x1s,x2s,x1r=x1r,x2r=x2r,section=seis)

io::IOStream=open("seis.su","w")
try
    isopen(io)
catch
   println("Error opening file seis.su file")
end
#@assert isopen(io)
println("ntr = ",length(x1r))
println("ns  = ",ns)
@inbounds for itrc=1:length(x1r)
        seis[itrc].hdr.tracl=Int32(itrc)
        seis[itrc].hdr.ns=Int16(ns)
        seis[itrc].hdr.dt=Int16(2000)
        SUPutTrace(io,seis[itrc])
        #=
        @inbounds for is=1:ns
                write(io,Float32(seis[itrc].dat[is]))
        end
        =#
end
close(io)
nothing
end
end


text:
function ComputeFCFV(mesh, se, Tdir, tau)
    ae = zeros(mesh.nel)
    be = zeros(mesh.nel)
    ze = zeros(mesh.nel,2)

    # Assemble FCFV elements
    @avx for iel=1:mesh.nel

        be[iel] = be[iel]   + mesh.vole[iel]*se[iel]

        for ifac=1:mesh.nf_el

            nodei = mesh.e2f[iel,ifac]
            bc    = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau*mesh.ke[iel],dAi,mesh.vole[iel],mesh.type)                              # Stabilisation parameter for the face

            # Assemble
            ze[iel,1] += (bc==1) * dAi*ni_x*Tdir[nodei]  # Dirichlet
            ze[iel,2] += (bc==1) * dAi*ni_y*Tdir[nodei]  # Dirichlet
            be[iel]   += (bc==1) * dAi*taui*Tdir[nodei]  # Dirichlet
            ae[iel]   +=           dAi*taui

        end
    end
    return ae, be, ze
end

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

function ComputeElementValues(mesh, uh, ae, be, ze, Tdir, tau)
    ue          = zeros(mesh.nel);
    qx          = zeros(mesh.nel);
    qy          = zeros(mesh.nel);

    @avx for iel=1:mesh.nel

        ue[iel]  =  be[iel]/ae[iel]
        qx[iel]  = -mesh.ke[iel]/mesh.vole[iel]*ze[iel,1]
        qy[iel]  = -mesh.ke[iel]/mesh.vole[iel]*ze[iel,2]

        for ifac=1:mesh.nf_el

            # Face
            nodei = mesh.e2f[iel,ifac]
            bc    = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau*mesh.ke[iel],dAi,mesh.vole[iel],mesh.type)      # Stabilisation parameter for the face

            # Assemble
            ue[iel] += (bc!=1) *  dAi*taui*uh[mesh.e2f[iel, ifac]]/ae[iel]
            qx[iel] -= (bc!=1) *  mesh.ke[iel]/mesh.vole[iel]*dAi*ni_x*uh[mesh.e2f[iel, ifac]]
            qy[iel] -= (bc!=1) *  mesh.ke[iel]/mesh.vole[iel]*dAi*ni_y*uh[mesh.e2f[iel, ifac]]
         end
    end
    return ue, qx, qy
end

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

function ElementAssemblyLoop(mesh, ae, be, ze, Tdir, Tneu, tau)
    # Assemble element matrices and rhs
    f    = zeros(mesh.nf);
    Kv   = zeros(mesh.nel, mesh.nf_el, mesh.nf_el)
    fv   = zeros(mesh.nel, mesh.nf_el)

    @avx for iel=1:mesh.nel

        for ifac=1:mesh.nf_el

            nodei = mesh.e2f[iel,ifac]
            bci   = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau*mesh.ke[iel],dAi,mesh.vole[iel], mesh.type)

            for jfac=1:mesh.nf_el

                nodej = mesh.e2f[iel,jfac]
                bcj   = mesh.bc[nodej]
                dAj   = mesh.dA[iel,jfac]
                nj_x  = mesh.n_x[iel,jfac]
                nj_y  = mesh.n_y[iel,jfac]
                tauj  = StabParam(tau*mesh.ke[iel],dAj,mesh.vole[iel], mesh.type)

                # Delta
                del = 0.0 + (ifac==jfac)*1.0

                # Element matrix
                nitnj = ni_x*nj_x + ni_y*nj_y;
                Ke_ij =-dAi * (1.0/ae[iel] * dAj * taui*tauj - mesh.ke[iel]/mesh.vole[iel]*dAj*nitnj - taui*del);
                yes   = (bci!=1) & (bcj!=1)
                Kv[iel,ifac,jfac] = yes * Ke_ij
            end
            # RHS
            Xi     = 0.0 + (bci==2)*1.0;
            ti     = Tneu[nodei]
            nitze  = ni_x*ze[iel,1] + ni_y*ze[iel,2]
            fe_i   = (bci!=1) * -dAi * (mesh.ke[iel]/mesh.vole[iel]*nitze - ti*Xi - 1.0/ae[iel]*be[iel]*taui)
            fv[iel,ifac] += (bci!=1) * fe_i
            # Dirichlet nodes
            Kv[iel,ifac,ifac] += (bci==1) * 1.0
            fv[iel,ifac]      += (bci==1) * Tdir[nodei]
        end
    end
    return Kv, fv
end

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

function CreateTripletsSparse(mesh, Kv, fv)
    # Create triplets and assemble sparse matrix
    idof = 1:mesh.nf_el
    ii   = repeat(idof, 1, length(idof))'
    ij   = repeat(idof, 1, length(idof))
    Ki   = mesh.e2f[:,ii]
    Kif  = mesh.e2f[:,ii[1,:]]
    Kj   = mesh.e2f[:,ij]
    K    = sparse(Ki[:], Kj[:], Kv[:], mesh.nf, mesh.nf)
    f    = sparse(Kif[:], ones(size(Kif[:])), fv[:], mesh.nf, 1)
    f    = Array(f)
    droptol!(K, 1e-6)
    return K, f
end

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


text:
using DifferentiableStateSpaceModels, BenchmarkTools
using Test, Parameters, LinearAlgebra, Random

using DifferentiableStateSpaceModels
using Test, LinearAlgebra

# The BLAS threads is still an issue in Julia 1.7
# This has no effect with MKL
DifferentiableStateSpaceModels.set_blas_threads()

println("Running Testsuite with Threads.nthreads() = $(Threads.nthreads()) BLAS.vendor = $(BLAS.vendor()), and BLAS.num_threads = $(BLAS.get_num_threads()) \n")

# Setting miniumum number of evalations to avoid compilation
BenchmarkTools.DEFAULT_PARAMETERS.evals = 3  # at least 3

# Benchmark groups
const SUITE = BenchmarkGroup()
SUITE["utilities"] = include(pkgdir(DifferentiableStateSpaceModels) *
                             "/benchmark/utilities.jl")
SUITE["rbc"] = include(pkgdir(DifferentiableStateSpaceModels) * "/benchmark/rbc.jl")
SUITE["sgu"] = include(pkgdir(DifferentiableStateSpaceModels) * "/benchmark/sgu.jl")
SUITE["fvgq"] = include(pkgdir(DifferentiableStateSpaceModels) * "/benchmark/fvgq.jl")
#SUITE["rbc_likelihoods"] =
# include(pkgdir(DifferentiableStateSpaceModels) * "/benchmark/rbc_likelihoods.jl")

# See README.md
# Can execute with:
# results = run(SUITE; verbose = true)

#Then extract with things like:
# median(results["rbc_likelihoods"])
# results["rbc_likelihoods"]["joint"]["likelihood"]
# Or save stuff, make a change, and then use `judge(m1, m2)`