text:
module Weighting
#
# This module has routines and target constants to re-weight our main dataset so it hits current or future totals for employment,
# population by age and so on. See:
#
# Creedy, John. 2003. ‘Survey Reweighting for Tax Microsimulation Modelling’. Treasury Working Paper Series 03/17.
# New Zealand Treasury. http://ideas.repec.org/p/nzt/nztwps/03-17.html.
#
# for an overview of how this works, and the `SurveyDataWeighting` module for implementation details.
#
# The targets we currently weight for are:
#
# * Employment and Unemployment (from NOMIS)
# * Tenure Type (from Scotgov)
# * Household Type (from NRS)
# * Receipts of disability/caring benefits (Stat-Explore)
# * population in 5- year age bands
# * household totals by local authority.
#
# 82 targets in all, presently.
#
using DataFrames

using SurveyDataWeighting:
    DistanceFunctionType,
    ITERATIONS_EXCEEDED,
    chi_square,
    constrained_chi_square,
    d_and_s_constrained,
    d_and_s_type_a,
    d_and_s_type_b,
    do_reweighting

using ScottishTaxBenefitModel
using .ModelHousehold
using .Definitions

export
    DEFAULT_TARGETS,
    generate_weights,
    initialise_target_dataframe,
    make_target_dataset

# FIXME rewrite this to load from a file.

const NUM_HOUSEHOLDS = 2_477_000.0 # sum of all hhld types below

#
# This is our default set of targets for 2021/2.
#
# See `weighting_target_set_creation.md`
# and `target_generation_worksheet-aug-22-2021.ods`
#
const DEFAULT_TARGETS_2021 = [
    # Employment totals from NOMIS
    # We drop 'inactive' because of collinearity with NRS population
    1.02361127107591*1_330_149,  #	1	M- Total in employment- aged 16+ - the constant here scales the NOMIS over 16 populations to the NRS totals
    1.02361127107591*68_308,	 #	2	M- Total unemployed- aged 16+
    1.01039881362443*1_304_812,  #	3	F- Total in employment- aged 16+
    1.01039881362443*50_481,     #	4	F- Total unemployed- aged 16+
    # Tenure - the constant scales SGov occupied households counts to NRS total household counts.
    # OOs dropped for colinearity with NRS household counts.
    1.00987684414087*370_845,   #	5	private rented+rent free
    1.00987684414087*280_715,   #	6	housing association
    1.00987684414087*314_433,   #	7	las etc rented
    # NRS popn by age 2021/2
	135_959,	#	8	M – 0 - 4
	152_847,	#	9	5 - 9
	151_875,	#	10	10 - 14
	144_207,	#	11	15 - 19
	173_302,	#	12	20 - 24
	189_139,	#	13	25 - 29
	185_637,	#	14	 30 - 34
	174_079,	#	15	 35 - 39
	159_586,	#	16	40 - 44
	169_376,	#	17	 45 - 49
	189_355,	#	18	 50 - 54
	193_348,	#	19	 55 - 59
	170_701,	#	20	 60 - 64
	144_529,	#	21	 65 - 69
	135_910,	#	22	 70 - 74
	89_206,  #	23	 75 - 79
	106_156,	#	24	80+
	127_847,	#	25	F – 0 - 4
	145_056,	#	26	5 - 9
	146_206,	#	27	10 - 14
	137_913,	#	28	15 - 19
	168_453,	#	29	20 - 24
	188_065,	#	30	25 - 29
	188_432,	#	31	 30 - 34
	181_587,	#	32	 35 - 39
	164_780,	#	33	40 - 44
	180_548,	#	34	 45 - 49
	203_758,	#	35	 50 - 54
	205_996,	#	36	 55 - 59
	181_868,	#	37	 60 - 64
	155_904,	#	38	 65 - 69
	149_920,	#	39	 70 - 74
	109_004,	#	40	 75 - 79
	165_451,	#	41	80+
    # NRS households
	468_147,	#	42	 # 42 - 1 adult: male
	453_675,	#	43	 # 43 - 1 adult: female
	795_465,	#	44	 # 44 - 2 adults
	90_276,	#	45	 # 45 - 1 adult 1 child
	67_324,	#	46	 # 46 - 1 adult 2+ children
	440_062,	#	47	 # 47 - 2+ adults 1+ children
	208_147,	#	48	 # 48 - 3+ adults
    # Disability benefit receipts from Stat-Explore
	82_031,	#	49	CARERS
	124_192,	#	50	AA
	432_744,	#	51	PIP/DLA
    # Household (not popn) from NRS by LA. Aberdeen city dropped for collinearity.
	113_217,	#	52	 # S12000034 - 52 Aberdeenshire
	54_378,	#	53	 # S12000041 - Angus
	41_635,	#	54	 # S12000035 - Argyll and Bute
	243_954,	#	55	 # S12000036 - City of Edinburgh
	24_039,	#	56	 # S12000005 - Clackmannanshire
	69_824,	#	57	 # S12000006 - Dumfries and Galloway
	71_077,	#	58	 # S12000042 - Dundee City
	55_642,	#	59	 # S12000008 - East Ayrshire
	46_627,	#	60	 # S12000045 - East Dunbartonshire
	47_707,	#	61	 # S12000010 - East Lothian
	39_978,	#	62	 # S12000011 - East Renfrewshire
	73_621,	#	63	 # S12000014 - Falkirk
	170_311,	#	64	 # S12000047 - Fife
	299_004,	#	65	 # S12000049 - Glasgow City
	110_436,	#	66	 # S12000017 - Highland
	37_462,	#	67	 # S12000018 - Inverclyde
	40_993,	#	68	 # S12000019 - Midlothian
	43_406,	#	69	 # S12000020 - Moray
	12_798,	#	70	 # S12000013 - Na h-Eileanan Siar
	64_201,	#	71	 # S12000021 - North Ayrshire
	153_918,	#	72	 # S12000050 - North Lanarkshire
	10_701,	#	73	 # S12000023 - Orkney Islands
	69_706,	#	74	 # S12000048 - Perth and KinroS
	87_944,	#	75	 # S12000038 - Renfrewshire
	55_090,	#	76	 # S12000026 - Scottish Borders
	10_503,	#	77	 # S12000027 - Shetland Islands
	52_804,	#	78	 # S12000028 - South Ayrshire
	149_064,	#	79	 # S12000029 - South Lanarkshire
	40_383,	#	80	 # S12000030 - Stirling
	43_215,	#	81	 # S12000039 - West Dunbartonshire
	80_562	#	82	 # S12000040 - West Lothian
    ]

#
#
# see data/targets/target_set.ods for derivation of these targets.
#
const DEFAULT_TARGETS_2020 = [
    # 1.01337880471403 is nrs 16+ popn estimate/nomis count
    # this scales to employment and popn numbers are consistent
    # we drop 'inactive' because of collinearity
    1_340_609.0*1.01337880471403, # 1 - M- Total in employment- aged 16+
    60_635*1.01337880471403, # 2 - M- Total unemployed- aged 16+
    # 745_379, # 3 - M- Total economically inactive- aged 16+
    1_301_248*1.01337880471403, # 3 - F- Total in employment- aged 16+
    59_302*1.01337880471403, # 4 - F- Total unemployed- aged 16+
    # 974_767, # 6 - F- Total economically inactive- aged 16+
    # 1_540_890, # 7 - owner occupied
    # 0.987518476103907 HERE is ratio of total hhls
    # from tenure data and nrs hhld type counts
    370_502*0.987518476103907, # 5 - private rented+rent free
    282_482*0.987518476103907, # 6 - housing association
    314_433*0.987518476103907, # 7 - las etc rented
    139_982, # 8 - M- 0 - 4
    153_297, # 9 - M- 5 - 9
    150_487, # 10 - M- 10 – 14
    144_172, # 11 - M- 15 - 19
    176_066, # 12 - M- 20 - 24
    191_145, # 13 - M- 25 - 29
    182_635, # 14 -  M- 30 - 34
    172_624, # 15 -  M- 35 - 39
    156_790, # 16 - M- 40 - 44
    174_812, # 17 -  M- 45 - 49
    193_940, # 18 -  M- 50 - 54
    190_775, # 19 -  M- 55 - 59
    166_852, # 20 -  M- 60 - 64
    144_460, # 21 -  M- 65 - 69
    132_339, # 22 -  M- 70 - 74
    87_886, # 23 -  M- 75 - 79
    104_741, # 24 - M- 80’+
    131_733, # 25 - F- 0 - 4
    146_019, # 26 - F- 5 - 9
    144_187, # 27 - F- 10 - 14
    137_786, # 28 - F- 15 - 19
    171_390, # 29 - F- 20 - 24
    191_110, # 30 - F- 25 - 29
    186_828, # 31 -  F- 30 - 34
    179_898, # 32 -  F- 35 - 39
    162_642, # 33 - F- 40 - 44
    186_646, # 34 -  F- 45 - 49
    207_150, # 35 -  F- 50 - 54
    202_348, # 36 -  F- 55 - 59
    177_841, # 37 -  F- 60 - 64
    154_984, # 38 -  F- 65 - 69
    146_517, # 39 -  F- 70 - 74
    108_065, # 40 -  F- 75 - 79
    165_153, # 41 -  F- 80+
    439_000, # 42 - 1 adult: male
    467_000, # 43 - 1 adult: female
    797_000, # 44 - 2 adults
    70_000, # 45 - 1 adult 1 child
    66_000, # 46 - 1 adult 2+ children
    448_000, # 47 - 2+ adults 1+ children
    190_000, # 48 - 3+ adults
    77_842, # 49 - CARER’S ALLOWANCE
    127_307, # 50 - AA
    431_461, # 51 pip/dla
    # council areas
    # FIXME note these sum to 0.992538132778121 of the household type totals above: Check Again
    # 108381, # S12000033 - Aberdeen City  # removed for collinearity
    112114, # S12000034 - 52 Aberdeenshire
    54221, # S12000041 - Angus
    41789, # S12000035 - Argyll and Bute
    238269, # S12000036 - City of Edinburgh
    23890, # S12000005 - Clackmannanshire
    69699, # S12000006 - Dumfries and Galloway
    70685, # S12000042 - Dundee City
    55387, # S12000008 - East Ayrshire
    46228, # S12000045 - East Dunbartonshire
    46771, # S12000010 - East Lothian
    39345, # S12000011 - East Renfrewshire
    72672, # S12000014 - Falkirk
    169239, # S12000047 - Fife
    294622, # S12000049 - Glasgow City
    109514, # S12000017 - Highland
    37614, # S12000018 - Inverclyde
    39733, # S12000019 - Midlothian
    42932, # S12000020 - Moray
    12833, # S12000013 - Na h-Eileanan Siar
    64140, # S12000021 - North Ayrshire
    152443, # S12000050 - North Lanarkshire
    10589, # S12000023 - Orkney Islands
    69003, # S12000048 - Perth and KinroS
    86683, # S12000038 - Renfrewshire
    54715, # S12000026 - Scottish Borders
    10439, # S12000027 - Shetland Islands
    52588, # S12000028 - South Ayrshire
    147434, # S12000029 - South Lanarkshire
    39654, # S12000030 - Stirling
    43030, # S12000039 - West Dunbartonshire
    78966 ] # S12000040 - West Lothian

const DEFAULT_TARGETS = DEFAULT_TARGETS_2021

function initialise_target_dataframe( n :: Integer ) :: DataFrame
    df = DataFrame(
        m_total_in_employment = zeros(n),
        m_total_unemployed = zeros(n),
        # m_total_economically_inactive = zeros(n),
        f_total_in_employment = zeros(n),
        f_total_unemployed = zeros(n),
        # f_total_economically_inactive = zeros(n),
        # owner_occupied = zeros(n),
        private_rented_plus_rent_free = zeros(n),
        housing_association = zeros(n),
        las_etc_rented = zeros(n),
        m_0_4 = zeros(n),
        m_5_9 = zeros(n),
        m_10_14 = zeros(n),
        m_15_19 = zeros(n),
        m_20_24 = zeros(n),
        m_25_29 = zeros(n),
        m_30_34 = zeros(n),
        m_35_39 = zeros(n),
        m_40_44 = zeros(n),
        m_45_49 = zeros(n),
        m_50_54 = zeros(n),
        m_55_59 = zeros(n),
        m_60_64 = zeros(n),
        m_65_69 = zeros(n),
        m_70_74 = zeros(n),
        m_75_79 = zeros(n),
        m_80_plus = zeros(n),
        f_0_4 = zeros(n),
        f_5_9 = zeros(n),
        f_10_14 = zeros(n),
        f_15_19 = zeros(n),
        f_20_24 = zeros(n),
        f_25_29 = zeros(n),
        f_30_34 = zeros(n),
        f_35_39 = zeros(n),
        f_40_44 = zeros(n),
        f_45_49 = zeros(n),
        f_50_54 = zeros(n),
        f_55_59 = zeros(n),
        f_60_64 = zeros(n),
        f_65_69 = zeros(n),
        f_70_74 = zeros(n),
        f_75_79 = zeros(n),
        f_80_plus = zeros(n),
        v_1_adult_male = zeros(n),
        v_1_adult_female = zeros(n),
        v_2_adults = zeros(n),
        v_1_adult_1_child = zeros(n),
        v_1_adult_2_plus_children = zeros(n),
        v_2_plus_adults_1_plus_children = zeros(n),
        v_3_plus_adults = zeros(n),
        ca = zeros(n),
        aa = zeros(n),
        pip_or_dla = zeros(n),
        # councils
        # S12000033 = zeros(n), #  Aberdeen City  # removed for collinearity
        S12000034 = zeros(n), #  Aberdeenshire
        S12000041 = zeros(n), #  Angus
        S12000035 = zeros(n), #  Argyll and Bute
        S12000036 = zeros(n), #  City of Edinburgh
        S12000005 = zeros(n), #  Clackmannanshire
        S12000006 = zeros(n), #  Dumfries and Galloway
        S12000042 = zeros(n), #  Dundee City
        S12000008 = zeros(n), #  East Ayrshire
        S12000045 = zeros(n), #  East Dunbartonshire
        S12000010 = zeros(n), #  East Lothian
        S12000011 = zeros(n), #  East Renfrewshire
        S12000014 = zeros(n), #  Falkirk
        S12000047 = zeros(n), #  Fife
        S12000049 = zeros(n), #  Glasgow City
        S12000017 = zeros(n), #  Highland
        S12000018 = zeros(n), #  Inverclyde
        S12000019 = zeros(n), #  Midlothian
        S12000020 = zeros(n), #  Moray
        S12000013 = zeros(n), #  Na h-Eileanan Siar
        S12000021 = zeros(n), #  North Ayrshire
        S12000050 = zeros(n), #  North Lanarkshire
        S12000023 = zeros(n), #  Orkney Islands
        S12000048 = zeros(n), #  Perth and KinroS
        S12000038 = zeros(n), #  Renfrewshire
        S12000026 = zeros(n), #  Scottish Borders
        S12000027 = zeros(n), #  Shetland Islands
        S12000028 = zeros(n), #  South Ayrshire
        S12000029 = zeros(n), #  South Lanarkshire
        S12000030 = zeros(n), #  Stirling
        S12000039 = zeros(n), #  West Dunbartonshire
        S12000040 = zeros(n) #  West Lothian
    )
    return df
end

function make_target_row!( row :: DataFrameRow, hh :: Household )
    num_male_ads = 0
    num_female_ads = 0
    num_u_16s = 0
    for (pid,pers) in hh.people
        if( pers.age < 16 )
            num_u_16s += 1;
        end
        if pers.sex == Male

            if( pers.age >= 16 )
                num_male_ads += 1;
            end
            if pers.employment_status in [
                Full_time_Employee,
                Part_time_Employee,
                Full_time_Self_Employed,
                Part_time_Self_Employed
                ]
                row.m_total_in_employment += 1
            elseif pers.employment_status in [Unemployed]
                row.m_total_unemployed += 1
            else
                # row.m_total_economically_inactive += 1 # delete!
            end
            if pers.age <= 4
                row.m_0_4 += 1
            elseif pers.age <= 9
                row.m_5_9 += 1
            elseif pers.age <= 14
                row.m_10_14 += 1
            elseif pers.age <= 19
                row.m_15_19 += 1
            elseif pers.age <= 24
                row.m_20_24 += 1
            elseif pers.age <= 29
                row.m_25_29 += 1
            elseif pers.age <= 34
                row.m_30_34 += 1
            elseif pers.age <= 39
                row.m_35_39 += 1
            elseif pers.age <= 44
                row.m_40_44 += 1
            elseif pers.age <= 49
                row.m_45_49 += 1
            elseif pers.age <= 54
                row.m_50_54 += 1
            elseif pers.age <= 59
                row.m_55_59 += 1
            elseif pers.age <= 64
                row.m_60_64 += 1
            elseif pers.age <= 69
                row.m_65_69 += 1
            elseif pers.age <= 74
                row.m_70_74 += 1
            elseif pers.age <= 79
                row.m_75_79 += 1
            else
                row.m_80_plus += 1
            end
        else  # female
            if( pers.age >= 16 )
                num_female_ads += 1;
            end
            if pers.employment_status in [
                Full_time_Employee,
                Part_time_Employee,
                Full_time_Self_Employed,
                Part_time_Self_Employed
                ]
                row.f_total_in_employment += 1
            elseif pers.employment_status in [Unemployed]
                row.f_total_unemployed += 1
            else
                # row.f_total_economically_inactive += 1
            end
            if pers.age <= 4
                row.f_0_4 += 1
            elseif pers.age <= 9
                row.f_5_9 += 1
            elseif pers.age <= 14
                row.f_10_14 += 1
            elseif pers.age <= 19
                row.f_15_19 += 1
            elseif pers.age <= 24
                row.f_20_24 += 1
            elseif pers.age <= 29
                row.f_25_29 += 1
            elseif pers.age <= 34
                row.f_30_34 += 1
            elseif pers.age <= 39
                row.f_35_39 += 1
            elseif pers.age <= 44
                row.f_40_44 += 1
            elseif pers.age <= 49
                row.f_45_49 += 1
            elseif pers.age <= 54
                row.f_50_54 += 1
            elseif pers.age <= 59
                row.f_55_59 += 1
            elseif pers.age <= 64
                row.f_60_64 += 1
            elseif pers.age <= 69
                row.f_65_69 += 1
            elseif pers.age <= 74
                row.f_70_74 += 1
            elseif pers.age <= 79
                row.f_75_79 += 1
            else
                row.f_80_plus += 1
            end

        end # female
        if get(pers.income,attendance_allowance,0.0) > 0 ### sp!!!!!
            row.aa += 1
        end
        if get(pers.income,carers_allowance,0.0) > 0
            row.ca += 1
        end
        if get( pers.income, personal_independence_payment_daily_living, 0.0 ) > 0 ||
           get( pers.income, personal_independence_payment_mobility, 0.0 ) > 0 ||
           get( pers.income, dlaself_care, 0.0 ) > 0 ||
           get( pers.income, dlamobility, 0.0 ) > 0
           row.pip_or_dla += 1
       end
    end # people
    if is_owner_occupier(hh.tenure)
        # row.owner_occupied = 1
    elseif hh.tenure == Council_Rented
        row.las_etc_rented = 1
    elseif hh.tenure == Housing_Association
        row.housing_association = 1
    else
        row.private_rented_plus_rent_free = 1
    end
    num_people = num_u_16s+num_male_ads+num_female_ads
    num_adults = num_male_ads+num_female_ads
    if num_people == 1
        if num_male_ads == 1
            row.v_1_adult_male = 1
        else
            row.v_1_adult_female = 1
        end
    elseif num_adults == 1
        if num_u_16s == 1
            row.v_1_adult_1_child = 1
        else
            row.v_1_adult_2_plus_children = 1
        end
    elseif num_adults == 2 && num_u_16s == 0
        row.v_2_adults = 1
    elseif num_adults > 2 && num_u_16s == 0
        row.v_3_plus_adults = 1
    elseif num_adults >= 2 && num_u_16s > 0
        row.v_2_plus_adults_1_plus_children = 1
    else
        @assert false "should never get here num_male_ads=$num_male_ads num_female_ads=$num_female_ads num_u_16s=$num_u_16s"
    end
    ## las
    if hh.council !== :S12000033
        # aberdeen city, dropped to avoid collinearity, otherwise ...
        row[hh.council] = 1
    end
end

function make_target_dataset( nhhlds :: Integer ) :: Matrix
    df :: DataFrame = initialise_target_dataframe( nhhlds )
    for hno in 1:nhhlds
        hh = FRSHouseholdGetter.get_household( hno )
        make_target_row!( df[hno,:], hh )
    end
    return Matrix{Float64}(df) # convert( Matrix, df )
end

#
# generate weights for the dataset and
#
#
function generate_weights(
    nhhlds :: Integer;
    weight_type :: DistanceFunctionType = constrained_chi_square,
    lower_multiple :: Real = 0.20, # these values can be narrowed somewhat, to around 0.25-4.7
    upper_multiple :: Real = 5,
    targets :: Vector = DEFAULT_TARGETS ) :: Vector

    data :: Matrix = make_target_dataset( nhhlds )
    nrows = size( data )[1]
    ncols = size( data )[2]
    ## FIXME parameterise this
    initial_weights = ones(nhhlds)*NUM_HOUSEHOLDS/nhhlds
    println( "initial_weights $(initial_weights[1])")

     # any smaller min and d_and_s_constrained fails on this dataset
    rw = do_reweighting(
         data               = data,
         initial_weights    = initial_weights,
         target_populations = targets,
         functiontype       = weight_type,
         lower_multiple     = lower_multiple,
         upper_multiple     = upper_multiple,
         tolx               = 0.000001,
         tolf               = 0.000001 )
    # println( "results for method $weight_type = $(rw.rc)" )
    @assert rw.rc[:error] == 0 "non zero return code from weights gen $(rw.rc)"
    weights = rw.weights
    weighted_popn = (weights' * data)'
    # println( "weighted_popn = $weighted_popn" )
    @assert weighted_popn ≈ targets

    if weight_type in [constrained_chi_square, d_and_s_constrained ]
      # check the constrainted methods keep things inside ll and ul
        for r in 1:nrows
            @assert weights[r] <= initial_weights[r]*upper_multiple
            @assert weights[r] >= initial_weights[r]*lower_multiple
        end
    end
    for hno in 1:nhhlds
        hh = FRSHouseholdGetter.get_household( hno )
        hh.weight = weights[hno]
    end
    return weights
end

end # package

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# Ordinary metrics

###########################################################
#
#   Metric types
#
###########################################################

type Euclidean <: Metric end
type SqEuclidean <: SemiMetric end
type Chebyshev <: Metric end
type Cityblock <: Metric end
type Jaccard <: Metric end
type RogersTanimoto <: Metric end

immutable Minkowski{T <: Real} <: Metric
    p::T
end

type Hamming <: Metric end

type CosineDist <: SemiMetric end
type CorrDist <: SemiMetric end

type ChiSqDist <: SemiMetric end
type KLDivergence <: PreMetric end
type JSDivergence <: SemiMetric end

type SpanNormDist <: SemiMetric end

typealias UnionMetrics Union{Euclidean, SqEuclidean, Chebyshev, Cityblock, Minkowski, Hamming, Jaccard, RogersTanimoto, CosineDist, CorrDist, ChiSqDist, KLDivergence, JSDivergence, SpanNormDist}

###########################################################
#
#  Define Evaluate
#
###########################################################

function evaluate(d::UnionMetrics, a::AbstractArray, b::AbstractArray)
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    if length(a) == 0
        return zero(result_type(d, a, b))
    end
    s = eval_start(d, a, b)
    if size(a) == size(b)
        @simd for I in eachindex(a, b)
            @inbounds ai = a[I]
            @inbounds bi = b[I]
            s = eval_reduce(d, s, eval_op(d, ai, bi))
        end
    else
        for (Ia, Ib) in zip(eachindex(a), eachindex(b))
            @inbounds ai = a[Ia]
            @inbounds bi = b[Ib]
            s = eval_reduce(d, s, eval_op(d, ai, bi))
        end
    end
    return eval_end(d, s)
end
result_type{T1, T2}(dist::UnionMetrics, ::AbstractArray{T1}, ::AbstractArray{T2}) =
    typeof(eval_end(dist, eval_op(dist, one(T1), one(T2))))
eval_start(d::UnionMetrics, a::AbstractArray, b::AbstractArray) =
    zero(result_type(d, a, b))
eval_end(d::UnionMetrics, s) = s

evaluate{T <: Number}(dist::UnionMetrics, a::T, b::T) = eval_end(dist, eval_op(dist, a, b))

# SqEuclidean
@inline eval_op(::SqEuclidean, ai, bi) = abs2(ai - bi)
@inline eval_reduce(::SqEuclidean, s1, s2) = s1 + s2

sqeuclidean(a::AbstractArray, b::AbstractArray) = evaluate(SqEuclidean(), a, b)
sqeuclidean{T <: Number}(a::T, b::T) = evaluate(SqEuclidean(), a, b)

# Euclidean
@inline eval_op(::Euclidean, ai, bi) = abs2(ai - bi)
@inline eval_reduce(::Euclidean, s1, s2) = s1 + s2
eval_end(::Euclidean, s) = sqrt(s)
euclidean(a::AbstractArray, b::AbstractArray) = evaluate(Euclidean(), a, b)
euclidean(a::Number, b::Number) = evaluate(Euclidean(), a, b)

# Cityblock
@inline eval_op(::Cityblock, ai, bi) = abs(ai - bi)
@inline eval_reduce(::Cityblock, s1, s2) = s1 + s2
cityblock(a::AbstractArray, b::AbstractArray) = evaluate(Cityblock(), a, b)
cityblock{T <: Number}(a::T, b::T) = evaluate(Cityblock(), a, b)

# Chebyshev
@inline eval_op(::Chebyshev, ai, bi) = abs(ai - bi)
@inline eval_reduce(::Chebyshev, s1, s2) = max(s1, s2)
# if only NaN, will output NaN
@inline eval_start(::Chebyshev, a::AbstractArray, b::AbstractArray) = abs(a[1] - b[1])
chebyshev(a::AbstractArray, b::AbstractArray) = evaluate(Chebyshev(), a, b)
chebyshev{T <: Number}(a::T, b::T) = evaluate(Chebyshev(), a, b)

# Minkowski
@inline eval_op(dist::Minkowski, ai, bi) = abs(ai - bi) .^ dist.p
@inline eval_reduce(::Minkowski, s1, s2) = s1 + s2
eval_end(dist::Minkowski, s) = s .^ (1/dist.p)
minkowski(a::AbstractArray, b::AbstractArray, p::Real) = evaluate(Minkowski(p), a, b)
minkowski{T <: Number}(a::T, b::T, p::Real) = evaluate(Minkowski(p), a, b)

# Hamming
@inline eval_op(::Hamming, ai, bi) = ai != bi ? 1 : 0
@inline eval_reduce(::Hamming, s1, s2) = s1 + s2
hamming(a::AbstractArray, b::AbstractArray) = evaluate(Hamming(), a, b)
hamming{T <: Number}(a::T, b::T) = evaluate(Hamming(), a, b)

# Cosine dist
function eval_start{T<:AbstractFloat}(::CosineDist, a::AbstractArray{T}, b::AbstractArray{T})
    zero(T), zero(T), zero(T)
end
@inline eval_op(::CosineDist, ai, bi) = ai * bi, ai * ai, bi * bi
@inline function eval_reduce(::CosineDist, s1, s2)
    a1, b1, c1 = s1
    a2, b2, c2 = s2
    return a1 + a2, b1 + b2, c1 + c2
end
function eval_end(::CosineDist, s)
    ab, a2, b2 = s
    max(1 - ab / (sqrt(a2) * sqrt(b2)), zero(eltype(ab)))
end
cosine_dist(a::AbstractArray, b::AbstractArray) = evaluate(CosineDist(), a, b)

# Correlation Dist
_centralize(x::AbstractArray) = x .- mean(x)
evaluate(::CorrDist, a::AbstractArray, b::AbstractArray) = cosine_dist(_centralize(a), _centralize(b))
corr_dist(a::AbstractArray, b::AbstractArray) = evaluate(CorrDist(), a, b)
result_type(::CorrDist, a::AbstractArray, b::AbstractArray) = result_type(CosineDist(), a, b)

# ChiSqDist
@inline eval_op(::ChiSqDist, ai, bi) = abs2(ai - bi) / (ai + bi)
@inline eval_reduce(::ChiSqDist, s1, s2) = s1 + s2
chisq_dist(a::AbstractArray, b::AbstractArray) = evaluate(ChiSqDist(), a, b)

# KLDivergence
@inline eval_op(::KLDivergence, ai, bi) = ai > 0 ? ai * log(ai / bi) : zero(ai)
@inline eval_reduce(::KLDivergence, s1, s2) = s1 + s2
kl_divergence(a::AbstractArray, b::AbstractArray) = evaluate(KLDivergence(), a, b)

# JSDivergence
@inline function eval_op{T}(::JSDivergence, ai::T, bi::T)
    u = (ai + bi) / 2
    ta = ai > 0 ? ai * log(ai) / 2 : zero(log(one(T)))
    tb = bi > 0 ? bi * log(bi) / 2 : zero(log(one(T)))
    tu = u > 0 ? u * log(u) : zero(log(one(T)))
    ta + tb - tu
end
@inline eval_reduce(::JSDivergence, s1, s2) = s1 + s2
js_divergence(a::AbstractArray, b::AbstractArray) = evaluate(JSDivergence(), a, b)

# SpanNormDist
function eval_start(::SpanNormDist, a::AbstractArray, b::AbstractArray)
    a[1] - b[1], a[1] - b[1]
end
@inline eval_op(::SpanNormDist, ai, bi)  = ai - bi
@inline function eval_reduce(::SpanNormDist, s1, s2)
    min_d, max_d = s1
    if s2 > max_d
        max_d = s2
    elseif s2 < min_d
        min_d = s2
    end
    return min_d, max_d
end

eval_end(::SpanNormDist, s) = s[2] - s[1]
spannorm_dist(a::AbstractArray, b::AbstractArray) = evaluate(SpanNormDist(), a, b)
function result_type{T1, T2}(dist::SpanNormDist, ::AbstractArray{T1}, ::AbstractArray{T2})
    typeof(eval_op(dist, one(T1), one(T2)))
end

# Jaccard

@inline eval_start(::Jaccard, a::AbstractArray, b::AbstractArray) = 0, 0
@inline function eval_op(::Jaccard, s1, s2)
    denominator = max(s1, s2)
    numerator = min(s1, s2)
    numerator, denominator
end
@inline function eval_reduce(::Jaccard, s1, s2)
    a = s1[1] + s2[1]
    b = s1[2] + s2[2]
    a, b
end
@inline eval_end(::Jaccard, a) = 1 - (a[1]/a[2])
jaccard(a::AbstractArray, b::AbstractArray) = evaluate(Jaccard(), a, b)

# Tanimoto

@inline eval_start(::RogersTanimoto, a::AbstractArray, b::AbstractArray) = 0, 0, 0, 0
@inline function eval_op(::RogersTanimoto, s1, s2)
  tt = s1 && s2
  tf = s1 && !s2
  ft = !s1 && s2
  ff = !s1 && !s2
  tt, tf, ft, ff
end
@inline function eval_reduce(::RogersTanimoto, s1, s2)
    a = s1[1] + s2[1]
    b = s1[2] + s2[2]
    c = s1[3] + s2[3]
    d = s1[4] + s1[4]
    a, b, c, d
end
@inline function eval_end(::RogersTanimoto, a)
    numerator = 2(a[2] + a[3])
    denominator = a[1] + a[4] + 2(a[2] + a[3])
    numerator / denominator
end
rogerstanimoto{T <: Bool}(a::AbstractArray{T}, b::AbstractArray{T}) = evaluate(RogersTanimoto(), a, b)

###########################################################
#
#  Special method
#
###########################################################

# SqEuclidean
function pairwise!(r::AbstractMatrix, dist::SqEuclidean, a::AbstractMatrix, b::AbstractMatrix)
    At_mul_B!(r, a, b)
    sa2 = sumabs2(a, 1)
    sb2 = sumabs2(b, 1)
    pdist!(r, sa2, sb2)
end
function pdist!(r, sa2, sb2)
    for j = 1 : size(r,2)
        sb = sb2[j]
        @simd for i = 1 : size(r,1)
            @inbounds r[i,j] = sa2[i] + sb - 2 * r[i,j]
        end
    end
    r
end
function pairwise!(r::AbstractMatrix, dist::SqEuclidean, a::AbstractMatrix)
    m, n = get_pairwise_dims(r, a)
    At_mul_B!(r, a, a)
    sa2 = sumsq_percol(a)
    for j = 1 : n
        for i = 1 : j-1
            @inbounds r[i,j] = r[j,i]
        end
        @inbounds r[j,j] = 0
        for i = j+1 : n
            @inbounds r[i,j] = sa2[i] + sa2[j] - 2 * r[i,j]
        end
    end
    r
end

# Euclidean
function pairwise!(r::AbstractMatrix, dist::Euclidean, a::AbstractMatrix, b::AbstractMatrix)
    m, na, nb = get_pairwise_dims(r, a, b)
    At_mul_B!(r, a, b)
    sa2 = sumsq_percol(a)
    sb2 = sumsq_percol(b)
    for j = 1 : nb
        for i = 1 : na
            @inbounds v = sa2[i] + sb2[j] - 2 * r[i,j]
            @inbounds r[i,j] = isnan(v) ? NaN : sqrt(max(v, 0.))
        end
    end
    r
end

function pairwise!(r::AbstractMatrix, dist::Euclidean, a::AbstractMatrix)
    m, n = get_pairwise_dims(r, a)
    At_mul_B!(r, a, a)
    sa2 = sumsq_percol(a)
    for j = 1 : n
        for i = 1 : j-1
            @inbounds r[i,j] = r[j,i]
        end
        @inbounds r[j,j] = 0
        for i = j+1 : n
            @inbounds v = sa2[i] + sa2[j] - 2 * r[i,j]
            @inbounds r[i,j] = isnan(v) ? NaN : sqrt(max(v, 0.))
        end
    end
    r
end

# CosineDist

function pairwise!(r::AbstractMatrix, dist::CosineDist, a::AbstractMatrix, b::AbstractMatrix)
    m, na, nb = get_pairwise_dims(r, a, b)
    At_mul_B!(r, a, b)
    ra = sqrt!(sumsq_percol(a))
    rb = sqrt!(sumsq_percol(b))
    for j = 1 : nb
        @simd for i = 1 : na
            @inbounds r[i,j] = max(1 - r[i,j] / (ra[i] * rb[j]), 0)
        end
    end
    r
end
function pairwise!(r::AbstractMatrix, dist::CosineDist, a::AbstractMatrix)
    m, n = get_pairwise_dims(r, a)
    At_mul_B!(r, a, a)
    ra = sqrt!(sumsq_percol(a))
    for j = 1 : n
        @simd for i = j+1 : n
            @inbounds r[i,j] = max(1 - r[i,j] / (ra[i] * ra[j]), 0)
        end
        @inbounds r[j,j] = 0
        for i = 1 : j-1
            @inbounds r[i,j] = r[j,i]
        end
    end
    r
end

# CorrDist
_centralize_colwise(x::AbstractVector) = x .- mean(x)
_centralize_colwise(x::AbstractMatrix) = x .- mean(x, 1)
function colwise!(r::AbstractVector, dist::CorrDist, a::AbstractMatrix, b::AbstractMatrix)
    colwise!(r, CosineDist(), _centralize_colwise(a), _centralize_colwise(b))
end
function colwise!(r::AbstractVector, dist::CorrDist, a::AbstractVector, b::AbstractMatrix)
    colwise!(r, CosineDist(), _centralize_colwise(a), _centralize_colwise(b))
end
function pairwise!(r::AbstractMatrix, dist::CorrDist, a::AbstractMatrix, b::AbstractMatrix)
    pairwise!(r, CosineDist(), _centralize_colwise(a), _centralize_colwise(b))
end
function pairwise!(r::AbstractMatrix, dist::CorrDist, a::AbstractMatrix)
    pairwise!(r, CosineDist(), _centralize_colwise(a))
end

ext: jl


text:
"""Compute the entropy of the specified dataset."""
function entropy(dataset::AbstractArray{T, 2}, target::Integer,
                c::Integer=2) where {T}
    class = dataset[:, target]
    n = length(class)
    uniques = unique(class)

    return -mapreduce(u -> begin
        subset = dataset[class .== u, :]
        p = size(subset, 1) / n
        return p * log2(p)
    end, +, uniques) / log2(c)
end

entropy(dataset::AbstractArray{T, 2}, target::String,
        header::AbstractArray{String}, c::Integer=2) where {T} =
    entropy(dataset, something(findfirst(isequal(target), header), 0), c)

"""Compute the information gain of the specified attribute."""
function information_gain(dataset::AbstractArray{T, 2}, attribute::Integer,
                        target::Integer, c::Integer=2) where {T}
    n = size(dataset, 1)
    attr = dataset[:, attribute]
    uniques = unique(attr)

    attr_sub_entropy = mapreduce(u -> begin
        subset = dataset[attr .== u, :]
        return size(subset, 1) / n * entropy(subset, target, c)
    end, +, uniques)

    return entropy(dataset, target, c) - attr_sub_entropy
end

information_gain(dataset::AbstractArray{T, 2}, attribute::String,
                target::String, header::AbstractArray{String},
                c::Integer=2) where {T} =
    information_gain(dataset,
            something(findfirst(isequal(attribute), header), 0),
            something(findfirst(isequal(target), header), 0), c)

information_gain(dataset::AbstractArray{T, 2}, attribute::String,
                target::Integer, header::AbstractArray{String},
                c::Integer=2) where {T} =
    information_gain(dataset,
            something(findfirst(isequal(attribute), header), 0), target, c)

information_gain(dataset::AbstractArray{T, 2}, attribute::Integer,
                target::String, header::AbstractArray{String},
                c::Integer=2) where {T} =
    information_gain(dataset, attribute,
            something(findfirst(isequal(target), header), 0), c)


text:
@testset "Update Type Tree" begin
    @test MonteCarlo.AbstractLocalUpdate <: MonteCarlo.AbstractUpdate
    @test MonteCarlo.AbstractGlobalUpdate <: MonteCarlo.AbstractUpdate
    @test MonteCarlo.Adaptive <: MonteCarlo.AbstractUpdate

    @test LocalSweep <: MonteCarlo.AbstractLocalUpdate

    @test MonteCarlo.AbstractParallelUpdate <: MonteCarlo.AbstractGlobalUpdate
    @test GlobalShuffle <: MonteCarlo.AbstractGlobalUpdate
    @test GlobalFlip <: MonteCarlo.AbstractGlobalUpdate
    @test SpatialShuffle <: MonteCarlo.AbstractGlobalUpdate
    @test TemporalShuffle <: MonteCarlo.AbstractGlobalUpdate
    @test Denoise <: MonteCarlo.AbstractGlobalUpdate
    @test DenoiseFlip <: MonteCarlo.AbstractGlobalUpdate
    @test StaggeredDenoise <: MonteCarlo.AbstractGlobalUpdate

    @test ReplicaExchange <: MonteCarlo.AbstractParallelUpdate
    @test ReplicaPull <: MonteCarlo.AbstractParallelUpdate
end

struct BadUpdate <: MonteCarlo.AbstractLocalUpdate end
MonteCarlo.name(::BadUpdate) = "BadUpdate"
MonteCarlo.update(::BadUpdate, args...) = 0.0

struct GoodUpdate <: MonteCarlo.AbstractLocalUpdate end
MonteCarlo.name(::GoodUpdate) = "GoodUpdate"
MonteCarlo.update(::GoodUpdate, args...) = 1.0

@testset "Scheduler" begin
    model = HubbardModelAttractive(2,2)

    # These do not advance sweeps and should therefore error
    @test_throws ErrorException SimpleScheduler()
    @test_throws ErrorException SimpleScheduler(Adaptive())
    @test_throws ErrorException SimpleScheduler(GlobalFlip())

    @test_throws MethodError AdaptiveScheduler()
    @test_throws ErrorException AdaptiveScheduler(tuple(), tuple())
    @test_throws ErrorException AdaptiveScheduler((Adaptive(),), tuple())
    @test_throws ErrorException AdaptiveScheduler((GlobalFlip(),), tuple())

    # Check constructor
    schedulers = (
        SimpleScheduler(LocalSweep(), GlobalFlip(), LocalSweep(2)),
        AdaptiveScheduler((LocalSweep(), GlobalFlip(), LocalSweep(2)), tuple()),
    )
    for scheduler in schedulers
        Random.seed!(123)
        @test scheduler.sequence == MonteCarlo.AcceptanceStatistics.(
            (LocalSweep(), GlobalFlip(), LocalSweep(), LocalSweep())
        )
        @test scheduler.idx == 0

        mc = DQMC(model, beta=1.0, scheduler = scheduler)
        MonteCarlo.reverse_build_stack(mc, mc.stack)
        MonteCarlo.propagate(mc)

        @test mc.scheduler.idx == 0
        @test mc.last_sweep == 0

        MonteCarlo.update(mc.scheduler, mc, mc.model)
        @test mc.scheduler.sequence[1].total == 1
        @test mc.scheduler.idx == 1
        @test mc.last_sweep == 1

        MonteCarlo.update(mc.scheduler, mc, mc.model)
        @test mc.scheduler.sequence[2].total == 1
        @test mc.scheduler.idx == 2
        @test mc.last_sweep == 2

        MonteCarlo.update(mc.scheduler, mc, mc.model)
        @test mc.scheduler.sequence[3].total == 1
        @test mc.scheduler.idx == 3
        @test mc.last_sweep == 3
    end

    schedulers = (
        SimpleScheduler(GoodUpdate(), GoodUpdate(), BadUpdate()),
        AdaptiveScheduler((GoodUpdate(), GoodUpdate(), BadUpdate()), tuple())
    )
    for scheduler in schedulers
        mc = DQMC(model, beta=1.0, scheduler = scheduler)
        MonteCarlo.reverse_build_stack(mc, mc.stack)
        MonteCarlo.propagate(mc)

        for _ in 1:300
            MonteCarlo.update(mc.scheduler, mc, mc.model)
        end

        io = IOBuffer()
        MonteCarlo.show_statistics(io, mc.scheduler)
        @test String(take!(io)) == "Update statistics (since start):\n\tBadUpdate              0.0% accepted   (  0 / 100)\n\tGoodUpdate           100.0% accepted   (200 / 200)\n\t--------------------------------------------------\n\tTotal                 66.7% accepted   (200 / 300)\n"
    end
end

@testset "AdaptiveScheduler" begin
    scheduler = AdaptiveScheduler(
        (BadUpdate(), Adaptive()),
        (GoodUpdate(), BadUpdate())
    )

    # check defaults
    @test scheduler.grace_period == 99
    @test scheduler.minimum_sampling_rate == 0.01
    @test scheduler.adaptive_rate == 9.0

    model = HubbardModelAttractive(2,2)
    mc = DQMC(model, beta=1.0, scheduler = scheduler)

    # Checks without adaptive corrections
    # should be save to run, 1 / 2^50 chance of already adapting
    for _ in 1:200
        MonteCarlo.update(scheduler, mc, model)
    end

    @test scheduler.sampling_rates == [0.5, 0.5, 1e-10]
    @test scheduler.sequence[1].accepted == 0.0
    @test scheduler.sequence[1].total == 100
    # no adaptive stuff yet, both updates should be picked at 50%
    @test scheduler.adaptive_pool[1].accepted == scheduler.adaptive_pool[1].total
    @test 30 <= scheduler.adaptive_pool[1].total <= 70
    @test scheduler.adaptive_pool[2].accepted == 0
    @test 30 <= scheduler.adaptive_pool[2].total <= 70

    # Checks with adaptive corrections
    for _ in 1:800
        MonteCarlo.update(scheduler, mc, model)
    end
    @test scheduler.sampling_rates[1] > 0.925
    @test scheduler.sampling_rates[2] < 0.075
    @test scheduler.sampling_rates[3] == 1e-10
    @test scheduler.sequence[1].accepted == 0.0
    @test scheduler.sequence[1].total == 500

    # With 0 acceptance rate we should have the hard limit for total:
    i_min = ceil(Int, log(
        scheduler.adaptive_rate / (scheduler.adaptive_rate+1),
        scheduler.minimum_sampling_rate / 0.5
    )) + 99 + 1 # +1 to be save
    @test scheduler.adaptive_pool[1].accepted == scheduler.adaptive_pool[1].total
    @test scheduler.adaptive_pool[1].total > 500 - i_min
    @test scheduler.adaptive_pool[2].accepted == 0
    @test scheduler.adaptive_pool[2].total < i_min
end

using MonteCarlo: conf, current_slice, nslices

@testset "Global update" begin
    # Verify that split greens calculation matches the normal one
    pivot = Vector{Int}(undef, 16)
    tempv = Vector{Float64}(undef, 16)

    Ul = Matrix{Float64}(undef, 16, 16)
    Dl = Vector{Float64}(undef, 16)
    Tl = Matrix{Float64}(undef, 16, 16)

    Ur = Matrix{Float64}(undef, 16, 16)
    Dr = Vector{Float64}(undef, 16)
    Tr = Matrix{Float64}(undef, 16, 16)

    for i in 1:10
        Bl = rand(16, 16)
        copyto!(Tl, Bl)
        MonteCarlo.udt_AVX_pivot!(Ul, Dl, Tl, pivot, tempv)
        Br = rand(16, 16)
        copyto!(Tr, Br)
        MonteCarlo.udt_AVX_pivot!(Ur, Dr, Tr, pivot, tempv)
        G1 = Matrix{Float64}(undef, 16, 16)
        MonteCarlo.calculate_greens_AVX!(Ul, Dl, Tl, Ur, Dr, Tr, G1, pivot, tempv)

        copyto!(Tl, Bl)
        MonteCarlo.udt_AVX_pivot!(Ul, Dl, Tl)
        copyto!(Tr, Br)
        MonteCarlo.udt_AVX_pivot!(Ur, Dr, Tr)
        G2 = Matrix{Float64}(undef, 16, 16)
        MonteCarlo.calculate_inv_greens_udt(Ul, Dl, Tl, Ur, Dr, Tr, G2, pivot, tempv)
        MonteCarlo.finish_calculate_greens(Ul, Dl, Tl, Ur, Dr, Tr, G2, pivot, tempv)

        @test G1 ≈ G2
    end

    models = (HubbardModelAttractive(2,2,mu=0.5), HubbardModelRepulsive(2,2))
    for model in models
        @testset "$(typeof(model))" begin
            mc1 = DQMC(model, beta=2.0)
            mc2 = DQMC(model, beta=2.0)

            # Verify probabilities and greens in global updates
            # This is the backbone check for all global and parallel updates
            for _ in 1:10
                # re-initialize everything with a random conf
                mc1.conf .= rand(DQMC, model, nslices(mc1))
                mc2.conf .= mc1.conf
                for mc in (mc1, mc2)
                    MonteCarlo.init!(mc)
                    MonteCarlo.reverse_build_stack(mc, mc.stack)
                    MonteCarlo.propagate(mc)
                end

                # global update
                temp_conf = shuffle(deepcopy(conf(mc1)))
                detratio, ΔE_boson, passthrough = MonteCarlo.propose_global_from_conf(mc1, model, temp_conf)
                global_p = exp(- ΔE_boson) * detratio
                MonteCarlo.accept_global!(mc1, model, temp_conf, passthrough)

                # global update through successive local updates
                local_p = 1.0
                for t in 1:nslices(mc2)
                    for i in 1:length(lattice(mc2))
                        if mc2.conf[i, current_slice(mc2)] != temp_conf[i, current_slice(mc2)]
                            detratio, ΔE_boson, passthrough = MonteCarlo.propose_local(
                                mc2, model, i, current_slice(mc2), conf(mc2)
                            )
                            local_p *= real(exp(- ΔE_boson) * detratio)
                            MonteCarlo.accept_local!(
                                mc2, model, i, current_slice(mc2), conf(mc2), detratio, ΔE_boson, passthrough
                            )
                        end
                    end
                    MonteCarlo.propagate(mc2)
                end

                # Move to correct time slice
                for t in 1:nslices(mc2)
                    MonteCarlo.propagate(mc2)
                end

                # Verify
                @test local_p ≈ global_p
                @test mc1.conf == mc2.conf
                @test current_slice(mc1) == current_slice(mc2)
                @test mc1.stack.greens ≈ mc2.stack.greens
            end
        end
    end

    function setup()
        model = HubbardModelAttractive(8,2)
        mc = DQMC(model, beta=1.0)
        MonteCarlo.init!(mc)
        MonteCarlo.reverse_build_stack(mc, mc.stack)
        MonteCarlo.propagate(mc)
        c = deepcopy(conf(mc))
        return mc, model, c
    end

    # Check config adjustments for global updates
    mc, model, c = setup()
    MonteCarlo.update(GlobalFlip(), mc, model)
    @test mc.temp_conf == -c

    # make this unlikely to fail randomly via big conf size
    mc, model, c = setup()
    MonteCarlo.update(GlobalShuffle(), mc, model)
    @test sum(mc.temp_conf) == sum(c)
    @test mc.temp_conf != c

    mc, model, c = setup()
    u = SpatialShuffle()
    MonteCarlo.init!(mc, u)
    MonteCarlo.update(u, mc, model)
    @test mc.temp_conf == c[u.indices, :]

    mc, model, c = setup()
    u = TemporalShuffle()
    MonteCarlo.init!(mc, u)
    MonteCarlo.update(u, mc, model)
    @test mc.temp_conf == c[:, u.indices]

    #=
    Lattice with conf:

        7 8 9                    + - -
        -----                    -----
    3 | 1 2 3 | 1            + | - + + | -
    6 | 4 5 6 | 4            - | + + - | +
    9 | 7 8 9 | 7            - | + - - | +
        -----                    -----
        1 2 3                    - + +

    Denoise                   DenoiseFlip             StaggeredDenoise
    1 = - -> +                1 = - -> -              1 = - -> -
    2 = + -> +                2 = + -> -              2 = + -> +
    3 = + -> -                3 = + -> +              3 = + -> +
    4 = + -> +                4 = + -> -              4 = + -> +
    5 = + -> +                5 = + -> -              5 = + -> -
    6 = - -> +                6 = - -> -              6 = - -> +
    7 = + -> -                7 = + -> +              7 = + -> +
    8 = - -> +                8 = - -> -              8 = - -> +
    9 = - -> -                9 = - -> +              9 = - -> +
    to surrounding            to -surrounding         mult by +1 for even, -1 for odd sites
    =#
    spatial = [-1, +1, +1, +1, +1, -1, +1, -1, -1]
    _conf = [spatial[i] for i in 1:9, slice in 1:10]

    model = HubbardModelAttractive(3, 2)
    mc = DQMC(model, beta=1.0)
    MonteCarlo.init!(mc)
    mc.conf .= _conf
    MonteCarlo.reverse_build_stack(mc, mc.stack)
    MonteCarlo.propagate(mc)
    MonteCarlo.update(Denoise(), mc, model)
    @test mc.temp_conf == [[+1, +1, -1, +1, +1, +1, -1, +1, -1][i] for i in 1:9, slice in 1:10]

    model = HubbardModelAttractive(3, 2)
    mc = DQMC(model, beta=1.0)
    MonteCarlo.init!(mc)
    mc.conf .= _conf
    MonteCarlo.reverse_build_stack(mc, mc.stack)
    MonteCarlo.propagate(mc)
    MonteCarlo.update(DenoiseFlip(), mc, model)
    @test mc.temp_conf == [[-1, -1, +1, -1, -1, -1, +1, -1, +1][i] for i in 1:9, slice in 1:10]

    model = HubbardModelAttractive(3, 2)
    mc = DQMC(model, beta=1.0)
    MonteCarlo.init!(mc)
    mc.conf .= _conf
    MonteCarlo.reverse_build_stack(mc, mc.stack)
    MonteCarlo.propagate(mc)
    MonteCarlo.update(StaggeredDenoise(), mc, model)
    @test mc.temp_conf == [[-1, +1, +1, +1, -1, +1, +1, +1, +1][i] for i in 1:9, slice in 1:10]
end

using Distributed
@testset "Parallel" begin
    @test isempty(MonteCarlo.connected_ids)

    # backend
    MonteCarlo.add_worker!(myid())
    @test MonteCarlo.connected_ids == [myid()]
    MonteCarlo.add_worker!(myid())
    @test MonteCarlo.connected_ids == [myid()]
    MonteCarlo.add_worker!(-1)
    @test MonteCarlo.connected_ids == [myid(), -1]
    MonteCarlo.remove_worker!(-1)
    @test MonteCarlo.connected_ids == [myid()]
    MonteCarlo.remove_worker!(myid())
    @test MonteCarlo.connected_ids == Int[]

    connect(myid())
    yield()
    @test MonteCarlo.connected_ids == [myid()]
    disconnect(myid())
    yield()
    @test MonteCarlo.connected_ids == Int[]
    connect([myid()])
    yield()
    @test MonteCarlo.connected_ids == [myid()]
    disconnect([myid()])
    yield()
    @test MonteCarlo.connected_ids == Int[]

    @test !isready(MonteCarlo.weight_probability)
    MonteCarlo.put_weight_prob!(myid(), 0.9, 0.6)
    yield()
    @test isready(MonteCarlo.weight_probability)
    @test take!(MonteCarlo.weight_probability) == (myid(), 0.9, 0.6)
end

@testset "Utilities/Printing" begin
    @test MonteCarlo.name(NoUpdate()) == "NoUpdate"
    @test MonteCarlo.name(LocalSweep()) == "LocalSweep"
    @test MonteCarlo.name(GlobalFlip()) == "GlobalFlip"
    @test MonteCarlo.name(GlobalShuffle()) == "GlobalShuffle"
    @test MonteCarlo.name(ReplicaExchange(1)) == "ReplicaExchange"
    @test MonteCarlo.name(ReplicaPull()) == "ReplicaPull"
    @test MonteCarlo.name(Adaptive()) == "Adaptive"

    # Idk if I want to keep these methods...
    model = HubbardModelAttractive(2,2)
    mc = DQMC(model, beta=1.0)
    for T in (LocalSweep, GlobalFlip, GlobalShuffle, NoUpdate, ReplicaPull)
        @test T(mc, model) == T()
    end
    @test ReplicaExchange(mc, model, 1) == ReplicaExchange(1)
end


text:
## Unary operations
-(x::C) where C <: CasadiSymbolicObject = pycall(casadi.minus, C, 0, x)

## Binary operations
+(x::C, y::Real) where C <: CasadiSymbolicObject = pycall(casadi.plus, C, x, y)
-(x::C, y::Real) where C <: CasadiSymbolicObject = pycall(casadi.minus, C, x, y)
/(x::C, y::Real) where C <: CasadiSymbolicObject = pycall(casadi.mrdivide, C, x, y)
^(x::C, y::Real) where C <: CasadiSymbolicObject = pycall(casadi.power, C, x, y)
^(x::C, y::Integer) where C <: CasadiSymbolicObject = pycall(casadi.power, C, x, y)
\(x::C, y::Real) where C <: CasadiSymbolicObject = pycall(casadi.solve, C, x, y)

function *(x::C, y::Real) where C <: CasadiSymbolicObject
    if size(x,2) == size(y,1)
        pycall(casadi.mtimes, C, x, y)
    else
        pycall(casadi.times, C, x, y)
    end
end

## Comparisons
>=(x::C, y::Real) where C <: CasadiSymbolicObject = pycall(casadi.ge, C, x, y)
>(x::C, y::Real)  where C <: CasadiSymbolicObject = pycall(casadi.gt, C, x, y)
<=(x::C, y::Real) where C <: CasadiSymbolicObject = pycall(casadi.le, C, x, y)
<(x::C, y::Real)  where C <: CasadiSymbolicObject = pycall(casadi.lt, C, x, y)
==(x::C, y::Real) where C <: CasadiSymbolicObject = pycall(casadi.eq, C, x, y)

## Linear algebra
×(x::C, y::Real) where C <: CasadiSymbolicObject = pycall(casadi.cross, C, x, y)


text:
############################################################################################
# Include files
include("parameter.jl")

############################################################################################
# Export


text:
# operation on MATLAB engine sessions

###########################################################
#
#   Session open & close
#
###########################################################
const default_startcmd = matlab_startcmd() * " -nosplash"

# 64 K buffer should be sufficient to store the output text in most cases
const default_output_buffer_size = 64 * 1024

type MSession
    ptr::Ptr{Void}
    buffer::Vector{UInt8}
    bufptr::Ptr{UInt8}

    function MSession(bufsize::Integer = default_output_buffer_size)
        ep = ccall(eng_open[], Ptr{Void}, (Ptr{UInt8},), default_startcmd)
        ep == C_NULL && throw(MEngineError("failed to open a MATLAB engine session"))
        # hide the MATLAB command window on Windows
        is_windows() && ccall(eng_set_visible[], Cint, (Ptr{Void}, Cint), ep, 0)

        buf = Array(UInt8, bufsize)
        if bufsize > 0
            bufptr = pointer(buf)
            ccall(eng_output_buffer[], Cint, (Ptr{Void}, Ptr{UInt8}, Cint),
                ep, bufptr, bufsize)
        else
            bufptr = convert(Ptr{UInt8}, C_NULL)
        end

        self = new(ep, buf, bufptr)
        finalizer(self, release)
        return self
    end
end

function unsafe_convert(::Type{Ptr{Void}}, m::MSession)
    ptr = m.ptr
    ptr == C_NULL && throw(UndefRefError())
    return ptr
end

function release(session::MSession)
    ptr = session.ptr
    if ptr != C_NULL
        ccall(eng_close[], Cint, (Ptr{Void},), ptr)
    end
    session.ptr = C_NULL
    return nothing
end

function close(session::MSession)
    # close a MATLAB Engine session
    ret = ccall(eng_close[], Cint, (Ptr{Void},), session)
    ret != 0 && throw(MEngineError("failed to close a MATLAB engine session (err = $ret)"))
    session.ptr = C_NULL
    return nothing
end

# default session

default_msession = nothing

function restart_default_msession(bufsize::Integer = default_output_buffer_size)
    global default_msession
    if default_msession !== nothing && default_msession.ptr != C_NULL
        close(default_msession)
    end
    default_msession = MSession(bufsize)
    return nothing
end

function get_default_msession()
    global default_msession
    if default_msession === nothing
        default_msession = MSession()
    end
    return default_msession::MSession
end

function close_default_msession()
    global default_msession
    if default_msession !== nothing
        close(default_msession)
        default_msession = nothing
    end
    return nothing
end

function show_msession(m::MSession = get_default_msession())
    ret = ccall(eng_set_visible[], Cint, (Ptr{Void}, Cint), m, 1)
    ret != 0 && throw(MEngineError("failed to show MATLAB engine session (err = $ret)"))
    return nothing
end

function hide_msession(m::MSession = get_default_msession())
    ret = ccall(eng_set_visible[], Cint, (Ptr{Void}, Cint), m, 0)
    ret != 0 && throw(MEngineError("failed to hide MATLAB engine session (err = $ret)"))
    return nothing
end

###########################################################
#
#   communication with MATLAB session
#
###########################################################

function eval_string(session::MSession, stmt::String)
    # evaluate a MATLAB statement in a given MATLAB session
    ret = ccall(eng_eval_string[], Cint, (Ptr{Void}, Ptr{UInt8}), session, stmt)
    ret != 0 && throw(MEngineError("invalid engine session (err = $ret)"))

    bufptr = session.bufptr
    if bufptr != C_NULL
        bs = unsafe_string(bufptr)
        if ~isempty(bs)
            print(bs)
        end
    end
    return nothing
end

eval_string(stmt::String) = eval_string(get_default_msession(), stmt)

function put_variable(session::MSession, name::Symbol, v::MxArray)
    # put a variable into a MATLAB engine session
    ret = ccall(eng_put_variable[], Cint, (Ptr{Void}, Ptr{UInt8}, Ptr{Void}), session, string(name), v)
    ret != 0 && throw(MEngineError("failed to put the variable $(name) into a MATLAB session (err = $ret)"))
    return nothing
end

put_variable(session::MSession, name::Symbol, v) = put_variable(session, name, mxarray(v))

put_variable(name::Symbol, v) = put_variable(get_default_msession(), name, v)

function get_mvariable(session::MSession, name::Symbol)
    pv = ccall(eng_get_variable[], Ptr{Void}, (Ptr{Void}, Ptr{UInt8}), session, string(name))
    pv == C_NULL && throw(MEngineError("failed to get the variable $(name) from a MATLAB session"))
    return MxArray(pv)
end

get_mvariable(name::Symbol) = get_mvariable(get_default_msession(), name)

get_variable(name::Symbol) = jvalue(get_mvariable(name))
get_variable(name::Symbol, kind) = jvalue(get_mvariable(name), kind)

###########################################################
#
#   macro to simplify syntax
#
###########################################################

function _mput_multi(vs::Symbol...)
    nv = length(vs)
    if nv == 1
        v = vs[1]
        :( MATLAB.put_variable($(Meta.quot(v)), $(v)) )
    else
        stmts = Array(Expr, nv)
        for i = 1 : nv
            v = vs[i]
            stmts[i] = :( MATLAB.put_variable($(Meta.quot(v)), $(v)) )
        end
        Expr(:block, stmts...)
    end
end

macro mput(vs...)
    esc( _mput_multi(vs...) )
end

function make_getvar_statement(v::Symbol)
    :( $(v) = MATLAB.get_variable($(Meta.quot(v))) )
end

function make_getvar_statement(ex::Expr)
    if !(ex.head == :(::))
        error("Invalid expression for @mget.")
    end
    v::Symbol = ex.args[1]
    k::Symbol = ex.args[2]

    :( $(v) = MATLAB.get_variable($(Meta.quot(v)), $(k)) )
end

function _mget_multi(vs::Union{Symbol, Expr}...)
    nv = length(vs)
    if nv == 1
        make_getvar_statement(vs[1])
    else
        stmts = Array(Expr, nv)
        for i = 1:nv
            stmts[i] = make_getvar_statement(vs[i])
        end
        Expr(:block, stmts...)
    end
end

macro mget(vs...)
    esc( _mget_multi(vs...) )
end

macro matlab(ex)
    :( MATLAB.eval_string($(mstatement(ex))) )
end

###########################################################
#
#   mxcall
#
###########################################################

# MATLAB does not allow underscore as prefix of a variable name
_gen_marg_name(mfun::Symbol, prefix::String, i::Int) = "jx_$(mfun)_arg_$(prefix)_$(i)"

function mxcall(session::MSession, mfun::Symbol, nout::Integer, in_args...)
    nin = length(in_args)

    # generate temporary variable names

    in_arg_names = Array(String, nin)
    out_arg_names = Array(String, nout)

    for i = 1:nin
        in_arg_names[i] = _gen_marg_name(mfun, "in", i)
    end

    for i = 1:nout
        out_arg_names[i] = _gen_marg_name(mfun, "out", i)
    end

    # generate MATLAB statement

    buf = IOBuffer()
    if nout > 0
        if nout > 1
            print(buf, "[")
        end
        join(buf, out_arg_names, ", ")
        if nout > 1
            print(buf, "]")
        end
        print(buf, " = ")
    end

    print(buf, string(mfun))
    print(buf, "(")
    if nin > 0
        join(buf, in_arg_names, ", ")
    end
    print(buf, ");")

    stmt = takebuf_string(buf)

    # put variables to MATLAB

    for i = 1:nin
        put_variable(session, Symbol(in_arg_names[i]), in_args[i])
    end

    # execute MATLAB statement

    eval_string(session, stmt)

    # get results from MATLAB

    ret = if nout == 1
        jvalue(get_mvariable(session, Symbol(out_arg_names[1])))
    elseif nout >= 2
        results = Array(Any, nout)
        for i = 1 : nout
            results[i] = jvalue(get_mvariable(session, Symbol(out_arg_names[i])))
        end
        tuple(results...)
    else
        nothing
    end

    # clear temporaries from MATLAB workspace

    for i = 1:nin
        eval_string(session, string("clear ", in_arg_names[i], ";"))
    end

    for i = 1:nout
        eval_string(session, string("clear ", out_arg_names[i], ";"))
    end

    return ret
end

mxcall(mfun::Symbol, nout::Integer, in_args...) = mxcall(get_default_msession(), mfun, nout, in_args...)


text:
@sampler SimpleSampler = begin
    N::MvNormal

    SimpleSampler(μ, Σ) = new(MvNormal(μ, Σ))

    @sample Scenario begin
        x = rand(sampler.N)
        return Scenario(q₁ = x[1], q₂ = x[2], d₁ = x[3], d₂ = x[4], probability = pdf(sampler.N, x))
    end
end

μ = [28, 32, 300, 300]
Σ = [2 0.5 0 0
     0.5 1 0 0
     0 0 50 20
     0 0 20 30]

sampler = SimpleSampler(μ, Σ)


text:
mutable struct PerftTree
    nodes :: Array{Int64,1}
    divide :: Dict{String,Int64}
end

function perft(board::Bitboard, depth::Int64, color::String)

    pt = PerftTree(zeros(depth), Dict{String,Int64}())

    return explore(pt, board, depth, 1, color)
end

function print_perftree(pt::PerftTree)
    println("Nodes      ", pt.nodes)
    for x in sort(collect(pt.divide))
        println(x)
    end
end

function sumdivide(pt::PerftTree)
    c = 0
    for k in keys(pt.divide)
        c += pt.divide[k]
    end
    println(c)
end

function explore(pt::PerftTree, board::Bitboard,
    max_depth::Int64, depth::Int64, color::String="white",
    move_name::String="")

    moves = get_all_valid_moves(board, color)

    if length(moves) == 0
        return pt
    end

    if depth == 1
        for m in moves
            push!(pt.divide,
                m.piece_type*"-"*UINT2PGN[m.source]*UINT2PGN[m.target]=>0)
        end
    end

    if depth > max_depth
        return pt
    end

    pt.nodes[depth] += length(moves)

    if depth >= max_depth && depth > 1
        pt.divide[move_name] += length(moves)
        return pt
    end

    new_color = change_color(color)
    c = 1
    for m in moves
        # newb = deepcopy(board)
        board = move_piece(board, m, color)
        # board = update_attacked(board)
        board = update_castling_rights(board)

        if depth == 1
            move_name = m.piece_type*"-"*UINT2PGN[m.source]*UINT2PGN[m.target]
        else
            pt.divide[move_name] += 1
        end

        pt = explore(pt, board, max_depth, depth+1, new_color, move_name)
        board = unmove_piece(board, m, color)
        # board = update_attacked(board)
        # board = update_castling_rights(board)
    end

    return pt
end


text:
img2array(img) = reshape(
    convert(Array{Float32}, img), size(img) ..., 1, 1
)

array2img(x) = convert(
    Array{Gray{Normed{UInt8,8}},2},
    reshape(x.data, size(x)[1], size(x)[2])
)

function unet_dimensions(input)
    output = input - 4
    for depth in 1:4
        output = output ÷ 2 - 4
    end

    for expand in 1:4
        output = 2 * output -4
    end

    return output
end

function required_padding(img)
    # calculate how much padding is needed around img
    height, width = size(img)
    output_height = unet_dimensions(height)

    return (height - output_height) ÷ 2
end

function unet_tiling(img)
    # given an input image, tiles it for the unet
    padding = required_padding(img)
    return padarray(img, Pad(:symmetric, padding, padding))
end


text:
#White Noise kernel

@doc """
# Description
Constructor for the Noise kernel (covariance)

k(x,x') = σ²δ(x-x'), where δ is a Kronecker delta function and equals 1 iff x=x' and zero otherwise
# Arguments:
* `lσ::Float64`: Log of the signal standard deviation σ
""" ->
type Noise <: Kernel
    lσ::Float64      # Log of Signal std
    Noise(lσ::Float64) = new(lσ)
end

function cov(noise::Noise, x::Vector{Float64}, y::Vector{Float64})
    sigma2 = exp(2*noise.lσ)
    prec   = eps()            #machine precision

    K =  sigma2*(norm(x-y)<prec)
    return K
end

get_params(noise::Noise) = Float64[noise.lσ]
get_param_names(noise::Noise) = [:lσ]
num_params(noise::Noise) = 1

function set_params!(noise::Noise, hyp::Vector{Float64})
    length(hyp) == 1 || throw(ArgumentError("Noise kernel only has one parameter"))
    noise.lσ = hyp[1]
end

function grad_kern(noise::Noise, x::Vector{Float64}, y::Vector{Float64})
    sigma2 = exp(2*noise.lσ)
    prec   = eps()            #machine precision

    dK_sigma = 2.0*sigma2*(norm(x-y)<prec)

    dK_theta = [dK_sigma]
    return dK_theta
end


text:
# Fourier method
# First published by Milstein, 1994

struct Milstein <: AbstractIteratedIntegralAlgorithm end

convorder(::Milstein) = 1//2
errcoeff(m, h, ::Milstein, ::MaxL2) = h/(√2*π)
norv(m, n, ::Milstein) = 2*m*n+m

"""
    levyarea(W, n, alg::Milstein)

Simulates an approximation of the iterated Itô-integrals ``\\int_0^1W_s\\otimes dW_s``
of the given ``m``-dimensional increment of a Wiener process with step size 1.
The parameter ``n`` specifies the number of terms in the approximation and thus determines the accuracy.
This is an efficient implementation of the algorithm proposed in [Milstein, 1994](@ref milstein1994).
It is based on a Fourier expansion of the Wiener process.
The algorithm needs approximately ``m^2+2\\cdot m\\cdot n`` Float's.
The time complexity is ``\\mathcal{O}(m^2\\cdot n)``.
"""
function levyarea(W::AbstractVector{T}, n::Integer, alg::Milstein) where {T}
    rng = default_rng()
    m = length(W)
    X = randn(rng, T, n, m) # allocates m*n Floats
    Y = randn(rng, T, m, n) # allocates m*n Floats
    Y .= (Y .- √(T(2)).*W) ./ (1:n)'
    A = Y * X # allocates m*m Floats
    # Add simple rest approximation (a₀)
    # M = randn(rng, T, m) # allocates m Floats
    M = randn!(rng, view(Y, :, 1))
    a = T(√(2*trigamma(n+1)))
    A .+= a .* W .* M'
    # Antisymmetrize
    # G = inv(2*T(π)).*(A .- A') # allocates m*m Floats
    for i = 1:m
        @inbounds A[i,i] = zero(T)
        for j = (i+1):m
            @inbounds A[i,j] = (A[i,j] - A[j,i])/(2*T(π))
            @inbounds A[j,i] = -A[i,j]
        end
    end
    return A
end

max_stars_repo_name: mrbauff/MasterThesis

max_issues_repo_name: mrbauff/MasterThesis


text:
include("MPC.jl")
# using RigidBodyDynamics
# using StaticArrays
# using RobotDynamics
# using Altro
# using TrajectoryOptimization

# model = Arthur()
# n,m = size(model)

# tf = 4.1
# dt = 0.1
# N = Int(round(tf/dt) + 1)

qref = typeof(zeros(7))[]
push!(qref,  [4.9268856741946365e-05, 0.26016423071338934, 3.140069344584098, -2.27008100955155, -2.40769461674617e-05, 0.9598460188468598, 1.569995694393751])
push!(qref,  [-0.003768053224579125, 0.26079787208394173, 3.1395327288847428, -2.269096150780195, -0.004224076946167462, 0.9616197667116058, 1.5721591274301712])
push!(qref,  [-0.01522001946854234, 0.2626987961955988, 3.137922881786677, -2.2661415744661304, -0.016824076946167464, 0.9669410103058437, 1.578649426539432])
push!(qref,  [-0.03430662987514771, 0.26586700304836064, 3.135239803289901, -2.2612172806093556, -0.03782407694616747, 0.9758097496295736, 1.5894665917215332])
push!(qref,  [-0.0610278844443952, 0.2703024926422273, 3.131483493394414, -2.254323269209871, -0.06722407694616747, 0.9882259846827953, 1.6046106229764747])
push!(qref,  [-0.0953296828066516, 0.2759962847990445, 3.126661557197997, -2.2454734980167035, -0.10496455313664366, 1.0041645773154986, 1.6240508594052145])
push!(qref,  [-0.13259401741002397, 0.2821818315115795, 3.121423165847148, -2.23585940048681, -0.14596455313664364, 1.0214797350427807, 1.6451700866655072])
push!(qref,  [-0.1698583520133963, 0.2883673782241145, 3.116184774496299, -2.226245302956917, -0.18696455313664362, 1.0387948927700625, 1.6662893139258002])
push!(qref,  [-0.20712268661676866, 0.2945529249366496, 3.1109463831454507, -2.2166312054270234, -0.2279645531366436, 1.0561100504973449, 1.6874085411860928])
push!(qref,  [-0.24438702122014103, 0.3007384716491846, 3.1057079917946018, -2.20701710789713, -0.2689645531366436, 1.0734252082246267, 1.7085277684463855])
push!(qref,  [-0.2816513558235134, 0.3069240183617196, 3.100469600443753, -2.1974030103672364, -0.3099645531366436, 1.090740365951909, 1.7296469957066785])
push!(qref,  [-0.3189156904268857, 0.31310956507425464, 3.095231209092904, -2.1877889128373433, -0.3509645531366436, 1.108055523679191, 1.7507662229669712])
push!(qref,  [-0.3561800250302581, 0.3192951117867897, 3.089992817742055, -2.17817481530745, -0.39196455313664363, 1.1253706814064732, 1.7718854502272638])
push!(qref,  [-0.3934443596336305, 0.3254806584993247, 3.0847544263912066, -2.1685607177775563, -0.43296455313664367, 1.1426858391337553, 1.7930046774875568])
push!(qref,  [-0.43070869423700286, 0.33166620521185974, 3.0795160350403576, -2.158946620247663, -0.47396455313664365, 1.1600009968610374, 1.8141239047478495])
push!(qref,  [-0.4679730288403753, 0.3378517519243948, 3.0742776436895087, -2.1493325227177693, -0.5149645531366437, 1.1773161545883195, 1.8352431320081424])
push!(qref,  [-0.5052373634437477, 0.34403729863692983, 3.06903925233866, -2.1397184251878762, -0.5559645531366437, 1.1946313123156016, 1.856362359268435])
push!(qref,  [-0.54250169804712, 0.3502228453494649, 3.063800860987811, -2.1301043276579827, -0.5969645531366438, 1.2119464700428837, 1.8774815865287278])
push!(qref,  [-0.5797660326504925, 0.3564083920619999, 3.0585624696369624, -2.1204902301280892, -0.6379645531366438, 1.2292616277701658, 1.8986008137890207])
push!(qref,  [-0.6170303672538648, 0.36259393877453494, 3.0533240782861135, -2.1108761325981957, -0.6789645531366438, 1.2465767854974479, 1.9197200410493134])
push!(qref,  [-0.6542947018572371, 0.36877948548706996, 3.0480856869352646, -2.1012620350683022, -0.7199645531366438, 1.26389194322473, 1.940839268309606])
push!(qref,  [-0.6915590364606096, 0.37496503219960503, 3.0428472955844157, -2.091647937538409, -0.7609645531366438, 1.281207100952012, 1.961958495569899])
push!(qref,  [-0.728823371063982, 0.38115057891214005, 3.0376089042335668, -2.0820338400085157, -0.8019645531366438, 1.2985222586792942, 1.983077722830192])
push!(qref,  [-0.7660877056673544, 0.3873361256246751, 3.0323705128827183, -2.072419742478622, -0.8429645531366439, 1.3158374164065763, 2.0041969500904844])
push!(qref,  [-0.8033520402707267, 0.39352167233721014, 3.0271321215318694, -2.0628056449487286, -0.8839645531366439, 1.3331525741338583, 2.0253161773507773])
push!(qref,  [-0.8406163748740991, 0.39970721904974515, 3.0218937301810205, -2.053191547418835, -0.9249645531366439, 1.3504677318611407, 2.04643540461107])
push!(qref,  [-0.8778807094774715, 0.4058927657622802, 3.0166553388301716, -2.043577449888942, -0.965964553136644, 1.3677828895884225, 2.067554631871363])
push!(qref,  [-0.915145044080844, 0.41207831247481524, 3.0114169474793226, -2.0339633523590486, -1.006964553136644, 1.3850980473157049, 2.0886738591316556])
push!(qref,  [-0.9524093786842164, 0.4182638591873503, 3.006178556128474, -2.024349254829155, -1.0479645531366442, 1.402413205042987, 2.1097930863919485])
push!(qref,  [-0.9896737132875887, 0.42444940589988533, 3.0009401647776253, -2.0147351572992616, -1.088964553136644, 1.419728362770269, 2.1309123136522414])
push!(qref,  [-1.026938047890961, 0.43063495261242035, 2.9957017734267763, -2.005121059769368, -1.1299645531366442, 1.4370435204975511, 2.1520315409125343])
push!(qref,  [-1.0642023824943334, 0.43682049932495537, 2.9904633820759274, -1.9955069622394748, -1.1709645531366442, 1.4543586782248332, 2.173150768172827])
push!(qref,  [-1.101466717097706, 0.44300604603749044, 2.9852249907250785, -1.9858928647095815, -1.2119645531366443, 1.4716738359521153, 2.1942699954331197])
push!(qref,  [-1.1387310517010785, 0.4491915927500255, 2.97998659937423, -1.976278767179688, -1.2529645531366442, 1.4889889936793974, 2.215389222693412])
push!(qref,  [-1.1759953863044506, 0.45537713946256053, 2.974748208023381, -1.9666646696497945, -1.2939645531366444, 1.5063041514066795, 2.236508449953705])
push!(qref,  [-1.213259720907823, 0.46156268617509555, 2.969509816672532, -1.9570505721199012, -1.3349645531366443, 1.5236193091339616, 2.257627677213998])
push!(qref,  [-1.2478021278200226, 0.46729641720217435, 2.9646540572035307, -1.9481387245716801, -1.3729697583408722, 1.5396697023673502, 2.277204276379416])
push!(qref,  [-1.2748022884966306, 0.4717782027187738, 2.9608585403999896, -1.9411727561509207, -1.402676624155692, 1.5522155332738579, 2.292506375168055])
push!(qref,  [-1.2941678050105965, 0.47499270549426853, 2.9581362549951593, -1.9361765052728712, -1.423983489970512, 1.5612138684508734, 2.3034816078838536])
push!(qref,  [-1.3058986773619201, 0.4769399255286585, 2.9564872009890393, -1.9331499719375316, -1.4368903557853316, 1.566664707898397, 2.310129974526812])
push!(qref,  [-1.3099949055506017, 0.47761986282194374, 2.95591137838163, -1.932093156144902, -1.4413972216001514, 1.5685680516164289, 2.3124514750969296])
push!(qref,  [-1.3100000000000003, 0.47762070845504206, 2.955910662235214, -1.9320918417907373, -1.4414028267565981, 1.568570418791257, 2.312454362330413]);

q̇ref = typeof(zeros(7))[]
push!(q̇ref,  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
push!(q̇ref,  [-0.07634644162642143, 0.012672827411047383, -0.010732313987104905, 0.019697175427098694, -0.084, 0.03547495729491942, 0.04326866072840471])
push!(q̇ref,  [-0.15269288325284286, 0.025345654822094766, -0.02146462797420981, 0.03939435085419739, -0.168, 0.07094991458983883, 0.08653732145680942])
push!(q̇ref,  [-0.22903932487926432, 0.03801848223314215, -0.032196941961314715, 0.05909152628129609, -0.252, 0.10642487188475827, 0.12980598218521414])
push!(q̇ref,  [-0.3053857665056857, 0.05069130964418953, -0.04292925594841962, 0.07878870170839478, -0.336, 0.14189982917967767, 0.17307464291361885])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3726433460337236, 0.061855467125350325, -0.05238391350848822, 0.09614097529893409, -0.41, 0.17315157727282096, 0.21119227260292772])
push!(q̇ref,  [-0.3081748275792907, 0.05115426887151827, -0.04332132502896183, 0.07950827192114293, -0.3390686581481977, 0.14319578771253505, 0.17465531825059288])
push!(q̇ref,  [-0.2318283859528692, 0.038481441460470875, -0.03258901104185691, 0.05981109649404422, -0.2550686581481976, 0.1077208304176156, 0.13138665752218814])
push!(q̇ref,  [-0.1554819443264477, 0.02580861404942348, -0.021856697054752, 0.04011392106694551, -0.17106865814819755, 0.07224587312269615, 0.08811799679378339])
push!(q̇ref,  [-0.07913550270002623, 0.013135786638376092, -0.011124383067647084, 0.0204167456398468, -0.08706865814819748, 0.036770915827776696, 0.04484933606537864])
push!(q̇ref,  [-0.002789061073605048, 0.00046295922732875194, -0.00039206908054222056, 0.0007195702127481829, -0.0030686581481977893, 0.0012959585328574108, 0.0015806753369740911])
push!(q̇ref,  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]);

params = MPC_Params()

Xref = typeof(zeros(14))[]
for k = 1:length(qref)
    push!(Xref, [qref[k]; q̇ref[k]])
end

for k = 1:params.H
    push!(Xref, [qref[end]; q̇ref[end]])
end

# prob = ArthurProblem(Xref, params=params)
# altro = ALTROSolver(prob, params.opts)
# solve!(altro)
# Z_track = TrajectoryOptimization.get_trajectory(prob)
global x0 = [-0.012784051833563126, 0.2722778684209022, 3.170027245975045, -2.2336492735211575, 0.03408274619129692, 1.226885242471706, 1.5562949837009983, -0.0005575828038824746, 3.2918256112210627e-6, 0.0011079314450628355, -7.003196207031308e-6, -0.0003270387146453352, 6.465030966857781e-7, -5.557828320653166e-5]
global prob_mpc = ArthurHorizonProblem(Xref, x0, params.H, start=1)
global altro_mpc = ALTROSolver(prob_mpc, params.opts)
solve!(altro_mpc)
println(states(altro_mpc))
println(controls(altro_mpc))

# num_iters = length(Z_track) - prob_mpc.N
# X_traj = [zero(state(Z_track[1])) for i = 1:num_iters+1]
# X_traj[1] = state(Z_track[1])
# U_traj = [zero(control(Z_track[1])) for i = 1:num_iters]

# t0 = 0
# k_mpc = 1
# for i = 1:num_iters
#     global t0 += params.dt
#     global k_mpc += 1

#     # Update the ALTRO solution, advancing forward by 1 time step
#     U_traj[i] = control(prob_mpc.Z[1])
#     # mpc_update(altro_mpc, prob_mpc, Z_track, t0, k_mpc)
#     x0 = rk4(prob_mpc.model, state(prob_mpc.Z[1]), control(prob_mpc.Z[1]), prob_mpc.Z[1].dt)
#     println(x0)
#     mpc_update(altro_mpc, prob_mpc, Z_track, x0, t0)
#     X_traj[i+1] = prob_mpc.x0

#     solve!(altro_mpc)
#     # println(i)
# end

# X_traj = typeof(zero(state(Z_track[1])))[]
# push!(X_traj, state(Z_track[1]))
# U_traj = typeof(zero(control(Z_track[1])))[]

X_traj = typeof(Xref[1])[]
push!(X_traj, x0)
U_traj = typeof(zeros(7))[]

t0 = 0
global iter = 1
max_iters = length(Xref) + params.H
# errors = Float64[]
# push!(errors, Inf)
# push!(errors, norm(X_traj[end] - Xref[end]))
while norm(X_traj[end] - Xref[end]) > 0.1# && iter < max_iters#&& norm(errors[iter+1] - errors[iter]) > 0.003
    # println(iter)
    global iter += 1
    global t0 += params.dt
    # Update the ALTRO solution, advancing forward by 1 time step
    push!(U_traj, control(prob_mpc.Z[1]))
    # mpc_update(altro_mpc, prob_mpc, Z_track, t0, k_mpc)
    x0 = rk4(prob_mpc.model, state(prob_mpc.Z[1]), control(prob_mpc.Z[1]), params.dt)
    # x0 = [-0.012784051833563126, 0.2722778684209022, 3.170027245975045, -2.2336492735211575, 0.03408274619129692, 1.226885242471706, 1.5562949837009983, -0.0005575828038824746, 3.2918256112210627e-6, 0.0011079314450628355, -7.003196207031308e-6, -0.0003270387146453352, 6.465030966857781e-7, -5.557828320653166e-5]
    # x0 = Xref[1]
    # mpc_update(altro_mpc, prob_mpc, Z_track, cons, x0, t0)
    k_mpc = argmin(norm.([(Xref[k][1:7] - x0[1:7]) for k=1:length(Xref)]))
    # k_mpc = 1
    println(k_mpc)
    # println(norm(X_traj[end] - Xref[end]))
    # println(X_traj[end])
    # println(prob_mpc.Z[1].dt)
    global prob_mpc = ArthurHorizonProblem(Xref, x0, params.H, start=k_mpc)
    global altro_mpc = ALTROSolver(prob_mpc, params.opts)
    solve!(altro_mpc)
    # println(states(altro_mpc))
    # println(controls(altro_mpc))

    push!(X_traj, prob_mpc.x0)
    # println(U_traj[end])
    # push!(errors, norm(X_traj[end] - Xref[end]))
    # solve!(altro_mpc)
end
# for traj in X_traj
#     println(traj)
# end
# for traj in U_traj
#     println(traj)
# end

X_sim = typeof(Xref[1])[]
push!(X_sim, X_traj[1])
for k = 1:length(X_traj)-1
    # X_sim[k+1] = discrete_dynamics(TrajectoryOptimization.integration(prob_mpc), params.model, X_sim[k], U_traj[k], 0.0, params.dt)
    push!(X_sim, rk4(params.model, X_sim[k], U_traj[k], params.dt))
end
for traj in X_sim
    println(traj)
end
println(argmin(norm.([(X_sim[k] - Xref[end]) for k = 1:length(X_sim)])))
println(X_sim[argmin(norm.([(X_sim[k] - Xref[end]) for k = 1:length(X_sim)]))])
println(X_sim[argmin(norm.([(X_sim[k] - Xref[end]) for k = 1:length(X_sim)]))] - Xref[end])
println(norm(X_sim[argmin(norm.([(X_sim[k] - Xref[end]) for k = 1:length(X_sim)]))][1:7] - Xref[end][1:7]))
# println([(errors[k+1] - errors[k]) for k = 1:(length(errors)-1)])
# # TODO:
# # Uref could be just U equilibrium (just gravity compensation),
# # Xref could be the desired joint and cartesian positions and velocities, and then 0 for the forces

# state = MechanismState(model.mechanism)
# zero!(state)
# v̇ = similar(velocity(state))
# for joint in joints(state.mechanism)
#     if joint == joints(state.mechanism)[2] || joint == joints(state.mechanism)[3] || joint == joints(state.mechanism)[4] || joint == joints(state.mechanism)[5] || joint == joints(state.mechanism)[6] || joint == joints(state.mechanism)[7] || joint == joints(state.mechanism)[8]
#         v̇[joint][1] = 0.0
#     end
# end

# Qref = Vector{Float64}[]
# Uref = Vector{Float64}[]
# push!(Qref, copy(Xref[1][1:14]))
# for k = 1:N-1
#     set_configuration!(state, Qref[k][1:7])
#     set_velocity!(state, Qref[k][8:end])
#     push!(Uref, inverse_dynamics(state, v̇))
#     push!(Qref, rk4(model, Qref[k], Uref[k], dt))
# end

# dtref = [dt for k=1:N]
# traj = RobotDynamics.Traj(Xref, Uref, dtref, cumsum(dtref) .- dtref[1]);

# # Set up
# # TODO: Tune weights
# Q = 100.0*Diagonal(@SVector ones(n))
# Qf = 100*Diagonal(@SVector ones(n))
# R = 1.0e-1*Diagonal(@SVector ones(m))
# obj = TrajectoryOptimization.TrackingObjective(Q, R, traj, Qf=Qf)

# # Create Empty ConstraintList
# conSet = ConstraintList(n,m,N)

# # Control Bounds based on Robot Specs (Joint torque limits)
# u_bnd = [39.0, 39.0, 39.0, 39.0, 9.0, 9.0, 9.0]
# control_bnd = BoundConstraint(n,m, u_min=-u_bnd, u_max=u_bnd)
# add_constraint!(conSet, control_bnd, 1:N-1)

# # State Bounds based on Robot Specs (Joint velocity and speed limits)
# x_bnd = zeros(n)
# x_bnd[1:7] = [Inf, deg2rad(128.9), Inf, deg2rad(147.8), Inf, deg2rad(120.3), Inf] # rad
# x_bnd[8:14] = [1.39, 1.39, 1.39, 1.39, 1.22, 1.22, 1.22] # rad/sec
# # x_bnd[8:14] = [0.25, 0.05, 0.05, 0.05, 0.25, 0.15, 0.15] # rad/sec
# # x_bnd[15:end] = [Inf for k=1:(n-14)] # Constraints on force elsewhere
# state_bnd = BoundConstraint(n,m, x_min=-x_bnd, x_max=x_bnd)
# add_constraint!(conSet, state_bnd, 1:N)

# # # Cartesian Velocity Bound
# # ẋ_max = 0.0005 # m/s
# # vel_bnd = NormConstraint(n, m, ẋ_max, Inequality(), 21:23)
# # add_constraint!(conSet, vel_bnd, 1:N)

# # # Force Bound (Fx Fy Fz)
# # F_max = 20 # Newtons
# # F_bnd = NormConstraint(n, m, F_max, Inequality(), 27:29)
# # add_constraint!(conSet, F_bnd, 1:N)

# # Other possible constraints:
# # Cartesian angular velocity bounds?
# # External moment bounds?

# # Goal Constraint - only if you want the final state to be the desired state
# # goal = GoalConstraint(xf)
# # add_constraint!(conSet, goal, N)

# prob = Problem(model, obj, Xref[end], tf, x0=Xref[1], constraints=conSet, X0=Xref, U0=Uref)

# opts = SolverOptions(
#     cost_tolerance_intermediate=1e-2,
#     penalty_scaling=10.,
#     penalty_initial=1.0,
#     projected_newton=false,
# )

# altro = ALTROSolver(prob, opts)
# set_options!(altro, show_summary=true)
# for k = 1:10
#     solve!(altro);
# end

# # Get some info on the solve
# max_violation(altro)
# cost(altro)
# iterations(altro)

# # Extract the solution
# X = states(altro)
# U = controls(altro)

# # Extract the solver statistics
# stats = Altro.stats(altro)   # alternatively, solver.stats
# stats.iterations             # equivalent to iterations(solver)
# stats.iterations_outer       # (number of Augmented Lagrangian iterations)
# stats.iterations_pn          # (number of projected newton iterations)
# stats.cost[end]              # terminal cost
# stats.c_max[end]             # terminal constraint satisfaction
# stats.gradient[end]          # terminal gradient of the Lagrangian
# dstats = Dict(stats)         # get the per-iteration stats as a dictionary (can be converted to DataFrame)

# using Printf
# for i in 1:N
#     print("x")
#     print(i-1)
#     print(":\t")
#     println(round.(X[i][8:14]; digits=5))
# end
# println()

# for i in 1:N-1
#     print("u")
#     print(i)
#     print(": ")
#     println(U[i])
# end
# println()

# for i in 1:N
#     print("x")
#     print(i-1)
#     print(": ")
#     println(norm(X[i] - Xref[i]))
# end
# println()

# print("Terminal Error: ")
# println(norm(X[end] - Xref[end]))
# println(X[end] - Xref[end])

# for i in 1:N-1
#     print("u")
#     print(i)
#     print(": ")
#     println(U[i] - Uref[i])
# end

# using MeshCat, MeshCatMechanisms, Blink
# urdf = "../../../kortex_description/arms/gen3/7dof/urdf/GEN3_URDF_V12 copy.urdf"
# body = findbody(model.mechanism, "Bracelet_Link")
# point = Point3D(default_frame(body), 0., 0, -0.07)
# # Create the visualizer
# vis = MechanismVisualizer(model.mechanism, URDFVisuals(urdf))

# # Render our target point attached to the robot as a sphere with radius 0.07
# setelement!(vis, point, 0.01)

# qs = Vector{Float64}[]
# for x in X
#     push!(qs, copy(x[1:7]))
# end
# ts = collect(0:dt:(N-1)*dt);

# setanimation!(vis, Animation(vis, ts, qs))
# render(vis)


text:
module RDates

include("abstracts.jl")
# Rounding and Compounds defined upfront as we need them for the grammar
include("rounding.jl")
include("compounds.jl")

include("grammar.jl")
include("monthinc.jl")
include("invalidday.jl")

# The various basic implementations (along with shows and grammar registrations)
include("primitives.jl")

include("calendars.jl")
include("ranges.jl")
# include("io.jl")
# Export the macro and non-macro parsers.
export @rd_str
export rdate
export is_holiday
export calendar
export apply
export with_cal_mgr

include("build.jl")

end # module RDates


text:
mutable struct History{A,O}
    actions::Vector{A}
    observations::Vector{O}
    History{A,O}() where {A,O} = new(Array{A}(0), Array{O}(0))
end

function add{A,O}(history::History{A,O},
            action::A,
            obs::O)

    push!(history.actions, action)
    push!(history.observations, obs)
end

function remove_last{A,O}(history::History{A,O})
    pop!(history.actions)
    pop!(history.observations)
end

function history_size{A,O}(history::History{A,O})
    return length(history.actions)
end

function truncate{A,O}(history::History{A,O}, d::Int64)
    history.actions = history.actions[1:d]
    history.observations = history.observations[1:d]
end


text:
"""
managing  uniform values - global values in shaders
"""
module Uniforms
using StaticArrays,ModernGL,Rocket,  ..ForDisplayStructs, Dictionaries, Parameters, ColorTypes

export changeMainTextureContribution,changeTextureContribution,coontrolMinMaxUniformVals,createStructsDict, setCTWindow,setMaskColor,setTextureVisibility, setTypeOfMainSampler!
export @uniforms
export @uniforms!

########## this from https://github.com/jorge-brito/Glutils.jl/blob/4b7da5895a3e927792994ad0b6643498a8125362/src/utils.jl
struct TypeRef{T} end

Base.adjoint(x::Ref) = x[]
Base.adjoint(::Type{T}) where {T} = TypeRef{T}()
Base.getindex(::TypeRef{T}) where {T} = Ref{T}()
Base.getindex(::TypeRef{T}, x::S) where {T, S} = Ref{T}(x)

const SymString = Union{AbstractString, Symbol}

macro map(expr...)
    T, ex =
        if length(expr) == 1
            nothing, first(expr)
        else
            first(expr), last(expr)
    end

    args = map(ex.args) do arg
        if Meta.isexpr(arg, :(=))
            key, value = arg.args
            return :( $(QuoteNode(key)) => $value )
        elseif Meta.isexpr(arg, :(:=))
            key, value = arg.args
            return :( $key => $value )
        else
            error("Invalid syntax for @map, expected expression of type `foo = bar` or `foo := bar`")
        end
    end

    result =
        if isnothing(T)
            Expr(:call, :Dict, args...)
        else
            Expr(:call, Expr(:curly, :Dict, T.args...), args...)
        end

    return esc(result)
end

const GL_ENUM_TYPE = @map Symbol, GLenum {
    GLboolean = GL_BOOL,
    GLchar    = GL_BYTE,
    GLdouble  = GL_DOUBLE,
    GLfloat   = GL_FLOAT,
    GLint     = GL_INT,
    GLshort   = GL_SHORT,
    GLubyte   = GL_UNSIGNED_BYTE,
    GLuint    = GL_UNSIGNED_INT,
    GLushort  = GL_UNSIGNED_SHORT,
}

macro genumof(type)
    :( GL_ENUM_TYPE[ $(QuoteNode(type)) ] )
end

genumof(::Type{Cint})    = GL_INT
genumof(::Type{Cchar})   = GL_BYTE
genumof(::Type{Cfloat})  = GL_FLOAT
genumof(::Type{Cshort})  = GL_SHORT
genumof(::Type{Int64})   = GL_DOUBLE
genumof(::Type{Cuint})   = GL_UNSIGNED_INT
genumof(::Type{Cuchar})  = GL_UNSIGNED_BYTE
genumof(::Type{Cushort}) = GL_UNSIGNED_SHORT

gbool(x::GLenum) = x == GL_TRUE
gbool(x::Bool) = x ? GL_TRUE : GL_FALSE

"""
        @g_str(name) -> GLenum
Returns a OpenGL constant.
# Examples
```julia
g"bool" == GL_BOOL
g"Texture 2D" == GL_TEXTURE_2D
g"Clamp.to.edge" == GL_CLAMP_TO_EDGE
g"texture-wrap-t" == GL_TEXTURE_WRAP_T
```
"""
macro g_str(ex)
    name = uppercase(replace(ex, r"(\s|\.|\-)+" => "_"))
    esc( Symbol("GL_$name") )
end

############# part below directly copied from https://github.com/jorge-brito/Glutils.jl/blob/master/src/uniforms.jl

# this function is used to get the
# suffix of the glUniform function
# for the correct type
_type_suffix(::Type{<:Bool}) = "i"
_type_suffix(::Type{<:Signed}) = "i"
_type_suffix(::Type{<:Unsigned}) = "ui"
_type_suffix(::Type{<:AbstractFloat}) = "f"

"""
        uniform!(location, values...) -> Nothing
Set the uniform `vec` variable at `location`.
# Examples
```julia
myuniform = getuniform(program, "someUniform")
# set float vec3
uniform!(myuniform, 1.0, 1.0, 0.5)
# set float vec4
uniform!(myuniform, Cfloat[1, 2, 3, 4])
# set a 4x4 matrix
uniform!(myuniform, rand(Cfloat, 4, 4))
```
"""
function uniform! end

@generated function uniform!(location, values::Vararg{T, N}) where {N, T <: Real}
    suffix = _type_suffix(T)
    glFunc = Symbol("glUniform$(N)$suffix")
    return :( $(glFunc)(location, values...) )
end

@generated function uniform!(location, vector::SVector{N, T}) where {N, T <: Real}
    suffix = _type_suffix(T)
    glFunc = Symbol("glUniform$(N)$(suffix)v")
    return :( $(glFunc)(location, 1, vector) )
end

@generated function uniform!(location, matrix::SMatrix{N, M, T}, transposed::Bool = false) where {N, M, T <: Real}
    glFunc = Symbol(
        if N == M
            "glUniformMatrix$(N)fv"
        else
            "glUniformMatrix$(N)x$(M)fv"
        end
    )
    return :( $(glFunc)(location, 1, transposed, Cfloat[matrix...]) )
end

function uniform!(location, vector::AbstractVector{T}) where {T <: Real}
    N = length(vector)
    uniform!(location, SVector{N, T}(vector...))
end

function uniform!(location, matrix::AbstractMatrix{T}, transposed::Bool = false) where {T <: Real}
    N, M = size(matrix)
    uniform!(location, SMatrix{N, M, T}(matrix), transposed)
end
"""
        getuniform(program, name) -> GLint
Gets the location of the uniform variable
identified by `name`.
"""
getuniform(program, name::SymString) = glGetUniformLocation(program, string(name))
"""
        @uniforms foo, bar, ... = program
Get the location of `foo`, `bar` and `etc` from `program`,
and assign each uniform to the corresponding variable.
This macro transform the following expression:
```julia
@uniforms foo, bar, qux: baz = program
```
Into:
```julia
foo = getuniform(program, :foo)
bar = getuniform(program, :bar)
qux = getuniform(program, :baz)
```
You can also rename the variables as following:
```julia
@uniforms a: foo, b: bar = program
```
The variable `a` contains the location of `foo`,
and `b` contains the loctation of `bar`.
"""
macro uniforms(ex)
    @assert ex.head == :(=)
    program = esc(ex.args[2])
    if !Meta.isexpr(ex.args[1], :tuple)
        args = [ex.args[1]]
    else
        args = ex.args[1].args
    end
    return Expr(:block,
        map(args) do arg
            if arg isa Symbol
                alias, name = arg, arg
            elseif Meta.isexpr(arg, :(=))
                alias, name = arg.args
            elseif Meta.isexpr(arg, :call)
                alias, name = arg.args[2:3]
            end
            return Expr(:(=), esc(alias), Expr(:call, :getuniform, program, QuoteNode(name)))
        end...
    )
end

macro uniforms!(ex)
    args = ex.args
    for i in eachindex(args)
        if Meta.isexpr(args[i], :(:=))
            name, value = args[i].args
            args[i] = Expr(:call, :uniform!, name, value)
        end
    end
    return esc(Expr(:block, args...))
end

############### my original part

"""
function cotrolling the window  for displaying CT scan  - min white and max max_shown_black
    uniformsStore - instantiated object holding references to uniforms controlling displayed window
 """
function setCTWindow(min_shown_whiteInner::Int32, max_shown_blackInner::Int32, uniformsStore ::MainImageUniforms)
    @uniforms! begin
    uniformsStore.min_shown_white:=min_shown_whiteInner
    uniformsStore.max_shown_black:=max_shown_blackInner
    uniformsStore.displayRange:= Float32(min_shown_whiteInner-max_shown_blackInner )
    end
end

"""
sets color of the mask

"""
function setMaskColor(color::RGB, uniformsStore ::MaskTextureUniforms)
    @uniforms! begin
    uniformsStore.colorsMaskRef:=Cfloat[color.r, color.g, color.b, 0.8]
    end

end#setMaskColor

"""
sets visibility of the texture
"""
function setTextureVisibility(isvisible::Bool, uniformsStore ::TextureUniforms)
    @uniforms! begin
    uniformsStore.isVisibleRef:= isvisible ? 1 : 0
    end

end#setTextureVisibility

"""
sets minimum and maximum value for display -
    in case of continuus colors it will clamp values - so all above max will be equaled to max ; and min if smallert than min
    in case of main CT mask - it will controll min shown white and max shown black
    in case of maks with single color associated we will step data so if data is outside the rande it will return 0 - so will not affect display
"""
function coontrolMinMaxUniformVals(textur::TextureSpec)
    newMin=textur.minAndMaxValue[1]
    newMax=textur.minAndMaxValue[2]
    uniformsStore= textur.uniforms
    range = newMax-newMin
    if(range<1) range=1 end#if

    @uniforms! begin
    uniformsStore.maskMinValue:= newMin
    uniformsStore.maskMAxValue:= newMax
    uniformsStore.maskRangeValue:= range
    end

end#coontrolMinMaxUniformVals

"""
controlls contribution  of given mask to the overall image - maximum value is 1 minimum 0 if we have 3 masks and all control contribution is set to 1 and all are visible their corresponding influence to pixel color is 33%
      if plus is pressed it will increse contribution by 0.1
      if minus is pressed it will decrease contribution by 0.1
it also modifies given TextureSpec
change - how should the texture spec be modified
"""
function changeTextureContribution(textur::TextureSpec, change::Float32)
    newValue= textur.maskContribution + change
    if(newValue>=0 &&newValue <=1)
        textur.maskContribution =  newValue
        @uniforms! begin
        textur.uniforms.maskContribution:= newValue

        end#@uniforms!
    end#if

end#changeTextureContribution

function changeMainTextureContribution(textur::TextureSpec, change::Float32,actor::SyncActor{Any, ActorWithOpenGlObjects} )
    newValue= textur.maskContribution + change
    if(newValue>=0 &&newValue <=1)
        textur.maskContribution =  newValue
        @uniforms! begin
        actor.actor.mainForDisplayObjects.mainImageUniforms.mainImageContribution:= newValue
        end#@uniforms!
    end#if

end#changeTextureContribution

end #module