text:
module Definitions
#
# This module contains enumerated types and constants that are used throughout the model.
#
using ScottishTaxBenefitModel
using ScottishTaxBenefitModel.Utils

export
   Employment_Status,  # mapped from empstat
   Employee,
   Missing_Employment_Status,
   Self_employed,
   FRS_DIR,
   HBAI_DIR,
   MATCHING_DIR,
   MODEL_DATA_DIR,
   MODEL_DATA_DIR,
   MODEL_PARAMS_DIR,
   PRICES_DIR,
   PROJECT_DIR,
   TEST_DIR,
   DOCS_DIR,
   SCRIPTS_DIR

const MODEL_NAME       = "Scottish Tax Benefit Model"
# FIXME these might be better at the top of ScottishTaxBenefitModel.jl itself
const PROJECT_DIR      = joinpath(dirname(pathof(ScottishTaxBenefitModel)),".." )
const MODEL_DATA_DIR   = joinpath( PROJECT_DIR, "data" )
const PRICES_DIR       = joinpath( MODEL_DATA_DIR, "prices", "indexes" )
const MATCHING_DIR     = joinpath( MODEL_DATA_DIR, "merging" )
const MODEL_PARAMS_DIR = joinpath( PROJECT_DIR, "params" )
const SCRIPTS_DIR      = joinpath( PROJECT_DIR, "scripts" )
const TEST_DIR         = joinpath( PROJECT_DIR, "test" )
const DOCS_DIR         = joinpath( PROJECT_DIR, "docs" )

# FIXME these should be in some config file probably.
const RAW_DATA = "/mnt/data/"
const FRS_DIR = "$RAW_DATA/frs/"
const HBAI_DIR = "$RAW_DATA/hbai/"

export NUM_REPEATS
const NUM_REPEATS = 30 # simulates a longer calculation

export LegacyOrUC, legacy_bens, uc_bens
@enum LegacyOrUC legacy_bens uc_bens

export old_style, new_style, OldOrNew
@enum OldOrNew begin
   missing_style = -1
   old_style = 1
   new_style = 2
end

export BereavementType,
   missing_bereave,
   bereavement_allowance,
   widowed_parents,
   bereavement_support

@enum BereavementType begin
   missing_bereave = -1
   bereavement_allowance = 1
   widowed_parents = 2
   bereavement_support = 3 # old style 1 off payment. Only recorded in 2015 FRS I think
end

export Illness_Length  # mapped from limitl
export Less_than_six_months, Between_six_months_and_12_months, v_12_months_or_more
export Missing_Illness_Length, RT

@enum Illness_Length begin  # mapped from limitl
   Missing_Illness_Length = -1
   Less_than_six_months = 1
   Between_six_months_and_12_months = 2
   v_12_months_or_more = 3
end

@enum PIPType begin
    no_pip = -1
    standard_pip = 1
    enhanced_pip = 2
end
export PIPType, no_pip, standard_pip, enhanced_pip

@enum LowMiddleHigh begin
    missing_lmh = -1 # fixme global change to "no_lmh"
    low = 1
    mid = 2
    high = 3
end
export LowMiddleHigh, missing_lmh, low, mid, high

export ADLS_Inhibited  # mapped from condit
export reduced_a_lot, reduced_a_little, not_reduced
export Missing_ADLS_Inhibited

@enum ADLS_Inhibited begin  # mapped from condit
   Missing_ADLS_Inhibited = -1
   reduced_a_lot = 1
   reduced_a_little = 2
   not_reduced = 3
end

@enum Employment_Status begin  # mapped from empstat
   Missing_Employment_Status = -1
   Employee = 1
   Self_employed = 2
end

export ILO_Employment  # mapped from empstati
export Full_time_Employee,
       Part_time_Employee,
       Full_time_Self_Employed,
       Part_time_Self_Employed,
       Unemployed,
       Retired,
       Student,
       Looking_after_family_or_home,
       Permanently_sick_or_disabled,
       Temporarily_sick_or_injured,
       Other_Inactive
export Missing_ILO_Employment

@enum ILO_Employment begin  # mapped from empstati
   Missing_ILO_Employment = -1
   Full_time_Employee = 1
   Part_time_Employee = 2
   Full_time_Self_Employed = 3
   Part_time_Self_Employed = 4
   Unemployed = 5
   Retired = 6
   Student = 7
   Looking_after_family_or_home = 8
   Permanently_sick_or_disabled = 9
   Temporarily_sick_or_injured = 10
   Other_Inactive = 11
end

export Ethnic_Group  # mapped from ethgr3
export White,
       Mixed_or_Multiple_ethnic_groups,
       Asian_or_Asian_British,
       Black_or_African_or_Caribbean_or_Black_British,
       Other_ethnic_group
export Missing_Ethnic_Group

@enum Ethnic_Group begin  # mapped from ethgr3
   Missing_Ethnic_Group = -1
   White = 1
   Mixed_or_Multiple_ethnic_groups = 2
   Asian_or_Asian_British = 3
   Black_or_African_or_Caribbean_or_Black_British = 4
   Other_ethnic_group = 5
end

export Ethnic_Scotland  # mapped from ethgrps
export Scottish,
       Other_British,
       Irish,
       Gypsy_or_Traveller,
       Polish,
       Any_other_white_ethnic_group,
       Any_mixed_or_multiple_ethnic_group,
       Pakistani_Pakistani_Scottish_or_Pakistani_British,
       Indian_Indian_Scottish_or_Indian_British,
       Bangladeshi_Bangladeshi_Scottish_or_Bangladeshi_British,
       Chinese_Chinese_Scottish_or_Chinese_British,
       Any_other_Asian,
       African_African_Scottish_or_African_British,
       Any_other_African,
       Caribbean_Caribbean_Scottish_or_Caribbean_British,
       Black_Black_Scottish_or_Black_British,
       Any_other_Caribbean_or_Black,
       Arab_Arab_Scottish_or_Arab_British,
       Any_other_ethnic_group
export Missing_Ethnic_Scotland

@enum Ethnic_Scotland begin  # mapped from ethgrps
   Missing_Ethnic_Scotland = -1
   Scottish = 1
   Other_British = 2
   Irish = 3
   Gypsy_or_Traveller = 4
   Polish = 5
   Any_other_white_ethnic_group = 6
   Any_mixed_or_multiple_ethnic_group = 7
   Pakistani_Pakistani_Scottish_or_Pakistani_British = 8
   Indian_Indian_Scottish_or_Indian_British = 9
   Bangladeshi_Bangladeshi_Scottish_or_Bangladeshi_British = 10
   Chinese_Chinese_Scottish_or_Chinese_British = 11
   Any_other_Asian = 12
   African_African_Scottish_or_African_British = 13
   Any_other_African = 14
   Caribbean_Caribbean_Scottish_or_Caribbean_British = 15
   Black_Black_Scottish_or_Black_British = 16
   Any_other_Caribbean_or_Black = 17
   Arab_Arab_Scottish_or_Arab_British = 18
   Any_other_ethnic_group = 19
end

export Sex  # mapped from sex
export Male, Female
export Missing_Sex

@enum Sex begin  # mapped from sex
   Missing_Sex = -1
   Male = 1
   Female = 2
end

export Health_Status  # mapped from heathad
export Very_Good, Good, Fair, Bad, Very_Bad
export Missing_Health_Status

@enum Health_Status begin  # mapped from heathad
   Missing_Health_Status = -1
   Very_Good = 1
   Good = 2
   Fair = 3
   Bad = 4
   Very_Bad = 5
end

export Marital_Status  # mapped from marital
export Married_or_Civil_Partnership,
       Cohabiting,
       Single,
       Widowed,
       Separated,
       Divorced_or_Civil_Partnership_dissolved
export Missing_Marital_Status

@enum Marital_Status begin  # mapped from marital
   Missing_Marital_Status = -1
   Married_or_Civil_Partnership = 1
   Cohabiting = 2
   Single = 3
   Widowed = 4
   Separated = 5
   Divorced_or_Civil_Partnership_dissolved = 6
end

export Qualification_Type  # mapped from dvhiqual
export Doctorate_or_MPhil,
       Masters_PGCE_or_some_other_post_graduate_degree,
       Degree_inc_foundation_degree_or_professional_institute_member,
       Teaching_qualification_excluding_PGCE,
       Any_other_foreign_qualification_at_degree_level,
       Any_other_work_related_or_vocational_qualification_at_degr,
       Any_other_professional_qualification_at_degree_level,
       Other_Higher_Education_qualification_below_degree_level,
       Nursing_or_other_medical_qualification_not_yet_mentioned,
       Diploma_in_higher_education,
       HNC_or_HND,
       BTEC_BEC_TEC_EdExcel_or_LQL_at_higher_level_Higher_Lev,
       SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_higher_level_level,
       NVQ_or_SVQ_Level_4,
       NVQ_or_SVQ_Level_5,
       RSA_or_OCR_a_higher_diploma_or_OCR_Level_4,
       A_Level_or_GCE_in_Applied_Subjects_or_equivalent,
       Welsh_Baccalaureate_at_Advanced_level,
       Scottish_Baccalaureate,
       International_Baccalaureate,
       AS_level_or_equivalent,
       Certificate_of_6th_Year_Studies_CSYS_Scotland,
       Access_to_Higher_Education,
       Advanced_Higher_or_Higher_or_Intermediate_or_Access_Qualification_Sc,
       Skills_for_work_Higher,
       ONC_or_OND,
       BTEC_BEC_TEC_EdExcel_or_LQL_at_Nat_Cert_or_Nat_Dipl_leve,
       SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_full_Nat_Cert_Level,
       New_Diploma_an_Advanced_diploma_level3,
       New_Diploma_a_Progression_diploma_level_3,
       NVQ_or_SVQ_Level_3,
       GNVQ_or_GSVQ_Advanced,
       RSA_or_OCR_an_advanced_diploma_or_advanced_certificate_or_OCR_L,
       City_and_Guilds_advanced_craft_or_part_3,
       Welsh_Baccalaureate_at_the_Intermediate_level,
       O_Level_or_equivalent_five_or_more,
       Standard_Grade_or_Ordinary_Grade_or_Lower_Scotland_five_or,
       GCSE_or_equivalent_five_or_more,
       CSE_five_or_more,
       Scottish_National_level_5,
       Skills_for_work_National_level_5,
       BTEC_BEC_TEC_EdExcel_or_LQL_first_diploma_or_general_diplo,
       SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_first_diploma_or_genera,
       New_Diploma_a_Higher_diploma_level_2,
       NVQ_or_SVQ_Level_2,
       GNVQ_or_GSVQ_Full_Intermediate,
       RSA_or_OCR_a_diploma_or_OCR_Level_2,
       City_and_Guilds_craft_or_part_2,
       Any_other_qualification_high_school_leavers_qualification,
       BTEC_BEC_TEC_EdExcel_or_LQL,
       BTEC_BEC_TEC_EdExcel_or_LQL_first_cert_or_general_cert_Le,
       SCOTVEC_SCOTEC_or_SCOTBEC_Scotland,
       SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_first_cert_or_general_c,
       SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_modules_towards_Nat,
       New_Diploma,
       New_Diploma_Foundation_Diploma_level_1,
       Welsh_Baccalaureate,
       Welsh_Baccalaureate_at_the_foundation_level,
       NVQ_or_SVQ,
       NVQ_or_SVQ_Level_1,
       GNVQ_or_GSVQ,
       GNVQ_or_GSVQ_Part_One_Intermediate,
       GNVQ_or_GSVQ_Full_Foundation,
       GNVQ_or_GSVQ_Part_One_Foundation,
       O_Level_or_equivalent,
       O_Level_or_equivalent_fewer_than_5,
       Standard_Grade_or_Ordinary_Grade_or_Lower_Scotland,
       Standard_Grade_or_Ordinary_Grade_or_Lower_Scotland_fewer_t,
       GCSE_or_equivalent,
       GCSE_or_equivalent_fewer_than_5,
       Scottish_National_level_1_to_4,
       Scottish_National_level,
       Skills_for_work_National_level_3_and_4,
       Skills_for_work,
       CSE,
       CSE_fewer_than_5,
       RSA_or_OCR,
       RSA_or_OCR_some_other_RSA_including_Stage_I_II_and_III_or_OC,
       City_and_Guilds,
       City_and_Guilds_foundation_or_part_1,
       YT_Certificate_or_YTP,
       Key_Skills_or_Core_Skills_Scotland_or_Essential_Skills_Wales,
       Basic_Skills_skills_for_life_or_literacy_or_numeracy_or_language,
       Entry_Level_Qualifications,
       Award_Certificate_or_Diploma_at_entry_level_levels_1_to_8,
       Any_other_professional_or_vocational_or_foreign_qualifications
export Missing_Highest_Qualification

@enum Qualification_Type begin  # mapped from dvhiqual
   Missing_Highest_Qualification = -1
   Doctorate_or_MPhil = 1
   Masters_PGCE_or_some_other_post_graduate_degree = 2
   Degree_inc_foundation_degree_or_professional_institute_member = 3
   Teaching_qualification_excluding_PGCE = 4
   Any_other_foreign_qualification_at_degree_level = 5
   Any_other_work_related_or_vocational_qualification_at_degr = 6
   Any_other_professional_qualification_at_degree_level = 7
   Other_Higher_Education_qualification_below_degree_level = 8
   Nursing_or_other_medical_qualification_not_yet_mentioned = 9
   Diploma_in_higher_education = 10
   HNC_or_HND = 11
   BTEC_BEC_TEC_EdExcel_or_LQL_at_higher_level_Higher_Lev = 12
   SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_higher_level_level = 13
   NVQ_or_SVQ_Level_4 = 14
   NVQ_or_SVQ_Level_5 = 15
   RSA_or_OCR_a_higher_diploma_or_OCR_Level_4 = 16
   A_Level_or_GCE_in_Applied_Subjects_or_equivalent = 17
   Welsh_Baccalaureate_at_Advanced_level = 18
   Scottish_Baccalaureate = 19
   International_Baccalaureate = 20
   AS_level_or_equivalent = 21
   Certificate_of_6th_Year_Studies_CSYS_Scotland = 22
   Access_to_Higher_Education = 23
   Advanced_Higher_or_Higher_or_Intermediate_or_Access_Qualification_Sc = 24
   Skills_for_work_Higher = 25
   ONC_or_OND = 26
   BTEC_BEC_TEC_EdExcel_or_LQL_at_Nat_Cert_or_Nat_Dipl_leve = 27
   SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_full_Nat_Cert_Level = 28
   New_Diploma_an_Advanced_diploma_level3 = 29
   New_Diploma_a_Progression_diploma_level_3 = 30
   NVQ_or_SVQ_Level_3 = 31
   GNVQ_or_GSVQ_Advanced = 32
   RSA_or_OCR_an_advanced_diploma_or_advanced_certificate_or_OCR_L = 33
   City_and_Guilds_advanced_craft_or_part_3 = 34
   Welsh_Baccalaureate_at_the_Intermediate_level = 35
   O_Level_or_equivalent_five_or_more = 36
   Standard_Grade_or_Ordinary_Grade_or_Lower_Scotland_five_or = 37
   GCSE_or_equivalent_five_or_more = 38
   CSE_five_or_more = 39
   Scottish_National_level_5 = 40
   Skills_for_work_National_level_5 = 41
   BTEC_BEC_TEC_EdExcel_or_LQL_first_diploma_or_general_diplo = 42
   SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_first_diploma_or_genera = 43
   New_Diploma_a_Higher_diploma_level_2 = 44
   NVQ_or_SVQ_Level_2 = 45
   GNVQ_or_GSVQ_Full_Intermediate = 46
   RSA_or_OCR_a_diploma_or_OCR_Level_2 = 47
   City_and_Guilds_craft_or_part_2 = 48
   Any_other_qualification_high_school_leavers_qualification = 49
   BTEC_BEC_TEC_EdExcel_or_LQL = 50
   BTEC_BEC_TEC_EdExcel_or_LQL_first_cert_or_general_cert_Le = 51
   SCOTVEC_SCOTEC_or_SCOTBEC_Scotland = 52
   SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_first_cert_or_general_c = 53
   SCOTVEC_SCOTEC_or_SCOTBEC_Scotland_modules_towards_Nat = 54
   New_Diploma = 55
   New_Diploma_Foundation_Diploma_level_1 = 56
   Welsh_Baccalaureate = 57
   Welsh_Baccalaureate_at_the_foundation_level = 58
   NVQ_or_SVQ = 59
   NVQ_or_SVQ_Level_1 = 60
   GNVQ_or_GSVQ = 61
   GNVQ_or_GSVQ_Part_One_Intermediate = 62
   GNVQ_or_GSVQ_Full_Foundation = 63
   GNVQ_or_GSVQ_Part_One_Foundation = 64
   O_Level_or_equivalent = 65
   O_Level_or_equivalent_fewer_than_5 = 66
   Standard_Grade_or_Ordinary_Grade_or_Lower_Scotland = 67
   Standard_Grade_or_Ordinary_Grade_or_Lower_Scotland_fewer_t = 68
   GCSE_or_equivalent = 69
   GCSE_or_equivalent_fewer_than_5 = 70
   Scottish_National_level_1_to_4 = 71
   Scottish_National_level = 72
   Skills_for_work_National_level_3_and_4 = 73
   Skills_for_work = 74
   CSE = 75
   CSE_fewer_than_5 = 76
   RSA_or_OCR = 77
   RSA_or_OCR_some_other_RSA_including_Stage_I_II_and_III_or_OC = 78
   City_and_Guilds = 79
   City_and_Guilds_foundation_or_part_1 = 80
   YT_Certificate_or_YTP = 81
   Key_Skills_or_Core_Skills_Scotland_or_Essential_Skills_Wales = 82
   Basic_Skills_skills_for_life_or_literacy_or_numeracy_or_language = 83
   Entry_Level_Qualifications = 84
   Award_Certificate_or_Diploma_at_entry_level_levels_1_to_8 = 85
   Any_other_professional_or_vocational_or_foreign_qualifications = 86
end

export Standard_Region  # mapped from gvtregn
export North_East,
       North_West,
       Yorks_and_the_Humber,
       East_Midlands,
       West_Midlands,
       East_of_England,
       London,
       South_East,
       South_West,
       Scotland,
       Wales,
       Northern_Ireland
export Missing_Standard_Region

@enum Standard_Region begin  # mapped from gvtregn
   Missing_Standard_Region = -1
   North_East = 112000001
   North_West = 112000002
   Yorks_and_the_Humber = 112000003
   East_Midlands = 112000004
   West_Midlands = 112000005
   East_of_England = 112000006
   London = 112000007
   South_East = 112000008
   South_West = 112000009
   Scotland = 299999999
   Wales = 399999999
   Northern_Ireland = 499999999
end

export Socio_Economic_Group  # mapped from nssec
export Employers_in_large_organisations,
       Higher_managerial_occupations,
       Higher_professional_occupations_New_self_employed,
       Lower_prof_and_higher_technical_Traditional_employee,
       Lower_managerial_occupations,
       Higher_supervisory_occupations,
       Intermediate_clerical_and_administrative,
       Employers_in_small_organisations_non_professional,
       Own_account_workers_non_professional,
       Lower_supervisory_occupations,
       Lower_technical_craft,
       Semi_routine_sales,
       Routine_sales_and_service,
       Never_worked,
       Full_time_student,
       Not_classified_or_inadequately_stated,
       Not_classifiable_for_other_reasons
export Missing_Socio_Economic_Group

@enum Socio_Economic_Group begin  # mapped from nssec
   Missing_Socio_Economic_Group = -1
   Employers_in_large_organisations = 1
   Higher_managerial_occupations = 2
   Higher_professional_occupations_New_self_employed = 3
   Lower_prof_and_higher_technical_Traditional_employee = 4
   Lower_managerial_occupations = 5
   Higher_supervisory_occupations = 6
   Intermediate_clerical_and_administrative = 7
   Employers_in_small_organisations_non_professional = 8
   Own_account_workers_non_professional = 9
   Lower_supervisory_occupations = 10
   Lower_technical_craft = 11
   Semi_routine_sales = 12
   Routine_sales_and_service = 13
   Never_worked = 14
   Full_time_student = 15
   Not_classified_or_inadequately_stated = 16
   Not_classifiable_for_other_reasons = 17
end

export SIC_2007  # mapped from sic
export Undefined_SIC,
       Crop_and_animal_production_hunting_and_related_serviceX,
       Forestry_and_logging,
       Fishing_and_aquaculture,
       Mining_of_coal_and_lignite,
       Extraction_of_crude_petroleum_and_natural_gas,
       Mining_of_metal_ores,
       Other_mining_and_quarrying,
       Mining_support_service_activities,
       Manufacture_of_food_products,
       Manufacture_of_beverages,
       Manufacture_of_tobacco_products,
       Manufacture_of_textiles,
       Manufacture_of_wearing_apparel,
       Manufacture_of_leather_and_related_products,
       Manufacture_of_wood_and_of_products_of_wood_and_corkX,
       Manufacture_of_paper_and_paper_products,
       Printing_and_reproduction_of_recorded_media,
       Manufacture_of_coke_and_refined_petroleum_products,
       Manufacture_of_chemicals_and_chemical_products,
       Manufacture_of_basic_pharmaceutical_products_andX,
       Manufacture_of_rubber_and_plastic_products,
       Manufacture_of_other_non_metallic_mineral_products,
       Manufacture_of_basic_metals,
       Manufacture_of_fabricated_metal_products_except_machineryX,
       Manufacture_of_computer_electronic_and_optical_products,
       Manufacture_of_electrical_equipment,
       Manufacture_of_machinery_and_equipment_nec,
       Manufacture_of_motor_vehicles_trailers_and_semi_trailers,
       Manufacture_of_other_transport_equipment,
       Manufacture_of_furniture,
       Other_manufacturing,
       Repair_and_installation_of_machinery_and_equipment,
       Electricity_gas_steam_and_air_conditioning_supply,
       Water_collection_treatment_and_supply,
       Sewerage,
       Waste_collection_treatment_and_disposal_activitiesX,
       Remediation_activities_and_other_waste_management_services,
       Construction_of_buildings,
       Civil_engineering,
       Specialised_construction_activities,
       Wholesale_and_retail_trade_and_repair_of_motor_vehiclesX,
       Wholesale_trade_except_of_motor_vehicles_and_motorcycles,
       Retail_trade_except_of_motor_vehicles_and_motorcycles,
       Land_transport_and_transport_via_pipelines,
       Water_transport,
       Air_transport,
       Warehousing_and_support_activities_for_transportation,
       Postal_and_courier_activities,
       Accommodation,
       Food_and_beverage_service_activities,
       Publishing_activities,
       Motion_picture_video_and_television_programme_productionX,
       Programming_and_broadcasting_activities,
       Telecommunications,
       Computer_programming_consultancy_and_related_activities,
       Information_service_activities,
       Financial_service_activities_except_insurance_and_pensionX,
       Insurance_reinsurance_and_pension_funding_exceptX,
       Activities_auxiliary_to_financial_services_and_insuranceX,
       Real_estate_activities,
       Legal_and_accounting_activities,
       Activities_of_head_offices_management_consultancyX,
       Architectural_and_engineering_activities_technical_testingX,
       Scientific_research_and_development,
       Advertising_and_market_research,
       Other_professional_scientific_and_technical_activities,
       Veterinary_activities,
       Rental_and_leasing_activities,
       Employment_activities,
       Travel_agency_tour_operator_and_other_reservation_serviceX,
       Security_and_investigation_activities,
       Services_to_buildings_and_landscape_activities,
       Office_administrative_office_support_and_other_businessX,
       Public_administration_and_defence_compulsory_socialX,
       Education,
       Human_health_activities,
       Residential_care_activities,
       Social_work_activities_without_accommodation,
       Creative_arts_and_entertainment_activities,
       Libraries_archives_museums_and_other_cultural_activities,
       Gambling_and_betting_activities,
       Sports_activities_and_amusement_and_recreation_activities,
       Activities_of_membership_organisations,
       Repair_of_computers_and_personal_and_household_goods,
       Other_personal_service_activities,
       Activities_of_households_as_employers_of_domestic_personnel,
       Undifferentiated_goods_and_services_producing_activitiesX,
       Activities_of_extraterritorial_organisations_and_bodies
export Missing_SIC_2007

@enum SIC_2007 begin  # mapped from sic
   Missing_SIC_2007 = -1
   Undefined_SIC = 0
   Crop_and_animal_production_hunting_and_related_serviceX = 1
   Forestry_and_logging = 2
   Fishing_and_aquaculture = 3
   Mining_of_coal_and_lignite = 5
   Extraction_of_crude_petroleum_and_natural_gas = 6
   Mining_of_metal_ores = 7
   Other_mining_and_quarrying = 8
   Mining_support_service_activities = 9
   Manufacture_of_food_products = 10
   Manufacture_of_beverages = 11
   Manufacture_of_tobacco_products = 12
   Manufacture_of_textiles = 13
   Manufacture_of_wearing_apparel = 14
   Manufacture_of_leather_and_related_products = 15
   Manufacture_of_wood_and_of_products_of_wood_and_corkX = 16
   Manufacture_of_paper_and_paper_products = 17
   Printing_and_reproduction_of_recorded_media = 18
   Manufacture_of_coke_and_refined_petroleum_products = 19
   Manufacture_of_chemicals_and_chemical_products = 20
   Manufacture_of_basic_pharmaceutical_products_andX = 21
   Manufacture_of_rubber_and_plastic_products = 22
   Manufacture_of_other_non_metallic_mineral_products = 23
   Manufacture_of_basic_metals = 24
   Manufacture_of_fabricated_metal_products_except_machineryX = 25
   Manufacture_of_computer_electronic_and_optical_products = 26
   Manufacture_of_electrical_equipment = 27
   Manufacture_of_machinery_and_equipment_nec = 28
   Manufacture_of_motor_vehicles_trailers_and_semi_trailers = 29
   Manufacture_of_other_transport_equipment = 30
   Manufacture_of_furniture = 31
   Other_manufacturing = 32
   Repair_and_installation_of_machinery_and_equipment = 33
   Electricity_gas_steam_and_air_conditioning_supply = 35
   Water_collection_treatment_and_supply = 36
   Sewerage = 37
   Waste_collection_treatment_and_disposal_activitiesX = 38
   Remediation_activities_and_other_waste_management_services = 39
   Construction_of_buildings = 41
   Civil_engineering = 42
   Specialised_construction_activities = 43
   Wholesale_and_retail_trade_and_repair_of_motor_vehiclesX = 45
   Wholesale_trade_except_of_motor_vehicles_and_motorcycles = 46
   Retail_trade_except_of_motor_vehicles_and_motorcycles = 47
   Land_transport_and_transport_via_pipelines = 49
   Water_transport = 50
   Air_transport = 51
   Warehousing_and_support_activities_for_transportation = 52
   Postal_and_courier_activities = 53
   Accommodation = 55
   Food_and_beverage_service_activities = 56
   Publishing_activities = 58
   Motion_picture_video_and_television_programme_productionX = 59
   Programming_and_broadcasting_activities = 60
   Telecommunications = 61
   Computer_programming_consultancy_and_related_activities = 62
   Information_service_activities = 63
   Financial_service_activities_except_insurance_and_pensionX = 64
   Insurance_reinsurance_and_pension_funding_exceptX = 65
   Activities_auxiliary_to_financial_services_and_insuranceX = 66
   Real_estate_activities = 68
   Legal_and_accounting_activities = 69
   Activities_of_head_offices_management_consultancyX = 70
   Architectural_and_engineering_activities_technical_testingX = 71
   Scientific_research_and_development = 72
   Advertising_and_market_research = 73
   Other_professional_scientific_and_technical_activities = 74
   Veterinary_activities = 75
   Rental_and_leasing_activities = 77
   Employment_activities = 78
   Travel_agency_tour_operator_and_other_reservation_serviceX = 79
   Security_and_investigation_activities = 80
   Services_to_buildings_and_landscape_activities = 81
   Office_administrative_office_support_and_other_businessX = 82
   Public_administration_and_defence_compulsory_socialX = 84
   Education = 85
   Human_health_activities = 86
   Residential_care_activities = 87
   Social_work_activities_without_accommodation = 88
   Creative_arts_and_entertainment_activities = 90
   Libraries_archives_museums_and_other_cultural_activities = 91
   Gambling_and_betting_activities = 92
   Sports_activities_and_amusement_and_recreation_activities = 93
   Activities_of_membership_organisations = 94
   Repair_of_computers_and_personal_and_household_goods = 95
   Other_personal_service_activities = 96
   Activities_of_households_as_employers_of_domestic_personnel = 97
   Undifferentiated_goods_and_services_producing_activitiesX = 98
   Activities_of_extraterritorial_organisations_and_bodies = 99
end

export Standard_Occupational_Classification  # mapped from soc2010
export Undefined_SOC,
       Managers_Directors_and_Senior_Officials,
       Professional_Occupations,
       Associate_Prof_and_Technical_Occupations,
       Admin_and_Secretarial_Occupations,
       Skilled_Trades_Occupations,
       Caring_leisure_and_other_service_occupations,
       Sales_and_Customer_Service,
       Process_Plant_and_Machine_Operatives,
       Elementary_Occupations
export Missing_Standard_Occupational_Classification

@enum Standard_Occupational_Classification begin  # mapped from soc2010
   Missing_Standard_Occupational_Classification = -1
   Undefined_SOC = 0
   Managers_Directors_and_Senior_Officials = 1000
   Professional_Occupations = 2000
   Associate_Prof_and_Technical_Occupations = 3000
   Admin_and_Secretarial_Occupations = 4000
   Skilled_Trades_Occupations = 5000
   Caring_leisure_and_other_service_occupations = 6000
   Sales_and_Customer_Service = 7000
   Process_Plant_and_Machine_Operatives = 8000
   Elementary_Occupations = 9000
end

export Tenure_Type  # mapped from tentyp2
export Council_Rented,
       Housing_Association,
       Private_Rented_Unfurnished,
       Private_Rented_Furnished,
       Mortgaged_Or_Shared,
       Owned_outright,
       Rent_free,
       Squats
export Missing_Tenure_Type,is_owner_occupier,renter, is_social_renter, is_private_renter

@enum Tenure_Type begin  # mapped from tentyp2
   Missing_Tenure_Type = -1
   Council_Rented = 1
   Housing_Association = 2
   Private_Rented_Unfurnished = 3
   Private_Rented_Furnished = 4
   Mortgaged_Or_Shared = 5
   Owned_outright = 6
   Rent_free = 7
   Squats = 8
end

function renter(tt::Tenure_Type)::Bool
   tt < Mortgaged_Or_Shared
end

function is_owner_occupier(tt::Tenure_Type)::Bool
   tt in [Mortgaged_Or_Shared, Owned_outright]
end

function is_social_renter( tt::Tenure_Type)::Bool
   tt in [Council_Rented, Housing_Association]
end

function is_private_renter( tt::Tenure_Type)::Bool
   tt in [Private_Rented_Furnished, Private_Rented_Unfurnished]
end

export CT_Band  # mapped from ctband
export Band_A,
       Band_B,
       Band_C,
       Band_D,
       Band_E,
       Band_F,
       Band_G,
       Band_H,
       Band_I,
       Household_not_valued_separately
export Missing_CT_Band

@enum CT_Band begin  # mapped from ctband
   Missing_CT_Band = -1
   Band_A = 1
   Band_B = 2
   Band_C = 3
   Band_D = 4
   Band_E = 5
   Band_F = 6
   Band_G = 7
   Band_H = 8
   Band_I = 9
   Household_not_valued_separately = 10
end

export Disability_Type,
       Disability_Dict,
       vision,
       hearing,
       mobility,
       dexterity,
       learning,
       memory,
       mental_health,
       stamina,
       socially,
       other_difficulty

@enum Disability_Type begin  # mapped from "disd01..10" "cdisd01..10"
   vision = 1
   hearing = 2
   mobility = 3
   dexterity = 4
   learning = 5
   memory = 6
   mental_health = 7
   stamina = 8
   socially = 9
   other_difficulty = 10
end
Disability_Dict = Dict{Disability_Type,T} where T<:Real

export Incomes_Type, Incomes_Dict, Incomes_Set, Expenses

export wages,
       self_employment_income,
       self_employment_expenses,
       self_employment_losses,
       odd_jobs,
       private_pensions,
       national_savings,
       bank_interest,
       stocks_shares,
       individual_savings_account,
       # dividends,
       property,
       royalties,
       bonds_and_gilts,
       other_investment_income,
       other_income,
       alimony_and_child_support_received,
       health_insurance,
       alimony_and_child_support_paid,
       care_insurance,
       trade_unions_etc,
       friendly_societies,
       work_expenses,
       avcs,
       other_deductions,
       loan_repayments,
       student_loan_repayments,
       pension_contributions_employee,
       pension_contributions_employer,
       education_allowances,
       foster_care_payments,
       student_grants,
       student_loans,
       income_tax,
       national_insurance,
       local_taxes,
       free_school_meals,
       dlaself_care,
       dlamobility,
       child_benefit,
       pension_credit,
       state_pension,
       bereavement_allowance_or_widowed_parents_allowance_or_bereavement,
       armed_forces_compensation_scheme,
       war_widows_or_widowers_pension,
       severe_disability_allowance,
       attendance_allowance, ## FIXME SP!
       carers_allowance,
       jobseekers_allowance,
       industrial_injury_disablement_benefit,
       employment_and_support_allowance,
       incapacity_benefit,
       income_support,
       maternity_allowance,
       maternity_grant_from_social_fund,
       funeral_grant_from_social_fund,
       any_other_ni_or_state_benefit,
       trade_union_sick_or_strike_pay,
       friendly_society_benefits,
       private_sickness_scheme_benefits,
       accident_insurance_scheme_benefits,
       hospital_savings_scheme_benefits,
       government_training_allowances,
       guardians_allowance,
       widows_payment,
       unemployment_or_redundancy_insurance,
       winter_fuel_payments,
       dwp_third_party_payments_is_or_pc,
       dwp_third_party_payments_jsa_or_esa,
       social_fund_loan_repayment_from_is_or_pc,
       social_fund_loan_repayment_from_jsa_or_esa,
       extended_hb,
       permanent_health_insurance,
       any_other_sickness_insurance,
       critical_illness_cover,
       working_tax_credit,
       child_tax_credit,
       working_tax_credit_lump_sum,
       child_tax_credit_lump_sum,
       housing_benefit,
       universal_credit,
       personal_independence_payment_daily_living,
       personal_independence_payment_mobility,
       other_benefits

@enum Incomes_Type begin
   wages = 1
   self_employment_income = 2
   self_employment_expenses = 3
   self_employment_losses = 4
   odd_jobs = 5
   private_pensions = 6
   national_savings = 7
   bank_interest = 8
   stocks_shares = 9
   individual_savings_account = 10
   # dividends = 11 ### FIXME NOT USED NEEDS DELETED. Use stocks_shares instead
   property = 12
   royalties = 13
   bonds_and_gilts = 14
   other_investment_income = 15
   other_income = 16
   alimony_and_child_support_received = 17
   health_insurance = 18
   alimony_and_child_support_paid = 19
   care_insurance = 20
   trade_unions_etc = 21
   friendly_societies = 22
   work_expenses = 23
   avcs = 24
   other_deductions = 25
   loan_repayments = 26
   student_loan_repayments = 27
   pension_contributions_employee = 28
   pension_contributions_employer = 29
   education_allowances = 30
   foster_care_payments = 31
   student_grants = 32
   student_loans = 33
   income_tax = 34
   national_insurance = 35
   local_taxes = 36
   free_school_meals = 37

   dlaself_care = 2001
   dlamobility = 2002
   child_benefit = 2003
   pension_credit = 2004
   state_pension = 2005
   bereavement_allowance_or_widowed_parents_allowance_or_bereavement = 2006
   armed_forces_compensation_scheme = 2008
   war_widows_or_widowers_pension = 2009
   severe_disability_allowance = 2010
   attendance_allowance = 2012
   carers_allowance = 2013
   jobseekers_allowance = 2014
   industrial_injury_disablement_benefit = 2015
   employment_and_support_allowance = 2016
   incapacity_benefit = 2017
   income_support = 2019
   maternity_allowance = 2021
   maternity_grant_from_social_fund = 2022
   funeral_grant_from_social_fund = 2024
   any_other_ni_or_state_benefit = 2030
   trade_union_sick_or_strike_pay = 2031
   friendly_society_benefits = 2032
   private_sickness_scheme_benefits = 2033
   accident_insurance_scheme_benefits = 2034
   hospital_savings_scheme_benefits = 2035
   government_training_allowances = 2036
   guardians_allowance = 2037
   widows_payment = 2060
   unemployment_or_redundancy_insurance = 2061
   winter_fuel_payments = 2062
   dwp_third_party_payments_is_or_pc = 2065
   dwp_third_party_payments_jsa_or_esa = 2066
   social_fund_loan_repayment_from_is_or_pc = 2069
   social_fund_loan_repayment_from_jsa_or_esa = 2070
   extended_hb = 2078
   permanent_health_insurance = 2081
   any_other_sickness_insurance = 2082
   critical_illness_cover = 2083
   working_tax_credit = 2090
   child_tax_credit = 2091
   working_tax_credit_lump_sum = 2092
   child_tax_credit_lump_sum = 2093
   housing_benefit = 2094
   universal_credit = 2095
   personal_independence_payment_daily_living = 2096
   personal_independence_payment_mobility = 2097

   other_benefits = 3000
end

Incomes_Dict = Dict{Incomes_Type,T} where T<:Real
Incomes_Set = Set{Incomes_Type}

const Expenses = Incomes_Set([
   permanent_health_insurance,
   health_insurance,
   alimony_and_child_support_paid,
   care_insurance,
   trade_unions_etc,
   friendly_societies,
   self_employment_expenses,
   work_expenses,
   pension_contributions_employee,
   pension_contributions_employer, # kinda sorta
   avcs,
   other_deductions,
   loan_repayments,
   student_loan_repayments
])

export Asset_Type, Asset_Dict, Asset_Set

export A_Current_account,
       A_NSB_Ordinary_account,
       A_NSB_Investment_account,
       A_Not_Used,
       A_Savings_investments_etc,
       A_Government_Gilt_Edged_Stock,
       A_Unit_or_Investment_Trusts,
       A_Stocks_Shares_Bonds_etc,
       A_PEP,
       A_National_Savings_capital_bonds,
       A_Index_Linked_National_Savings_Certificates,
       A_Fixed_Interest_National_Savings_Certificates,
       A_Pensioners_Guaranteed_Bonds,
       A_SAYE,
       A_Premium_bonds,
       A_National_Savings_income_bonds,
       A_National_Savings_deposit_bonds,
       A_First_Option_bonds,
       A_Yearly_Plan,
       A_ISA,
       A_Fixd_Rate_Svngs_Bonds_or_Grntd_Incm_Bonds_or_Grntd_Growth_Bonds,
       A_GEB,
       A_Basic_Account,
       A_Credit_Unions,
       A_Endowment_Policy_Not_Linked
export Missing_Asset_Type

@enum Asset_Type begin  # mapped from assetype
   Missing_Asset_Type = -1
   A_Current_account = 1
   A_NSB_Ordinary_account = 2
   A_NSB_Investment_account = 3
   A_Not_Used = 4
   A_Savings_investments_etc = 5
   A_Government_Gilt_Edged_Stock = 6
   A_Unit_or_Investment_Trusts = 7
   A_Stocks_Shares_Bonds_etc = 8
   A_PEP = 9
   A_National_Savings_capital_bonds = 10
   A_Index_Linked_National_Savings_Certificates = 11
   A_Fixed_Interest_National_Savings_Certificates = 12
   A_Pensioners_Guaranteed_Bonds = 13
   A_SAYE = 14
   A_Premium_bonds = 15
   A_National_Savings_income_bonds = 16
   A_National_Savings_deposit_bonds = 17
   A_First_Option_bonds = 18
   A_Yearly_Plan = 19
   A_ISA = 21
   A_Fixd_Rate_Svngs_Bonds_or_Grntd_Incm_Bonds_or_Grntd_Growth_Bonds = 25
   A_GEB = 26
   A_Basic_Account = 27
   A_Credit_Unions = 28
   A_Endowment_Policy_Not_Linked = 29
end

@enum Included_In_Pay_Type begin  # mapped from assetype
   ssp = 1
   smp = 2
   mileage = 4
   motoring_expenses = 5
   spp = 6
   sap = 7
end

Included_In_Pay_Dict = Dict{Included_In_Pay_Type,T} where T<:Real

export Included_In_Pay_Type,
   ssp,
   smp,
   mileage,
   motoring_expenses,
   spp,
   sap,
   Included_In_Pay_Dict

Asset_Dict = Dict{Asset_Type,T} where T<:Real
Asset_Set = Set{Asset_Type}

export Relationship  # mapped from relhrp
export
       This_Person,
       Spouse,
       Cohabitee,
       Son_or_daughter_incl_adopted,
       Step_son_or_daughter,
       Foster_child,
       Son_in_law_or_daughter_in_law,
       Parent,
       Step_parent,
       Foster_parent,
       Parent_in_law,
       Brother_or_sister_incl_adopted,
       Step_brother_or_sister,
       Foster_brother_or_sister,
       Brother_or_sister_in_law,
       Grand_child,
       Grand_parent,
       Other_relative,
       Other_non_relative,
       Civil_Partner,
       Missing_Relationship

@enum Relationship begin  # mapped from relhrp
   Missing_Relationship = -1
   This_Person = 0 # GKS added
   Spouse = 1
   Cohabitee = 2
   Son_or_daughter_incl_adopted = 3
   Step_son_or_daughter = 4
   Foster_child = 5
   Son_in_law_or_daughter_in_law = 6
   Parent = 7
   Step_parent = 8
   Foster_parent = 9
   Parent_in_law = 10
   Brother_or_sister_incl_adopted = 11
   Step_brother_or_sister = 12
   Foster_brother_or_sister = 13
   Brother_or_sister_in_law = 14
   Grand_child = 15
   Grand_parent = 16
   Other_relative = 17
   Other_non_relative = 18
   Civil_Partner = 20
end

export Relationship_Dict
Relationship_Dict = Dict{BigInt,Relationship}

export Employment_Type  # mapped from etype
export An_Employee,
       Running_a_business_or_prof_practice,
       Partner_in_a_business_or_practice,
       Working_for_myself,
       Sub_Contractor,
       Doing_freelance_work,
       Self_employed_in_some_other_way
export Missing_Employment_Type

@enum Employment_Type begin  # mapped from etype
   Missing_Employment_Type = -1
   An_Employee = 1
   Running_a_business_or_prof_practice = 2
   Partner_in_a_business_or_practice = 3
   Working_for_myself = 4
   Sub_Contractor = 5
   Doing_freelance_work = 6
   Self_employed_in_some_other_way = 7
end

export Employment_Sector  # mapped from jobsect
export Private, Public
export Missing_Employment_Sector

@enum Employment_Sector begin  # mapped from jobsect
   Missing_Employment_Sector = -1
   Private = 1
   Public = 2
end

export Account_Type  # mapped from account
export Current_account,
       NSB_Direct_Saver,
       NSB_Investment_account,
       Not_Used,
       Savings_investments_etc,
       Government_Gilt_Edged_Stock,
       Unit_or_Investment_Trusts,
       Stocks_Shares_Bonds_etc,
       PEP,
       National_Savings_capital_bonds,
       Index_Linked_National_Savings_Certificates,
       Fixed_Interest_National_Savings_Certificates,
       Pensioners_Guaranteed_Bonds,
       SAYE,
       Premium_bonds,
       National_Savings_income_bonds,
       National_Savings_deposit_bonds,
       First_Option_bonds,
       Yearly_Plan,
       ISA,
       Profit_sharing,
       Company_Share_Option_Plans,
       Member_of_Share_Club,
       Fixed_Rate_Savings_or_Guaranteed_Income_or_Guaranteed_Growth_Bonds,
       Guaranteed_Equity_Bond,
       Basic_Account,
       Credit_Unions,
       Endowment_Policy_Not_Linked,
       Post_Office_Card_Account
export Missing_Account_Type

@enum Account_Type begin  # mapped from account
   Missing_Account_Type = -1
   Current_account = 1
   NSB_Direct_Saver = 2
   NSB_Investment_account = 3
   Not_Used = 4
   Savings_investments_etc = 5
   Government_Gilt_Edged_Stock = 6
   Unit_or_Investment_Trusts = 7
   Stocks_Shares_Bonds_etc = 8
   PEP = 9
   National_Savings_capital_bonds = 10
   Index_Linked_National_Savings_Certificates = 11
   Fixed_Interest_National_Savings_Certificates = 12
   Pensioners_Guaranteed_Bonds = 13
   SAYE = 14
   Premium_bonds = 15
   National_Savings_income_bonds = 16
   National_Savings_deposit_bonds = 17
   First_Option_bonds = 18
   Yearly_Plan = 19
   ISA = 21
   Profit_sharing = 22
   Company_Share_Option_Plans = 23
   Member_of_Share_Club = 24
   Fixed_Rate_Savings_or_Guaranteed_Income_or_Guaranteed_Growth_Bonds = 25
   Guaranteed_Equity_Bond = 26# Guaranteed/Fixed rate E bond
   Basic_Account = 27
   Credit_Unions = 28
   Endowment_Policy_Not_Linked = 29
   Post_Office_Card_Account = 30
end

export Account_Tax_Status  # mapped from invtax
export After_tax, Before_tax
export Missing_Account_Tax_Status

@enum Account_Tax_Status begin  # mapped from invtax
   Missing_Account_Tax_Status = -1
   After_tax = 1
   Before_tax = 2
end

export Benefit_Type  # mapped from benefit
export DLAself_care,
       DLAmobility,
       Child_Benefit,
       Pension_credit,
       State_Pension,
       Bereavement_Allowance_or_Widowed_Parents_Allowance_or_Bereavement,
       Armed_Forces_Compensation_Scheme,
       War_Widows_or_Widowers_Pension,
       Severe_Disability_Allowance,
       Attendance_Allowance,
       Carers_Allowance,
       Jobseekers_Allowance,
       Industrial_Injury_Disablement_Benefit,
       Employment_and_Support_Allowance,
       Incapacity_Benefit,
       Income_Support,
       Maternity_Allowance,
       Maternity_Grant_from_Social_Fund,
       Funeral_Grant_from_Social_Fund,
       Any_other_NI_or_State_benefit,
       Trade_Union_sick_or_strike_pay,
       Friendly_Society_Benefits,
       Private_sickness_scheme_benefits,
       Accident_insurance_scheme_benefits,
       Hospital_savings_scheme_benefits,
       Government_training_allowances,
       Guardians_Allowance,
       Widows_Payment,
       Unemployment_or_Redundancy_Insurance,
       Winter_Fuel_Payments,
       DWP_third_party_payments_IS_or_PC,
       DWP_third_party_payments_JSA_or_ESA,
       Social_Fund_Loan_Repayment_from_IS_or_PC,
       Social_Fund_Loan_Repayment_from_JSA_or_ESA,
       Extended_HB,
       Permanent_health_insurance,
       Any_other_sickness_insurance,
       Critical_Illness_Cover,
       Working_Tax_Credit,
       Child_Tax_Credit,
       Working_Tax_Credit_Lump_Sum,
       Child_Tax_Credit_Lump_Sum,
       Housing_Benefit,
       Universal_Credit,
       Personal_Independence_Payment_Daily_Living,
       Personal_Independence_Payment_Mobility,
       A_loan_from_the_DWP_and_DfC,
       A_loan_or_grant_from_local_authority,
       Future_Pension_Credit,
       Future_Universal_Credit,
       Future_Housing_Benefit,
       Future_Working_Tax_Credit,
       Future_Child_Tax_Credit,
       Future_Income_Support,
       Future_Jobseekers_Allowance,
       Future_Employment_and_Support_Allowance,
       DWP_third_party_payments_UC,
       Social_Fund_Loan_UC
export Missing_Benefit_Type

@enum Benefit_Type begin  # mapped from benefit
   Missing_Benefit_Type = -1
   DLAself_care = 1
   DLAmobility = 2
   Child_Benefit = 3
   Pension_credit = 4
   State_Pension = 5
   Bereavement_Allowance_or_Widowed_Parents_Allowance_or_Bereavement = 6
   Armed_Forces_Compensation_Scheme = 8
   War_Widows_or_Widowers_Pension = 9
   Severe_Disability_Allowance = 10
   Attendance_Allowance = 12
   Carers_Allowance = 13
   Jobseekers_Allowance = 14
   Industrial_Injury_Disablement_Benefit = 15
   Employment_and_Support_Allowance = 16
   Incapacity_Benefit = 17
   Income_Support = 19
   Maternity_Allowance = 21
   Maternity_Grant_from_Social_Fund = 22
   Funeral_Grant_from_Social_Fund = 24
   Any_other_NI_or_State_benefit = 30
   Trade_Union_sick_or_strike_pay = 31
   Friendly_Society_Benefits = 32
   Private_sickness_scheme_benefits = 33
   Accident_insurance_scheme_benefits = 34
   Hospital_savings_scheme_benefits = 35
   Government_training_allowances = 36
   Guardians_Allowance = 37
   Widows_Payment = 60
   Unemployment_or_Redundancy_Insurance = 61
   Winter_Fuel_Payments = 62
   DWP_third_party_payments_IS_or_PC = 65
   DWP_third_party_payments_JSA_or_ESA = 66
   Social_Fund_Loan_Repayment_from_IS_or_PC = 69
   Social_Fund_Loan_Repayment_from_JSA_or_ESA = 70
   Extended_HB = 78
   Permanent_health_insurance = 81
   Any_other_sickness_insurance = 82
   Critical_Illness_Cover = 83
   Working_Tax_Credit = 90
   Child_Tax_Credit = 91
   Working_Tax_Credit_Lump_Sum = 92
   Child_Tax_Credit_Lump_Sum = 93
   Housing_Benefit = 94
   Universal_Credit = 95
   Personal_Independence_Payment_Daily_Living = 96
   Personal_Independence_Payment_Mobility = 97
   A_loan_from_the_DWP_and_DfC = 98
   A_loan_or_grant_from_local_authority = 99
   Future_Pension_Credit = 102
   Future_Universal_Credit = 103
   Future_Housing_Benefit = 104
   Future_Working_Tax_Credit = 105
   Future_Child_Tax_Credit = 106
   Future_Income_Support = 107
   Future_Jobseekers_Allowance = 108
   Future_Employment_and_Support_Allowance = 109
   DWP_third_party_payments_UC = 110
   Social_Fund_Loan_UC = 111
end

export Child_Care_Type  # mapped from chlook
export Playgroup_or_pre_school,
       Day_nursery_or_crèche,
       Nursery_school,
       Nursery_class_attached_to_Primary_or_Infants_School,
       Reception_class_at_Primary_or_Infants_School,
       Breakfast_club,
       After_school_club_or_activities,
       Holiday_scheme_or_club,
       Special_day_school_or_nursery_unit_for_children_with_special_ed,
       Other_formal,
       Grand_parents,
       Non_resident_parent_or_ex_spouse_or_ex_partner,
       Childd_brother_or_sister,
       Other_relatives,
       Childminder,
       Nanny_or_Au_Pair,
       Friends_or_neighbours,
       Other_non_relatives_includes_babysitters
export Missing_Child_Care_Type

@enum Child_Care_Type begin  # mapped from chlook
   Missing_Child_Care_Type = -1
   Playgroup_or_pre_school = 1
   Day_nursery_or_crèche = 2
   Nursery_school = 3
   Nursery_class_attached_to_Primary_or_Infants_School = 4
   Reception_class_at_Primary_or_Infants_School = 5
   Breakfast_club = 6
   After_school_club_or_activities = 7
   Holiday_scheme_or_club = 8
   Special_day_school_or_nursery_unit_for_children_with_special_ed = 9
   Other_formal = 10
   Grand_parents = 11
   Non_resident_parent_or_ex_spouse_or_ex_partner = 12
   Childs_brother_or_sister = 13
   Other_relatives = 14
   Childminder = 15
   Nanny_or_Au_Pair = 16
   Friends_or_neighbours = 17
   Other_non_relatives_includes_babysitters = 18
end

export Fuel_Type  # mapped from fueltyp
export Petrol, Diesel, Biofuel_eg_E85_fuel, Hybrid_use_a_combination_of_petrol_and_electricity, Electric, LPG, Other, Dont_know
export No_Fuel, Missing_Fuel_Type, Fuel_Type_Dict

@enum Fuel_Type begin  # mapped from fueltyp
   Missing_Fuel_Type = -1
   No_Fuel = 0
   Petrol = 1
   Diesel = 2
   Biofuel_eg_E85_fuel = 3
   Hybrid_use_a_combination_of_petrol_and_electricity = 4
   Electric = 5
   LPG = 6
   Other = 7
   Dont_know = 8
end

const Fuel_Type_Dict = Dict{Fuel_Type,T} where T <: Number

export
   Uprate_Item_Type,
   upr_no_uprate,
   upr_earnings,
   upr_housing_rents,
   upr_housing_oo,
   upr_unearned,
   upr_costs,
   upr_cpi,
   upr_gdp_deflator,
   upr_nominal_gdp,
   upr_shares

@enum Uprate_Item_Type begin
    upr_no_uprate
    upr_earnings
    upr_housing_rents
    upr_housing_oo
    upr_unearned
    upr_costs
    upr_cpi
    upr_gdp_deflator
    upr_nominal_gdp
    upr_shares
end

@enum JSAType begin
    no_jsa = -1
    contributory_jsa = 1
    income_related_jsa =2
    both_jsa = 3
end
# FIXME we also use this for ESA, which makes the names confusing
export JSAType, contributory_jsa, income_related_jsa, both_jsa, no_jsa

@enum LMTBenefitType esa hb ctr is jsa pc wtc ctc sc
const LMTBenefitSet = Set{LMTBenefitType}
export LMTBenefitType, esa, hb, ctr, is, jsa, pc, wtc, ctc, sc, LMTBenefitSet

@enum LMTPremia  begin
   family_premium = 1
   disabled_child = 2
   severe_disability_single = 3
   severe_disability_couple = 4
   carer_single = 5
   carer_couple = 6
   enhanced_disability_child = 7
   enhanced_disability_single = 8
   enhanced_disability_couple = 9
   disability_single = 10
   disability_couple = 11
   pensioner_is = 12
end

const LMTPremiaDict = Dict{LMTPremia,T} where T <:Real
const LMTPremiaSet = Set{LMTPremia}

export LMTPremia,LMTPremiaDict, LMTPremiaSet, disabled_child,
    severe_disability_single,severe_disability_couple,
    carer_single,carer_couple,enhanced_disability_child,
    enhanced_disability_single,enhanced_disability_couple,disability_single,disability_couple,pensioner_is

export BIG_NOTHING
const BIG_NOTHING = BigInt(0)

export DataSource, FRS, OtherSource, ExampleSource

@enum DataSource begin  # mapped from relhrp
   FRS = 1
   OtherSource = 2
   ExampleSource = 3
end

export get_pid, safe_assign, safe_inc, from_pid
"""
  get a unique ID for a person from (say) a certain year of the FRS, given hhld number and number inside the household
"""
function get_pid(source::DataSource, year::Integer, hid::Integer, pno::Integer)::BigInt
   (Int(source) * 10^11) + (year * 10^7) + (hid * 10^2) + pno
end

function from_pid( pid :: Integer ) :: NamedTuple
   ids = extract_digits( pid, 1:1 )
   source = DataSource(ids)
   year = extract_digits(pid, 2:5)
   hid = extract_digits( pid, 6:10)
   pno = extract_digits( pid, 11:12)
   return (datasource=source, year=year, hid=hid, pno=pno )
end

export DEFAULT_MISSING_VALUES, safe_inc

const DEFAULT_MISSING_VALUES = [-9.0, -8.0, -7.0, -6.0, -5.0, -4.0, -3.0, -2.0, -1.0]

"""
a = a+b if b is not missing or (-9 to -1), else a
"""
function safe_inc(a::Real, b::Union{Real,Missing})::Real
   if ismissing(b) || (b in DEFAULT_MISSING_VALUES)
      return a
   end
   a + b
end

"""
map missing values or (-9 to -2) to -1, else a
used for mapping values to enums, where we add a missing enum with value -1
"""
function safe_assign(a::Union{Number,Missing,AbstractString},deflt=-1)
   if ismissing(a)
      return deflt
   elseif typeof(a) <: AbstractString
      a = tryparse( typeof( deflt ), a )
      if a === nothing
         return deflt
      end
   end
   if a < -1
      return deflt
   end
   return a
end

export make_sym_for_frame, make_sym_for_asset, make_sym_from_frame
"""
"income", :fred => :income_fred
"""
function make_sym_for_frame(prefix::AbstractString, enum::Enum)::Symbol
    sym = Symbol(enum)
    Symbol(lowercase(prefix * "_" * String(Symbol(sym))))
end

"""
"income", :a_fred => :income_fred
"""
function make_sym_for_asset(enum::Enum)::Symbol
    s = String(Symbol(enum))[3:end]
   # println( "s=$s")
    Symbol(lowercase("asset_" * s))
end

"""
"income", :income_fred => :fred
"""
function make_sym_from_frame(prefix::AbstractString, sym::Symbol)::Symbol
   # FIXME got to be a simpler way
    matchstr = "$(prefix)(.*)"
    re = Regex(matchstr)
    rm = match(re, String(sym))
    Symbol(rm[1])
end

export safe_to_bool, m2z

function safe_to_bool( thing ) :: Bool
   ty = typeof( thing )
   if ty == Bool
      return thing
   end
   if ismissing( thing )
      return false
   end
   if ty <: Number
      return thing == 1
   end
   return false
end

function m2z( s )
   if ismissing( s )
      return 0
   end
   s
end

export
   age_range,
   age_ends,
   age_str

#
# FIXME this is a dup of the thing in `weights`
#
const age_ends = [4,9,14,19,24,29,34,39,44,49,54,59,64,69,74,79,99999]

function age_str( band )
   a2 = "$(age_ends[band])"
   a1 = ""
   if band==1
      a1 = "0"
   elseif band==lastindex(age_ends)
      return "$(age_ends[band-1]+1) and above"
   else
      a1 = "$(age_ends[band-1]+1)"
   end
   return "$(a1)-$(a2)"
end

function age_range( age :: Integer ) :: Integer
   for a in eachindex(age_ends)
      if age <= age_ends[a]
         return a
      end
   end
end

#
# This is for quick lookups of hhld/people by hid/pid
# you need the `isequal`/`hash` to make sure `Dict`s can use these
# as keys.

export OneIndex, in_indexes
struct OneIndex
   id :: BigInt
   data_year :: Int
end

function in_indexes( s :: Set{OneIndex}, id :: BigInt, data_year:: Int)
   return OneIndex( id, data_year ) in s
end

function in_indexes( s :: Set{OneIndex}, id :: BigInt )
   # this is silly ...
   bits = from_pid( id )
   return in_indexes( s, id, bits.year )
end

Base.isequal( a :: OneIndex, b :: OneIndex ) = (a.id==b.id)&&(a.data_year == b.data_year )

Base.hash( a :: OneIndex ) = UInt64(a.id*1000)+UInt64(a.data_year)

function Base.isless( a :: OneIndex, b :: OneIndex )
   if b.data_year > a.data_year
       return true
   elseif b.data_year < a.data_year
       return false
   else
       return b.id > a.id
   end
end

export UBIMTBenTreatment,
   ub_as_is,
   ub_abolish,
   ub_keep_housing

@enum UBIMTBenTreatment ub_as_is ub_abolish ub_keep_housing

end # module

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using PyCall
using HDF5

@pyimport keras.models as models
@pyimport keras.layers.core as core
@pyimport keras.layers.convolutional as kconv
@pyimport keras.callbacks as kcallbacks
@pyimport keras.utils.io_utils as k_io

"""
Temporary utility function until next release of PyCall (as of 1.4.0).
  Adds `revdims` flag
"""
function to_python{T<:PyCall.NPY_TYPES}(a::StridedArray{T}, revdims::Bool=false)
    @PyCall.npyinitialize
    size_a = revdims ? reverse(size(a)) : size(a)
    strides_a = revdims ? reverse(strides(a)) : strides(a)
    p = ccall(PyCall.npy_api[:PyArray_New], PyPtr,
              (PyPtr,Cint,Ptr{Int},Cint, Ptr{Int},Ptr{T}, Cint,Cint,PyPtr),
              PyCall.npy_api[:PyArray_Type],
              ndims(a), Int[size_a...], PyCall.npy_type(T),
              Int[strides_a...] * sizeof(eltype(a)), a, sizeof(eltype(a)),
              PyCall.NPY_ARRAY_ALIGNED | PyCall.NPY_ARRAY_WRITEABLE,
              C_NULL)
    return PyObject(p, a)
end

"""
Convert array to row major, C style arrays.
Does not modify memory during conversion
"""
function to_python_array(arr::AbstractArray)
    to_python_array(arr, size(arr)...)
end

"""
Convert array to row major, C style arrays.
Does not modify buffer during conversion

Specifically for bit arrays -> UInt8

Specify dimensions IN JULIA COORDINATES
Given Julia array with dims: (1, 2, 3)
Return Python array with dims: (3, 2, 1)
"""
function to_python_array(arr::AbstractArray, dims...)
    reshaped = reshape(arr, dims)
    rounded = round(UInt8, reshaped)
    to_python(rounded, true)
end

"Simple softmax classifier"
function CLF_LINEAR(features::Vector{Function})
    input_shape = get_input_size(features)
    KerasNetwork(models.Sequential([
                                    core.Flatten(input_shape=input_shape),
                                    core.Dense(N*N, input_dim=(N*N)),
                                    core.Activation("softmax")
                                    ]),
                 features)
end

"""
    Network from `Teaching Deep Convolutional Neural Networks to Play Go`
    Designed for 19x19, but should be fine for other sizes as well
    """
function CLF_DCNN(features::Vector{Function})
    # Reverse because julia <--> python
    input_shape = reverse(get_input_size(features))

    KerasNetwork(models.Sequential([
                                    kconv.Convolution2D(64, 7, 7, activation="relu", border_mode="same", input_shape=input_shape),
                                    kconv.Convolution2D(64, 5, 5, activation="relu", border_mode="same", input_shape=input_shape),
                                    kconv.Convolution2D(64, 5, 5, activation="relu", border_mode="same", input_shape=input_shape),
                                    kconv.Convolution2D(48, 5, 5, activation="relu", border_mode="same", input_shape=input_shape),
                                    kconv.Convolution2D(48, 5, 5, activation="relu", border_mode="same", input_shape=input_shape),
                                    kconv.Convolution2D(32, 5, 5, activation="relu", border_mode="same", input_shape=input_shape),
                                    kconv.Convolution2D(32, 5, 5, activation="relu", border_mode="same", input_shape=input_shape),
                                    core.Flatten(),
                                    core.Dense(N*N, activation="softmax")
                                    ]),
                 features)
end

function compile(network::KerasNetwork)
    println(STDERR, "Compiling model...")
    network.model.compile(loss="categorical_crossentropy",
                          optimizer="adadelta",
                          metrics=["accuracy"])
end

function _callbacks(hf5path=nothing)
    callbacks = Vector()
    if hf5path != nothing
        checkpointer = kcallbacks.ModelCheckpoint(filepath=hf5path, verbose=1, save_best_only=true)
        push!(callbacks, checkpointer)
    end
    push!(callbacks, kcallbacks.ProgbarLogger())
    callbacks
end

# TODO: More advanced buffering and shuffling across multiple files
# Doesn't have a fixed batch size right now
@pydef type SGFGenerator
    __init__(self, filenames, features) = begin
        self[:filenames] = Vector{AbstractString}(filenames)
        self[:features] = features
        self[:buffer] = Vector()
        self[:fi] = 1
        self[:bi] = 1
        self[:tm] = 0.0
    end
    __next__(self) = self[:next]()
    next(self) = begin
        features = Vector{Function}(self[:features])
        # Load in the next file and return as one batch
        examples = Vector()
        while length(examples) == 0
            filename = AbstractString(self[:filenames][self[:fi]])
            examples = generate_training_data(filename, features=features)
            # Loop back to the beginning if we run out of files
            self[:fi] = (self[:fi] % length(self[:filenames])) + 1
        end
        X = go.batch_training_examples([x[1] for x in examples])
        Y = go.batch_training_examples([x[2][:] for x in examples])
        return (to_python_array(X),
                to_python_array(Y))
    end
end

# TODO: Add an in memory buffer to allow better shuffling
@pydef type HF5Generator
    __init__(self, hf5path, features, batch_size, start, stop) = begin
        h5 = h5open(AbstractString(hf5path))
        self[:data] = h5["data"]
        self[:label] = h5["label"]
        self[:features] = features
        self[:batch_size] = batch_size
        self[:start] = start
        self[:stop] = stop
        self[:i] = 0
    end
    __next__(self) = self[:next]()
    next(self) = begin
        features = Vector{Function}(self[:features])

        batch_size = self[:batch_size]
        start = self[:start]
        stop = self[:stop]
        i = self[:i]
        if start + i * batch_size > stop - batch_size
            i = 0
        end
        self[:i] = i + 1

        bstart = start + i * batch_size + 1
        bstop = bstart + batch_size - 1

        X = self[:data][:,:,:, bstart:bstop]
        Y = self[:label][:, bstart:bstop]

        return (to_python_array(X),
                to_python_array(Y))
    end
end

"""
Train a network that loads data from
"""
function keras_train_generator(network::KerasNetwork, generator; epochs=20, samples_per_epoch=1000000, recompile=true, validation_data=nothing, checkpoint_path=nothing)
    if recompile
        compile(network)
    end

    println(STDERR, "Fitting model...")
    network.model.fit_generator(generator, samples_per_epoch=samples_per_epoch, nb_epoch=epochs,
                                verbose=2,
                                validation_data=validation_data,
                                callbacks=_callbacks(checkpoint_path))
end

" Train from an HDF5 file "
function keras_train_h5(network::KerasNetwork, hdf5_path::AbstractString; epochs=20, batch_size=32, recompile=true, n_validation=10000, n_train=nothing, checkpoint_path=nothing)
    if recompile
        compile(network)
    end

    #Open the file and see how many examples there are
    h5 = h5open(hdf5_path)
    X = h5["data"]
    Y = h5["label"]
    n_examples = size(X)[end]
    close(h5)

    max_n_train = n_examples - n_validation
    n_train = min(n_train, max_n_train)

    # These act like matrices
    trainX = k_io.HDF5Matrix(hdf5_path, "data", 0, n_train)
    trainY = k_io.HDF5Matrix(hdf5_path, "label", 0, n_train)

    valX = k_io.HDF5Matrix(hdf5_path, "data", n_train, n_train + n_validation)
    valY = k_io.HDF5Matrix(hdf5_path, "label", n_train, n_train + n_validation)

    network.model.fit(trainX, trainY, nb_epoch=epochs, batch_size=batch_size, verbose=2,
                      validation_data=(valX, valY),
                      shuffle="batch", # Don't want to shuffle data too large for memory
                      callbacks=_callbacks(checkpoint_path))
end

function train_model(network::KerasNetwork, X, Y; epochs=20, batch_size=32, recompile=true, validation_split=0.25, checkpoint_path=nothing)
    if recompile
        compile(network)
    end

    # Handle the julia <--> python array conversion
    println(STDERR, "Doing python <--> Julia conversion...")
    X = to_python_array(X)
    Y = to_python_array(Y)
    println(STDERR, "Fitting model...")
    network.model.fit(X, Y, nb_epoch=epochs, batch_size=batch_size, verbose=2,
                      validation_split=validation_split,
                      callbacks=_callbacks(checkpoint_path))
end

function save_model(network::KerasNetwork, folder::AbstractString, name::AbstractString; save_yaml=true, save_weights=true, overwrite=false)
    if save_weights
        hf5path = joinpath(folder, string(name, ".hf5"))
        isfile(hf5path) && !overwrite && (println(STDERR, "File exists: $(hf5path)"); return)
        network.model.save_weights(hf5path)
    end

    if save_yaml
        ymlpath = joinpath(folder, string(name, ".yml"))
        isfile(ymlpath) && !overwrite && (println(STDERR, "File exists: $(ymlpath)"); return)
        yaml = network.model.to_yaml()

        # Add feature names
        # e.g. -> player_colors, liberties, ...
        # the split is to remove "go." - caused problems when reading in
        features = join(map(x -> split(string(x), ".")[end], network.features), ", ")
        yaml = string(yaml, "\nfeatures: [$(features)]")

        open(joinpath(folder, string(name, ".yml")), "w") do file
            write(file, yaml)
        end
    end
end

ext: jl


text:
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using BSeries
using Latexify
using OrdinaryDiffEq
using RootedTrees
using StaticArrays
using Symbolics
import SymPy; sp=SymPy;
using PyCall, LaTeXStrings; import PyPlot; plt=PyPlot

inset_locator = pyimport("mpl_toolkits.axes_grid.inset_locator")

cycler = pyimport("cycler").cycler
line_cycler   = (cycler(color=["#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442","#E69F00"]) +
                 cycler(linestyle=["-", "--", "-.", ":", "-", "--", "-."]))
marker_cycler = (cycler(color=["#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442","#E69F00"]) +
                 cycler(linestyle=["none", "none", "none", "none", "none", "none", "none"]) +
                 cycler(marker=["4", "2", "3", "1", "+", "x", "."]))

plt.rc("axes", prop_cycle=line_cycler)
plt.rc("font", family="serif", size=18.)
plt.rc("savefig", dpi=100)
plt.rc("legend", loc="best", fontsize="medium", fancybox=true, framealpha=0.5)
plt.rc("lines", linewidth=2.5, markersize=5, markeredgewidth=2.5)

max_stars_repo_name: ketch/bseries-package-notebooks

max_issues_repo_name: ketch/bseries-package-notebooks

max_forks_repo_name: ketch/bseries-package-notebooks


text:
"""
Compute spectral geopotential `ϕ` from spectral temperature `Tabs`
and spectral topography `ϕ0`.
"""
function geopotential!( ϕ::Array{Complex{NF},3},      # geopotential
                        ϕ0::Array{Complex{NF},2},     # surface geopotential
                        Tabs::Array{Complex{NF},3},   # absolute Temperature
                        G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    mx,nx,nlev = size(ϕ)

    @boundscheck size(ϕ) == size(Tabs) || throw(BoundsError())
    @boundscheck (mx,nx) == size(ϕ0)   || throw(BoundsError())

    @unpack xgeop1, xgeop2, lapserate_correction = G.geometry

    # BOTTOM LAYER
    # is last index k=end, integration over half a layer
    for j in 1:nx
        for i in 1:mx
            ϕ[i,j,end] = ϕ0[i,j] + xgeop1[end]*Tabs[i,j,end]
        end
    end

    # OTHER LAYERS
    # integrate two half-layers from bottom to top
    for k in nlev-1:-1:1
        for j in 1:nx
            for i in 1:mx
                ϕ[i,j,k] = ϕ[i,j,k+1] + xgeop2[k+1]*Tabs[i,j,k+1] + xgeop1[k]*Tabs[i,j,k]
            end
        end
    end

    # LAPSERATE CORRECTION IN THE FREE TROPOSPHERE (>nlev)
    # TODO only for spectral coefficients 1,: ?
    for k in 2:nlev-1
        for j in 1:nx
            ϕ[1,j,k] = ϕ[1,j,k] +
                    lapserate_correction[k-1]*(Tabs[1,j,k+1] - Tabs[1,j,k-1])
        end
    end
end


text:
# equality and speed tests for ZPoly, GF2Poly, and Poly

module PTests

using Polynomials, IntModN, Compat, Base.Test

import Base: promote_rule, convert, ==

export tests

convert{T}(::Type{ZPoly{T}}, p::Poly{T}) = ZP(reverse(p.a))
# cannot use promotion with poly as not a Number
=={T}(a::ZPoly{T}, b::Poly{T}) = a == convert(ZPoly{T}, b)
=={T}(a::Poly{T}, b::ZPoly{T}) = convert(ZPoly{T}, a) == b
=={U<:Unsigned,J<:Integer,I<:Integer}(a::GF2Poly{U}, b::Poly{ZField{J,2,I}}) = convert(ZPoly{ZField{J,2,I}}, a) == convert(ZPoly{ZField{J,2,I}}, b)
=={U<:Unsigned,J<:Integer,I<:Integer}(a::Poly{ZField{J,2,I}}, b::GF2Poly{U}) = convert(ZPoly{ZField{J,2,I}}, a) == convert(ZPoly{ZField{J,2,I}}, b)

function make_random(deg, modulus)
    T = ZField{Int,modulus,Int}
    a = rand!(T, Array(T, rand(0:deg+1)))
    p = ZP(a)
    p, Poly(reverse(a))
end

function make_randoms(n, deg, modulus)
    a = @compat Tuple{ZPoly{ZField{Int,modulus,Int}}, Poly{ZField{Int,modulus,Int}}}[]
    for _ in 1:n
        push!(a, make_random(deg, modulus))
    end
    (a, @compat Tuple{ZPoly{ZField{Int,modulus,Int}}, Poly{ZField{Int,modulus,Int}}})
end

function make_random2(deg)
    T = ZField{Int,2,Int}
    a = rand!(T, Array(T, rand(0:deg+1)))
    p = ZP(a)
    convert(GF2Poly{UInt}, p), p, Poly(reverse(a))
end

function make_randoms2(n, deg)
    a = @compat Tuple{GF2Poly{UInt}, ZPoly{ZField{Int,2,Int}}, Poly{ZField{Int,2,Int}}}[]
    for _ in 1:n
        push!(a, make_random2(deg))
    end
    (a, (GF2Poly{UInt}, ZPoly{ZField{Int,2,Int}}, Poly{ZField{Int,2,Int}}))
end

function test_op(a, idx, op, T)
    check_zero = in(op, (rem, div, %, /))
    ZERO = zero(T[1])
    for i in 1:length(a)
        for j in 1:length(a)
            if !check_zero || a[j][1] != ZERO
                op(a[i][idx], a[j][idx])
            end
        end
    end
end

function do_timing(n, deg)

    # warm up
    a, T = make_randoms2(10, deg)
    for op in (+, -, *, /, %)
        for idx in 1:3
            test_op(a, idx, op, T)
        end
    end

    srand(1)
    a, T = make_randoms2(n, deg)
    println("\nfirst poly: ", a[1][1])
    for idx in 1:3
        print("$(idx): $(T[idx])]\n")
    end
    for op in (+, -, *, /, %)
        println("\n$op")
        for idx in 1:3
            gc()
            if VERSION < v"0.4-"
                gc_disable()
            else
                gc_enable(false)
            end
            print("$(idx): ")
            @time test_op(a, idx, op, T)
            if VERSION < v"0.4-"
                gc_enable()
            else
                gc_enable(true)
            end
        end
    end
end

#do_timing(1000, 8)

function test_eq(a, T, op)
    for i1 in 1:length(a)
        for i2 in 1:length(a)
            for t1 in 1:length(a[1])  # all same length
                p1 = a[i1][t1]
                q1 = a[i2][t1]
                for x in op
                    r1 = nothing
                    try
                        r1 = x(p1, q1)
                    catch
                    end
                    for t2 in 1:length(a[1])
#                        println("test $i1 $i2 $(T[t1]) $(T[t2]) $x")
                        p2 = a[i1][t2]
                        q2 = a[i2][t2]
                        r2 = nothing
                        try
                            r2 = x(p2, q2)
                        catch
                        end
                        if r1 != r2
                            println("($p1) $x ($q1)")
                            println("r1  $(typeof(r1))  $(r1)")
                            println("($p2) $x ($q2)")
                            println("r2  $(typeof(r2))  $(r2)")
                        end
                        @test r1 == r2
                    end
                end
            end
        end
    end
end

function do_eq(n, deg)
    print("\ntesting eq (mod 2)...")
    a, T = make_randoms2(n, deg)
    test_eq(a, T, (+, -, *, /, %))
    println("ok")

    print("testing eq (mod 5)...")
    a, T = make_randoms(n, deg, 5)
    test_eq(a, T, (+, -, *, /, %))
    println("ok")
end

#do_eq(30, 8)

function ptests()
    do_timing(1000, 8)
    do_eq(30, 2)
end

end


text:
module ApproxFunFromFile

export ApproxFunVandermonde

using ApproxFun

#pts=points(c,100) #This is a built in function; Approxfun would then do eg. vals=cos(pts) to get the values
#show(pts) # Density looks like a 'U' with greater sampling density near the extremes of the range

# From: https://github.com/ApproxFun/ApproxFun.jl/issues/275 , courtesy of private communication with Sheehan Olver
# Least squares approximation of data on an evenly spaced grid with Chebyshev series
function vandermonde(S,n,x::AbstractVector)
    V=Array(Float64,length(x),n)
    for k=1:n
        V[:,k]=Fun([zeros(k-1);1],S)(x)
    end
    V
end

# For ...(this)... case, make sure `length(pts) >> n`.
function ApproxFunVandermonde(filename,n=20,  lower=0.0, upper=360.0)
    c=Chebyshev([lower,upper]) #Define Chebyshev domain in this range (to match data imported)

    # Standard two column data form
    df=readdlm(filename)

    pts=df[:,1] # Points
    vals=df[:,2] #Values at these points

    V=vandermonde(c,n,pts)
    # Are you ready for the magic?
    af=Fun(V\vals,c) # Approximate Function (af)
    # me is now an ApproxFun representation of the tabulated data.
    # As a Chebyshev polynomial fit we can do all sorts of differentiation + integration.
    return af,df
end

function graphFun(af,df)
    # Logscale version...
    plot(af,label="Chebyshev fn")
    plot!(df[:,1],df[:,2],label="Raw fn")
    yaxis!("Value fn(x)")
    xaxis!("x")
end

function basictest()
    U,df=ApproxFunVandermonde("test.dat")
    using Plots
    graphFun(U,df)
    printFun(U)

    println(U)
    println(df)
end

function printFun(af)
    for i=0:360
        @printf("%f %f\n",i,af(i))
    end
end

end

max_issues_repo_name: jarvist/Teclo

max_forks_repo_name: jarvist/Teclo


text:
# @doc """
# EM algorithm for HMM
#     phghm, ph1, pvgh, loglikelihood = HMMem(v, H, V; MaxIterations=10, PlotProgress=false)
#
# ##### Inputs:
# * `v`: cell array containing sequences, so v{3}(5) is the 3rd sequence, 5th timestep
# * `H`: number of hidden states
# * `V`: number of visible (observation) states
#
# ##### Outputs:
# * `phghm`: learned transition p(h(t)|h(t-1))
# * `ph1`: learned initial distribution p(h(1))
# * `pvgh`: learned emission p(v|h)
# * `loglikelihood`: log likelihood of the sequences
# """ ->
function HMMem(v::Array{Any,1}, H, V; MaxIterations=10, PlotProgress=false)
    N = length(v) # number of sequences

    # random initialisation:
    phghm = condp(rand(H,H))    # transition distribution p(h(t)|h(t-1))
    pvgh = condp(rand(V,H))     # emission distribution p(v(t)|h(t))
    ph1 = condp(rand(H,1))      # initial p(h(1))
    llik = zeros(MaxIterations, N)
    loglik = zeros(N)

    if PlotProgress
        PyPlot.figure()
        ax=PyPlot.subplot(1,1,1)
    end

    for emloop = 1:MaxIterations
        A = zeros(H, H)
        a = zeros(H, 1)
        B = zeros(V, H)
        for n = 1:N

        # Perform Inference tasks (E step):
            alpha, loglik = HMMforward(v[n], phghm, ph1, pvgh)
            llik[emloop, :] = llik[emloop, :] + loglik
            beta = HMMbackward(v[n], phghm, pvgh)
            phtgV1T, phthtpgV1T = HMMsmooth(v[n], phghm, ph1, pvgh, alpha, beta, ReturnPairwiseMarginals=true)

        # Collect the statistics for the M-step:
            A = A + sum(phthtpgV1T, 3)[:, :, 1]
            a = a + phtgV1T[:, 1]
            for t = 1:length(v[n])
                i = v[n][t]
                #B[i, :] = B[i, :] + phtgV1T[:, t]'
                B[i, :] = B[i, :] + phtgV1T[:, t]
            end
        end
        # Perform M-step
        ph1 = condp(a)
        phghm = condp(A')
        pvgh = condp(B)
    end
    totalLogLik=sum(llik,2)
    if PlotProgress
        PyPlot.plot(totalLogLik); PyPlot.suptitle("log likelihood")
        println(minimum(totalLogLik))
        ax[:set_ylim](totalLogLik[2],totalLogLik[end])
    end
    loglikelihood = loglik[end]
    return phghm, ph1, pvgh, loglikelihood
end


text:
function sensorloc_data(seed, n)

    rng = Random.MersenneTwister(seed)

    # n = number of sensors points
    # m = number of anchor points
    m = floor(Int, 0.1 * n)

    # Sensor true position (2 dimensional)
    x_true = Random.rand(rng, Float64, (2, n))
    # Distances from sensors to sensors
    d = Dict((i, j) => norm(x_true[:, i] - x_true[:, j]) for i in 1:n for j in 1:i)

    # Anchor positions
    a = Dict(i => Random.rand(rng, Float64, (2, 1)) for i in 1:m)
    # Distances from anchor to sensors
    d_bar = Dict((k, j) => norm(x_true[:, j] - a[k]) for k in 1:m for j in 1:n)

    return m, x_true, a, d, d_bar
end

function sensorloc_eval(n, m, x_true, XX)
    @show norm(x_true - XX[1:2, 3:n + 2])
    @show rank = length([eig for eig in eigen(XX).values if eig > 1e-7])
    return nothing
end


text:
using Revise
using SparseTimeSeries

y1 = EventSeries(range(0, step=0.1, length=5) |> collect, rand(1:100, 5)) # create an EventSeries from a timestamps and a value array

##

y2 = EventSeries(range(0, step=0.05, length=6), ['A','A','A','B','B','C']) # removes repeated values by default

##

y = TaggedEventSeries()
y[:Rand] = y1 # assign EventSeries by setindex!
y[:Char] = y2
push!(y, TaggedEvent(:Int, 3.0, 123)) # push value to new tag
push!(y, TaggedEvent(:Char, 10.0, 'D')) # push new value (more recent) to existing tag

tagged_events(y) # sorts all events according to time

##

fill_forward_event(y, 2.0)

##

fill_forward_value(y, 5.0)

##
EventSeries(y)


text:
using Statistics
using Test
using GeoPhyInv

@init_parallel_stencil(2, false, Float64, 4)

# testing attenuation
#=
np2 = nextpow(2, 2 * length(tgrid));
fnpow2grid = FFTW.fftfreq(np2, inv(step(tgrid)));
fill!(medium.m[:Q], 100)
medium.fc = NamedArray([20.0, 100 * 2 * pi], ([:freqmin, :freqmax],))
medium.ic = NamedArray([5], ([:nsls],))

vp0 = medium.ref[:vp]
rho0 = medium.ref[:rho]
Ku = rho0 * vp0 * vp0
tau_epsilon = medium[:tau_epsilon]
tau_sigma = medium[:tau_sigma]

Kc = GeoPhyInv.complexK(Ku, 2.0 * pi .* fnpow2grid, tau_sigma, tau_epsilon)
=#

function check()
    # least-squares misfit
    paerr = GeoPhyInv.VNamedD_misfit(rec1, pa.c.data[1])
    err = GeoPhyInv.func_grad!(paerr)

    # normalize error
    error = err[1] / paerr.ynorm

    # desired accuracy?
    @test error < 1e-2
end

# without attenuation
medium = Medium(:acou_homo2D, 5);
ageom = AGeom(medium.mgrid, :xwell);
tgrid = range(0.0, stop = 2.0, length = 2000)
wav = ricker(10.0, tgrid, tpeak = 0.25);
srcwav = SrcWav(tgrid, ageom, [:p])
update!(srcwav, [:p], wav)

vp0 = mean(medium[:vp])
rho0 = mean(medium[:rho])
rec1 = GeoPhyInv.Born.mod(
    medium,
    medium_pert = medium,
    ageom = ageom,
    srcwav = srcwav,
    tgridmod = tgrid,
    src_flag = 2,
)

pa = SeisForwExpt(
    FdtdAcou(),
    npw = 1,
    medium = medium,
    ageom = [ageom],
    srcwav = [srcwav],
    sflags = [2],
    rflags = [1],
    tgrid = tgrid,
    verbose = true,
);

@time update!(pa);

check()

# @info "testing with attenuation"
# medium = Medium(:acou_homo2D);
# medium = Medium(medium.mgrid, [:vp, :rho, :Q])
# update!(medium, [:vp, :rho, :Q], [[1500, 2500], [1500, 2500], [10, 10]])
# fill!(medium)

# ageom = AGeom(medium.mgrid, :xwell);
# tgrid = range(0.0, stop = 2.0, length = 1000)
# wav = ricker(10.0, tgrid, tpeak = 0.25);
# srcwav = SrcWav(tgrid, ageom, [:p])
# update!(srcwav, [:p], wav)

# pa = SeisForwExpt(
#     FdtdAcouVisco(),
#     npw = 1,
#     medium = medium,
#     ageom = [ageom],
#     srcwav = [srcwav],
#     sflags = [2],
#     rflags = [1],
#     tgrid = tgrid,
#     verbose = true,
# );

# @time update!(pa);

# vp0 = mean(medium[:vp])
# rho0 = mean(medium[:rho])
# rec1 = GeoPhyInv.Born.mod(
#     pa.c.exmedium,
#     medium_pert = pa.c.exmedium,
#     ageom = ageom,
#     srcwav = srcwav,
#     tgridmod = tgrid,
#     src_flag = 2,
# )

# check()


text:
"""
$(TYPEDSIGNATURES)
Create a context with the provided container.
"""
function Context(id::UInt, container::T) where {T}
    history = Stack{T}()
    push!(history, container)
    return Context(id, history, Symbol[])
end

"""
$(TYPEDSIGNATURES)
Save the current context to history.
"""
function save(c::Context)
    push!(getfield(c, :history), deepcopy(c.data))
end

"""
$(TYPEDSIGNATURES)
Restore to the last saved context.
"""
function restore(c::Context)
    pop!(getfield(c, :history))
end

# Extensions to Base functions

function Base.show(io::IO, c::Context)
    print(io, "Context(id=", c.hex_id, ",generations=", c.generations, ")")
end

# Standard context management

function Base.push!(c::Context, entry)
    push!(c.data, entry)
end

Base.getindex(c::Context, index) = c.data[index]

Base.empty!(c::Context) = empty!(c.data)

Base.length(c::Context) = length(c.data)

Base.iterate(c::Context) = iterate(c.data)
Base.iterate(c::Context, state) = iterate(c.data, state)

# Property interface

Base.propertynames(c::Context) = (:id, :data, :generations, :hex_id, :path)

function Base.getproperty(c::Context, s::Symbol)
    s === :id && return getfield(c, :id)
    s === :data && return first(getfield(c, :history))
    s === :generations && return length(getfield(c, :history))
    s === :hex_id && return string("0x", string(getfield(c, :id); base = 16))
    s === :path && return getfield(c, :path)
    throw(UndefVarError(s))
end


text:
# Bio.Structure
# =============
#
# Module for structural biology.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Structure

using Compat

import Bio.IO: FileFormat, AbstractParser

include("model.jl")
include("pdb.jl")
include("spatial.jl")

end # Structure


text:
using Enzyme, BenchmarkTools, FiniteDiff, LinearAlgebra

f(x) = x^2

f(x) = x[1]*x[2]

inp =  [2.0 1.0
        0.0 1.0]
dinp = [0.0 0.0
        0.0 0.0]

f2(x) = inv(x)

@benchmark autodiff(det, Duplicated(inp, dinp))

@benchmark FiniteDi


text:
using BranchAndPrune

struct Interval
    lo::Float64
    hi::Float64
end

# Search the zero of a monotonic function
struct ZeroSearch <: AbstractDepthFirstSearch{Interval}
    f::Function
    initial::Interval
    tol::Float64
end

function BranchAndPrune.process(search::ZeroSearch, interval)
    ylo = search.f(interval.lo)
    yhi = search.f(interval.hi)

    # If both have the same time they are on the same side of the zero
    if ylo*yhi > 0
        return :discard, interval
    elseif interval.hi - interval.lo < search.tol
        return :store, interval
    else
        return :bisect, interval
    end
end

function BranchAndPrune.bisect(::ZeroSearch, interval)
    m = (interval.hi + interval.lo)/2
    return Interval(interval.lo, m), Interval(m, interval.hi)
end

function find_zero(f, interval)
    search = ZeroSearch(f, interval, 1e-10)

    local endtree = nothing
    niter = 0

    for wt in search
        endtree = wt
        niter += 1
    end

    println("Search finished in $niter iterations.")
    d = data(endtree)  # Retrieve the data from the tree
    if length(d) == 0
        println("The function has no zero.")
    else
        # If there is a zero, the tree will have only one data available
        z = d[1]
        println("The function has a zero in the interval [$(z.lo), $(z.hi)].")
    end
end

find_zero(x -> x/3 + 5, Interval(-20, 20))  # Exact solution is -15
find_zero(x -> (x - 4)^3 - 8, Interval(-20, 20))  # Exact solution is 6


text:
"$(SIGNATURES) Generate the indices into the concatenated constraint vector for each constraint.
Determines the bandedness of the Jacobian"
function gen_con_inds(conSet::ConstraintList, structure=:by_knotpoint)
	n,m = conSet.n, conSet.m
    N = length(conSet.p)
    numcon = length(conSet.constraints)
    conLen = length.(conSet.constraints)

    cons = [[@SVector ones(Int,length(con)) for j in eachindex(conSet.inds[i])]
		for (i,con) in enumerate(conSet.constraints)]

    # Dynamics and general constraints
    idx = 0
	if structure == :by_constraint
	    for (i,con) in enumerate(conSet.constraints)
			for (j,k) in enumerate(conSet.inds[i])
				cons[i][TO._index(con,k)] = idx .+ (1:conLen[i])
				idx += conLen[i]
	        end
	    end
	elseif structure == :by_knotpoint
		for k = 1:N
			for (i,con) in enumerate(conSet.constraints)
				inds = conSet.inds[i]
				if k in inds
					j = k -  inds[1] + 1
					cons[i][j] = idx .+ (1:conLen[i])
					idx += conLen[i]
				end
			end
		end
	end
    return cons
end

"$(SIGNATURES)
Get the constraint Jacobian structure as a sparse array, and fill in the linear indices
used for filling a vector of the non-zero elements of the Jacobian"
function constraint_jacobian_structure(solver::ConstrainedSolver,
		structure=:by_knopoint)
    n,m,N = size(solver)
    conSet = get_constraints(solver)
    idx = 0.0
    linds = jacobian_linear_inds(solver)

    NN = num_primals(solver)
    NP = num_duals(solver)
    D = spzeros(Int,NP,NN)

    # Number of elements in each block
    blk_len = map(con->length(con.∇c[1]), conSet.constraints)

    # Number of knot points for each constraint
    con_len = map(con->length(con.∇c), conSet.constraints)

    # Linear indices
	if structure == :by_constraint
	    for (i,con) in enumerate(conSet.constraints)
	        for (j,k) in enumerate(con.inds)
	            inds = idx .+ (1:blk_len[i])
	            linds[i][j] = inds
	            con.∇c[j] = inds
	            idx += blk_len[i]
	        end
	    end
	elseif structure == :by_knotpoint
		for k = 1:N
			for (i,con) in enumerate(conSet.constraints)
				if k in con.inds
					inds = idx .+ (1:blk_len[i])
					j = TO._index(con,k)
					linds[i][j] = inds
					con.∇c[j] = inds
					idx += blk_len[i]
				end
			end
		end
	end

    copy_jacobians!(D, solver)
    return D
end


text:
immutable InverseGaussian <: ContinuousUnivariateDistribution
    μ::Float64
    λ::Float64

    function InverseGaussian(μ::Real, λ::Real)
        (μ > zero(μ) && λ > zero(λ)) || error("InverseGaussian's μ and λ must be positive")
        new(float64(μ), float64(λ))
    end

    InverseGaussian(μ::Real) = InverseGaussian(μ, 1.0)
    InverseGaussian() = new(1.0, 1.0)
end

@distr_support InverseGaussian 0.0 Inf

#### Parameters

shape(d::InverseGaussian) = d.λ
params(d::InverseGaussian) = (d.μ, d.λ)

#### Statistics

mean(d::InverseGaussian) = d.μ

var(d::InverseGaussian) = d.μ^3 / d.λ

skewness(d::InverseGaussian) = 3.0 * sqrt(d.μ / d.λ)

kurtosis(d::InverseGaussian) = 15.0 * d.μ / d.λ

function mode(d::InverseGaussian)
    μ, λ = params(d)
    r = μ / λ
    μ * (sqrt(1.0 + 2.25 * r^2) - 1.5 * r)
end

#### Evaluation

function pdf(d::InverseGaussian, x::Float64)
    if x > 0.0
        μ, λ = params(d)
        return sqrt(λ / (twoπ * x^3)) * exp(-λ * (x - μ)^2 / (2.0 * μ^2 * x))
    else
        return 0.0
    end
end

function logpdf(d::InverseGaussian, x::Float64)
    if x > 0.0
        μ, λ = params(d)
        return 0.5 * (log(λ) - (log2π + 3.0 * log(x)) - λ * (x - μ)^2 / (μ^2 * x))
    else
        return -Inf
    end
end

function cdf(d::InverseGaussian, x::Float64)
    if x > 0.0
        μ, λ = params(d)
        u = sqrt(λ / x)
        v = x / μ
        return Φ(u * (v - 1.0)) + exp(2.0 * λ / μ) * Φ(-u * (v + 1.0))
    else
        return 0.0
    end
end

function ccdf(d::InverseGaussian, x::Float64)
    if x > 0.0
        μ, λ = params(d)
        u = sqrt(λ / x)
        v = x / μ
        Φc(u * (v - 1.0)) - exp(2.0 * λ / μ) * Φ(-u * (v + 1.0))
    else
        return 1.0
    end
end

function logcdf(d::InverseGaussian, x::Float64)
    if x > 0.0
        μ, λ = params(d)
        u = sqrt(λ / x)
        v = x / μ
        a = logΦ(u * (v -1.0))
        b = 2.0 * λ / μ + logΦ(-u * (v + 1.0))
        a + log1pexp(b - a)
    else
        return -Inf
    end
end

function logccdf(d::InverseGaussian, x::Float64)
    if x > 0.0
        μ, λ = params(d)
        u = sqrt(λ / x)
        v = x / μ
        a = logΦc(u * (v - 1.0))
        b = 2.0 * λ / μ + logΦ(-u * (v + 1.0))
        a + log1mexp(b - a)
    else
        return 0.0
    end
end

# TODO: need a more accurate method
function quantile(d::InverseGaussian, p::Float64)
    if 0.0 < p < 1.0
        # Whitmore and Yalovsky (1978) approximation
        μ, λ = params(d)
        w = Φinv(p)
        φ = λ / μ
        x = μ * exp(w / sqrt(φ) - 0.5 / φ)
        # one multiplicitive Newton iteration
        xn = x * exp((p - cdf(d, x)) / (pdf(d, x) * exp(x)))
        # additive Newton iterations
        while !isapprox(xn,x)
            x = xn
            xn = x + (p - cdf(d, x)) / pdf(d, x)
        end
        return xn
    elseif p <= 0.0
        return 0.0
    else
        return Inf
    end
end

#### Sampling

# rand method from:
#   John R. Michael, William R. Schucany and Roy W. Haas (1976)
#   Generating Random Variates Using Transformations with Multiple Roots
#   The American Statistician , Vol. 30, No. 2, pp. 88-90
function rand(d::InverseGaussian)
    μ, λ = params(d)
    z = randn()
    v = z * z
    w = μ * v
    x1 = μ + μ / (2.0 * λ) * (w - sqrt(w * (4.0 * λ + w)))
    p1 = μ / (μ + x1)
    u = rand()
    u >= p1 ? μ^2 / x1 : x1
end


text:
using DifferentialEquations, Flux, DiffEqFlux
using Plots, NPZ
using IterTools: ncycle
using BSON: @save, @load
include("./Custom_DataLoader.jl")
# using .Custom_DataLoader

# Julia starts from 1
true_data = npzread("SWE_RNN_Training_Data.npy")
true_data = permutedims(true_data, (1,3,2))
data_len = size(true_data)[3]
state_data = true_data[:,1:8,:] # last index inclusive

tspan = (0.0,true_data[1,9,data_len])
tsteps = 0.0:tspan[2]/data_len:tspan[2]
tsteps = tsteps[1:100]
T_ = Array{Float64, 3}(undef, size(true_data)[1], 1, size(true_data)[3])

# Outer For-loop
for i in 1:size(true_data)[1]
    T_[i,1,:] = tsteps[:]
end

# print(size(T_),size(state_data))

test_data = npzread("SWE_RNN_Testing_Data.npy")
test_data = permutedims(test_data, (1,3,2))
u0 = test_data[1,1:8,1]

my_nn = Chain(
  Dense(8, 32, relu),
  Dense(32, 8)) |> f32

p, re = Flux.destructure(my_nn)

function right_hand_side(du,u,p,t)
  m = re(p)
  nn_output = m(u)
  du[1] = nn_output[1]
  du[2] = nn_output[2]
end

prob = ODEProblem(right_hand_side,u0,tspan,p)
sol = solve(prob, Tsit5())

# Plot the solution
plot(sol)
savefig("Initial_NN_ODE.png")

function predict_ode(state_batch,time_batch)
  tmp_prob = remake(prob; p = p,u=state_batch[1])
  pred = Array(solve(tmp_prob, Tsit5(), saveat = time_batch))
  return pred
end

function loss_ode(state_batch,time_batch)
    pred = predict_ode(state_batch,time_batch)
    loss = sum(abs2, state_batch-pred)
    return loss
end

k = 20
train_loader = MyDataLoader((state_data,T_), batchsize = k, shuffle=true)
numEpochs = 200
losses=[]

# test_batch = iterate(train_loader)
# print(size(test_batch[1][2]))
# loss_ode(test_batch[1][1],test_batch[1][2])
# loss_ode(test_data[1,1:8,:],tsteps)

cb() = begin

  l=loss_ode(test_data[1,1:8,:],tsteps)
  push!(losses, l)
  @show l

  # pred = predict_ode(state_data,tsteps)

  # plt = plot(x_data,hcat(transpose(pred),state_data),
  #   line=(4, [:solid :solid :dash :dash]), ylim = (-6, 6),
  #   label=["y1-pred" "y2-pred" "y1-true" "y2-true"])

  # display(plt)
  # Tell sciml_train to not halt the optimization. If return true, then
  # optimization stops.
  return false
end

opt=ADAM(0.01)
Flux.train!(loss_ode, Flux.params(p), ncycle(train_loader,numEpochs), opt, cb=cb) # cb=Flux.throttle(cb, 10)

my_nn_final = re(p)
@save "mymodel.bson" my_nn_final

# Inference
@load "mymodel.bson" my_nn_final
p, re = Flux.destructure(my_nn_final)

prob_final = ODEProblem(right_hand_side, u0, tspan, p)
sol_final = solve(prob_final,Tsit5())
plot(sol_final)
savefig("NN_ode_optimized.png") # Both state components close to 1

max_stars_repo_name: Romit-Maulik/Tutorials-Demos-Practice

max_issues_repo_name: Romit-Maulik/Tutorials-Demos-Practice

max_forks_repo_name: Romit-Maulik/Tutorials-Demos-Practice


text:
using Dolo

# model = yaml_import("examples/models/consumption_savings_iid.yaml")

model = yaml_import("examples/models/rbc.yaml")

Dolo.perturb(model)

exit()

@time sol = Dolo.time_iteration(model, verbose=true, ignore_constraints=false, maxit=5);

# simulate(model, sol.dr)

sol = Dolo.improved_time_iteration(model; verbose=true, ignore_constraints=false)

simulate(model, sol.dr)

exit()

using SimplePlots

tab = tabulate(model, sol.dr, :w)

plot(tab[:w], tab[:c])

function step(model)
    sol = Dolo.improved_time_iteration(model; verbose=false, ignore_constraints=false);
    dist = Dolo.ergodic_distribution(model, sol);
end

@time step(model);

using Dolo: Euler, MSM
using StaticArrays

function dFdp(F::Euler, x0::MSM, x1::MSM, p0::SVector, dp::SVector, ε=1e-8)

    no = maximum(abs, dp)
    println(no)
    p1 = p0 + (dp/no)*ε

    f0 = F(x0, x1; p=p0)
    f1 = F(x0, x1; p=p1)

    df = ((f1-f0)/ε) * no

    return df

end

function update_guess(F, x0, p0, dp, tol=1e-8, maxit=1000)
    df = dFdp(F, x0, x0, p0, dp)
    J = Dolo.df_A(F, x0, x0)
    L = Dolo.df_B(F, x0, x0)
    Dolo.prediv!(L, J)
    Dolo.mult!(L, -1.0)
    dπ = -J\df
    dx = dπ
    i = 0
    while i<=maxit
        i += 1
        err = Dolo.norm(dπ)
        if err<tol
            break
        end
        dπ = L*dπ
        dx += dπ
    end
    return dx, i

end

sol =  improved_time_iteration(model, verbose=false);

dp = SVector(0.01, 0, 0, 0, 0, 0)

@time res = Dolo.time_iteration(model, ignore_constraints=false, verbose=false);

xx = MSM([sol.dr(1, F.grid_endo.nodes)])

@time dx, i = update_guess(F, xx, F.p, dp);

F(xx+dx, xx+dx; p=p+dp, set_future=true);

@time J = Dolo.df_A(F, F.x0, F.x0, inplace=false);
@time J = Dolo.df_A(F, F.x0, F.x0, inplace=true);

Jii = Dolo.df_A(F, F.x0, F.x0, inplace=true);

f0 = z0*0
fi = z0*0
xi = z0*0

using StaticArrays
using Dolo: MSM
n_x = length(f0.data[1])
N = length(xi.data)
cata = zeros(SMatrix{n_x, n_x, Float64, n_x*n_x}, N)
J = MSM(cata, f0.sizes)
out = (f0, fi, xi, J)

function fun!(out, u)
    r = F(u,F.x0)
    out.data[:] .= r.data
end

fun!(f0, F.x0)

fp = f0*0

fp, Jp = Dolo.DiffFun!(fun!, F.x0, 1e-8, out);

fp.data[1]

F(F.x0, F.x0).data[1]


text:
m = rand(3, 3)
axs = (1:3, 1:3)
ty = Interpolations.BSplineInterpolation{Float64, ndims(m), typeof(m), BSpline{Linear{Throw{OnGrid}}}, typeof(axs)}
etty = Interpolations.Extrapolation{Float64, 2, Interpolations.BSplineInterpolation{Float64, 2, Matrix{Float64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Flat{Nothing}}
fieldsum(x) = mapreduce(y->getfield(x, y), +, ntuple(identity, nfields(x)))

@testset "Interpolations.BSplineInterpolation constructor gradient" begin
    zg = Zygote.gradient((x,y,z)->sum(ty(x,y,z)), m, axs, BSpline(Linear()))
    fd = grad(central_fdm(5,1),
              x->sum(ty(x, axs, BSpline(Linear()))),
              m)
    @test isapprox(zg[1], fd[1])
    @test zg[2:3] == (nothing, nothing)
end

@testset "Interpolations.BSplineInterpolation functor gradient" begin
    zg = Zygote.gradient((x,y,z,w)->sum(Interpolations.BSplineInterpolation(x,y,z,w)), Float64, m, BSpline(Linear()), axs)
    fd = grad(central_fdm(5,1),
              x->sum(Interpolations.BSplineInterpolation(Float64, x, BSpline(Linear()), axs)),
              m)
    @test isapprox(zg[2], fd[1])
    @test findall(==(nothing), zg) == [1, 3, 4]
end

@testset "Interpolations.Extrapolation constructor gradient" begin
    itp = ty(m, axs, BSpline(Linear()))
    zg = Zygote.gradient((x,y)->sum(etty(x,y)), itp, Flat())
    fd = grad(central_fdm(5,1),
              x->sum(etty(itp, Flat())),
              m)
    @test isapprox(zg[1], fd[1])
    @test findall(==(nothing), zg) == [2]
end

@testset "Interpolations.copy_with_padding function gradient" begin
    zg = Zygote.gradient((x,y,z)->sum(Interpolations.copy_with_padding(x,y,z)), Float64, m, BSpline(Linear()))
    fd = grad(central_fdm(5,1),
              x->sum(Interpolations.copy_with_padding(Float64, x, BSpline(Linear()))),
              m)
    @test isapprox(zg[2], fd[1])
    @test findall(==(nothing), zg) == [1, 3]
end

@testset "Interpolations.FilledExtrapolation constructor gradient" begin
    @test_broken begin
        itp = ty(m, axs, BSpline(Linear()))
        zg = Zygote.gradient((x,y)->sum(Interpolations.FilledExtrapolation(x,y)), itp, 0.0)
        fd = grad(central_fdm(5,1),
                x->sum(Interpolations.FilledExtrapolation(x, 0.0)),
                itp)
        @test_broken isapprox(zg[1], fd[1])
        @test_broken findall(==(nothing), zg) == [2]
    end # TODO: Remove the @test_broken after tests work.
end

@testset "SVector and SMatrix gradient" begin
    for t in (Float32, Float64, RGB{Float32}, RGB{Float64})
        inp = rand(t, 2)
        inp_mat = rand(t, 3, 3)
        _abs(c)  = mapreducec(v->abs(float(v)), +, 0, c)
        @test Zygote.gradient(x->_abs(sum(SVector{2, t}(x))), inp)[1] == ones(t, 2)
        @test Zygote.gradient(x->_abs(sum(SMatrix{3, 3, t, 9}(x))), inp_mat)[1] == ones(t, 3, 3)
    end
end

@testset "DiffImages.Homography method gradient" begin
    h = DiffImages.Homography{Float64}()
    v = rand(SVector{2, Float64})
    zy = Zygote.gradient((x,y)->sum(y(x)), v, h)
    fd = grad(central_fdm(5,1), (x,y)->sum(y(x)), v, h)
    @test (zy[1] .≈ fd[1]) == Bool[1, 1]
    @test zy[2].H ≈ fd[2].H
    @test typeof(zy[2].H) == typeof(fd[2].H)
end

@testset "CoordinateTransformations.LinearMap gradient" begin
    m = CoordinateTransformations.LinearMap(rand(3, 3))
    sm = CoordinateTransformations.LinearMap(rand(SMatrix{3, 3, Float64}))
    v = rand(3)
    sv = rand(SVector{3, Float64})
    zy = Zygote.gradient((x, y)->sum(x(y)), m, v)
    fd = grad(central_fdm(5,1), (x, y)->sum(x(y)), m, v)
    zy_s = Zygote.gradient((x, y)->sum(x(y)), sm, sv)
    fd_s = grad(central_fdm(5,1), (x, y)->sum(x(y)), sm, sv)

    @test zy[1].linear ≈ fd[1].linear
    @test zy[2] ≈ fd[2]
    @test typeof(zy[1].linear) == typeof(fd[1].linear)
    @test zy_s[1].linear ≈ fd_s[1].linear
    @test zy_s[2] ≈ fd_s[2]
    @test_broken typeof(zy_s[1].linear) == typeof(fd_s[1].linear) # TODO: fix in another issue.
end

@testset "Rotations.RotMatrix{N, T, L} gradient" begin
    rm = RotMatrix(π/4)
    v = rand(2)
    zy = Zygote.gradient((x, y)->sum(LinearMap(x)(y)), rm, v)
    fd = grad(central_fdm(5,1), (x, y)->sum(LinearMap(x)(y)), rm, v)
    @test zy[1] ≈ fd[1]
    @test zy[2] ≈ fd[2]
end

@testset "ImageTransformations._getindex gradient" begin
    _sep(x) = x.r + x.g + x.b
    for t in (Float32, Float64, RGB{Float32}, RGB{Float64})
        itp = extrapolate(interpolate(rand(t, 3, 3), BSpline(Linear())), zero(t))
        for ind in ((2.5, 2.5), (5, 5))
            if t <: Colorant
                zy = Zygote.gradient((x,y)->_sep(ImageTransformations._getindex(x,y)), itp, ind)
                @test zy[2] ≈ fieldsum.(Interpolations.gradient(itp, ind...))
            else
                zy = Zygote.gradient((x,y)->ImageTransformations._getindex(x,y), itp, ind)
                @test zy[2] ≈ Interpolations.gradient(itp, ind...)
            end
        end
    end
end