text:
using Test
using ScottishTaxBenefitModel
using .ModelHousehold:
    Household,
    Person,
    People_Dict,
    PeopleArray,
    default_bu_allocation,
    get_benefit_units,
    get_head,
    get_spouse,
    num_people,
    make_benefit_unit

using .ExampleHouseholdGetter
using .Definitions
using Dates: Date
using .IncomeTaxCalculations:
    old_enough_for_mca,
    apply_allowance,
    calc_income_tax!

using .STBParameters: IncomeTaxSys
using .Results:
    init_household_result,
    IndividualResult,
    BenefitUnitResult,
    ITResult,
    NIResult,
    map_incomes
using .STBIncomes
using .ExampleHelpers
using .RunSettings

function get_tax(; scotland = false ) :: IncomeTaxSys
    it = get_default_it_system( year=2019, scotland=scotland, weekly=false )
    it.non_savings_rates ./= 100.0
    it.savings_rates ./= 100.0
    it.dividend_rates ./= 100.0
    it.personal_allowance_withdrawal_rate /= 100.0
    it.mca_credit_rate /= 100.0
    it.mca_withdrawal_rate /= 100.0
    it.pension_contrib_withdrawal_rate /= 100.0
    it
end

@testset "Melville 2019 ch2 examples 1; basic calc Scotland vs RUK" begin
    # BASIC IT Calcaulation on
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )

    @time names = ExampleHouseholdGetter.initialise( Settings())
    income = [11_730,14_493,30_000,33_150.0,58_600,231_400]
    ntests = size(income)[1]
    for i in 1:ntests-1
        income[i] += itsys_scot.personal_allowance # weird way this is expessed in Melville
    end
    taxes_ruk = [2_346.0,2898.60,6_000,6_630.0,15_940.00,89_130.0]
    taxes_scotland = [2_325.51,2_898.60,6_155.07, 7260.57,17_695.07,92_613.07]
    @test size( income ) == size( taxes_ruk) == size( taxes_scotland )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    scottish = ExampleHouseholdGetter.get_household( "mel_c2_scot" )

    bus = default_bu_allocation( scottish )
    nbus = size(bus)[1]
    println( bus )
    # @test nbus == 1 == size( bus[1])[1]
    pers = bus[1][1]
    for i in 1:ntests
        prsc = IndividualResult{Float64}()
        prsc.income[WAGES] = income[i]
        calc_income_tax!( prsc, pers, itsys_scot )
        println( "Scotland $i : calculated $(prsc.income[INCOME_TAX]) expected $(taxes_scotland[i])")
        @test prsc.income[INCOME_TAX] ≈ taxes_scotland[i]
        pruk = IndividualResult{Float64}()
        pruk.income[WAGES] = income[i]

        calc_income_tax!( pruk, pers, itsys_ruk )
        println( "rUK $i : calculated $(prsc.income[INCOME_TAX]) expected $(taxes_ruk[i])")
        @test pruk.income[INCOME_TAX]  ≈ taxes_ruk[i]
        println( ruk.people[RUK_PERSON].income )
    end
end # example 1

@testset "Melville 2019 ch2 example 2; personal savings allowance" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    scottish = ExampleHouseholdGetter.get_household( "mel_c2_scot" )
    income = [20_000,37_501,64_000,375_000.0]
    psa = [1_000.0,500.0,500.0,0.0]
    @test size( income ) == size( psa ) # same Scotland and RUK
    pers = ruk.people[RUK_PERSON] # doesn't matter S/RUK

    for i in size(income)[1]
        pruk = IndividualResult{Float64}()
        pruk.income[WAGES] = income[i]
        println( "case $i income = $(income[i])")
        calc_income_tax!( pruk, pers, itsys_ruk )
        @test pruk.it.personal_savings_allowance == psa[i]
    end
end # example 2

@testset "ch2 example 3; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    tax_due_scotland = 5680.07
    pruk = IndividualResult{Float64}()
    pruk.income[SELF_EMPLOYMENT_INCOME] = 40_000.00
    pruk.income[BANK_INTEREST] = 1_250.00
    prsc = IndividualResult{Float64}()
    prsc.income[SELF_EMPLOYMENT_INCOME] = 40_000.00
    prsc.income[BANK_INTEREST] = 1_250.00

    calc_income_tax!( prsc, pers, itsys_scot )
    #
    @test prsc.income[INCOME_TAX] ≈ tax_due_scotland

    pruk = IndividualResult{Float64}()
    pruk.income[SELF_EMPLOYMENT_INCOME] = 40_000.00
    pruk.income[BANK_INTEREST] = 1_250.00
    calc_income_tax!( pruk, pers, itsys_ruk )
    tax_due_ruk = 5_550.00
    @test pruk.income[INCOME_TAX] ≈ tax_due_ruk
end # example 3

@testset "ch2 example 4; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )

    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    tax_due_ruk = 840.00
    tax_due_scotland = 819.51
    prsc = IndividualResult{Float64}()
    prsc.income[PROPERTY] = 16_700.00
    prsc.income[BANK_INTEREST] = 1_100.00
    calc_income_tax!( prsc, pers, itsys_scot )
    @test prsc.income[INCOME_TAX] ≈ tax_due_scotland
    pruk = IndividualResult{Float64}()
    pruk.income[PROPERTY] = 16_700.00
    pruk.income[BANK_INTEREST] = 1_100.00

    calc_income_tax!( pruk, pers, itsys_ruk )
    @test pruk.income[INCOME_TAX] ≈ tax_due_ruk
end # example 4

@testset "ch2 example 5; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    tax_due_ruk = 11_232.00
    tax_due_scotland = 12_864.57
    prsc = IndividualResult{Float64}()
    prsc.income[SELF_EMPLOYMENT_INCOME] = 58_850.00
    prsc.income[BANK_INTEREST] = 980.00
    calc_income_tax!( prsc, pers, itsys_scot )
    @test prsc.income[INCOME_TAX] ≈ tax_due_scotland
    pruk = IndividualResult{Float64}()
    pruk.income[SELF_EMPLOYMENT_INCOME] = 58_850.00
    pruk.income[BANK_INTEREST] = 980.00

    calc_income_tax!( pruk, pers, itsys_ruk )
    @test pruk.income[INCOME_TAX] ≈ tax_due_ruk
end # example 5

@testset "ch2 example 6; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]

    tax_due_ruk = 93_825.75
    tax_due_scotland = 97_397.17
    prsc = IndividualResult{Float64}()
    prsc.income[SELF_EMPLOYMENT_INCOME] = 240_235.00
    prsc.income[BANK_INTEREST] = 1_600.00
    calc_income_tax!( prsc, pers, itsys_scot )
    @test prsc.income[INCOME_TAX]  ≈ tax_due_scotland
    pruk = IndividualResult{Float64}()
    pruk.income[SELF_EMPLOYMENT_INCOME] = 240_235.00
    pruk.income[BANK_INTEREST] = 1_600.00
    calc_income_tax!( pruk, pers, itsys_ruk )
    @test pruk.income[INCOME_TAX]  ≈ tax_due_ruk
end # example 6

@testset "ch2 example 7; savings calc" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[self_employment_income] = 10_000.00
    pers.income[bank_interest] = 3_380.00
    pers.income[other_investment_income] = 36_680.00/0.8 # gross up at basic
    tax_due_ruk = 10_092.00 # inc already deducted at source
    tax_due_scotland = 10_092.00

    prsc = IndividualResult{Float64}()
    prsc.income = map_incomes( pers )
    calc_income_tax!( prsc, pers, itsys_scot )
    @test prsc.income[INCOME_TAX] ≈ tax_due_ruk

    pruk = IndividualResult{Float64}()
    pruk.income = map_incomes( pers )
    calc_income_tax!( pruk, pers, itsys_ruk )
    @test pruk.income[INCOME_TAX] ≈ tax_due_ruk

end # example 7

#
# stocks_shares
#

@testset "ch2 example 8; simple stocks_shares" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[property] = 28_590.00
    pers.income[bank_interest] = 1_050.00
    pers.income[stocks_shares] = 204_100.0 # gross up at basic
    tax_due_ruk = 74_834.94 # inc already deducted at source
    tax_due_scotland = 74_834.94+140.97

    prsc = IndividualResult{Float64}()
    prsc.income = map_incomes( pers )
    calc_income_tax!( prsc, pers, itsys_scot )
    @test prsc.income[INCOME_TAX] ≈ tax_due_scotland

    pruk = IndividualResult{Float64}()
    pruk.income = map_incomes( pers )
    calc_income_tax!( pruk, pers, itsys_ruk )
    @test pruk.income[INCOME_TAX] ≈ tax_due_ruk

end # example 8

@testset "ch2 example 9; simple stocks_shares" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[private_pensions] = 17_750.00
    pers.income[bank_interest] = 195.00
    pers.income[stocks_shares] = 1_600.0 # gross up at basic
    tax_due_ruk = 1_050.00 # inc already deducted at source
    tax_due_scotland = 1_050.00-20.49
    prsc = IndividualResult{Float64}()
    prsc.income = map_incomes( pers )
    calc_income_tax!( prsc, pers, itsys_scot )
    @test prsc.income[INCOME_TAX] ≈ tax_due_scotland

    pruk = IndividualResult{Float64}()
    pruk.income = map_incomes( pers )
    calc_income_tax!( pruk, pers, itsys_ruk )
    @test pruk.income[INCOME_TAX] ≈ tax_due_ruk
end # example 9

@testset "ch3 personal allowances ex 1 - hr allowance withdrawal" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )

    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.income[self_employment_income] = 110_520.00
    tax_due_ruk = 33_812.00
    pruk = IndividualResult{Float64}()
    pruk.income = map_incomes( pers )
    calc_income_tax!( pruk, pers, itsys_ruk )

    @test pruk.income[INCOME_TAX] ≈ tax_due_ruk
    pers.income[self_employment_income] += 100.0

    pruk = IndividualResult{Float64}()
    pruk.income = map_incomes( pers )
    calc_income_tax!( pruk, pers, itsys_ruk )
    tax_due_ruk = 33_812.00+60.0
    @test pruk.income[INCOME_TAX] ≈ tax_due_ruk

    # tax_due_scotland = 33_812.00+61.5 ## FIXME actually, check this by hand

end # example1 ch3

@testset "ch3 personal allowances ex 2 - marriage allowance" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )

    names = ExampleHouseholdGetter.initialise( Settings() )
    scot = ExampleHouseholdGetter.get_household( "mel_c2_scot" ) # scots are a married couple
    head = scot.people[SCOT_HEAD]
    spouse = scot.people[SCOT_SPOUSE]
    head.income[self_employment_income] = 11_290.0
    head_tax_due = 0.0
    spouse.income[self_employment_income] = 20_000.0
    spouse_tax_due_ruk = 1_258.0
    bu = make_benefit_unit( PeopleArray([head,spouse]), head.pid, spouse.pid )
    bruk = init_benefit_unit_result( bu )
    calc_income_tax!( bruk, head, spouse, itsys_ruk )
    println( bruk )
    println( "HEAD\n$(inctostr(bruk.pers[head.pid].income ))")
    println( bruk.pers[head.pid].it)
    println( "SPOUSE\n$(inctostr(bruk.pers[spouse.pid].income ))")
    println( bruk.pers[spouse.pid].it)
    @test bruk.pers[spouse.pid].income[INCOME_TAX] ≈ spouse_tax_due_ruk
end # example 2 ch3

@testset "ch3 blind person" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )

    names = ExampleHouseholdGetter.initialise( Settings() )
    ruk = ExampleHouseholdGetter.get_household( "mel_c2" )
    pers = ruk.people[RUK_PERSON]
    pers.registered_blind = true
    bu = make_benefit_unit( PeopleArray([pers]), pers.pid, BigInt(-1) )
    bruk = init_benefit_unit_result( bu )
    result = calc_income_tax!( bruk, pers, nothing, itsys_ruk )
    @test bruk.pers[pers.pid].it.allowance ≈
        itsys_ruk.personal_allowance + itsys_ruk.blind_persons_allowance
    # test that tax is 2450xmr
end

@testset "ch3 tax reducers" begin
    # check that mca is always 10% of amount
    # check marriage transfer is always basic rate tax credit
    # checl MCA only available if 1 spouse born before 6th April  1935
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )

    names = ExampleHouseholdGetter.initialise( Settings() )
    scot = ExampleHouseholdGetter.get_household( "mel_c2_scot" ) # scots are a married couple
    head = scot.people[SCOT_HEAD]
    spouse = scot.people[SCOT_SPOUSE]
    bu = make_benefit_unit( PeopleArray([head,spouse]), head.pid, spouse.pid )
    head_ages = [75,90,90,70] # after 1935
    spouse_ages = [90,70,70,90]
    head_incomes = [19_100.0, 29_710.0,41_080.0,0.0]
    spouse_incomes = [12_450.0,0,13_950.0,49_300.0]
    for i in 1:4
        head.income[private_pensions] = head_incomes[i]
        spouse.income[private_pensions] = spouse_incomes[i]
        # head.income[private_pension] = head_incomes[i]
        head.age = head_ages[i]
        spouse.age = spouse_ages[i]
        bruk = init_benefit_unit_result( bu )
        brscot = init_benefit_unit_result( bu )
        calc_income_tax!( bruk, head, spouse, itsys_ruk )
        calc_income_tax!( brscot, head, spouse, itsys_scot )
        if i == 1
            @test bruk.pers[head.pid].it.mca ≈ 891.50 ≈
                brscot.pers[head.pid].it.mca
            @test bruk.pers[spouse.pid].it.mca ≈ 0.0 ≈
                brscot.pers[spouse.pid].it.mca
        elseif i == 2
            @test bruk.pers[head.pid].it.mca ≈ 886.00 ≈
                brscot.pers[head.pid].it.mca
            @test bruk.pers[spouse.pid].it.mca ≈ 0.0 ≈
                brscot.pers[spouse.pid].it.mca
        elseif i == 3
            @test bruk.pers[head.pid].it.mca ≈ 345.00 ≈
                brscot.pers[head.pid].it.mca
            @test bruk.pers[spouse.pid].it.mca ≈ 0.0 ≈
                brscot.pers[spouse.pid].it.mca
        elseif i == 4
            @test bruk.pers[head.pid].it.mca ≈ 0.0 ≈
                brscot.pers[head.pid].it.mca
            @test bruk.pers[spouse.pid].it.mca ≈ 345.0 ≈
                brscot.pers[spouse.pid].it.mca
        end
    end
end

## TODO car tax; pension contributions

@testset "pension contributions tax relief standard case:  - Melville ch14 example 2(b)" begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = ExampleHouseholdGetter.initialise( Settings() )
    scot = ExampleHouseholdGetter.get_household( "mel_c2_scot" ) # scots are a married couple
    alana = scot.people[SCOT_HEAD]
    alana.income = Incomes_Dict{Float64}() # clear
    alana.income[self_employment_income] = 62_000.0
    alana.income[pension_contributions_employee] = (400.00*12)*0.8 # net contribs per month; expressed gross in example

    bu = make_benefit_unit( PeopleArray([alana]), alana.pid, BigInt(-1))
    bruk = init_benefit_unit_result( bu )
    brscot = init_benefit_unit_result( bu )

    calc_income_tax!( bruk, alana, nothing, itsys_ruk );
    calc_income_tax!( brscot, alana, nothing, itsys_scot );
    @test bruk.pers[alana.pid].it.pension_relief_at_source ≈ 960.0
    @test bruk.pers[alana.pid].it.pension_eligible_for_relief ≈ 400.0*12*0.8
    @test brscot.pers[alana.pid].it.non_savings_thresholds[1] ≈ itsys_scot.non_savings_thresholds[1]+400.0*12
    @test brscot.pers[alana.pid].it.non_savings_thresholds[2] ≈ itsys_scot.non_savings_thresholds[2]+400.0*12
    @test brscot.pers[alana.pid].it.pension_relief_at_source ≈ 960.0
    @test brscot.pers[alana.pid].it.pension_eligible_for_relief ≈ 400.0*12*0.8

end

@testset "pension: Tax Relief Minima - Melville ch14 ex1(a)"  begin
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    names = ExampleHouseholdGetter.initialise( Settings() )
    scot = ExampleHouseholdGetter.get_household( "mel_c2_scot" ) # scots are a married couple
    gordon = scot.people[SCOT_HEAD]
    gordon.income = Incomes_Dict{Float64}() # clear
    gordon.income[self_employment_income] = 27_800.0
    gordon.income[pension_contributions_employee] =  27_800.00 # net contribs per month; expressed gross in example
    bu = make_benefit_unit( PeopleArray([gordon]),gordon.pid, BigInt(-1))
    bruk = init_benefit_unit_result( bu )
    brscot = init_benefit_unit_result( bu )
    calc_income_tax!( bruk, gordon, nothing, itsys_ruk )
    calc_income_tax!( brscot, gordon, nothing, itsys_scot )
    @test bruk.pers[gordon.pid].it.pension_eligible_for_relief ≈ 27_800.0
    @test brscot.pers[gordon.pid].it.pension_eligible_for_relief ≈ 27_800.0
    gordon.income[self_employment_income] = 2_500.0
    gordon.income[pension_contributions_employee] =  27_800.00 # net contribs per month; expressed gross in example
    bruk = init_benefit_unit_result( bu )
    brscot = init_benefit_unit_result( bu )
    calc_income_tax!( bruk, gordon, nothing, itsys_ruk )
    calc_income_tax!( brscot, gordon, nothing, itsys_scot )
    @test bruk.pers[gordon.pid].it.pension_eligible_for_relief ≈ 3_600.0
    @test brscot.pers[gordon.pid].it.pension_eligible_for_relief ≈ 3_600.0
end

#
# no actual example for this and I'm not 100% sure this
#is exactly how it works, but should be near enough.
#
@testset "pension: Tax Relief Annual Allowance Charge"  begin
    itsys_ruk :: IncomeTaxSys = get_tax( scotland = false )
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    names = ExampleHouseholdGetter.initialise( Settings() )
    scot = ExampleHouseholdGetter.get_household( "mel_c2_scot" ) # scots are a married couple
    gordon = scot.people[SCOT_HEAD]
    gordon.income = Incomes_Dict{Float64}() # clear
    gordon.income[self_employment_income] = 60_000.0
    gordon.income[pension_contributions_employer] =  50_000.00 # net contribs per month; expressed gross in example
    bu = make_benefit_unit( PeopleArray([gordon]),gordon.pid, BigInt(-1))
    bruk = init_benefit_unit_result( bu )
    brscot = init_benefit_unit_result( bu )

    calc_income_tax!( bruk, gordon, nothing, itsys_ruk );
    @test bruk.pers[gordon.pid].it.pension_eligible_for_relief ≈ 40_000

    gordon.income[self_employment_income] = 300_000.0
    bruk = init_benefit_unit_result( bu )
    calc_income_tax!( bruk, gordon, nothing, itsys_ruk );
    @test bruk.pers[gordon.pid].it.pension_eligible_for_relief ≈ 10_000

    gordon.income[self_employment_income] = 150_000.0
    bruk = init_benefit_unit_result( Float64, bu )
    res_uk = calc_income_tax!( bruk, gordon, nothing, itsys_ruk );
    # tapering thing seems really weird - this is approximarely right
    @test bruk.pers[gordon.pid].it.pension_eligible_for_relief ≈ 40_000

    # complete-ish calc from combes, tutin&rowes, 2018 edn, p199, updated to 2019/20 rates
    gordon.income[self_employment_income] = 180_000.0
    gordon.income[pension_contributions_employer] =  14_400.00 # net contribs per month; expressed gross in example
    bruk = init_benefit_unit_result( bu )
    calc_income_tax!( bruk, gordon, nothing, itsys_ruk );
    @test bruk.pers[gordon.pid].income[INCOME_TAX] ≈ 61_500.0

end

@testset "Crude MCA Age Check" begin
    # cut-off for jan 2010 should be age 85

    d = Date( 2020, 1, 28 )
    itsys :: IncomeTaxSys = get_tax( scotland = false )
    @test old_enough_for_mca( itsys, 85, d )
    @test ! old_enough_for_mca( itsys, 84, d )
    @test old_enough_for_mca( itsys, 86, d )
end

@testset "Apply Allowance" begin
    allowance = 10_000
    allowance,t1 = apply_allowance( allowance, 5_000 )
    @test t1 == 0
    @test allowance == 5_000
    allowance,t2 = apply_allowance( allowance, 4_000 )
    @test t2 == 0
    @test allowance == 1_000
    allowance,t3 = apply_allowance( allowance, 4_000 )
    @test t3 == 3_000
    @test allowance == 0
end

@testset "Run on actual Data" begin
    nhhs,npeople = init_data()
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )
    it2 = deepcopy(itsys_scot)
    it2.non_savings_rates[1:3] .+= 0.5
    println( it2 )
    for hno in 1:nhhs
        hh = get_household(hno)
        if hno % 100 == 0
            println( "hhno $hno")
        end
        bus = get_benefit_units( hh )
        for bu in bus
            # income tax, with some nonsense for
            # what remains of joint taxation..
            head = get_head( bu )
            spouse = get_spouse( bu )
            bures = init_benefit_unit_result( bu )
            calc_income_tax!(
                bures,
                head,
                spouse,
                itsys_scot )
            calc_income_tax!(
                bures,
                head,
                spouse,
                it2 )
                for chno in bu.children
                child = bu.people[chno]
                calc_income_tax!(
                    bures.pers[child.pid],
                    child,
                    itsys_scot )
            end  # child loop
        end # bus loop
    end # hhld loop
end #

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
n = 3
m = 2

for elty in (Float32, Float64, ComplexF32, ComplexF64)

    A0 = randn(elty, n, n)
    B0 = randn(elty, n, m)
    B1 = randn(elty, n, m)

    A = copy(A0)

    # Full matrix
    linws = FastLapackInterface.LinSolveWs{elty,Int64}(n)

    FastLapackInterface.lu!(A, linws)
    @test A == A0
    F = lu(A)
    @test UpperTriangular(reshape(linws.lu, n, n)) ≈ F.U

    FastLapackInterface.lu!(A', linws)
    @test UpperTriangular(reshape(linws.lu, n, n)) ≈ F.U

    B = copy(B0)
    FastLapackInterface.linsolve_core!(A, B, linws)
    @test A == A0
    @test B ≈ A \ B0

    copy!(B, B1)
    FastLapackInterface.linsolve_core_no_lu!(A, B, linws)
    @test A == A0
    @test B ≈ A \ B1

    copy!(B, B0)
    FastLapackInterface.linsolve_core!(A', B, linws)
    @test B ≈ A' \ B0

    copy!(B, B1)
    FastLapackInterface.linsolve_core_no_lu!(A', B, linws)
    @test A == A0
    @test B ≈ A' \ B1

    # view of a matrix in upper left corner
    linws1 = FastLapackInterface.LinSolveWs{elty,Int64}(n - 1)
    C = view(A, 1:n-1, 1:n-1)
    D = view(B, 1:n-1, 1:m-1)
    D0 = copy(D)
    D1 = copy(D)
    FastLapackInterface.linsolve_core!(C, D, linws1)
    @test C == view(A, 1:n-1, 1:n-1)
    @test D ≈ C \ D0

    D = copy(D1)
    FastLapackInterface.linsolve_core_no_lu!(C, D, linws1)
    @test C == view(A, 1:n-1, 1:n-1)
    @test D ≈ C \ D1

    D = copy(D0)
    FastLapackInterface.linsolve_core!(C', D, linws1)
    @test C == view(A, 1:n-1, 1:n-1)
    @test D ≈ C' \ D0

    D = copy(D1)
    FastLapackInterface.linsolve_core_no_lu!(C', D, linws1)
    @test C == view(A, 1:n-1, 1:n-1)
    @test D ≈ C' \ D1

    # view of a matrix in lower left corner
    linws1 = FastLapackInterface.LinSolveWs{elty,Int64}(n - 1)
    C = view(A, 2:n, 1:n-1)
    C1 = copy(C)
    FastLapackInterface.lu!(C1, linws1)
    F = LinearAlgebra.lu!(C)
    @test triu(reshape(linws1.lu, n - 1, n - 1)) ≈ F.U
    @test tril(reshape(linws1.lu, n - 1, n - 1), -1) ≈ tril(F.L, -1)

    D0 = view(B0, 2:n, 1:m-1)
    D1 = view(B1, 2:n, 1:m-1)
    D = copy(D0)
    FastLapackInterface.linsolve_core!(C, D, linws1)
    @test C == view(A, 2:n, 1:n-1)
    @test D ≈ C \ D0

    D = copy(D1)
    FastLapackInterface.linsolve_core_no_lu!(C, D, linws1)
    @test C == view(A, 2:n, 1:n-1)
    @test D ≈ C \ D1

    D = copy(D0)
    FastLapackInterface.linsolve_core!(C', D, linws1)
    @test C == view(A, 2:n, 1:n-1)
    @test D ≈ C' \ D0

    D = copy(D1)
    FastLapackInterface.linsolve_core_no_lu!(C', D, linws1)
    @test C == view(A, 2:n, 1:n-1)
    @test D ≈ C' \ D1

    # using too big a work space
    linws1 = FastLapackInterface.LinSolveWs{elty,Int64}(n)
    C = view(A, 2:n, 1:n-1)
    C1 = copy(C)
    FastLapackInterface.lu!(C1, linws1)
    F = LinearAlgebra.lu!(C)
    @test triu(reshape(linws1.lu[1:(n-1)^2], n - 1, n - 1)) ≈ F.U
    @test tril(reshape(linws1.lu[1:(n-1)^2], n - 1, n - 1), -1) ≈ tril(F.L, -1)

    D0 = view(B0, 2:n, 1:m-1)
    D1 = view(B1, 2:n, 1:m-1)
    D = copy(D0)
    FastLapackInterface.linsolve_core!(C, D, linws1)
    @test C == view(A, 2:n, 1:n-1)
    @test D ≈ C \ D0

    D = copy(D1)
    FastLapackInterface.linsolve_core_no_lu!(C, D, linws1)
    @test C == view(A, 2:n, 1:n-1)
    @test D ≈ C \ D1

    D = copy(D0)
    FastLapackInterface.linsolve_core!(C', D, linws1)
    @test C == view(A, 2:n, 1:n-1)
    @test D ≈ C' \ D0

    D = copy(D1)
    FastLapackInterface.linsolve_core_no_lu!(C', D, linws1)
    @test C == view(A, 2:n, 1:n-1)
    @test D ≈ C' \ D1

end

ext: jl


text:
using VDPTag2
using Test
using POMDPSimulators
using MCTS
using POMDPs
using ParticleFilters
using ProgressMeter
using LinearAlgebra
using Random
using POMDPModelTools

Random.seed!(1)
pomdp = VDPTagPOMDP()
gen = NextMLFirst(mdp(pomdp), MersenneTwister(31))
s = TagState([1.0, 1.0], [-1.0, -1.0])

struct MyNode end
MCTS.n_children(::MyNode) = rand(1:10)

@inferred next_action(gen, pomdp, s, MyNode())
@inferred next_action(gen, pomdp, initial_state_distribution(pomdp), MyNode())

for a in linspace(0.0, 2*pi, 100)
    s = TagState(Vec2(0,0), Vec2(1,1))
    barriers = CardinalBarriers(0.2, 1.8)
    agent_speed = 1.0
    step_size = 0.5
    delta = agent_speed*step_size*Vec2(cos(a), sin(a))
    agent = VDPTag2.barrier_stop(barriers, s.agent, delta)
    @test agent == s.agent+delta
end

pomdp = VDPTagPOMDP()
for sao in stepthrough(pomdp, RandomPolicy(pomdp), "sao", max_steps=10)
    @show sao
end

dpomdp = AODiscreteVDPTagPOMDP(pomdp, 30, 0.5)
for sao in stepthrough(dpomdp, RandomPolicy(dpomdp), "sao", max_steps=10)
    @show sao
end

pomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.2, 1.8)))
filter = SIRParticleFilter(pomdp, 1000)
for sao in stepthrough(pomdp, ToNextML(pomdp), filter, "sao", max_steps=10)
    @show sao
end

# test to make sure it can't pass through any walls
pomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.0, 100.0)))
for quadrant in [Vec2(1,1), Vec2(-1,1), Vec2(1,-1), Vec2(-1,-1)]
    @showprogress for i in 1:100
        is = initial_state(pomdp, Base.GLOBAL_RNG)
        is = TagState(quadrant, is.target)
        for (s, sp) in stepthrough(pomdp, ToNextML(pomdp), filter, "s,sp", max_steps=100, initial_state=is)
            @test all(s.agent.*quadrant .>= 0.0)
            if s == sp
                println("did not move (this should not happen a bunch of times)")
            end
        end
        for (s, sp) in stepthrough(pomdp, RandomPolicy(pomdp), "s,sp", max_steps=100, initial_state=is)
            @test all(s.agent.*quadrant .>= 0.0)
            if s == sp
                println("did not move (this should not happen a bunch of times)")
            end
        end
    end
end

# make sure it does pass through walls sometimes without any barriers
pomdp = VDPTagPOMDP()
N = 100
for quadrant in [Vec2(1,1), Vec2(-1,1), Vec2(1,-1), Vec2(-1,-1)]
    in_other = falses(N)
    @showprogress for i in 1:N
        is = initial_state(pomdp, Base.GLOBAL_RNG)
        is = TagState(quadrant, is.target)
        hr = HistoryRecorder(max_steps=100)
        hist = simulate(hr, pomdp, ToNextML(pomdp), filter)
        in_other[i] = any(any(s.agent.*quadrant .< 0.0) for s in state_hist(hist))
    end
    @show sum(in_other)/length(in_other)
    println("Should be near one?")
    @test any(in_other)
end
for quadrant in [Vec2(1,1), Vec2(-1,1), Vec2(1,-1), Vec2(-1,-1)]
    in_other = falses(N)
    @showprogress for i in 1:N
        is = initial_state(pomdp, Base.GLOBAL_RNG)
        is = TagState(quadrant, is.target)
        hr = HistoryRecorder(max_steps=2)
        hist = simulate(hr, pomdp, RandomPolicy(pomdp))
        in_other[i] = any(any(s.agent.*quadrant .< 0.0) for s in state_hist(hist))
    end
    @show sum(in_other)/length(in_other)
    println("Should be near 3/4?")
    @test any(in_other)
end


text:
#==
This example shows a disconnected Pareto front.

Tanaka, Masahiro, et al.
"GA-based decision support system for multicriteria optimization."
Systems, Man and Cybernetics, 1995. Intelligent Systems for the 21st Century.,
EEE International Conference on. Vol. 2. IEEE, 1995.
==#

using MultiJuMP, JuMP
using Ipopt

m = multi_model(Ipopt.Optimizer)
@variable(m, 0 <= x[i=1:2] <= π)
@NLexpression(m, f1, x[1])
@NLexpression(m, f2, x[2])
@NLconstraint(m, x[1]^2+x[2]^2 -1 - 0.1cos(16atan(x[1]/x[2])) >= 0)
@NLconstraint(m, (x[1]-0.5)^2+(x[2]-0.5)^2 <= 0.5)

obj1 = SingleObjective(f1, sense = MOI.MIN_SENSE,
                       iv = Dict{String, Any}("x[$i]" => j for (i,j) in enumerate([0., 1.])))

# As the problem is nonconvex, we have to supply the
# initial value to get the global minimum of f2
obj2 = SingleObjective(f2, sense = MOI.MIN_SENSE,
                       iv = Dict{String, Any}("x[$i]" => j for (i,j) in enumerate([1., 0.])))

multim = get_multidata(m)
multim.objectives = [obj1, obj2]
multim.pointsperdim = 40
optimize!(m, method = NBI())

using Plots
plot(multim)


text:
"""
Struct holding the parameters needed at runtime in number format NF.
"""
@with_kw struct Constants{NF<:AbstractFloat}
    R_earth::NF      # Radius of Earth
    Ω::NF            # Angular frequency of Earth's rotation
    g::NF            # Gravitational acceleration
    akap::NF         # Ratio of gas constant to specific heat of dry air at constant pressure
    R::NF            # Gas constant
    γ::NF            # Reference temperature lapse rate (-dT/dz in deg/km)
    hscale::NF       # Reference scale height for pressure (in km)
    hshum::NF        # Reference scale height for specific humidity (in km)
    rh_ref::NF       # Reference relative humidity of near-surface air

    # TIME STEPPING
    Δt::NF                  # time step [s]
    robert_filter::NF       # Robert (1966) time filter coefficient to suppress comput. mode
    williams_filter::NF     # Williams time filter (Amezcua 2011) coefficient for 3rd order acc

    # DIFFUSION AND DRAG
    sdrag::NF               # drag [1/s] for zonal wind in the stratosphere
end

"""
Generator function for a Constants struct.
"""
function Constants{NF}(P::Params) where NF      # number format NF

    @unpack R_earth, Ω, g, akap, R, γ, hscale, hshum, rh_ref, Δt = P
    @unpack robert_filter, williams_filter = P
    @unpack tdrs = P

    # stratospheric drag [1/s] from drag time timescale tdrs [hrs]
    sdrag = 1/(tdrs*3600)

    # This implies conversion to NF
    return Constants{NF}(R_earth,Ω,g,akap,R,γ,hscale,hshum,rh_ref,
                            Δt,robert_filter,williams_filter,
                            sdrag)
end


text:
# common setup for Svalbard runs
#
# Todo

println("Time to load BITEModel.jl")
@time import BITEModel
const BM=BITEModel

update_cache = false
verbose = false
plotyes = false

# RGI data
all_svalbard = readdlm("data/rgi-6.0-huss/outlines/xyzn/svalbard/svalbard.dat")
cutoff_km2 = (5,10^10) # km^2
westfonna_ll = [18.5, 79.7] # long lat
austfonna_ll = [20, 79.2] # long lat
svalbard_nrs = Int[all_svalbard[i,1] for i=1:size(all_svalbard,1) if
                   cutoff_km2[1]<=all_svalbard[i,end]<=cutoff_km2[2] && # sizes
                   # remove Austfonna and Westfonna:
                   !(all_svalbard[i,2]>westfonna_ll[1] && all_svalbard[i,3]>westfonna_ll[2]) &&
                   !(all_svalbard[i,2]>austfonna_ll[1] && all_svalbard[i,3]>austfonna_ll[2])]
#with_thick = ["00027", "00408", "00409", "00427", "00428", "00429", "00434"][2:end] # 27 is on west or austfonna
with_thick = ["00409", "00428"] # the others are useless

brucebreen = BM.RGIGlacier("RGI60-07.01455", "Brucebreen") # area 6km^2, land, long lat: 17.291 78.485
rabotbreen = BM.RGIGlacier("RGI60-07.01478", "Rabotbreen") # area 70km^2, land
tunabreen = BM.RGIGlacier("RGI60-07.01458", "Tunabreen") # area 160km^2, tide-water
region_svalbard = 7

init_svalbard_gl(nr::Int; kws...) = init_svalbard_gl(BM.RGIGlacier(nr, region_svalbard); kws...)
init_svalbard_gl(rgi::Union{String,Symbol}; kws...) = init_svalbard_gl(BM.RGIGlacier(rgi); kws...)
function init_svalbard_gl(gid::BM.GlacierID; kws...)
    dt, pl = init_svalbard_gl_dt(gid; kws...)
    init_svalbard_gl(dt, pl)
end

function init_svalbard_gl(dt::BM.DataTable, pl)
    gl,pp,pm,pn,pl = BM.load_glacier(dt,pl)
    gb, pm = BM.make_bands(gl, pp, pm, pn);
    return gl,gb,pp,pm,pn,pl,dt
end
function init_svalbard_gl_dt(gid; kws...)
    pl = BM.LoadPara(gid; kws...)
    dt, pl = BM.make_datatable(gid,pl)

    # set temperature profile:
    dt.temp.opts[:temp] = [-10, -1.0] # colder at tongue than at top
    dt.temp.opts[:ele] = linspace(-30.0, 1717, 2) # a bit below sea to height of highest peak on Svalbard
    return dt, pl
end


text:
#!/usr/bin/env julia
#
# Example usage:
#
#   julia reproducibility_scanning_tools.jl libusb
#
# Alternatively, you can get a shell and interactively do things:
#
#   julia -i reproducibility_scanning_tools.jl
#
#      julia> check_project("libusb")
#      julia> do_full_test()

using Random, SHA, JSON, BinaryBuilder, Tar, TranscodingStreams, CodecZlib, SimpleBufferStream

skip_dirs = ("0_RootFS",)
skip_projects = ("OpenBLAS", "LLVM")

function decompress!(input::IO, output::IO; blocksize::Int = 2*1024*1024)
    output = TranscodingStream(GzipDecompressor(), output)

    # Write that first chunk, then all the rest
    while !eof(input)
        write(output, read(input, blocksize))
    end

    # Close the TranscodingStream (we expect the caller to close the real `output`)
    close(output)
end

function collect_projects(root::AbstractString = @__DIR__)
    projects = String[]
    for d in readdir(root; join=true)
        if !isdir(d) || basename(d) in skip_dirs
            continue
        end

        for proj in readdir(d; join=true)
            if !isdir(proj) || basename(proj) in skip_projects
                continue
            end
            if !isfile(joinpath(proj, "build_tarballs.jl"))
                continue
            end

            push!(projects, proj)
        end
    end
    return projects
end

function collect_platforms(project::AbstractString)
    platforms = String[]
    cd(project) do
        mktempdir() do tmpdir
            meta_json = joinpath(tmpdir, "meta.json")
            if !success(`$(Base.julia_cmd()) build_tarballs.jl --meta-json=$(meta_json)`)
                error("Collecting the platforms for $(project) failed")
            end

            json = String(read(meta_json))
            buff = IOBuffer(strip(json))
            while !eof(buff)
                obj = BinaryBuilder.JSON.parse(buff)
                append!(platforms, get(obj, "platforms", String[]))
            end
        end
    end
    return platforms
end

function inlogdir(hdr)
    path = hdr.path
    if startswith(path, "./")
        path = path[3:end]
    end
    # Log files can be either directly under `logs` or in its subdirectories
    return dirname(path) == "logs" || occursin(r"^logs/", dirname(path))
end

function collect_tar_tree(io::IO;
                          HashType = SHA.SHA256_CTX,
                          buf::Vector{UInt8} = Vector{UInt8}(undef, Tar.DEFAULT_BUFFER_SIZE),)
    tree = Dict{String,Any}()
    Tar.read_tarball(!inlogdir, io; buf=buf) do hdr, parts
        isempty(parts) && return
        name = pop!(parts)
        node = tree
        for part in parts
            node′ = get(node, part, nothing)
            if !(node′ isa Dict)
                node′ = node[part] = Dict{String,Any}()
            end
            node = node′
        end
        if hdr.type == :directory
            node[name] = Dict{String,Any}()
            return
        end
        if hdr.type == :symlink
            mode = "120000"
            hash = Tar.git_object_hash("blob", HashType) do io
                write(io, hdr.link)
            end
        elseif hdr.type == :file
            mode = iszero(hdr.mode & 0o100) ? "100644" : "100755"
            hash = Tar.git_file_hash(io, hdr.size, HashType, buf=buf)
        else
            error("unsupported type for git tree hashing: $(hdr.type)")
        end
        node[name] = (mode, hash)
    end
    return tree
end

function collect_compressed_tar_tree(filename)
    open(filename) do io_gz
        io = BufferStream()
        t_decomp = @async decompress!(io_gz, io)
        t_tree = @async collect_tar_tree(io)
        wait(t_decomp)
        return fetch(t_tree)
    end
end

function compare(a::Dict, a_filename::AbstractString, b::Dict, b_filename::AbstractString, prefix::AbstractString = "/")
    function key_equivalence(test_dict, check_dict, check_filename)
        all_ok = true
        for k in keys(test_dict)
            if k ∉ keys(check_dict)
                @error("$(prefix)/$(k) not present in $(check_filename)")
                all_ok = false
            end
            if isa(test_dict[k], Dict) && !isa(check_dict[k], Dict)
                @error("$(prefix)/$(k) is not a directory in $(check_filename)")
                all_ok = false
            end
        end
        return all_ok
    end
    all_ok = true
    all_ok &= key_equivalence(a, b, b_filename)
    all_ok &= key_equivalence(b, a, a_filename)

    for k in keys(a)
        if isa(a[k], Dict)
            # Recurse into dicts
            all_ok &= compare(a[k], a_filename, b[k], b_filename, joinpath(prefix, k))
        else
            # Check permissions and hashes of leaf nodes
            a_perms, a_hash = a[k]
            b_perms, b_hash = b[k]

            if a_perms != b_perms
                @error("Permissions mismatch!", path=joinpath(prefix,k), a_filename, b_filename, a_perms, b_perms)
                all_ok = false
            end
            if a_hash != b_hash
                @error("Content hash mismatch!", path=joinpath(prefix,k), a_filename, b_filename, a_hash, b_hash)
                all_ok = false
            end
        end
    end
    return all_ok
end

function check_project(project, platform)
    cd(project) do
        rm("products"; recursive=true, force=true)
        rm("build"; recursive=true, force=true)
        # Do first build
        @info("Building $(project) for $(platform) for the first time...")
        run(`$(Base.julia_cmd()) --color=yes build_tarballs.jl --verbose --debug $(platform)`)
        for f in readdir("products"; join=true)
            mv(f, "$(f).old")
        end
        # Then do second build
        @info("Building $(project) for $(platform) for the second time...")
        run(`$(Base.julia_cmd()) --color=yes build_tarballs.jl $(platform)`)

        # Now compare everything except for the `log/` directory:
        for f in readdir("products"; join=true)
            if !isfile("$(f).old")
                continue
            end
            new_tree = collect_compressed_tar_tree(f)
            old_tree = collect_compressed_tar_tree("$(f).old")
            if !compare(new_tree, f, old_tree, "$(f).old")
                error("$(basename(project)) is not reproducible for $(platform)!")
            else
                @info("$(basename(project)) passes for $(platform)!")
            end
        end
    end
end

function check_project(project)
    if basename(project) == project && !isdir(project)
        project = joinpath(uppercase(project[1:1]), project)
    end
    if !isdir(project)
        error("Unknown project $(project)")
    end

    @info("Investigating $(project)")
    platforms = collect_platforms(project)

    # Try to build for linux64 if we can, as that's most convenient
    linux64_platforms = sort!(filter(p -> occursin("x86_64-linux-gnu", p), platforms))
    local initial_platform
    if isempty(linux64_platforms)
        initial_platform = first(platforms)
    else
        initial_platform = first(linux64_platforms)
    end

    # Perform initial build with initial platform
    check_project(project, initial_platform)

    # If that worked, check all other platforms too!
    for other_platform in platforms
        if other_platform == initial_platform
            continue
        end
        check_project(project, other_platform)
    end
end

# For the truly insane
function do_full_test()
    projects = collect_projects()
    for project in shuffle(projects)
        check_project(project)
    end
end

# Auto-run projects if given on command line:
if !isempty(ARGS)
    for project in ARGS
        check_project(project)
    end
else
    println("Tools loaded, try check_project(\"libusb\")")
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
fBoundary1D(x,t) = x > 0 ? 1 : 0
fBoundary2D_x(x,y,t) = 0
fBoundary2D_y(x,y,t) = y > 0.5 ? 1 : 0
fBoundary3D(x,y,z,t) = 1

@testset "Medium" begin

    #Test we can call the variable and the function
    @test_throws ErrorException try @eval @agent(1, u::Medium) catch err; throw(err.error) end
    @test begin m = @agent(1, u::Medium, Boundary = BoundaryFlat(1)); m.declaredSymbols["Medium"] == [:u] end
    @test begin m = @agent(1, [u,v]::Medium, Boundary = BoundaryFlat(1)); m.declaredSymbols["Medium"] == [:u,:v] end
    @test begin m = @agent(1, UpdateMedium = ∂t(u) = Δ(u) + δx(0.)*δy(0.)*1.); "UpdateMedium" in keys(m.declaredUpdates) end

    #Medium access in Community
    m = @agent(3,
        [u,v]::Medium,
        Boundary = BoundaryFlat(3,Bounded(),Bounded(),Bounded())
        )
    mc = compile(m)
    # println(mc.program)
    com = Community(mc,N=10,mediumN=[10,10,10])

    #Check call and boundary
    @test_nowarn b = com.u #Read
    @test_nowarn com.v .= 1 #Assign scalar
    @test_nowarn com.v .= zeros(10,10,10) #Assign vector
    @test_nowarn com.v[1,2,3] = 1; #Assign single value

    @test begin
        com.u .= 0
        com.simulationBox .= [-10 10;-10 10;-10 10]
        comt = mc.evolve(com,dt=0.1,tMax=1)

        all(comt[end].u .== zeros(10,10,10))
    end

    for platform in testplatforms

        ########################################################################################################33
        #Dirichlet boundary
        ########################################################################################################33

        ## Default

        ###Dim 1
        m = @agent(1,
            [u,v]::Medium,

            UpdateMedium = ∂t(u) = Δx(u),

            Boundary = BoundaryFlat(1,
                            Bounded(medium=DirichletBoundaryCondition()),
                            )
            )
        mc = compile(m,platform=platform)
        # println(mc.program)

        @test begin #Test symmetric decay and loose of concentration
            com = Community(mc,N=10,mediumN=[10])
            com.u .= 1
            com.simulationBox .= [-10 10]
            comt = mc.evolve(com,dt=0.1,tMax=1000)
            all(comt[end].u .< 10E-8) && all(abs.(comt[100].u - reverse(comt[100].u)) .< 10E-8)
        end

        ###Dim 2
        m = @agent(2,
            [u,v]::Medium,

            UpdateMedium = ∂t(u) = Δx(u),

            Boundary = BoundaryFlat(2,
                            Bounded(medium=DirichletBoundaryCondition()),
                            Bounded(medium=DirichletBoundaryCondition()),
                            )
            )
        mc = compile(m,platform=platform)
        # println(mc.program)

        @test begin #Test symmetric decay and loose of concentration
            com = Community(mc,N=10,mediumN=[10,15])
            com.u .= 1
            com.simulationBox .= [-10 10;-20 20]
            comt = mc.evolve(com,dt=0.1,tMax=1000)
            all(comt[end].u .< 10E-8) && all(abs.(comt[100].u - reverse(comt[100].u)) .< 10E-8)
        end

        ###Dim 3
        m = @agent(3,
            [u,v]::Medium,

            UpdateMedium = ∂t(u) = Δx(u),

            Boundary = BoundaryFlat(3,
                            Bounded(medium=DirichletBoundaryCondition()),
                            Bounded(medium=DirichletBoundaryCondition()),
                            Bounded(medium=DirichletBoundaryCondition()),
                            )
            )
        mc = compile(m,platform=platform)
        # println(mc.program)

        @test begin #Test symmetric decay and loose of concentration
            com = Community(mc,N=10,mediumN=[5,10,15])
            com.u .= 1
            com.simulationBox .= [-10 10;-20 20;-30 30]
            comt = mc.evolve(com,dt=0.1,tMax=1000)
            all(comt[end].u .< 10E-8) && all(abs.(comt[100].u - reverse(comt[100].u)) .< 10E-8)
        end

        ## With a user defined function

        ###Dim 1
        m = @agent(1,
            [u,v]::Medium,

            UpdateMedium = ∂t(u) = Δx(u),

            Boundary = BoundaryFlat(1,
                            Bounded(medium=DirichletBoundaryCondition(fBoundary1D)),
                            )
            )
        mc = compile(m,platform=platform)
        # println(mc.program)

        @test begin #Test symmetric decay and loose of concentration
            com = Community(mc,N=10,mediumN=[10])
            com.u .= 1
            com.simulationBox .= [-10 10]
            comt = mc.evolve(com,dt=0.1,tMax=1000)
            all(abs.(comt[end].u .- range(0.1,.9, length=10)) .< 10E-3)
        end

        ###Dim 2
        m = @agent(2,
            [u,v]::Medium,

            UpdateMedium = ∂t(u) = Δx(u) + Δy(u),

            Boundary = BoundaryFlat(2,
                            Bounded(medium=DirichletBoundaryCondition(fBoundary2D_x)),
                            Bounded(medium=DirichletBoundaryCondition(fBoundary2D_y)),
                            )
            )
        mc = compile(m,platform=platform)
        # println(mc.program)

        @test begin #Test symmetric decay and loose of concentration
            Nx = 10; Ny = 10
            com = Community(mc,N=0,mediumN=[Nx,Ny])
            com.u .= 1
            com.simulationBox .= [0 1;0 1]
            comt = mc.evolve(com,dt=0.001,tMax=2000)
            m = zeros(Nx,Ny)
            for (i,x) in enumerate(range(1/Nx,(Nx-1)/Nx,length=Nx))
                for (j,y) in enumerate(range(1/Ny,(Ny-1)/Ny,length=Ny))
                    for n in 1:2:40
                        m[i,j] += 4/(π*n*sinh(π*n))*sinh(n*π*y)*sin(n*π*x)
                    end
                end
            end
            all(abs.(comt[end].u .- m) .< 10E-2)
        end

        ###Dim 3
        m = @agent(3,
            [u,v]::Medium,

            UpdateMedium = ∂t(u) = Δx(u),

            Boundary = BoundaryFlat(3,
                            Bounded(medium=DirichletBoundaryCondition(fBoundary3D)),
                            Bounded(medium=DirichletBoundaryCondition(fBoundary3D)),
                            Bounded(medium=DirichletBoundaryCondition(fBoundary3D)),
                            )
            )
        mc = compile(m,platform=platform)
        # println(mc.program)

        @test begin #Test symmetric decay and loose of concentration
            com = Community(mc,N=10,mediumN=[5,10,15])
            com.u .= 0
            com.simulationBox .= [-10 10;-20 20;-30 30]
            comt = mc.evolve(com,dt=0.1,tMax=1000)
            all(comt[end].u .- 1 .< 10E-8)
        end

        #######################################################################################################33
        #Newmann boundary
        #######################################################################################################33

        ## Default

        ###Dim 1
        m = @agent(1,
            [u,v]::Medium,

            UpdateMedium = ∂t(u) = Δx(u),

            Boundary = BoundaryFlat(1,
                            Bounded(medium=NewmannBoundaryCondition()),
                            )
            )
        mc = compile(m,platform=platform)
        # println(mc.program)

        @test begin #Test conservation of matter and homogeneization
            Nx = 10
            com = Community(mc,N=10,mediumN=[10])
            com.u .= exp.(- (range(0,Nx,length=Nx).-Nx/2).^2)
            com.simulationBox .= [-10 10]
            comt = mc.evolve(com,dt=0.1,tMax=1000)
            all( sum(comt[1].u) .- sum(comt[end].u) .< 0.02 ) && (abs(minimum(comt[end].u) - maximum(comt[end].u)) .< 10E-4)
        end

        ###Dim 2
        m = @agent(2,
            [u,v]::Medium,

            UpdateMedium = ∂t(u) = Δx(u) + Δy(u),

            Boundary = BoundaryFlat(2,
                            Bounded(medium=NewmannBoundaryCondition()),
                            Bounded(medium=NewmannBoundaryCondition())
                            )
            )
        mc = compile(m,platform=platform)
        # println(mc.program)

        @test begin #Test conservation of matter and homogeneization
            Nx = 10; Ny = 10
            com = Community(mc,N=10,mediumN=[Nx,Ny])
            for (i,x) in enumerate(range(0,Nx,length=Nx))
                for (j,y) in enumerate(range(0,Ny,length=Ny))
                    com.u[i,j] = exp.(- (i.-Nx/2).^2- (j.-Ny/2).^2)
                end
            end
            com.simulationBox .= [-10 10;-10 10]
            comt = mc.evolve(com,dt=0.1,tMax=1000)
            all( sum(comt[1].u) .- sum(comt[end].u) .< 10E-1 ) && (minimum(comt[end].u) - maximum(comt[end].u) .< 10E-2)
        end

        ###Dim 2
        m = @agent(3,
            [u,v]::Medium,

            UpdateMedium = ∂t(u) = Δx(u) + Δy(u) + Δz(u),

            Boundary = BoundaryFlat(3,
                            Bounded(medium=NewmannBoundaryCondition()),
                            Bounded(medium=NewmannBoundaryCondition()),
                            Bounded(medium=NewmannBoundaryCondition())
                            )
            )
        mc = compile(m,platform=platform)
        # println(mc.program)

        @test begin #Test conservation of matter and homogeneization
            Nx = 10; Ny = 10; Nz = 10
            com = Community(mc,N=10,mediumN=[Nx,Ny,Nz])
            for (i,x) in enumerate(range(0,Nx,length=Nx))
                for (j,y) in enumerate(range(0,Ny,length=Ny))
                    for (k,z) in enumerate(range(0,Nz,length=Nz))
                        com.u[i,j,k] = exp.(- (i.-Nx/2).^2 - (j.-Ny/2).^2 - (k.-Nz/2).^2)
                    end
                end
            end
            com.simulationBox .= [-10 10;-10 10;-10 10]
            comt = mc.evolve(com,dt=0.1,tMax=1000)
            all( sum(comt[1].u) .- sum(comt[end].u) .< 10E-1 ) && (minimum(comt[end].u) - maximum(comt[end].u) .< 10E-2)
        end

        ########################################################################################################33
        #Delta sources
        ########################################################################################################33

        ##1 dim
        m = @agent(1,
                    [u,v]::Medium,
                    # UpdateMedium = ∂t(u) = Δx(u) + δx(0.),
                    Boundary = BoundaryFlat(1,Bounded())
                    )
        mc = compile(m,platform=platform)
        # println(prettify(mc.program))
        com = Community(mc,N=10,mediumN=[30])
        com.u .= 0
        com.simulationBox .= [-10 10]
        @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
        @test begin
            comt = mc.evolve(com,dt=0.01,tMax=4000,dtSave=200)

            all(abs.(comt[end].u .- comt[10].u) .< 0.001)
        end

        ##2 dim
        m = @agent(2, [u,v]::Medium, UpdateMedium = ∂t(u) = Δx(u) + Δy(u) + δx(0.) + δy(0.), Boundary = BoundaryFlat(2,Bounded(),Bounded()))
        mc = compile(m,platform=platform)
        # println(mc.program)
        com = Community(mc,N=10,mediumN=[30,30])
        com.u .= 0
        com.simulationBox .= [-10 10; -10 10]
        @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
        @test begin
            comt = mc.evolve(com,dt=0.01,tMax=4000,dtSave=200)

            all(abs.(comt[end].u .- transpose(comt[10].u)) .< 0.001)
        end

        ##3 dim
        m = @agent(3,
                    [u,v]::Medium,
                    UpdateMedium = ∂t(u) = Δx(u) + Δy(u) + Δz(u) + δx(0.) + δy(0.) + δz(0.),
                    Boundary = BoundaryFlat(3,Bounded(),Bounded(),Bounded()))
        mc = compile(m,platform=platform)
        # println(mc.program)
        com = Community(mc,N=10,mediumN=[5,5,5])
        com.u .= 0
        com.simulationBox .= [-10 10; -10 10; -10 10]
        @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
        @test begin
            comt = mc.evolve(com,dt=0.01,tMax=4000,dtSave=200)

            all(abs.(comt[end].u .- comt[10].u) .< 0.001)
        end

        ########################################################################################################33
        #Coupling sources
        ########################################################################################################33

        ########################################################################################################33
        #Agents reading from medium
        ########################################################################################################33

    end

    # for i in testplatforms #Check coumpling agents with medium

    #     m = @agent(1,
    #             p::Local,
    #             [u,v]::Medium,
    #             UpdateMedium = begin
    #                 ∂t(u) = Δx(u)
    #                 ∂t_v = Δx(v) + δx(0.)
    #             end,
    #             UpdateMediumInteraction = u += p*dt
    #             )
    #     s = SimulationFree(m,
    #                     box=[(:x,-20.,20.)],
    #                     medium=[MediumFlat("Dirichlet",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=1)
    #     com.u .= 0.
    #     com.v .= 0.
    #     com.x .= 0.
    #     com.p .= 1.
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=4000,dtSave=200)

    #         all(abs.(comt[end].u .- comt[end].v) .< 0.1)
    #     end

    # end

    # for i in testplatforms #Check agents reading from medium

    #     m = @agent(1,
    #             p::Local,
    #             u::Medium,
    #             k::Global,

    #             Equation = begin
    #                 d_p = (u-k*p)*dt
    #             end
    #             )
    #     s = SimulationFree(m,
    #                     box=[(:x,0.,100.)],
    #                     medium=[MediumFlat("Dirichlet",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=100)
    #     com.x .= [Float64(i) for i in 0:1:99]
    #     com.u .= [Float64(i) for i in 0:99]
    #     com.k = 1.
    #     com.p .= 0.
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=4000,dtSave=200)
    #         all(comt[end].p .<= [Float64(i) for i in 0:1:99])
    #     end

    # end

end


text:
# # Approximations of the bilinear term

# Ultimately we want to evolve, in time, an equation of the form # $$
#   \omega_t = G(\omega).
# $$
# For the 2D NSE, $G$ is of the form
# $$
#  G(\omega) = g + \nu \Delta \omega + B(\omega),
# $$
# where
# $$
#   B(\omega) = \left(\partial_x^2 - \partial_y^2\right)(uv) + \partial_{xy}\left(v^2 - u^2\right),
# $$
# with $\mathbf{u} = (u(x, y), v(x, y))$ being the velocity field associated with the vorticity, i.e.  $\omega = \boldsymbol{\nabla} \times \mathbf{u} = v_x - u_y$.

# In spectral space, $\Delta \omega$ is just a diagonal operator, while $g$ is a given scalar field, which are both easy to compute. The most computationally demanding term is $B(\omega)$.

# Here, we test different ways of approximating $B(\omega)$.

# The packages we are gonna need.

using FFTW
using Plots
using LinearAlgebra: mul!
using Test
using BenchmarkTools
using Random

@info "Threads: $(FFTW.nthreads())"

# ## Operators

# We first define a method to build the various operators acting in spectral space:

function get_operators(N, κ₀)
    ## Differentiation in spectral space
    Dx_hat = im * κ₀ * [
        ifelse(k1 ≤ div(N, 2) + 1, k1 - 1, k1 - 1 - N) for k2 in 1:div(N, 2)+1, k1 in 1:N
    ]
    Dy_hat = im * κ₀ * [k2 - 1 for k2 in 1:div(N, 2)+1, k1 in 1:N]

    Delta_hat = - κ₀^2 * [
        ifelse(k1 ≤ div(N, 2) + 1, (k1 - 1)^2 + (k2 - 1)^2, (k1 - 1 - N)^2 + (k2 - 1)^2)
        for k2 in 1:div(N, 2)+1, k1 in 1:N
    ]

    ## For the Basdevant formulation
    DxsqDysq_hat = Dx_hat.^2 .- Dy_hat.^2
    Dxy_hat = Dx_hat .* Dy_hat

    ## Recovering of the velocity field from the vorticity
    Hu_hat = - Dy_hat ./ Delta_hat
    Hu_hat[1, 1] = 0.0
    Hv_hat = Dx_hat ./ Delta_hat
    Hv_hat[1, 1] = 0.0

    return Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat, DxsqDysq_hat, Dxy_hat
end

# ## Generation of a scalar field

# Method to generate a scalar field from a list of wavenumbers and amplitudes:

function field_from_modes(L, N, modes::Matrix{<:Integer}, amps::Matrix{<:Real})
    κ₀ = 2π/L
    x = y = (L/N):(L/N):L
    field = sum(
        [
            2κ₀^2 * (k[1]^2 + k[2]^2) * (
                a[1] * cos.(κ₀ * (k[1] * one.(y) * x' + k[2] * y * one.(x)'))
                - a[2] * sin.(κ₀ * (k[1] * one.(y) * x' + k[2] * y * one.(x)'))
            )
            for (k, a) in zip(eachrow(modes), eachrow(amps))
        ]
    )
    return field
end

# Method to randomly generate a scalar field from a given random number generator and a given number wavenumbers to be excited:

function field_from_modes(rng::AbstractRNG, L, N, num_modes::Int)

    modes = rand(rng, 1:div(N,10), num_modes, 2)
    amps = rand(rng, num_modes, 2)

    field = field_from_modes(L, N, modes, amps)

    return field
end

# Method to randomly generate a scalar field from a given number wavenumbers to be excited:

field_from_modes(L, N, num_modes::Int) = field_from_modes(Xoshiro(), L, N, num_modes)

# ## Methods for computing the bilinear term

function bilinear_naive(vort_hat, params)
    operators, vars = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat = operators
    N, Nsub = vars

    u_hat = Hu_hat .* vort_hat
    v_hat = Hv_hat .* vort_hat
    u = irfft(u_hat, N)
    v = irfft(v_hat, N)
    vort = irfft(vort_hat, N)
    wu_hat = rfft(vort .* u)
    wv_hat = rfft(vort .* v)

    bilin_hat = Dx_hat .* wu_hat .+ Dy_hat .* wv_hat
    ## dealiasing
    bilin_hat[div(Nsub,2) + 1:end, :] .= 0.0im
    bilin_hat[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im
    return bilin_hat
end

#

function bilinear_auxs!(bilin_hat, vort_hat, params)
    operators, vars, auxs = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat = operators
    N, Nsub = vars
    u_hat, v_hat, u, v, w, wu, wv, wu_hat, wv_hat = auxs

    u_hat .= Hu_hat .* vort_hat
    v_hat .= Hv_hat .* vort_hat
    u .= irfft(u_hat, N)
    v .= irfft(v_hat, N)
    w .= irfft(vort_hat, N)
    wu .= w .* u
    wv .= w .* v
    wu_hat .= rfft(wu)
    wv_hat .= rfft(wv)

    bilin_hat .= Dx_hat .* wu_hat .+ Dy_hat .* wv_hat
    ## dealiasing
    bilin_hat[div(Nsub,2) + 1:end, :] .= 0.0im
    bilin_hat[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im

    return bilin_hat
end

#

function bilinear_plan_brdcst!(bilin_hat, vort_hat, params)
    operators, vars, auxs, plans = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat = operators
    N, Nsub = vars
    u_hat, v_hat, u, v, w, wu, wv, wu_hat, wv_hat = auxs
    plan, plan_inv = plans

    u_hat .= Hu_hat .* vort_hat
    v_hat .= Hv_hat .* vort_hat
    mul!(u, plan_inv, u_hat)
    mul!(v, plan_inv, v_hat)
    ## mul! on inverse plan mutates last argument, but we need to preserve vort_hat
    ## so, reuse v_hat for that:
    v_hat .= vort_hat
    mul!(w, plan_inv, v_hat)
    wu .= w .* u
    wv .= w .* v

    mul!(wu_hat, plan, wu)
    mul!(wv_hat, plan, wv)

    bilin_hat .= Dx_hat .* wu_hat .+ Dy_hat .* wv_hat
    ## dealiasing
    bilin_hat[div(Nsub,2) + 1:end, :] .= 0.0im
    bilin_hat[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im

    return bilin_hat
end

#

function bilinear_Basdevant_brdcst(vort_hat, params)
    operators, vars = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat = operators
    N, Nsub = vars

    u_hat = Hu_hat .* vort_hat
    v_hat = Hv_hat .* vort_hat
    u = irfft(u_hat, N)
    v = irfft(v_hat, N)
    uv_hat = rfft(u .* v)
    v2u2_hat = rfft(v.^2 - u.^2)

    bilin_hat = DxsqDysq_hat .* uv_hat .+ Dxy_hat .* v2u2_hat
    ## dealiasing
    bilin_hat[div(Nsub,2) + 1:end, :] .= 0.0im
    bilin_hat[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im

    return bilin_hat
end

#

function bilinear_Basdevant_plan_brdcst!(bilin_hat, vort_hat, params)
    operators, vars, auxs, plans = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat = operators
    N, Nsub = vars
    u_hat, v_hat, u, v, w, uv, v2u2, uv_hat, v2u2_hat = auxs
    plan, plan_inv = plans

    u_hat .= Hu_hat .* vort_hat
    v_hat .= Hv_hat .* vort_hat

    mul!(u, plan_inv, u_hat)
    mul!(v, plan_inv, v_hat)

    uv .= u .* v
    v2u2 .= v.^2 .- u.^2
    mul!(uv_hat, plan, uv)
    mul!(v2u2_hat, plan, v2u2)

    bilin_hat .= DxsqDysq_hat .* uv_hat .+ Dxy_hat .* v2u2_hat

    ## dealiasing
    bilin_hat[div(Nsub,2) + 1:end, :] .= 0.0im
    bilin_hat[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im

    return bilin_hat
end

#

function bilinear_Basdevant_plan_loop!(bilin_hat, vort_hat, params)
    operators, vars, auxs, plans = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat = operators
    N, Nsub = vars
    u_hat, v_hat, u, v, w, uv, v2u2, uv_hat, v2u2_hat = auxs
    plan, plan_inv = plans

    @inbounds for k in eachindex(u_hat)
        u_hat[k] = Hu_hat[k] * vort_hat[k]
        v_hat[k] = Hv_hat[k] * vort_hat[k]
    end

    mul!(u, plan_inv, u_hat)
    mul!(v, plan_inv, v_hat)

    @inbounds for k in eachindex(u)
        uv[k] = u[k] * v[k]
        v2u2[k] = v[k]^2 - u[k]^2
    end

    mul!(uv_hat, plan, uv)
    mul!(v2u2_hat, plan, v2u2)

    @inbounds for k in eachindex(vort_hat)
        bilin_hat[k] = DxsqDysq_hat[k] * uv_hat[k] + Dxy_hat[k] * v2u2_hat[k]
    end
    ## dealiasing
    @inbounds for j in 1:N, i in div(Nsub,2) + 1:div(N, 2) + 1
        bilin_hat[i, j] = 0.0im
    end
    @inbounds for j in div(Nsub,2) + 1:div(N,2) + div(Nsub,2), i in 1:div(N, 2) + 1
        bilin_hat[i, j] = 0.0im
    end

    return bilin_hat
end

#

function bilinear_Basdevant_plan_doubleloop!(bilin_hat, vort_hat, params)
    operators, vars, auxs, plans = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat = operators
    N, Nsub = vars
    u_hat, v_hat, u, v, w, uv, v2u2, uv_hat, v2u2_hat = auxs
    plan, plan_inv = plans

    M, N = size(vort_hat)

    @inbounds for j in 1:N, i in 1:M
        u_hat[i, j] = Hu_hat[i, j] * vort_hat[i, j]
        v_hat[i, j] = Hv_hat[i, j] * vort_hat[i, j]
    end

    mul!(u, plan_inv, u_hat)
    mul!(v, plan_inv, v_hat)

    @inbounds for j in 1:N, i in 1:N
        uv[i, j] = u[i, j] * v[i, j]
        v2u2[i, j] = v[i, j]^2 - u[i, j]^2
    end

    mul!(uv_hat, plan, uv)
    mul!(v2u2_hat, plan, v2u2)

    @inbounds for j in 1:N, i in 1:M
        bilin_hat[i, j] = DxsqDysq_hat[i, j] * uv_hat[i, j] + Dxy_hat[i, j] * v2u2_hat[i, j]
    end
    ## dealiasing
    @inbounds for j in 1:N, i in div(Nsub,2) + 1:M
        bilin_hat[i, j] = 0.0im
    end
    @inbounds for j in div(Nsub,2) + 1:div(N,2) + div(Nsub,2), i in 1:M
        bilin_hat[i, j] = 0.0im
    end

    return bilin_hat
end

# ## The spatial domain and its discretization

L = 2π
κ₀ = 2π/L
N = 256
Nsub = 162
x = y = (L/N):(L/N):L

# ## Preparation for tests

Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat, DxsqDysq_hat, Dxy_hat = get_operators(N, κ₀)

vort = field_from_modes(L, N, 48)
vort_hat = rfft(vort)

u_hat = similar(vort_hat)
v_hat = similar(vort_hat)
u = similar(vort)
v = similar(vort)
w = similar(vort)
uv = similar(vort)
v2u2 = similar(vort)
uv_hat = similar(vort_hat)
v2u2_hat = similar(vort_hat)
bilin_hat = similar(vort_hat)

plan = plan_rfft(vort, flags=FFTW.MEASURE)
plan_inv = plan_irfft(vort_hat, N, flags=FFTW.MEASURE)

operators = Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat, DxsqDysq_hat, Dxy_hat
vars = N, Nsub
auxs = u_hat, v_hat, u, v, w, uv, v2u2, uv_hat, v2u2_hat
plans = plan, plan_inv

# ## Check results

# As a first test, we check that all implementations return about the same vector

@testset "same result bilinear terms" begin

    bilin_hat_naive = bilinear_naive(
        vort_hat,
        (operators, vars)
    )

    bilinear_hat_brdcst = similar(vort_hat)
    bilinear_auxs!(bilinear_hat_brdcst, vort_hat, (operators, vars, auxs))

    bilin_hat_plan_brdcst = similar(vort_hat)
    bilinear_plan_brdcst!(
        bilin_hat_plan_brdcst,
        vort_hat,
        (operators, vars, auxs, plans)
    )

    bilin_hat_Basdevant_brdcst = bilinear_Basdevant_brdcst(
        vort_hat,
        (operators, vars)
    )

    bilin_hat_Basdevant_plan_brdcst = similar(vort_hat)
    bilinear_Basdevant_plan_brdcst!(
        bilin_hat_Basdevant_plan_brdcst,
        vort_hat,
        (operators, vars, auxs, plans)
    )

    bilin_hat_Basdevant_plan_loop = similar(vort_hat)
    bilinear_Basdevant_plan_loop!(
        bilin_hat_Basdevant_plan_loop,
        vort_hat,
        (operators, vars, auxs, plans)
    )

    bilin_hat_Basdevant_plan_doubleloop = similar(vort_hat)
    bilinear_Basdevant_plan_doubleloop!(
        bilin_hat_Basdevant_plan_doubleloop,
        vort_hat,
        (operators, vars, auxs, plans)
    )

    @test bilin_hat_naive ≈ bilinear_hat_brdcst
    @test bilinear_hat_brdcst ≈ bilin_hat_plan_brdcst
    @test bilin_hat_plan_brdcst ≈ bilin_hat_Basdevant_brdcst
    @test bilin_hat_Basdevant_brdcst ≈ bilin_hat_Basdevant_plan_brdcst
    @test bilin_hat_Basdevant_plan_brdcst ≈ bilin_hat_Basdevant_plan_loop
    @test bilin_hat_Basdevant_plan_loop ≈ bilin_hat_Basdevant_plan_doubleloop
end
nothing

# ## Orthogonality

# As a second test, we check the orthogonality condition of the bilinear term.

# Due to the periodicity and the divergence-free conditions, we should have $(B(\omega), \omega)_{L^2} = 0$. We compute $B(\omega)$ in spectral space and transform it to physical space in order to check that.

bilin_hat_naive = bilinear_naive(vort_hat, (operators, vars))
bilin_naive = irfft(bilin_hat_naive, N)
nothing

# We compute the inner product by summing over the nodal values and multiplying that by the area element $\Delta x \Delta y = (L/N)^2$. Notice computing the inner product as such has an inherent numerical error, which is also of the order of the area element, so we don't expect the output to be of machine precision.

sum(bilin_naive .* vort) * (L / N)^2

# In spectral space, due to Parseval's identity, this translates to the following (except we need to correct for some indices...)

M = iseven(N) ? div(N, 2) : div(N, 2) + 1

(
    2 * sum(bilin_hat_naive[2:M,:] .* conj.(vort_hat[2:M,:]))
    + sum(bilin_hat_naive[1,:] .* conj.(vort_hat[1,:]))
    + sum(bilin_hat_naive[M+1:end,:] .* conj.(vort_hat[M+1:end,:]))
) / N^2 * (L / N)^2

# Hmm, something is wrong...

# ## Benchmark

@info "bilinear_naive"
@btime bilinear_naive(vh, p) setup = (vh = copy($vort_hat); p = $(operators, vars));

@info "bilinear_auxs!"
@btime bilinear_auxs!(bh, vh, p) setup = (bh = similar($vort_hat); vh = copy($vort_hat); p = $(operators, vars, auxs, plans));

@info "bilinear_Basdevant_brdcst"
@btime bilinear_Basdevant_brdcst(vh, p) setup = (vh = copy($vort_hat); p = $(operators, vars));

@info "bilinear_plan_brdcst!"
@btime bilinear_plan_brdcst!(bh, vh, p) setup = (bh = similar($vort_hat); vh = copy($vort_hat); p = $(operators, vars, auxs, plans));

@info "bilin_hat_Basdevant_plan_brdcst!"
@btime bilinear_Basdevant_plan_brdcst!(bh, vh, p) setup = (bh = similar($vort_hat); vh = copy($vort_hat); p = $(operators, vars, auxs, plans));

@info "bilinear_Basdevant_plan_loop!"
@btime bilinear_Basdevant_plan_loop!(bh, vh, p) setup = (bh = similar($vort_hat); vh = copy($vort_hat); p = $(operators, vars, auxs, plans));

@info "bilinear_Basdevant_plan_doubleloop!"
@btime bilinear_Basdevant_plan_doubleloop!(bh, vh, p) setup = (bh = similar($vort_hat); vh = copy($vort_hat); p = $(operators, vars, auxs, plans));

# ## Conclusions

# Clearly, `bilinear_naive` is pretty bad, in great part because of the use of the allocating transforms `rfft` and `irfft`.

# The method `bilinear_auxs!` saves some allocations but is pretty slow as well, due to also using `rfft` and `irfft`.

# The method `bilinear_Basdevant_brdcst` saves some time by using the Basdevant formulation, but is still quite slow because of `rfft` and `irfft`.

# All the others are much faster, due to using the plans for the transforms and auxiliary functions to avoid allocation of temporary variable.

# Of the nonallocating ones, the ones that use the Basdevant formulation have a clear advantage. Indeed, computing the direct and inverse Fourier transform are the most costly operation here, so reducing it from five to four makes a significant difference.

# Of these nonallocating ones that use the Basdevant formulation, the two that use explicit loops perform better than the broadcast one, with a slight advantage for the one that uses a single loop with `eachindex` instead of a double loop.

# If one does [@code_warntype](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#man-code-warntype-1) on them, one sees that most are clean, but there are some yellow warnings with type `Union{Nothing, Tuple{Int64, Int64}}` precisely on the two fastest ones. Yellow warnings, fortunately, tend to be harmless, as is the case here. In this case, the warning is due to the way the [iterator interface](https://docs.julialang.org/en/v1/manual/interfaces/#Interfaces) is implemented, which returns `nothing` when the end of the iteration is reached. One can get rid of it with `while` loops, but it won't improve the performance.


text:
# Data source for london tube https://www.whatdotheyknow.com/request/distance_between_adjacent_underg
using CSV, DataFrames

df = CSV.read("data/london/london.csv", DataFrame)

id = 1
stations = Dict()
edges = Dict()
for row in Tables.namedtupleiterator(df)
    station_a = getindex(row, Symbol("Station from (A)"))
    if !haskey(stations, station_a)
        push!(stations, station_a => id)
        global id += 1
    end
    station_b = getindex(row, Symbol("Station to (B)"))
    if !haskey(stations, station_b)
        push!(stations, station_b => id)
        global id += 1
    end

    (from,to) = (stations[station_a], stations[station_b])
    distance = getindex(row, Symbol("Distance (Kms)"))
    if !haskey(edges, (from,to)) && !haskey(edges, (to,from))
        push!(edges, (from,to) => round(Int, distance * 1000)) # in meters
    end
end

n, E, C = id - 1, collect(keys(edges)), collect(values(edges))

output = open("data/london/london.gr", "w")
write(output, "c London Tube Stations Graph with distances\n")
write(output, "c Source: https://www.whatdotheyknow.com/request/distance_between_adjacent_underg\n")
write(output, "c\n")
write(output, "p st ", string(n), " ", string(length(E)), "\n")
write(output, "c\n")
for (i, (a, b)) in enumerate(E)
    write(output, "a ", string(a), " ", string(b), " ", string(C[i]), "\n")
end
close(output)


text:
module ImageText

using Images
using Crayons

key=[" ", ".", ":", "-", "=", "+", "*", "#", "%", "@"]

color_mat = [
0  0  0  0
1  128  0  0
2  0  128  0
3  128  128  0
4  0  0  128
5  128  0  128
6  0  128  128
7  192  192  192
8  128  128  128
9  255  0  0
10  0  255  0
11  255  255  0
12  0  0  255
13  255  0  255
14  0  255  255
15  255  255  255
16  0  0  0
17  0  0  95
18  0  0  135
19  0  0  175
20  0  0  215
21  0  0  255
22  0  95  0
23  0  95  95
24  0  95  135
25  0  95  175
26  0  95  215
27  0  95  255
28  0  135  0
29  0  135  95
30  0  135  135
31  0  135  175
32  0  135  215
33  0  135  255
34  0  175  0
35  0  175  95
36  0  175  135
37  0  175  175
38  0  175  215
39  0  175  255
40  0  215  0
41  0  215  95
42  0  215  135
43  0  215  175
44  0  215  215
45  0  215  255
46  0  255  0
47  0  255  95
48  0  255  135
49  0  255  175
50  0  255  215
51  0  255  255
52  95  0  0
53  95  0  95
54  95  0  135
55  95  0  175
56  95  0  215
57  95  0  255
58  95  95  0
59  95  95  95
60  95  95  135
61  95  95  175
62  95  95  215
63  95  95  255
64  95  135  0
65  95  135  95
66  95  135  135
67  95  135  175
68  95  135  215
69  95  135  255
70  95  175  0
71  95  175  95
72  95  175  135
73  95  175  175
74  95  175  215
75  95  175  255
76  95  215  0
77  95  215  95
78  95  215  135
79  95  215  175
80  95  215  215
81  95  215  255
82  95  255  0
83  95  255  95
84  95  255  135
85  95  255  175
86  95  255  215
87  95  255  255
88  135  0  0
89  135  0  95
90  135  0  135
91  135  0  175
92  135  0  215
93  135  0  255
94  135  95  0
95  135  95  95
96  135  95  135
97  135  95  175
98  135  95  215
99  135  95  255
100  135  135  0
101  135  135  95
102  135  135  135
103  135  135  175
104  135  135  215
105  135  135  255
106  135  175  0
107  135  175  95
108  135  175  135
109  135  175  175
110  135  175  215
111  135  175  255
112  135  215  0
113  135  215  95
114  135  215  135
115  135  215  175
116  135  215  215
117  135  215  255
118  135  255  0
119  135  255  95
120  135  255  135
121  135  255  175
122  135  255  215
123  135  255  255
124  175  0  0
125  175  0  95
126  175  0  135
127  175  0  175
128  175  0  215
129  175  0  255
130  175  95  0
131  175  95  95
132  175  95  135
133  175  95  175
134  175  95  215
135  175  95  255
136  175  135  0
137  175  135  95
138  175  135  135
139  175  135  175
140  175  135  215
141  175  135  255
142  175  175  0
143  175  175  95
144  175  175  135
145  175  175  175
146  175  175  215
147  175  175  255
148  175  215  0
149  175  215  95
150  175  215  135
151  175  215  175
152  175  215  215
153  175  215  255
154  175  255  0
155  175  255  95
156  175  255  135
157  175  255  175
158  175  255  215
159  175  255  255
160  215  0  0
161  215  0  95
162  215  0  135
163  215  0  175
164  215  0  215
165  215  0  255
166  215  95  0
167  215  95  95
168  215  95  135
169  215  95  175
170  215  95  215
171  215  95  255
172  215  135  0
173  215  135  95
174  215  135  135
175  215  135  175
176  215  135  215
177  215  135  255
178  215  175  0
179  215  175  95
180  215  175  135
181  215  175  175
182  215  175  215
183  215  175  255
184  215  215  0
185  215  215  95
186  215  215  135
187  215  215  175
188  215  215  215
189  215  215  255
190  215  255  0
191  215  255  95
192  215  255  135
193  215  255  175
194  215  255  215
195  215  255  255
196  255  0  0
197  255  0  95
198  255  0  135
199  255  0  175
200  255  0  215
201  255  0  255
202  255  95  0
203  255  95  95
204  255  95  135
205  255  95  175
206  255  95  215
207  255  95  255
208  255  135  0
209  255  135  95
210  255  135  135
211  255  135  175
212  255  135  215
213  255  135  255
214  255  175  0
215  255  175  95
216  255  175  135
217  255  175  175
218  255  175  215
219  255  175  255
220  255  215  0
221  255  215  95
222  255  215  135
223  255  215  175
224  255  215  215
225  255  215  255
226  255  255  0
227  255  255  95
228  255  255  135
229  255  255  175
230  255  255  215
231  255  255  255
232  8  8  8
233  18  18  18
234  28  28  28
235  38  38  38
236  48  48  48
237  58  58  58
238  68  68  68
239  78  78  78
240  88  88  88
241  98  98  98
242  108  108  108
243  118  118  118
244  128  128  128
245  138  138  138
246  148  148  148
247  158  158  158
248  168  168  168
249  178  178  178
250  188  188  188
251  198  198  198
252  208  208  208
253  218  218  218
254  228  228  228
255  238  238  238
]

function float_to_char(n; reverse=true, a = 0.0, b = 1.0)
    if reverse
        n=b-n
    end

    n = (n - a) / (b - a)

    if n < 0.05
        return key[1]
    elseif n < 0.1
        return key[2]
    elseif n < 0.15
        return key[3]
    elseif n < 0.22
        return key[4]
    elseif n < 0.33
        return key[5]
    elseif n < 0.55
        return key[6]
    elseif n < 0.7
        return key[7]
    elseif n < 0.8
        return key[8]
    elseif n < 0.9
        return key[9]
    elseif n <= 1.0
        return key[10]
    else
        println("outofrange $n should be zero to 1")
        return "X"
    end
    return "X"
end

function get_color(rgb)

    r = Int64(round(rgb.r * 256))
    g = Int64(round(rgb.g * 256))
    b = Int64(round(rgb.b * 256))

    best = 256
    match = 100000
    for i = 1:256
        t = abs(r - color_mat[i,2]) + abs(g - color_mat[i,3])  + abs(b - color_mat[i,4])
        if t < match
            match = t
            best = i
        end
    end

    return best-1

end

function go(img; blackwhite=false, reverse=false, background=missing, text=:color, rescale=true)

    gray_image = Gray.(img)

    if rescale
        arr = Float64.(gray_image)
        a = minimum(arr)
        b = maximum(arr)
    else
        a = 0.0
        b = 1.0
    end

    nx,ny = size(gray_image)

    str = []

    for x = 1:nx
        vs = String[]
        for y = 1:ny
            t = float_to_char(gray_image[x,y], reverse=reverse, a=a, b = b )
            push!(vs, t)
            color = get_color(img[x,y])
            if !blackwhite
                #                printstyled(t, color=color)

                if text == :color
                    forecolor=color
                else
                    forecolor=text
                end

                if ismissing(background)
                    backcolor=missing
                elseif background == :color
                    backcolor=color
                else
                    backcolor=background
                end

#                println("xx ", [forecolor, backcolor])
                if ismissing(background)
                    printstyled(t, color=forecolor)
                else
                    print(Crayon(foreground = forecolor, background=backcolor), t)
                end
            end
        end
        push!(str, vs)
        if blackwhite
            println(join(vs))
        else
            if ismissing(background)
                println()
            else
                print(Crayon(background=:white), " \n")
            end
#            println()
        end

    end

    return str

end

end  #end module


text:
function add_constrs!(model::Model, cbegins::IVec, inds::IVec, coeffs::FVec, rel::CVec, rhs::FVec)
    nnz   = length(inds)
    ncons = length(rhs)
    (nnz == length(coeffs)) || error("Incompatible constraint argument dimensions.")

    for k in 1:length(rel)
      if rel[k] == '>'
        rel[k] = convert(Cchar, 'G')
      elseif rel[k] == '<'
        rel[k] = convert(Cchar, 'L')
      elseif rel[k] == '='
        rel[k] = convert(Cchar, 'E')
      end
    end

    if ncons > 0
        stat = @cpx_ccall(addrows, Cint, (
                          Ptr{Void},        # environment
                          Ptr{Void},        # problem
                          Cint,             # num new cols
                          Cint,             # num new rows
                          Cint,             # num non-zeros
                          Ptr{Cdouble},     # rhs
                          Ptr{Cchar},       # sense
                          Ptr{Cint},        # matrix start
                          Ptr{Cint},        # matrix index
                          Ptr{Cdouble},     # matrix values
                          Ptr{Ptr{Cchar}},  # col names
                          Ptr{Ptr{Cchar}}   # row names
                          ),
                          model.env.ptr, model.lp, 0, ncons, nnz, rhs, rel, cbegins.-1, inds.-1, coeffs, C_NULL, C_NULL)

        if stat != 0
           throw(CplexError(model.env, stat))
        end
    end
end

add_constr!(model::Model, coef::Vector, sense::Char, rhs) = add_constrs!(model, [1], collect(1:length(coef)), coef, [sense], [rhs])

function add_constrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, rel::Vector, rhs::Vector)
    add_constrs!(model, ivec(cbeg), ivec(inds), fvec(coeffs), cvecx(rel, length(cbeg)), fvec(rhs))
end

function add_constrs_t!(model::Model, At::SparseMatrixCSC{Float64}, rel::GCharOrVec, b::Vector)
    n, m = size(At)
    (m == length(b) && n == num_var(model)) || error("Incompatible argument dimensions.")
    add_constrs!(model, At.colptr[1:At.n], At.rowval, At.nzval, rel, b)
end

function add_constrs_t!(model::Model, At::Matrix{Float64}, rel::GCharOrVec, b::Vector)
    n, m = size(At)
    (m == length(b) && n == num_var(model)) || error("Incompatible argument dimensions.")
    add_constrs_t!(model, sparse(At), rel, b)
end

function add_constrs!(model::Model, A::CoeffMat, rel::GCharOrVec, b::Vector{Float64})
    m, n = size(A)
    (m == length(b) && n == num_var(model)) || error("Incompatible argument dimensions.")
    add_constrs_t!(model, transpose(A), rel, b)
end

add_constrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, rel::Char, rhs::Vector) = add_constrs!(model, cbeg, inds, coeffs, cvecx(rel, length(cbeg)), rhs)

function add_rangeconstrs!(model::Model, cbegins::IVec, inds::IVec, coeffs::FVec, lb::FVec, ub::FVec)
    nnz   = length(inds)
    ncons = length(lb)
    (ncons  == length(ub) && nnz == length(coeffs)) || error("Incompatible constraint argument dimensions.")

    sense = fill!(Array(Cchar, ncons), 'R')

    for i in 1:ncons
        if lb[i] == -Inf
            lb[i] = -ub[i]
            ub[i] = Inf
            coeffs[cbegins[i]+1:cbegins[i+1]] = -coeffs[cbegins[i]+1:cbegins[i+1]]
        end
    end

    if ncons > 0
        stat = @cpx_ccall(addrows, Cint, (
                          Ptr{Void},        # environment
                          Ptr{Void},        # problem
                          Cint,             # num new cols
                          Cint,             # num new rows
                          Cint,             # num non-zeros
                          Ptr{Cdouble},     # rhs
                          Ptr{Cchar},       # sense
                          Ptr{Cint},        # matrix start
                          Ptr{Cint},        # matrix index
                          Ptr{Cdouble},     # matrix values
                          Ptr{Void},        # col names
                          Ptr{Void}         # row names
                          ),
                          model.env.ptr, model.lp, 0, ncons, nnz, lb, sense, cbegins[1:end-1], inds, coeffs, C_NULL, C_NULL)
        if stat != 0
            throw(CplexError(model.env, stat))
        end
        for i = 1:ncons
            stat = @cpx_ccall(chgrngval, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Ptr{Cint},
                      Ptr{Cdouble}
                      ),
                      model.env.ptr, model.lp, 1, Cint[i-1], [ub[i]-lb[i]])
            if stat != 0
                throw(CplexError(model.env, stat))
            end
        end
    end
end

function add_rangeconstrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, lb::Vector, ub::Vector)
    add_rangeconstrs!(model, ivec(cbeg), ivec(inds), fvec(coeffs), fvec(lb), fvec(ub))
end

function add_rangeconstrs_t!(model::Model, At::SparseMatrixCSC{Float64}, lb::Vector, ub::Vector)
    add_rangeconstrs!(model, At.colptr-1, At.rowval-1, At.nzval, lb, ub)
end

function add_rangeconstrs_t!(model::Model, At::Matrix{Float64}, lb::Vector, ub::Vector)
    add_rangeconstrs_t!(model, sparse(At), lb, ub)
end

function add_rangeconstrs!(model::Model, A::CoeffMat, lb::Vector, ub::Vector)
    m, n = size(A)
    (m == length(lb) == length(ub) && n == num_var(model)) || error("Incompatible constraint argument dimensions.")
    add_rangeconstrs_t!(model, transpose(A), lb, ub)
end

function num_constr(model::Model)
    ncons = @cpx_ccall(getnumrows, Cint, (
                       Ptr{Void},
                       Ptr{Void}
                       ),
                       model.env.ptr, model.lp)
    return ncons
end

function get_constr_senses(model::Model)
    ncons = num_constr(model)
    senses = Array(Cchar, ncons)
    stat = @cpx_ccall(getsense, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Ptr{Cchar},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, senses, 0, ncons-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return senses
end

function set_constr_senses!(model::Model, senses::Vector)
    ncons = num_constr(model)
    stat = @cpx_ccall(chgsense, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Ptr{Cint},
                      Ptr{Cchar}
                      ),
                      model.env.ptr, model.lp, ncons, Cint[0:ncons-1;], convert(Vector{Cchar},senses))
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function get_rhs(model::Model)
    ncons = num_constr(model)
    rhs = Array(Cdouble, ncons)
    stat = @cpx_ccall(getrhs, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, rhs, 0, ncons-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return rhs
end

function set_rhs!(model::Model, rhs::Vector)
    ncons = num_constr(model)
    @assert ncons == length(rhs)
    stat = @cpx_ccall(chgrhs, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Ptr{Cint},
                      Ptr{Cdouble}
                      ),
                      model.env.ptr, model.lp, ncons, Cint[0:ncons-1;], float(rhs))
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function get_constrLB(model::Model)
    senses = get_constr_senses(model)
    ret    = get_rhs(model)
    for i = 1:num_constr(model)
        if senses[i] == 'G' || senses[i] == 'E'
            # Do nothing
        else
            # LEQ constraint so LB is -Inf
            ret[i] = -Inf
        end
    end
    return ret
end

function get_constrUB(model::Model)
    senses = get_constr_senses(model)
    ret    = get_rhs(model)
    for i = 1:num_constr(model)
        if senses[i] == 'L' || senses[i] == 'E'
            # Do nothing
        else
            # GEQ constraint so UB is Inf
            ret[i] = +Inf
        end
    end
    return ret
end

function set_constrLB!(model::Model, lb)
    senses = get_constr_senses(model)
    rhs    = get_rhs(model)
    sense_changed = false
    for i = 1:num_constr(model)
        if senses[i] == 'G' || senses[i] == 'E'
            # Do nothing
        elseif senses[i] == 'L'
            if lb[i] != -Inf
                # LEQ constraint with non-NegInf LB implies a range
                if isapprox(lb[i], rhs[i])
                    # seems to be an equality
                    senses[i] = 'E'
                    sense_changed = true
                else
                    error("Tried to set LB != -Inf on a LEQ constraint (index $i)")
                end
            else
                lb[i] = rhs[i]
            end
        end
    end
    if sense_changed
        set_constr_senses!(model, senses)
    end
    set_rhs!(model, lb)
end

function set_constrUB!(model::Model, ub)
    senses = get_constr_senses(model)
    rhs    = get_rhs(model)
    sense_changed = false
    for i = 1:num_constr(model)
        if senses[i] == 'L' || senses[i] == 'E'
            # Do nothing
        elseif senses[i] == 'G'
            if ub[i] != Inf
                # GEQ constraint with non-PosInf UB implies a range
                if isapprox(ub[i], rhs[i])
                    # seems to be an equality
                    senses[i] = 'E'
                    sense_changed = true
                else
                    error("Tried to set UB != +Inf on a GEQ constraint (index $i)")
                end
            else
              ub[i] = rhs[i]
            end
        end
    end
    if sense_changed
        set_constr_senses!(model, senses)
    end
    set_rhs!(model, ub)
end

function get_nnz(model::Model)
  ret = @cpx_ccall(getnumnz, Cint, (Ptr{Void}, Ptr{Void}), model.env.ptr, model.lp)
  ret == 0 && throw(error("Could not query number of nonzeros in model"))
  return ret
end

function get_constr_matrix(model::Model)
  nzcnt_p = Array(Cint, 1)
  m = num_constr(model)
  n = num_var(model)
  nnz = get_nnz(model)
  cmatbeg = Array(Cint, n+1)
  cmatind = Array(Cint, nnz)
  cmatval = Array(Cdouble, nnz)
  surplus_p = Array(Cint, 1)
  stat = @cpx_ccall(getcols, Cint, (
                    Ptr{Void},
                    Ptr{Void},
                    Ptr{Cint},
                    Ptr{Cint},
                    Ptr{Cint},
                    Ptr{Cdouble},
                    Cint,
                    Ptr{Cint},
                    Cint,
                    Cint
                    ),
                    model.env.ptr, model.lp, nzcnt_p, cmatbeg, cmatind, cmatval, nnz, surplus_p, 0, n-1)
  if stat != 0 || surplus_p[1] < 0
    throw(CplexError(model.env, stat))
  end
  cmatbeg[end] = nnz # add the last entry that Julia wants
  return SparseMatrixCSC(m, n, convert(Vector{Int64}, cmatbeg.+1), convert(Vector{Int64}, cmatind.+1), cmatval)
end

get_num_sos(model::Model) = @cpx_ccall(getnumsos, Cint, (Ptr{Void}, Ptr{Void}), model.env.ptr, model.lp)

# int CPXaddsos( CPXCENVptr env, CPXLPptr lp, int numsos, int numsosnz, char const * sostype, int const * sosbeg, int const * sosind, double const * soswt, char ** sosname )
function add_sos!(model::Model, sostype::Symbol, idx::Vector{Int}, weight::Vector{Cdouble})
    ((nelem = length(idx)) == length(weight)) || error("Index and weight vectors of unequal length")
    (sostype == :SOS1) ? (typ = CPX_TYPE_SOS1) : ((sostype == :SOS2) ? (typ = CPX_TYPE_SOS2) : error("Invalid SOS constraint type"))
    stat = @cpx_ccall(addsos, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Cint,
                      Ptr{Cchar},
                      Ptr{Cint},
                      Ptr{Cint},
                      Ptr{Cdouble},
                      Ptr{Ptr{Cchar}}
                      ),
                      model.env.ptr, model.lp, convert(Cint, 1), convert(Cint, nelem), [convert(Cchar, typ)], Cint[0], convert(Vector{Cint}, idx.-1), weight, C_NULL)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    model.has_int = true
    return nothing
end

add_indicator_constraint(model::Model, idx, coeff, sense, rhs, indicator) =
    add_indicator_constraint(model::Model, idx, coeff, sense, rhs, indicator, 0)
add_indicator_constraint(model::Model, idx, coeff, sense, rhs, indicator, comp) =
    add_indicator_constraint(model, convert(Vector{Cint},idx), convert(Vector{Cdouble},coeff),
                             convert(Cchar,sense), convert(Cdouble,rhs), convert(Cint,indicator), convert(Cint,comp))
function add_indicator_constraint(model::Model, idx::Vector{Cint}, coeff::Vector{Cdouble}, sense::Cchar, rhs::Cdouble, indicator::Cint, comp::Cint)
    (nzcnt = length(idx)) == length(coeff) || error("Incompatible lengths in constraint specification")
    stat = @cpx_ccall(addindconstr, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Cint,
                      Cint,
                      Cdouble,
                      Cchar,
                      Ptr{Cint},
                      Ptr{Cdouble},
                      Ptr{Cchar}),
                      model.env.ptr, model.lp, indicator, comp,
                      nzcnt, rhs, sense, idx.-1, coeff, C_NULL)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return nothing
end


text:
##check whether two ranges intersect
#function intersectrng(r1::UnitRange, r2::UnitRange)
#    i = intersect(r1)
#    return first(i)<=last(i)
#end
##check whether all corresponding ranges in two lists of ranges intersect
#intersectrng(t1::Tuple, t2::Tuple) = intersectrng(first(t1), first(t2)) && intersectrng(Base.tail(t1), Base.tail(t2))
#intersectrng(t1::Tuple{}, t2::Tuple{}) = true

crop_k_inds(img::AbstractArray, k::Int) = crop_k_inds(axes(img), k)
crop_k_inds(inds::Tuple, k::Int) = ((first(first(inds))+k):(last(first(inds))-k), crop_k_inds(Base.tail(inds), k)...)
crop_k_inds(inds::Tuple{}, k::Int) = ()

#negranges(t::Tuple) = (-first(t):first(t), negranges(Base.tail(t))...)
roiranges(r::Tuple) = roiranges(ntuple(i->zero(first(r)),length(r)), r) #assume centered at origin
roiranges(ctr::Tuple, r::Tuple) = ((first(ctr)-first(r)):(first(ctr)+first(r)), roiranges(Base.tail(ctr), Base.tail(r))...)
roiranges(ctr::Tuple{}, r::Tuple{}) = ()

#return iterator over all indices adjacent to I
neighbor_inds(I::CartesianIndex) = filter(x->x!=I, CartesianIndices(roiranges(I.I, ntuple(i->1, length(I)))))

function find_local_maxima(img::AbstractArray, thresh)
    inds = CartesianIndices(crop_k_inds(img, 1))
    maxima = typeof(first(inds).I)[]
    for I in inds
        v = img[I]
        if v >= thresh
            is_maximum = true
            for i in neighbor_inds(I)
                if v <= img[i]
                    is_maximum = false
                    break
                end
            end
            if is_maximum
                push!(maxima, I.I)
            end
        end
    end
    return maxima
end

#for each local intensity maximum assign an integer index and paint a canvas with it
#if while painting the canvas we encounter an already-painted pixel then this roi overlaps
#...so we exclude it from the list of found beads
function exclude_overlapping(roi_ctrs, roi_radii, imgsz)
    canvas = zeros(Int, imgsz)
    keep = trues(length(roi_ctrs))
    for (i,c) in enumerate(roi_ctrs)
        #printed = false
        for I in CartesianIndices(roiranges(c, roi_radii))
            if checkbounds(Bool, canvas, I) #could be more efficient
                v = canvas[I]
                if v!=0
                    #if !printed
                        #print("overlap between roi centered at $c and roi centered at $(roi_ctrs[v])\n")
                        #print("radius was $roi_radii \n")
                        #printed=true
                    #end
                    keep[i] = false
                else
                    canvas[I] = i
                end
            else #the roi extends outside of the image
                keep[i] = false
            end
        end
    end
    return keep
end

#find beads that do not overlap within radiii roi_radii and return a list of rois (views into parent image)
#along with a list of coordinate centers in parent image
function find_beads(fullimg::AbstractArray{T,N}, roi_radii::NTuple{N,Int}, bias; noise_stds=5, sat_val=typemax(T)) where {T,N}
    #estimate noise level in black pixels, assuming at least 80% of pixels are "black"
    psorted = sort(reshape(fullimg, length(fullimg)))
    noise = std(view(psorted, 1:round(Int, 0.8*length(psorted))))
    #blur the image to eliminate multiple local maxima from nearby pixels of same bead
    kern = KernelFactors.gaussian(div.(roi_radii,2))
    kernmax = maximum(first(kern))
    imgf = imfilter(fullimg, kern)
    #find all local intensity maxima that are above a noise threshold (set by estimated noise)
    thresh = bias+kernmax*noise*noise_stds
    maxima = find_local_maxima(imgf, thresh)
    keep = exclude_overlapping(maxima, roi_radii, size(fullimg))
    if all(keep.==false)
        @warn("No valid bead ROIs found.  Common reasons for this are that the bead sample was too dense, or that the image is too noisy.")
        return [],[]
    else
        mkeep = maxima[keep]
        firstroi = view(fullimg, roiranges(mkeep[1], roi_radii)...)
        firstroi = OffsetArray(firstroi, roiranges(roi_radii)...)
        rois = typeof(firstroi)[]
        keep2 = trues(length(mkeep))
        for (i,c) in enumerate(mkeep)
            roi = view(fullimg, roiranges(c, roi_radii)...)
            roi = OffsetArray(roi, roiranges(roi_radii)...)
            if all(roi.<sat_val)
                push!(rois, roi)
            else
                #print("Excluding ROI centered at $c due to saturation\n")
                keep2[i] = false
            end
        end
        if all(keep2.==false)
            @warn("After removing ROIs with saturated pixels no valid bead ROIs remained.  Use less intense illumination or a shorter exposure time.")
            return [], []
        end
        return rois, mkeep[keep2]
    end
end

max_forks_repo_name: HolyLab/ExtractPSF


text:
export hessian_operator
"""
Return the value of the hessian of a multivariate function as a LinearOperator.
For more information on LinearOperators, see the doc for the package of the
same name.
"""
function hessian_operator(nlp,x)
    n = nlp.meta.nvar
    # temp = Array(Float64, n)
    # temp = Array{Float64}(n)
    temp = Array{Float64}(undef, n)
    return hess_op!(nlp, x, temp)
    #return LinearOperator(nlp.meta.nvar, nlp.meta.nvar, true, true, v -> hprod(nlp,x,v))
end