text:
module IOModule

using NetCDF
using NetCDF: NcFile,NcVar,NC_NETCDF4
using GridModule: Grid, LatLonGrid, UnstructuredGrid, get_number_of_dimensions
using FieldModule: Field, LatLonField,UnstructuredField,LatLonDirectionIndicators
using FieldModule: UnstructuredDirectionIndicators
using FieldModule: round,convert,invert,add_offset,get_data,maximum,equals
using HDModule: RiverParameters,RiverPrognosticFields
using LakeModule: LakeParameters,LatLonLakeParameters, UnstructuredLakeParameters
using LakeModule: GridSpecificLakeParameters,LakeFields
using MergeTypesModule

import LakeModule: write_lake_numbers_field,write_lake_volumes_field
import LakeModule: write_diagnostic_lake_volumes_field
import HDModule: write_river_initial_values,write_river_flow_field

function get_ncells(file_name::AbstractString)
  river_directions::Array{Float64,1} = NetCDF.ncread(file_name,"FDIR")
  return size(river_directions,1)
end

function get_additional_grid_information(file_name::AbstractString)
  clat::Array{Float64,1} = NetCDF.ncread(file_name,"clat")
  clon::Array{Float64,1} = NetCDF.ncread(file_name,"clon")
  clat_bounds::Array{Float64,2} = NetCDF.ncread(file_name,"clat_bnds")
  clon_bounds::Array{Float64,2} = NetCDF.ncread(file_name,"clon_bnds")
  return clat,clon,clat_bounds,clon_bounds
end

function load_field(file_handle::NcFile,grid::Grid,variable_name::AbstractString,
                    field_type::DataType,;timestep::Int64=-1)
  variable::NcVar = file_handle[variable_name]
  values::Array{field_type,get_number_of_dimensions(grid)} = NetCDF.readvar(variable)
  if isa(grid,LatLonGrid)
    if timestep == -1
      values = permutedims(values, [2,1])
    else
      values = permutedims(values, [2,1])[:,:,timestep]
    end
  end
  return Field{field_type}(grid,values)
end

function load_3d_field(file_handle::NcFile,grid::Grid,variable_name::AbstractString,
                       field_type::DataType,layer=Int64,;timestep::Int64=-1)
  variable::NcVar = file_handle[variable_name]
  values::Array{field_type,(timestep == -1) ?
                            get_number_of_dimensions(grid)+1 :
                            get_number_of_dimensions(grid)+2 } = NetCDF.readvar(variable)
  local layer_values::Array{field_type,get_number_of_dimensions(grid)}
  if isa(grid,UnstructuredGrid)
    if timestep == -1
      layer_values = values[:,layer]
    else
      layer_values = values[:,layer,timestep]
    end
  end
  return Field{field_type}(grid,layer_values)
end

function load_array_of_fields(file_handle::NcFile,grid::Grid,variable_base_name::AbstractString,
                              field_type::DataType,number_of_fields::Int64)
  if number_of_fields == 1 && isa(grid,LatLonGrid)
    return Field{field_type}[load_field(file_handle,grid,variable_base_name,field_type)]
  else
    array_of_fields::Array{Field{field_type},1} = Array{Field{field_type}}(undef,number_of_fields)
    for i = 1:number_of_fields
      if isa(grid,LatLonGrid)
        array_of_fields[i] = load_field(file_handle,grid,variable_base_name*string(i),field_type)
      else
        array_of_fields[i] = load_3d_field(file_handle,grid,variable_base_name,field_type,i)
      end
    end
    return array_of_fields
  end
end

function prep_array_of_fields(dims::Array{NcDim},
                              variable_base_name::AbstractString,
                              array_of_fields::Array{Field{T},1},
                              fields_to_write::Vector{Pair}) where {T}
  number_of_fields::Int64 = size(array_of_fields,1)
  if number_of_fields == 1
    return prep_field(dims,variable_base_name,array_of_fields[1],fields_to_write)
  else
    for i = 1:number_of_fields
      prep_field(dims,variable_base_name*string(i),
                 array_of_fields[i],fields_to_write)
    end
  end
end

function prep_dims(grid::LatLonGrid)
  lat = NcDim("Lat",grid.nlat)
  lon = NcDim("Lon",grid.nlon)
  dims::Array{NcDim} = NcDim[ lon; lat ]
  return dims
end

function prep_dims(grid::UnstructuredGrid)
  ncells = NcDim("ncells",grid.ncells)
  dims::Array{NcDim} = NcDim[ ncells ]
  return dims
end

function prep_field(dims::Array{NcDim},variable_name::AbstractString,
                    field::Field,fields_to_write::Vector{Pair})
  local variable::NcVar
  if isa(field,LatLonField)
    variable = NcVar(variable_name,dims)
    push!(fields_to_write,Pair(variable,permutedims(get_data(field), [2,1])))
  elseif isa(field,UnstructuredField)
    variable = NcVar(variable_name,dims,atts=Dict("grid_type"=>"unstructured",
                                                  "coordinates"=>"clat clon"))
    push!(fields_to_write,Pair(variable,get_data(field)))
  end
end

function write_fields(grid::Grid,fields_to_write::Vector{Pair},filepath::AbstractString,
                      dims::Array{NcDim})
  variables::Array{NcVar} = NcVar[ pair.first for pair in fields_to_write ]
  if isa(grid,UnstructuredGrid)
    prepare_icon_file(grid,filepath,dims,variables)
  else
    NetCDF.create(filepath,variables,mode=NC_NETCDF4)
  end
  for (variable,data) in fields_to_write
    NetCDF.putvar(variable,data)
  end
end

function prepare_icon_file(grid::UnstructuredGrid,filepath::AbstractString,
                           dims::Array{NcDim},variables::Array{NcVar})
  vertices = NcDim("vertices",3)
  extended_dims::Array{NcDim} = vcat(NcDim[ vertices ], dims)
  clat = NcVar("clat",dims,atts=Dict("long_name"=>"center_latitude",
                                     "units"=>"radians",
                                     "bounds"=>"clat_bnds",
                                     "standard_name"=>"latitude"))
  clon = NcVar("clon",dims,atts=Dict("long_name"=>"center_longitude",
                                     "units"=>"radians",
                                     "bounds"=>"clon_bnds",
                                     "standard_name"=>"longitude"))
  clat_bounds = NcVar("clat_bnds",extended_dims)
  clon_bounds = NcVar("clon_bnds",extended_dims)
  extended_variables::Array{NcVar} = vcat(variables,NcVar[ clat, clon, clat_bounds, clon_bounds])
  NetCDF.create(filepath,extended_variables,mode=NC_NETCDF4)
  NetCDF.putvar(clat,grid.clat)
  NetCDF.putvar(clon,grid.clon)
  NetCDF.putvar(clat_bounds,grid.clat_bounds)
  NetCDF.putvar(clon_bounds,grid.clon_bounds)
end

function write_field(grid::Grid,variable_name::AbstractString,
                     field::Field,filepath::AbstractString)
  dims::Array{NcDim} = prep_dims(grid)
  local variable::NcVar
  if isa(grid,UnstructuredGrid)
    variable = NcVar(variable_name,dims,atts=Dict("grid_type"=>"unstructured",
                                                  "coordinates"=>"clat clon"))
    variables::Array{NcVar} = NcVar[ variable ]
    prepare_icon_file(grid,filepath,dims,variables)
  else
    variable = NcVar(variable_name,dims)
    NetCDF.create(filepath,variable,mode=NC_NETCDF4)
  end
  if isa(grid,LatLonGrid)
    NetCDF.putvar(variable,permutedims(get_data(field), [2,1]))
  else
    NetCDF.putvar(variable,get_data(field))
  end
end

function load_river_parameters(hd_para_filepath::AbstractString,grid::Grid;
                               day_length=86400.0,step_length=86400.0)
  println("Loading: " * hd_para_filepath)
  file_handle::NcFile = NetCDF.open(hd_para_filepath)
  local landsea_mask::Field{Bool}
  try
    if isa(grid,LatLonGrid)
      direction_based_river_directions = load_field(file_handle,grid,"FDIR",Float64)
      flow_directions = LatLonDirectionIndicators(round(Int64,direction_based_river_directions))
      landsea_mask =
        load_field(file_handle,grid,"FLAG",Bool)
    elseif isa(grid,UnstructuredGrid)
      index_based_river_directions = load_field(file_handle,grid,"FDIR",Float64)
      flow_directions = UnstructuredDirectionIndicators(round(Int64,index_based_river_directions))
      landsea_mask_int::Field{Int64} =
        load_field(file_handle,grid,"MASK",Int64)
      landsea_mask = equals(landsea_mask_int,1)
    end
    landsea_mask = invert(landsea_mask)
    river_reservoir_nums::Field{Int64} =
      round(Int64,load_field(file_handle,grid,"ARF_N",Float64))
    overland_reservoir_nums::Field{Int64} =
      round(Int64,load_field(file_handle,grid,"ALF_N",Float64))
    base_reservoir_nums::Field{Int64} = Field{Int64}(grid,1)
    river_retention_coefficients::Field{Float64} =
      load_field(file_handle,grid,"ARF_K",Float64)
    overland_retention_coefficients::Field{Float64} =
      load_field(file_handle,grid,"ALF_K",Float64)
    base_retention_coefficients::Field{Float64} =
      load_field(file_handle,grid,"AGF_K",Float64)
    return  RiverParameters(flow_directions,
                            river_reservoir_nums,
                            overland_reservoir_nums,
                            base_reservoir_nums,
                            river_retention_coefficients,
                            overland_retention_coefficients,
                            base_retention_coefficients,
                            landsea_mask,grid,day_length,
                            step_length)
  finally
    NetCDF.close(file_handle)
  end
end

function load_river_initial_values(hd_start_filepath::AbstractString,grid::Grid,
                                   river_parameters::RiverParameters)
  river_prognostic_fields::RiverPrognosticFields =
    RiverPrognosticFields(river_parameters)
  local reservoir_units_adjustment_factor
  println("Loading: " * hd_start_filepath)
  file_handle::NcFile = NetCDF.open(hd_start_filepath)
  try
    if ! isa(grid,UnstructuredGrid)
      reservoir_units_adjustment = 1.0/river_parameters.step_length
      river_prognostic_fields.river_inflow =
        load_field(file_handle,river_parameters.grid,"FINFL",Float64)
    else
      river_prognostic_fields.river_inflow = Field{Float64}(grid,0.0)
      reservoir_units_adjustment = 1.0
    end
    river_prognostic_fields.base_flow_reservoirs =
      load_array_of_fields(file_handle,river_parameters.grid,"FGMEM",Float64,
                            maximum(river_parameters.base_reservoir_nums))
    river_prognostic_fields.base_flow_reservoirs *= reservoir_units_adjustment
    river_prognostic_fields.overland_flow_reservoirs =
      load_array_of_fields(file_handle,river_parameters.grid,"FLFMEM",Float64,
                            maximum(river_parameters.overland_reservoir_nums))
    river_prognostic_fields.overland_flow_reservoirs *= reservoir_units_adjustment
    river_prognostic_fields.river_flow_reservoirs =
      load_array_of_fields(file_handle,river_parameters.grid,"FRFMEM",Float64,
                            maximum(river_parameters.river_reservoir_nums))
    river_prognostic_fields.river_flow_reservoirs *= reservoir_units_adjustment
  finally
    NetCDF.close(file_handle)
  end
  return river_prognostic_fields
end

function write_river_initial_values(hd_start_filepath::AbstractString,
                                    river_parameters::RiverParameters,
                                    river_prognostic_fields::RiverPrognosticFields)
  println("Writing: " * hd_start_filepath)
  if ! isa(river_parameters.grid,UnstructuredGrid)
    reservoir_units_adjustment = river_parameters.step_length
  else
    reservoir_units_adjustment = 1.0
  end
  river_prognostic_fields.base_flow_reservoirs *= reservoir_units_adjustment
  river_prognostic_fields.overland_flow_reservoirs *= reservoir_units_adjustment
  river_prognostic_fields.river_flow_reservoirs *= reservoir_units_adjustment
  fields_to_write::Vector{Pair} = Pair[]
  dims::Array{NcDim} = prep_dims(river_parameters.grid)
  prep_field(dims,"FINFL",
             river_prognostic_fields.river_inflow,fields_to_write)
  prep_array_of_fields(dims,"FGMEM",
                       river_prognostic_fields.base_flow_reservoirs,
                       fields_to_write)
  prep_array_of_fields(dims,"FLFMEM",
                       river_prognostic_fields.overland_flow_reservoirs,
                       fields_to_write)
  prep_array_of_fields(dims,"FRFMEM",
                       river_prognostic_fields.river_flow_reservoirs,
                       fields_to_write)
  write_fields(river_parameters.grid,fields_to_write,hd_start_filepath,dims)
end

function load_lake_parameters(lake_para_filepath::AbstractString,grid::Grid,hd_grid::Grid,
                              surface_model_grid::Grid)
  println("Loading: " * lake_para_filepath)
  file_handle::NcFile = NetCDF.open(lake_para_filepath)
  try
    grid_specific_lake_parameters::GridSpecificLakeParameters =
      load_grid_specific_lake_parameters(file_handle,grid,surface_model_grid)
    lake_centers::Field{Bool} =
      load_field(file_handle,grid,"lake_centers",Bool)
    connection_volume_thresholds::Field{Float64} =
      load_field(file_handle,grid,"connection_volume_thresholds",Float64)
    flood_volume_thresholds::Field{Float64} =
      load_field(file_handle,grid,"flood_volume_thresholds",Float64)
    flood_local_redirect::Field{Bool} =
      load_field(file_handle,grid,"flood_local_redirect",Bool)
    connect_local_redirect::Field{Bool} =
      load_field(file_handle,grid,"connect_local_redirect",Bool)
    additional_flood_local_redirect::Field{Bool} =
      load_field(file_handle,grid,"additional_flood_local_redirect",Bool)
    additional_connect_local_redirect::Field{Bool} =
      load_field(file_handle,grid,"additional_connect_local_redirect",Bool)
    merge_points::Field{MergeTypes} =
      load_field(file_handle,grid,"merge_points",Int64)
    return LakeParameters(lake_centers,
                          connection_volume_thresholds,
                          flood_volume_thresholds,
                          flood_local_redirect,
                          connect_local_redirect,
                          additional_flood_local_redirect,
                          additional_connect_local_redirect,
                          merge_points,grid,hd_grid,
                          surface_model_grid,
                          grid_specific_lake_parameters)
  finally
    NetCDF.close(file_handle)
  end
end

function write_river_flow_field(river_parameters::RiverParameters,river_flow_field::Field{Float64};timestep::Int64=-1)
  variable_name::String = "river_flow"
  filepath::String = timestep == -1 ? "/Users/thomasriddick/Documents/data/temp/river_model_results.nc" : "/Users/thomasriddick/Documents/data/temp/river_model_results_$(timestep).nc"
  write_field(river_parameters.grid,variable_name,river_flow_field,filepath)
end

function load_grid_specific_lake_parameters(file_handle::NcFile,grid::LatLonGrid,
                                            surface_model_grid::LatLonGrid)
  flood_next_cell_lat_index::Field{Int64} =
      load_field(file_handle,grid,"flood_next_cell_lat_index",Int64)
  flood_next_cell_lon_index::Field{Int64} =
      load_field(file_handle,grid,"flood_next_cell_lon_index",Int64)
  connect_next_cell_lat_index::Field{Int64} =
      load_field(file_handle,grid,"connect_next_cell_lat_index",Int64)
  connect_next_cell_lon_index::Field{Int64} =
      load_field(file_handle,grid,"connect_next_cell_lon_index",Int64)
  flood_force_merge_lat_index::Field{Int64} =
      load_field(file_handle,grid,"flood_force_merge_lat_index",Int64)
  flood_force_merge_lon_index::Field{Int64} =
      load_field(file_handle,grid,"flood_force_merge_lon_index",Int64)
  connect_force_merge_lat_index::Field{Int64} =
      load_field(file_handle,grid,"connect_force_merge_lat_index",Int64)
  connect_force_merge_lon_index::Field{Int64} =
      load_field(file_handle,grid,"connect_force_merge_lon_index",Int64)
  flood_redirect_lat_index::Field{Int64} =
      load_field(file_handle,grid,"flood_redirect_lat_index",Int64)
  flood_redirect_lon_index::Field{Int64} =
      load_field(file_handle,grid,"flood_redirect_lon_index",Int64)
  connect_redirect_lat_index::Field{Int64} =
      load_field(file_handle,grid,"connect_redirect_lat_index",Int64)
  connect_redirect_lon_index::Field{Int64} =
      load_field(file_handle,grid,"connect_redirect_lon_index",Int64)
  additional_flood_redirect_lat_index::Field{Int64} =
      load_field(file_handle,grid,"additional_flood_redirect_lat_index",Int64)
  additional_flood_redirect_lon_index::Field{Int64} =
      load_field(file_handle,grid,"additional_flood_redirect_lon_index",Int64)
  additional_connect_redirect_lat_index::Field{Int64} =
      load_field(file_handle,grid,"additional_connect_redirect_lat_index",Int64)
  additional_connect_redirect_lon_index::Field{Int64} =
      load_field(file_handle,grid,"additional_connect_redirect_lon_index",Int64)
  corresponding_surface_cell_lat_index::Field{Int64} =
      load_field(file_handle,surface_model_grid,"corresponding_surface_cell_lat_index",Int64)
  corresponding_surface_cell_lon_index::Field{Int64} =
      load_field(file_handle,surface_model_grid,"corresponding_surface_cell_lon_index",Int64)
  add_offset(flood_next_cell_lat_index,1,Int64[-1])
  add_offset(flood_next_cell_lon_index,1,Int64[-1])
  add_offset(connect_next_cell_lat_index,1,Int64[-1])
  add_offset(connect_next_cell_lon_index,1,Int64[-1])
  add_offset(flood_force_merge_lat_index,1,Int64[-1])
  add_offset(flood_force_merge_lon_index,1,Int64[-1])
  add_offset(connect_force_merge_lat_index,1,Int64[-1])
  add_offset(connect_force_merge_lon_index,1,Int64[-1])
  add_offset(flood_redirect_lat_index,1,Int64[-1])
  add_offset(flood_redirect_lon_index,1,Int64[-1])
  add_offset(connect_redirect_lat_index,1,Int64[-1])
  add_offset(connect_redirect_lon_index,1,Int64[-1])
  add_offset(additional_flood_redirect_lat_index,1,Int64[-1])
  add_offset(additional_flood_redirect_lon_index,1,Int64[-1])
  add_offset(additional_connect_redirect_lat_index,1,Int64[-1])
  add_offset(additional_connect_redirect_lon_index,1,Int64[-1])
  add_offset(corresponding_surface_cell_lat_index,1,Int64[-1])
  add_offset(corresponding_surface_cell_lon_index,1,Int64[-1])
  return LatLonLakeParameters(flood_next_cell_lat_index,
                              flood_next_cell_lon_index,
                              connect_next_cell_lat_index,
                              connect_next_cell_lon_index,
                              flood_force_merge_lat_index,
                              flood_force_merge_lon_index,
                              connect_force_merge_lat_index,
                              connect_force_merge_lon_index,
                              flood_redirect_lat_index,
                              flood_redirect_lon_index,
                              connect_redirect_lat_index,
                              connect_redirect_lon_index,
                              additional_flood_redirect_lat_index,
                              additional_flood_redirect_lon_index,
                              additional_connect_redirect_lat_index,
                              additional_connect_redirect_lon_index,
                              corresponding_surface_cell_lat_index,
                              corresponding_surface_cell_lon_index)
end

function load_grid_specific_lake_parameters(file_handle::NcFile,grid::UnstructuredGrid,
                                            surface_model_grid::UnstructuredGrid)
  flood_next_cell_index::Field{Int64} =
      load_field(file_handle,grid,"flood_next_cell_index",Int64)
  connect_next_cell_index::Field{Int64} =
      load_field(file_handle,grid,"connect_next_cell_index",Int64)
  flood_force_merge_index::Field{Int64} =
      load_field(file_handle,grid,"flood_force_merge_index",Int64)
  connect_force_merge_index::Field{Int64} =
      load_field(file_handle,grid,"connect_force_merge_index",Int64)
  flood_redirect_index::Field{Int64} =
      load_field(file_handle,grid,"flood_redirect_index",Int64)
  connect_redirect_index::Field{Int64} =
      load_field(file_handle,grid,"connect_redirect_index",Int64)
  additional_flood_redirect_index::Field{Int64} =
      load_field(file_handle,grid,"additional_flood_redirect_index",Int64)
  additional_connect_redirect_index::Field{Int64} =
      load_field(file_handle,grid,"additional_connect_redirect_index",Int64)
  corresponding_surface_cell_index::Field{Int64} =
      load_field(file_handle,surface_model_grid,"corresponding_surface_cell_index",Int64)
  return UnstructuredLakeParameters(flood_next_cell_index,
                                    connect_next_cell_index,
                                    flood_force_merge_index,
                                    connect_force_merge_index,
                                    flood_redirect_index,
                                    connect_redirect_index,
                                    additional_flood_redirect_index,
                                    additional_connect_redirect_index,
                                    corresponding_surface_cell_index)
end

function load_lake_initial_values(lake_start_filepath::AbstractString,
                                  grid::Grid,hd_grid::Grid)
  local initial_water_to_lake_centers::Field{Float64}
  local initial_spillover_to_rivers::Field{Float64}
  println("Loading: " * lake_start_filepath)
  file_handle::NcFile = NetCDF.open(lake_start_filepath)
  try
    initial_water_to_lake_centers =
      load_field(file_handle,grid,"water_redistributed_to_lakes",Float64)
    initial_spillover_to_rivers =
      load_field(file_handle,hd_grid,"water_redistributed_to_rivers",Float64)
  finally
    NetCDF.close(file_handle)
  end
  return initial_water_to_lake_centers,initial_spillover_to_rivers
end

function write_lake_numbers_field(lake_parameters::LakeParameters,lake_fields::LakeFields;
                                  timestep::Int64=-1)
  variable_name::String = "lake_field"
  filepath::String = timestep == -1 ? "/Users/thomasriddick/Documents/data/temp/lake_model_results.nc" : "/Users/thomasriddick/Documents/data/temp/lake_model_results_$(timestep).nc"
  write_field(lake_parameters.grid,variable_name,
              lake_fields.lake_numbers,filepath)
end

function write_lake_volumes_field(lake_parameters::LakeParameters,lake_volumes::Field{Float64})
  variable_name::String = "lake_field"
  filepath::String = "/Users/thomasriddick/Documents/data/temp/lake_model_out.nc"
  write_field(lake_parameters.grid,variable_name,lake_volumes,filepath)
end

function write_diagnostic_lake_volumes_field(lake_parameters::LakeParameters,
                                             diagnostic_lake_volumes::Field{Float64};
                                             timestep::Int64=-1)
  variable_name::String = "diagnostic_lake_volumes"
  filepath::String = timestep == -1 ?
                     "/Users/thomasriddick/Documents/data/temp/lake_model_volume_results.nc" :
                     "/Users/thomasriddick/Documents/data/temp/lake_model_volume_results_$(timestep).nc"
  write_field(lake_parameters.grid,variable_name,diagnostic_lake_volumes,filepath)
end

function load_drainage_fields(drainages_filename::AbstractString,grid::Grid;
                              first_timestep::Int64=1,last_timestep::Int64=1)
  println("Loading: " * drainages_filename)
  file_handle::NcFile = NetCDF.open(drainages_filename)
  drainages::Array{Field{Float64},1} = Field{Float64}[]
  try
    for i::Int64 = first_timestep:last_timestep
      drainages[i] =  load_field(file_handle,grid,"drainages",
                                 Float64;timestep=i)
    end
  finally
    NetCDF.close(file_handle)
  end
  return drainages
end

function load_runoff_fields(runoffs_filename::AbstractString,grid::Grid;
                            first_timestep::Int64=1,last_timestep::Int64=1)
  println("Loading: " * runoffs_filename)
  file_handle::NcFile = NetCDF.open(runoffs_filename)
  runoffs::Array{Field{Float64},1} = Field{Float64}[]
  try
    for i::Int64 = first_timestep:last_timestep
      runoffs[i] = load_field(file_handle,grid,"runoffs",
                              Float64;timestep=i)
    end
  finally
    NetCDF.close(file_handle)
  end
  return runoffs
end

function load_lake_evaporation_fields(lake_evaporations_filename::AbstractString,
                                      grid::LatLonGrid;first_timestep::Int64=1,
                                      last_timestep::Int64=1)
  println("Loading: " * lake_evaporations_filename)
  file_handle::NcFile = NetCDF.open(lake_evaporations_filename)
  lake_evaporations::Array{Field{Float64},1} = Field{Float64}[]
  try
    for i::Int64 = first_timestep:last_timestep
      lake_evaporations[i] = load_field(file_handle,grid,"evaporation",
                              Float64;timestep=i)
    end
  finally
    NetCDF.close(file_handle)
  end
  return lake_evaporations
end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
function seat_to_bin(seat_str::String, letters::Tuple{String,String})::String
    for (i, letter) in enumerate(letters)
        seat_str = replace(seat_str, letter => i - 1)
    end
    return seat_str
end

function parse_ticket(ticket_str::String)::Tuple{Int,Int}
    row = parse(Int, seat_to_bin(ticket_str[1:7], ("F", "B")), base = 2)
    seat = parse(Int, seat_to_bin(ticket_str[8:10], ("L", "R")), base = 2)
    return (row, seat)
end

function seat_id(ticket_str::String)::Int
    (row, seat) = parse_ticket(ticket_str)
    return 8 * row + seat
end

function solve_one()
    return maximum(map(seat_id, readlines("data/day_five.txt")))
end

function solve_two()
    ids = map(seat_id, readlines("data/day_five.txt"))
    expected_tot = (length(ids) + 1) * (maximum(ids) + minimum(ids)) / 2
    return expected_tot - sum(ids)
end

println(solve_one())
println(solve_two())

ext: jl

max_stars_repo_name: ColCarroll/adventofcode

max_issues_repo_name: ColCarroll/adventofcode

max_forks_repo_name: ColCarroll/adventofcode


text:
module CustomSch

using DagScheduler
using Dagger

DagScheduler.META_IMPL[:node] = "DagScheduler.ShmemMeta.ShmemExecutorMeta"
DagScheduler.META_IMPL[:cluster] = "DagScheduler.ShmemMeta.ShmemExecutorMeta"

_drunenv = nothing
function compute_dag(ctx, d::Dagger.Thunk)
    global _drunenv
    if _drunenv === nothing
        isdir(".mempool") && rm(".mempool"; recursive=true)
        node1 = DagScheduler.NodeEnv(1, getipaddr(), [2,3])
        _drunenv = DagScheduler.RunEnv(; nodes=[node1])
    end
    DagScheduler.rundag(_drunenv, d)
end

function cleanup(ctx)
    global _drunenv
    (_drunenv === nothing) || DagScheduler.cleanup(_drunenv)
    _drunenv = nothing
end

Dagger.PLUGIN_CONFIGS[:scheduler] = "CustomSch"

end # module CustomSch


text:
using HDF5
using CRC32c
using Test

@testset "plain" begin

# Create a new file
fn = tempname()
f = h5open(fn, "w")
@test isopen(f)
# Write scalars
f["Float64"] = 3.2
f["Int16"] = Int16(4)
# compression of empty array (issue #246)
f["compressedempty", shuffle=(), compress=4] = Int64[]
# compression of zero-dimensional array (pull request #445)
f["compressed_zerodim", shuffle=(), compress=4] = fill(Int32(42), ())
f["bloscempty", blosc=4] = Int64[]
# Create arrays of different types
A = randn(3, 5)
write(f, "Afloat64", convert(Matrix{Float64}, A))
write(f, "Afloat32", convert(Matrix{Float32}, A))
Ai = rand(1:20, 2, 4)
write(f, "Aint8", convert(Matrix{Int8}, Ai))
f["Aint16"] = convert(Matrix{Int16}, Ai)
write(f, "Aint32", convert(Matrix{Int32}, Ai))
write(f, "Aint64", convert(Matrix{Int64}, Ai))
write(f, "Auint8", convert(Matrix{UInt8}, Ai))
write(f, "Auint16", convert(Matrix{UInt16}, Ai))

# test writing multiple variable (issue #599)
write(f, "Auint32", convert(Matrix{UInt32}, Ai), "Auint64", convert(Matrix{UInt64}, Ai))

# Arrays of bools (pull request #540)
Abool = [false, true, false]
write(f, "Abool", Abool)

salut = "Hi there"
ucode = "uniçº∂e"
write(f, "salut", salut)
write(f, "ucode", ucode)
# Manually write a variable-length string (issue #187)
let
    dtype = HDF5.Datatype(HDF5.h5t_copy(HDF5.H5T_C_S1))
    HDF5.h5t_set_size(dtype.id, HDF5.H5T_VARIABLE)
    HDF5.h5t_set_cset(dtype.id, HDF5.cset(typeof(salut)))
    dspace = HDF5.dataspace(salut)
    dset = HDF5.d_create(f, "salut-vlen", dtype, dspace)
    GC.@preserve salut begin
        HDF5.h5d_write(dset, dtype, HDF5.H5S_ALL, HDF5.H5S_ALL, HDF5.H5P_DEFAULT, [pointer(salut)])
    end
end
# Arrays of strings
salut_split = ["Hi", "there"]
write(f, "salut_split", salut_split)
salut_2d = ["Hi" "there"; "Salut" "friend"]
write(f, "salut_2d", salut_2d)
# Arrays of strings as vlen
vlen = HDF5.VLen(salut_split)
d_write(f, "salut_vlen", vlen)
# Empty arrays
empty = UInt32[]
write(f, "empty", empty)
# Empty strings
empty_string = ""
write(f, "empty_string", empty_string)
# Empty array of strings
empty_string_array = String[]
write(f, "empty_string_array", empty_string_array)
# Array of empty string
empty_array_of_strings = [""]
write(f, "empty_array_of_strings", empty_array_of_strings)
# Attributes
species = [["N", "C"]; ["A", "B"]]
attrs(f)["species"] = species
C∞ = 42
attrs(f)["C∞"] = C∞
dset = f["salut"]
@test !isempty(dset)
label = "This is a string"
attrs(dset)["typeinfo"] = label
close(dset)
# Scalar reference values in attributes
attrs(f)["ref_test"] = HDF5.Reference(f, "empty_array_of_strings")
@test read(attrs(f)["ref_test"]) === HDF5.Reference(f, "empty_array_of_strings")
# Group
g = g_create(f, "mygroup")
# Test dataset with compression
R = rand(1:20, 20, 40);
g["CompressedA", chunk=(5, 6), shuffle=(), compress=9] = R
g["BloscA", chunk=(5, 6), shuffle=(), blosc=9] = R
close(g)
# Copy group containing dataset
o_copy(f, "mygroup", f, "mygroup2")
# Copy dataset
g = g_create(f, "mygroup3")
o_copy(f["mygroup/CompressedA"], g, "CompressedA")
o_copy(f["mygroup/BloscA"], g, "BloscA")
close(g)
# Writing hyperslabs
dset = d_create(f, "slab", datatype(Float64), dataspace(20, 20, 5), chunk=(5, 5, 1))
Xslab = randn(20, 20, 5)
for i = 1:5
    dset[:,:,i] = Xslab[:,:,i]
end
# More complex hyperslab and assignment with "incorrect" types (issue #34)
d = d_create(f, "slab2", datatype(Float64), ((10, 20), (100, 200)), chunk=(1, 1))
d[:,:] = 5
d[1,1] = 4
# 1d indexing
d = d_create(f, "slab3", datatype(Int), ((10,), (-1,)), chunk=(5,))
@test d[:] == zeros(Int, 10)
d[3:5] = 3:5
# Create a dataset designed to be deleted
f["deleteme"] = 17.2
close(f)
@test !isopen(f)
# Test the h5read/write interface, with attributes
W = copy(reshape(1:120, 15, 8))
Wa = Dict("a" => 1, "b" => 2)
h5write(fn, "newgroup/W", W)
h5writeattr(fn, "newgroup/W", Wa)

# Read the file back in
fr = h5open(fn)
x = read(fr, "Float64")
@test x == 3.2 && isa(x, Float64)
y = read(fr, "Int16")
@test y == 4 && isa(y, Int16)
zerodim = read(fr, "compressed_zerodim")
@test zerodim == 42 && isa(zerodim, Int32)
bloscempty = read(fr, "bloscempty")
@test bloscempty == Int64[] && isa(bloscempty, Vector{Int64})
Af32 = read(fr, "Afloat32")
@test convert(Matrix{Float32}, A) == Af32
@test eltype(Af32) == Float32
Af64 = read(fr, "Afloat64")
@test convert(Matrix{Float64}, A) == Af64
@test eltype(Af64) == Float64
@test eltype(fr["Afloat64"]) == Float64  # issue 167
Ai8 = read(fr, "Aint8")
@test Ai == Ai8
@test eltype(Ai8) == Int8
Ai16 = read(fr, "Aint16")
@test Ai == Ai16
@test eltype(Ai16) == Int16
Ai32 = read(fr, "Aint32")
@test Ai == Ai32
@test eltype(Ai32) == Int32
Ai64 = read(fr, "Aint64")
@test Ai == Ai64
@test eltype(Ai64) == Int64
Ai8 = read(fr, "Auint8")
@test Ai == Ai8
@test eltype(Ai8) == UInt8
Ai16 = read(fr, "Auint16")
@test Ai == Ai16
@test eltype(Ai16) == UInt16
Ai32 = read(fr, "Auint32")
@test Ai == Ai32
@test eltype(Ai32) == UInt32
Ai64 = read(fr, "Auint64")
@test Ai == Ai64
@test eltype(Ai64) == UInt64

Abool_read = read(fr, "Abool")
@test Abool_read == Abool
@test eltype(Abool_read) == Bool

salutr = read(fr, "salut")
@test salut == salutr
salutr = read(fr, "salut-vlen")
@test salut == salutr
ucoder = read(fr, "ucode")
@test ucode == ucoder
salut_splitr = read(fr, "salut_split")
@test salut_splitr == salut_split
salut_2dr = read(fr, "salut_2d")
@test salut_2d == salut_2dr
salut_vlenr = read(fr, "salut_vlen")
#@test salut_vlenr == salut_split
Rr = read(fr, "mygroup/CompressedA")
@test Rr == R
Rr2 = read(fr, "mygroup2/CompressedA")
@test Rr2 == R
Rr3 = read(fr, "mygroup3/CompressedA")
@test Rr3 == R
Rr4 = read(fr, "mygroup/BloscA")
@test Rr4 == R
Rr5 = read(fr, "mygroup2/BloscA")
@test Rr5 == R
Rr6 = read(fr, "mygroup3/BloscA")
@test Rr6 == R
dset = fr["mygroup/CompressedA"]
@test get_chunk(dset) == (5, 6)
@test name(dset) == "/mygroup/CompressedA"
dset2 = fr["mygroup/BloscA"]
@test get_chunk(dset2) == (5, 6)
@test name(dset2) == "/mygroup/BloscA"
Xslabr = read(fr, "slab")
@test Xslabr == Xslab
Xslabr = h5read(fn, "slab", (:, :, :))  # issue #87
@test Xslabr == Xslab
Xslab2r = read(fr, "slab2")
target = fill(5, 10, 20)
target[1] = 4
@test Xslab2r == target
dset = fr["slab3"]
@test dset[3:5] == [3:5;]
emptyr = read(fr, "empty")
@test isempty(emptyr)
empty_stringr = read(fr, "empty_string")
@test empty_stringr == empty_string
empty_string_arrayr = read(fr, "empty_string_array")
@test empty_string_arrayr == empty_string_array
empty_array_of_stringsr = read(fr, "empty_array_of_strings")
@test empty_array_of_stringsr == empty_array_of_strings
@test a_read(fr, "species") == species
@test a_read(fr, "C∞") == C∞
dset = fr["salut"]
@test a_read(dset, "typeinfo") == label
close(dset)
# Test ref-based reading
Aref = fr["Afloat64"]
sel = (2:3, 1:2:5)
Asub = Aref[sel...]
@test Asub == A[sel...]
close(Aref)
# Test iteration, name, and parent
for obj in fr
    @test filename(obj) == fn
    n = name(obj)
    p = parent(obj)
end
# Test reading multiple vars at once
z = read(fr, "Float64", "Int16")
@test z == (3.2, 4)
@test typeof(z) == Tuple{Float64,Int16}
# Test function syntax
read(fr, "Float64") do x
    @test x == 3.2
end
read(fr, "Float64", "Int16") do x, y
    @test x == 3.2
    @test y == 4
end
# Test reading entire file at once
z = read(fr)
@test z["Float64"] == 3.2
close(fr)

# Test object deletion
fr = h5open(fn, "r+")
@test haskey(fr, "deleteme")
o_delete(fr, "deleteme")
@test !haskey(fr, "deleteme")
close(fr)

# Test the h5read interface
Wr = h5read(fn, "newgroup/W")
@test Wr == W
rng = (2:3:15, 3:5)
Wr = h5read(fn, "newgroup/W", rng)
@test Wr == W[rng...]
War = h5readattr(fn, "newgroup/W")
@test War == Wa

# issue #618
# Test that invalid writes treat implicit creation as a transaction, cleaning up the partial
# operation
hid = h5open(fn, "w")
A = rand(3, 3)'
@test !haskey(hid, "A")
@test_throws ArgumentError write(hid, "A", A)
@test !haskey(hid, "A")
dset = d_create(hid, "attr", datatype(Int), dataspace(0))
@test !haskey(attrs(dset), "attr")
# broken test - writing attributes does not check that the stride is correct
@test_skip @test_throws ArgumentError write(dset, "attr", A)
@test !haskey(attrs(dset), "attr")
close(hid)

# more do syntax
h5open(fn, "w") do fid
    g_create(fid, "mygroup") do g
        write(g, "x", 3.2)
    end
end
fid = h5open(fn, "r")
@test keys(fid) == ["mygroup"]
g = fid["mygroup"]
@test keys(g) == ["x"]
close(g)
close(fid)
rm(fn)

# more do syntax: atomic rename version
tmpdir = mktempdir()
outfile = joinpath(tmpdir, "test.h5")

# create a new file
h5rewrite(outfile) do fid
    g_create(fid, "mygroup") do g
        write(g, "x", 3.3)
    end
end
@test length(readdir(tmpdir)) == 1
h5open(outfile, "r") do fid
    @test keys(fid) == ["mygroup"]
    @test keys(fid["mygroup"]) == ["x"]
end

# fail to overwrite
@test_throws ErrorException h5rewrite(outfile) do fid
    g_create(fid, "mygroup") do g
        write(g, "oops", 3.3)
    end
    error("failed")
end
@test length(readdir(tmpdir)) == 1
h5open(outfile, "r") do fid
    @test keys(fid) == ["mygroup"]
    @test keys(fid["mygroup"]) == ["x"]
end

# overwrite
h5rewrite(outfile) do fid
    g_create(fid, "mygroup") do g
        write(g, "y", 3.3)
    end
end
@test length(readdir(tmpdir)) == 1
h5open(outfile, "r") do fid
    @test keys(fid) == ["mygroup"]
    @test keys(fid["mygroup"]) == ["y"]
end
rm(tmpdir, recursive=true)

test_files = joinpath(@__DIR__, "test_files")

d = h5read(joinpath(test_files, "compound.h5"), "/data")
@test typeof(d[1]) == NamedTuple{(:wgt, :xyz, :uvw, :E), Tuple{Float64, Array{Float64, 1}, Array{Float64, 1}, Float64}}

# get-datasets
fn = tempname()
fd = h5open(fn, "w")
fd["level_0"] = [1,2,3]
grp = g_create(fd, "mygroup")
fd["mygroup/level_1"] = [4, 5]
grp2 = g_create(grp, "deep_group")
fd["mygroup/deep_group/level_2"] = [6.0, 7.0]
datasets = get_datasets(fd)
@test sort(map(name, datasets)) ==  sort(["/level_0", "/mygroup/deep_group/level_2", "/mygroup/level_1"])
close(fd)
rm(fn)

# File creation and access property lists
cpl = p_create(HDF5.H5P_FILE_CREATE)
cpl[:userblock] = 1024
apl = p_create(HDF5.H5P_FILE_ACCESS)
apl[:libver_bounds] = (HDF5.H5F_LIBVER_EARLIEST, HDF5.H5F_LIBVER_LATEST)
h5open(fn, false, true, true, true, false, cpl, apl) do fid
    write(fid, "intarray", [1, 2, 3])
end
h5open(fn, "r", libver_bounds=(HDF5.H5F_LIBVER_EARLIEST, HDF5.H5F_LIBVER_LATEST)) do fid
    intarray = read(fid, "intarray")
    @test intarray == [1, 2, 3]
end

# Test null terminated ASCII string (e.g. exported by h5py) #332
h5open(joinpath(test_files, "nullterm_ascii.h5"), "r") do fid
    str = read(fid["test"])
    @test str == "Hello World"
end

@test HDF5.unpad(UInt8[0x43, 0x43, 0x41], 1) == "CCA"

# Test the h5read/write interface with a filename as a first argument, when
# the file does not exist
rm(fn)
h5write(fn, "newgroup/W", W)
Wr = h5read(fn, "newgroup/W")
@test Wr == W
close(f)
rm(fn)

if !isempty(HDF5.libhdf5_hl)
    # Test direct chunk writing
    h5open(fn, "w") do f
      d = d_create(f, "dataset", datatype(Int), dataspace(4, 4), chunk=(2, 2))
      raw = HDF5.ChunkStorage(d)
      raw[1,1] = 0, collect(reinterpret(UInt8, [1,2,5,6]))
      raw[3,1] = 0, collect(reinterpret(UInt8, [3,4,7,8]))
      raw[1,3] = 0, collect(reinterpret(UInt8, [9,10,13,14]))
      raw[3,3] = 0, collect(reinterpret(UInt8, [11,12,15,16]))
    end

    @test h5open(fn, "r") do f
      vec(f["dataset"][:,:])
    end == collect(1:16)

    close(f)
    rm(fn)
end

# Test that switching time tracking off results in identical files
fn1 = tempname(); fn2 = tempname()
h5open(fn1, "w") do f
    f["x", track_times=false] = [1, 2, 3]
end
sleep(1)
h5open(fn2, "w") do f
    f["x", track_times=false] = [1, 2, 3]
end
@test open(crc32c, fn1) == open(crc32c, fn2)
rm(fn1); rm(fn2)

end # testset plain

@testset "complex" begin
  HDF5.enable_complex_support()

  fn = tempname()
  f = h5open(fn, "w")

  f["ComplexF64"] = 1.0 + 2.0im
  attrs(f["ComplexF64"])["ComplexInt64"] = 1im

  Acmplx = rand(ComplexF64, 3, 5)
  write(f, "Acmplx64", convert(Matrix{ComplexF64}, Acmplx))
  write(f, "Acmplx32", convert(Matrix{ComplexF32}, Acmplx))

  dset = d_create(f, "Acmplx64_hyperslab", datatype(Complex{Float64}), dataspace(Acmplx))
  for i in 1:size(Acmplx, 2)
    dset[:, i] = Acmplx[:,i]
  end

  HDF5.disable_complex_support()
  @test_throws ErrorException f["_ComplexF64"] = 1.0 + 2.0im
  @test_throws ErrorException write(f, "_Acmplx64", convert(Matrix{ComplexF64}, Acmplx))
  @test_throws ErrorException write(f, "_Acmplx32", convert(Matrix{ComplexF32}, Acmplx))
  HDF5.enable_complex_support()

  close(f)

  fr = h5open(fn)
  z = read(fr, "ComplexF64")
  @test z == 1.0 + 2.0im && isa(z, ComplexF64)
  z_attrs = attrs(fr["ComplexF64"])
  @test read(z_attrs["ComplexInt64"]) == 1im

  Acmplx32 = read(fr, "Acmplx32")
  @test convert(Matrix{ComplexF32}, Acmplx) == Acmplx32
  @test eltype(Acmplx32) == ComplexF32
  Acmplx64 = read(fr, "Acmplx64")
  @test convert(Matrix{ComplexF64}, Acmplx) == Acmplx64
  @test eltype(Acmplx64) == ComplexF64

  dset = fr["Acmplx64_hyperslab"]
  Acmplx64_hyperslab = zeros(eltype(dset), size(dset))
  for i in 1:size(dset, 2)
    Acmplx64_hyperslab[:,i] = dset[:,i]
  end
  @test convert(Matrix{ComplexF64}, Acmplx) == Acmplx64_hyperslab

  HDF5.disable_complex_support()
  z = read(fr, "ComplexF64")
  @test isa(z, NamedTuple{(:r, :i), Tuple{Float64, Float64}})

  Acmplx32 = read(fr, "Acmplx32")
  @test eltype(Acmplx32) == NamedTuple{(:r, :i), Tuple{Float32, Float32}}
  Acmplx64 = read(fr, "Acmplx64")
  @test eltype(Acmplx64) == NamedTuple{(:r, :i), Tuple{Float64, Float64}}

  close(fr)

  HDF5.enable_complex_support()
end

# test strings with null and undefined references
@testset "undefined and null" begin
fn = tempname()
f = h5open(fn, "w")

# don't silently truncate data
@test_throws ArgumentError write(f, "test", ["hello","there","\0"])
@test_throws ArgumentError write(f, "trunc1", "\0")
@test_throws ArgumentError write(f, "trunc2", "trunc\0ateme")

# test writing uninitialized string arrays
undefstrarr = similar(Vector(1:3), String) # strs = String[#undef, #undef, #undef]
@test_throws UndefRefError write(f, "undef", undefstrarr)

close(f)
rm(fn)

end # testset null and undefined

# test writing abstract arrays
@testset "abstract arrays" begin

# test writing reinterpreted data
fn = tempname()
try
    h5open(fn, "w") do f
        data = reinterpret(UInt8, [true, false, false])
        write(f, "reinterpret array", data)
    end

    @test h5open(fn, "r") do f
        read(f, "reinterpret array")
    end == UInt8[0x01, 0x00, 0x00]
finally
    rm(fn)
end

# don't silently fail for arrays with a different stride
fn = tempname()
try
    data = rand(UInt16, 2, 3);
    pdv_data = PermutedDimsArray(data, (2, 1))

    @test_throws ArgumentError h5write(fn, "pdv_data", pdv_data)
finally
    rm(fn)
end

# test alignment
fn = tempname()
h5open(fn, "w", alignment=(0, 8)) do fid
    fid["x"] = zeros(10, 10)
end

end # writing abstract arrays

# issue #705
@testset "empty and 0-size arrays" begin
fn = tempname()
hfile = h5open(fn, "w")

# Write datasets with various 0-sizes
write(hfile, "empty", HDF5.EmptyArray{Int64}()) # HDF5 empty
write(hfile, "zerodim", fill(1.0π))   # 0-dimensional
write(hfile, "zerovec", zeros(0))     # 1-dimensional, size 0
write(hfile, "zeromat", zeros(0, 0))  # 2-dimensional, size 0
write(hfile, "zeromat2", zeros(0, 1)) # 2-dimensional, size 0 with non-zero axis
dempty = hfile["empty"]
dzerodim = hfile["zerodim"]
dzerovec = hfile["zerovec"]
dzeromat = hfile["zeromat"]
dzeromat2 = hfile["zeromat2"]

# Test that eltype is preserved (especially for EmptyArray)
@test eltype(dempty) == Int64
@test eltype(dzerodim) == Float64
@test eltype(dzerovec) == Float64
@test eltype(dzeromat) == Float64
@test eltype(dzeromat2) == Float64
# Test sizes are as expected
@test size(dempty) == ()
@test size(dzerovec) == (0,)
@test size(dzeromat) == (0, 0)
@test size(dzeromat2) == (0, 1)
@test HDF5.isnull(dempty)
@test !HDF5.isnull(dzerovec)
@test !HDF5.isnull(dzeromat)
@test !HDF5.isnull(dzeromat2)
# Reading back must preserve emptiness
@test read(dempty) isa HDF5.EmptyArray
# but 0-dimensional Array{T,0} are stored as HDF5 scalar
@test size(dzerodim) == ()
@test !HDF5.isnull(dzerodim)
@test read(dzerodim) == 1.0π

# Similar tests for writing to attributes
write(dempty, "attr", HDF5.EmptyArray{Float64}())
write(dzerodim, "attr", fill(1.0ℯ))
write(dzerovec, "attr", zeros(Int64, 0))
write(dzeromat, "attr", zeros(Int64, 0, 0))
write(dzeromat2, "attr", zeros(Int64, 0, 1))
aempty = dempty["attr"]
azerodim = dzerodim["attr"]
azerovec = dzerovec["attr"]
azeromat = dzeromat["attr"]
azeromat2 = dzeromat2["attr"]
# Test that eltype is preserved (especially for EmptyArray)
@test eltype(aempty) == Float64
@test eltype(azerodim) == Float64
@test eltype(azerovec) == Int64
@test eltype(azeromat) == Int64
@test eltype(azeromat2) == Int64
# Test sizes are as expected
@test size(aempty) == ()
@test size(azerovec) == (0,)
@test size(azeromat) == (0, 0)
@test size(azeromat2) == (0, 1)
@test HDF5.isnull(aempty)
@test !HDF5.isnull(azerovec)
@test !HDF5.isnull(azeromat)
@test !HDF5.isnull(azeromat2)
# Reading back must preserve emptiness
@test read(aempty) isa HDF5.EmptyArray
# but 0-dimensional Array{T,0} are stored as HDF5 scalar
@test size(azerodim) == ()
@test !HDF5.isnull(azerodim)
@test read(azerodim) == 1.0ℯ

# Concatenation of EmptyArrays is not supported
x = HDF5.EmptyArray{Float64}()
@test_throws ErrorException [x x]
@test_throws ErrorException [x; x]
@test_throws ErrorException [x x; x x]

close(hfile)
rm(fn)

# check that printing EmptyArray doesn't error
buf = IOBuffer()
show(buf, HDF5.EmptyArray{Int64}())
@test String(take!(buf)) == "HDF5.EmptyArray{Int64}()"
show(buf, MIME"text/plain"(), HDF5.EmptyArray{Int64}())
@test String(take!(buf)) == "HDF5.EmptyArray{Int64}()"
end # empty and 0-size arrays

@testset "generic read of native types" begin
fn = tempname()
hfile = h5open(fn, "w")

dtype_varstring = HDF5.Datatype(HDF5.h5t_copy(HDF5.H5T_C_S1))
HDF5.h5t_set_size(dtype_varstring, HDF5.H5T_VARIABLE)

write(hfile, "uint8_array", UInt8[(1:8)...])
write(hfile, "bool_scalar", true)

fixstring = "fix"
varstring = "var"
write(hfile, "fixed_string", fixstring)
vardset = d_create(hfile, "variable_string", dtype_varstring, dataspace(varstring))
GC.@preserve varstring begin
    HDF5.h5d_write(vardset, dtype_varstring, HDF5.H5S_ALL, HDF5.H5S_ALL, HDF5.H5P_DEFAULT, [pointer(varstring)])
end
flush(hfile)
close(dtype_varstring)

# generic read() handles concrete types with definite sizes transparently
d = read(hfile["uint8_array"], UInt8)
@test d isa Vector{UInt8}
@test d == 1:8
d = read(hfile["bool_scalar"], Bool)
@test d isa Bool
@test d == true
d = read(hfile["fixed_string"], HDF5.FixedString{length(fixstring),0})
@test d isa String
@test d == fixstring
d = read(hfile["variable_string"], Cstring)
@test d isa String
@test d == varstring
# will also accept memory-compatible reinterpretations
d = read(hfile["uint8_array"], Int8)
@test d isa Vector{Int8}
@test d == 1:8
d = read(hfile["bool_scalar"], UInt8)
@test d isa UInt8
@test d == 0x1
# but should throw on non-compatible types
@test_throws ErrorException("""
                            Type size mismatch
                            sizeof(UInt16) = 2
                            sizeof($(sprint(show, datatype(UInt8)))) = 1
                            """) read(hfile["uint8_array"], UInt16)

# Strings are not fixed size, but generic read still handles them if given the correct
# underlying FixedString or Cstring type; a method overload makes String work, too.
d = read(hfile["fixed_string"], String)
@test d isa String
@test d == fixstring
d = read(hfile["variable_string"], String)
@test d isa String
@test d == varstring

close(hfile)
rm(fn)
end # generic read of native types

@testset "show" begin
fn = tempname()

# First create data objects and sure they print useful outputs

hfile = h5open(fn, "w")
@test sprint(show, hfile) == "HDF5 data file: $fn"

group = g_create(hfile, "group")
@test sprint(show, group) == "HDF5 group: /group (file: $fn)"

dset = d_create(group, "dset", datatype(Int), dataspace((1,)))
@test sprint(show, dset) == "HDF5 dataset: /group/dset (file: $fn xfer_mode: 0)"

meta = a_create(dset, "meta", datatype(Bool), dataspace((1,)))
@test sprint(show, meta) == "HDF5 attribute: meta"

prop = p_create(HDF5.H5P_DATASET_CREATE)
@test sprint(show, prop) == "HDF5 property: dataset create class"

dtype = HDF5.Datatype(HDF5.h5t_copy(HDF5.H5T_IEEE_F64LE))
@test sprint(show, dtype) == "HDF5 datatype: H5T_IEEE_F64LE"

dspace = dataspace((1,))
@test occursin(r"^HDF5.Dataspace\(\d+\)", sprint(show, dspace))

# Now test printing after closing each object

close(dspace)
@test sprint(show, dspace) == "HDF5.Dataspace(-1)"

close(dtype)
@test sprint(show, dtype) == "HDF5 datatype: (invalid)"

close(prop)
@test sprint(show, prop) == "HDF5 property (invalid)"

close(meta)
@test sprint(show, meta) == "HDF5 attribute (invalid)"

close(dset)
@test sprint(show, dset) == "HDF5 dataset (invalid)"

close(group)
@test sprint(show, group) == "HDF5 group (invalid)"

close(hfile)
@test sprint(show, hfile) == "HDF5 data file (closed): $fn"

rm(fn)
end # show tests

@testset "split1" begin

@test HDF5.split1("a") == ("a", nothing)
@test HDF5.split1("/a/b/c") == ("/", "a/b/c")
@test HDF5.split1("a/b/c") == ("a", "b/c")
@test HDF5.split1(GenericString("a")) == ("a", nothing)
@test HDF5.split1(GenericString("/a/b/c")) == ("/", "a/b/c")
@test HDF5.split1(GenericString("a/b/c")) == ("a", "b/c")

# The following two paths have the same graphemes but different code unit structures:
# the first one is
#     <latin small letter a with circumflex> "/" <greek small leter alpha>
# while the second one is
#     "a" <combining circumflex accent> "/" <greek small letter alpha>
circa = "â" # <latin small leter a with circumflex>
acomb = "â" # "a" + <combining circumflex accent>
path1 = circa * "/α"
path2 = acomb * "/α"
# Sanity checks that the two strings are different but equivalent under normalization
@test path1 != path2
@test Base.Unicode.normalize(path1, :NFC) == Base.Unicode.normalize(path2, :NFC)
# Check split1 operates correctly
@test HDF5.split1(path1) == (circa, "α")
@test HDF5.split1(path2) == (acomb, "α")
@test HDF5.split1("/" * path1) == ("/", path1)
@test HDF5.split1("/" * path2) == ("/", path2)

end # split1 tests

# Also tests AbstractString interface
@testset "haskey" begin
fn = tempname()
hfile = h5open(fn, "w")

group1 = g_create(hfile, "group1")
group2 = g_create(group1, "group2")

@test haskey(hfile, GenericString("group1"))
@test !haskey(hfile, GenericString("groupna"))
@test haskey(hfile, "group1/group2")
@test !haskey(hfile, "group1/groupna")

dset1 = d_create(hfile, "dset1", datatype(Int), dataspace((1,)))
dset2 = d_create(group1, "dset2", datatype(Int), dataspace((1,)))

@test haskey(hfile, "dset1")
@test !haskey(hfile, "dsetna")
@test haskey(hfile, "group1/dset2")
@test !haskey(hfile, "group1/dsetna")

meta1 = a_create(dset1, "meta1", datatype(Bool), dataspace((1,)))
@test haskey(dset1, "meta1")
@test !haskey(dset1, "metana")

attribs = attrs(hfile)
attribs["test1"] = true
attribs["test2"] = "foo"

haskey(attribs, "test1")
haskey(attribs, "test2")
!haskey(attribs, "testna")

attribs = attrs(dset2)
attribs["attr"] = "foo"
haskey(attribs, GenericString("attr"))

close(hfile)
rm(fn)
end # haskey tests

@testset "AbstractString" begin

fn = GenericString(tempname())
hfile = h5open(fn, "w")
close(hfile)
hfile = h5open(fn); close(hfile)
hfile = h5open(fn, "w")

@test_nowarn g_create(hfile, GenericString("group1"))
@test_nowarn d_create(hfile, GenericString("dset1"), datatype(Int), dataspace((1,)))
@test_nowarn d_create(hfile, GenericString("dset2"), 1)

@test_nowarn hfile[GenericString("group1")]
@test_nowarn hfile[GenericString("dset1")]

dset1 = hfile["dset1"]
@test_nowarn a_create(dset1, GenericString("meta1"), datatype(Bool), dataspace((1,)))
@test_nowarn a_create(dset1, GenericString("meta2"), 1)
@test_nowarn dset1[GenericString("meta1")]
@test_nowarn dset1[GenericString("x")] = 2

array_of_strings = ["test",]
write(hfile, "array_of_strings", array_of_strings)
@test_nowarn attrs(hfile)[GenericString("ref_test")] = HDF5.Reference(hfile, GenericString("array_of_strings"))
@test read(attrs(hfile)[GenericString("ref_test")]) === HDF5.Reference(hfile, "array_of_strings")

hfile[GenericString("test")] = 17.2
@test_nowarn o_delete(hfile, GenericString("test"))
@test_nowarn a_delete(dset1, GenericString("meta1"))

# transient types
memtype_id = HDF5.h5t_copy(HDF5.H5T_NATIVE_DOUBLE)
dt = HDF5.Datatype(memtype_id)
@test !HDF5.h5t_committed(dt)
t_commit(hfile, GenericString("dt"), dt)
@test HDF5.h5t_committed(dt)

dt = datatype(Int)
ds = dataspace(0)
d = d_create(hfile, GenericString("d"), dt, ds)
g = g_create(hfile, GenericString("g"))
a = a_create(hfile, GenericString("a"), dt, ds)

for obj in (d, g)
   @test_nowarn a_write(obj, GenericString("a"), 1)
   @test_nowarn a_read(obj, GenericString("a"))
   @test_nowarn write(obj, GenericString("aa"), 1)
   @test_nowarn attrs(obj)["attr1"] = GenericString("b")
end
@test_nowarn write(d, "attr2", GenericString("c"))
@test_nowarn d_write(g, GenericString("ag"), GenericString("gg"))
@test_nowarn d_write(g, GenericString("ag_array"), [GenericString("a1"), GenericString("a2")])

genstrs = GenericString["fee", "fi", "foo"]
@test_nowarn a_write(d, GenericString("myattr"), genstrs)
@test genstrs == read(d["myattr"])

for obj in (hfile,)
    @test_nowarn d_open(obj, GenericString("d"))
    @test_nowarn d_write(obj, GenericString("dd"), 1)
    @test_nowarn d_read(obj, GenericString("dd"))
    @test_nowarn read(obj, GenericString("dd"))
    @test_nowarn read(obj, GenericString("dd")=>Int)
end
read(attrs(hfile), GenericString("a"))

write(hfile, GenericString("ASD"), GenericString("Aa"))
write(g, GenericString("ASD"), GenericString("Aa"))
write(g, GenericString("ASD1"), [GenericString("Aa")])

# test writing multiple variable
@test_nowarn write(hfile, GenericString("a1"), rand(2,2), GenericString("a2"), rand(2,2))

# copy methods
d1 = d_create(hfile, GenericString("d1"), dt, ds)
d1["x"] = 32
@test_nowarn o_copy(hfile, GenericString("d1"), hfile, GenericString("d1copy1"))
@test_nowarn o_copy(d1, hfile, GenericString("d1copy2"))

fn = GenericString(tempname())
A = Matrix(reshape(1:120, 15, 8))
@test_nowarn h5write(fn, GenericString("A"), A)
@test_nowarn h5read(fn, GenericString("A"))
@test_nowarn h5read(fn, GenericString("A"), (2:3:15, 3:5))

@test_nowarn h5write(fn, GenericString("x"), 1)
@test_nowarn h5read(fn, GenericString("x") => Int)

@test_nowarn h5rewrite(fn) do fid
    g_create(fid, "mygroup") do g
        write(g, "x", 3.3)
    end
end
@test_nowarn h5rewrite(fn) do fid
    g_create(fid, "mygroup") do g
        write(g, "y", 3.3)
    end
end

@test_nowarn h5write(fn, "W", [1 2; 3 4])
@test_nowarn h5writeattr(fn, GenericString("W"), Dict("a" => 1, "b" => 2))
@test_nowarn h5readattr(fn, GenericString("W"))

fn_external = GenericString(tempname())
dset = d_create_external(hfile, "ext", fn_external, Int, (10,20))

end

# length for FixedString
fix = HDF5.FixedString{4,0}((b"test"...,))
@test length(fix) == 4


text:
import KernelDensity: UnivariateKDE, default_bandwidth, kde
import Distributions: Distribution

"""
    bin(x::AbstractUncertainIndexValueDataset, binning::BinnedResampling{UncertainScalarKDE}) -> UncertainIndexValueDataset
    bin(x::AbstractUncertainIndexValueDataset, binning::BinnedResampling{UncertainScalarPopulation}) -> UncertainIndexValueDataset

Resample every element of `x` the number of times given by `binning.n`. After resampling,
distribute the values according to their indices, into the bins given by `binning.left_bin_edges`.

## Returns

Returns an `UncertainIndexValueDataset`. Indices are assumed to be uniformly distributed within each
bin, and are represented as `CertainValue`s at the bin centers. Values of the dataset have different
representations depending on what `binning` is:

- If `binning isa BinnedResampling{UncertainScalarKDE}`, then values in each bin are represented by a
    kernel density estimate to the distribution of the resampled values whose resampled indices
    fall in that bin.
- If `binning isa BinnedResampling{UncertainScalarPopulation}`, then values in each bin are
    represented by equiprobable populations consisting of the resampled values whose resampled
    indices fall in the bins.
"""
function bin(x::AbstractUncertainIndexValueDataset, binning::BinnedResampling); end

"""
    bin(x::AbstractUncertainIndexValueDataset, binning::BinnedResampling{RawValues}) -> Tuple(Vector, Vector{Vector})

Resample every element of `x` the number of times given by `binning.n`. After resampling,
distribute the values according to their indices, into the `N` bins given by `binning.left_bin_edges`.

## Returns

Return a tuple containing the `N` different bin centers and a `N`-length vector of
resampled values whose resampled indices fall in the `N` different bins.

## Example

```julia
# Some example data with unevenly spaced time indices
npts = 300
time, vals = sort(rand(1:1000, npts)), rand(npts)

# Add uncertainties to indices and values, and represent as
# UncertainIndexValueDataset
utime = [UncertainValue(Normal, t, 10) for t in time]
uvals = [UncertainValue(Normal, v, 0.1) for v in vals]

udata = UncertainIndexValueDataset(utime, uvals)

# Bin data into fall in 25 time step wide time bins ranging
# from time indices 100 to 900 and return a vector of raw
# values for each bin. Do this by resampling each uncertain
# data point 10000 times and distributing those draws among
# the bins.
left_bin_edges = 100:25:900
n_draws = 10000
binning = BinnedResampling(RawValues, left_bin_edges, n_draws)

bin_centers, bin_draws = bin(udata, binning)
```
"""
function bin(x::AbstractUncertainIndexValueDataset, binning::BinnedResampling{RawValues};
        nan_threshold = 0.0)
    # Pre-allocate some arrays into which we resample the values of the
    # index and value populations.
    idxs = fill(NaN, binning.n)
    vals = fill(NaN, binning.n)

    perminds = zeros(Int, binning.n)
    sorted_idxs = fill(NaN, binning.n)
    sorted_vals = fill(NaN, binning.n)

    bin_edges = binning.left_bin_edges
    n_bins = length(bin_edges) - 1

    # Used to compute the index of the bin into which a draw belongs
    mini = minimum(binning.left_bin_edges)
    s = step(binning.left_bin_edges)

    # Empty vectors that will contain draws.
    binvecs = [Vector{Float64}(undef, 0) for i = 1:n_bins]
    #[sizehint!(bv, resampling.n*n_bins) for bv in binvecs]

    @inbounds for (j, (idx, val)) in enumerate(zip(x.indices, x.values))
        # Resample the j-th index and j-th value
        resample!(idxs, idx)
        resample!(vals, val)

        # Get the vector that sorts the index vector, and use that to
        # sort the draws.
        sortperm!(perminds, idxs)
        sorted_idxs .= idxs[perminds]
        sorted_vals .= vals[perminds]

        # The vectors above are sorted sorted, so this can be done faster
        for i in 1:n_bins
            inbin = findall(bin_edges[i] .<= sorted_idxs .<= bin_edges[i+1])
            if length(inbin) > nan_threshold
                append!(binvecs[i], sorted_vals[inbin])
            end
        end
    end

    bin_centers = bin_edges[1:end-1] .+ step(bin_edges)/2

    return bin_centers, binvecs
end

function bin(x::AbstractUncertainIndexValueDataset, binning::BinnedResampling{UncertainScalarPopulation};
        nan_threshold = 0)

    # Get bin center and a vector for each bin containing the values falling in that bin.
    left_bin_edges = binning.left_bin_edges
    n = binning.n
    n_bins = length(left_bin_edges) - 1

    bin_centers, binvecs = bin(x, BinnedResampling(RawValues, left_bin_edges, n))

    # Estimate distributions in each bin by kernel density estimation
    estimated_value_dists = Vector{Union{CertainValue, UncertainScalarPopulation}}(undef, n_bins)
    binvec_lengths = length.(binvecs)

    for i in 1:n_bins
        L = binvec_lengths[i]
        # If bin contains enough values, represent as population. Otherwise,
        # set to NaN.
        if L > nan_threshold
            probs = Weights(repeat([1/L], L))
            estimated_value_dists[i] = UncertainScalarPopulation(binvecs[i], probs)
        else
            estimated_value_dists[i] = UncertainValue(NaN)
        end
    end

    new_inds = UncertainIndexDataset(bin_centers)
    new_vals = UncertainValueDataset(estimated_value_dists)

    UncertainIndexValueDataset(new_inds, new_vals)
end

function bin(x::AbstractUncertainIndexValueDataset, binning::BinnedResampling{UncertainScalarKDE};
        nan_threshold = 0,
        kernel::Type{D} = Normal,
        bw_factor = 4,
        npoints::Int = 2048) where {K <: UnivariateKDE, D <: Distribution}

    # Get bin center and a vector for each bin containing the values falling in that bin.
    left_bin_edges = binning.left_bin_edges
    n = binning.n
    n_bins = length(left_bin_edges) - 1

    bin_centers, binvecs = bin(x, BinnedResampling(RawValues, left_bin_edges, n))

    # Estimate distributions in each bin by kernel density estimation
    estimated_value_dists = Vector{Union{CertainValue, UncertainScalarKDE}}(undef, n_bins)
    binvec_lengths = length.(binvecs)

    for i in 1:n_bins
        L = binvec_lengths[i]
        # If bin contains enough values, represent as KDE estimate. Otherwise,
        # set to NaN.
        if L > nan_threshold
            bw = default_bandwidth(binvecs[i]) / bw_factor
            # Kernel density estimation
            KDE = kde(binvecs[i], npoints = npoints, kernel = kernel, bandwidth = bw)

            # Get the x value for which the density is estimated.
            xrange = KDE.x

            # Normalise estimated density
            density = KDE.density ./ sum(KDE.density)

            estimated_value_dists[i] = UncertainScalarKDE(KDE, binvecs[i], xrange, Weights(density))
        else
            estimated_value_dists[i] = UncertainValue(NaN)
        end
    end

    new_inds = UncertainIndexDataset(bin_centers)
    new_vals = UncertainValueDataset(estimated_value_dists)

    UncertainIndexValueDataset(new_inds, new_vals)
end

export bin


text:
using Parameters

"""
Configurable Parameters for the equality saturation process.
"""
@with_kw mutable struct SaturationParams
    timeout::Int = 7
    # default sizeout. TODO make this bytes
    # sizeout::Int = 2^14
    matchlimit::Int = 5000
    eclasslimit::Int = 5000
    stopwhen::Function = ()->false
    scheduler::Type{<:AbstractScheduler} = BackoffScheduler
    schedulerparams::Tuple=()
end


text:
# This file is a part of LegendUpROOTIO.jl, licensed under the MIT License (MIT).

module LegendUpROOTIO

using ArgCheck
using ArraysOfArrays
using ElasticArrays
using LegendDataTypes
using RadiationDetectorSignals
using StaticArrays
using StructArrays
using Tables
using Unitful
using UpROOT

import TypedTables

using LegendDataTypes: readdata, writedata, getunits, setunits!, units_from_string, units_to_string

using RadiationDetectorSignals: RealQuantity, ArrayOfDims, AosAOfDims, SArrayOfDims,
    recursive_ndims

using UpROOT: ROOTIOBuffer, ROOTClassHeader, ROOTClassVersion, TObjectContent, TNamedContent,
    TClonesArrayHeader, kByteCountMask

include("mgdo.jl")

end # module


text:
# WARNING: This should not be used as-is.  It's unfinished experimentation.

using Transformations
using StochasticOptimization
using PenaltyFunctions

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

#=
See http://www.breloff.com/DeepRL-OnlineGAE/ for details/explanation of OnlineGAE

Samples from a distribution D(ϕ), where ϕ contains all sufficient statistics.
In the case of a MultivariateNormal N(μ, Σ):
    Σ = U'U   (we output U to ensure Σ is positive definite)
    ϕ := vcat(μ, vec(U))
    ϕ(s,θ) is a learnable transformation from states to sufficient statistics

We then sample from our multivariate distribution:
    z ~ D(ϕ)
and (deterministically) map that sample to actions (i.e. project squashed samples `z` into the action space `A`):
    a = π(z)    or    a ~ π(s,θ)

---

Notes:

The gradient wrt transformation params θ can be broken into the components of ϕ which map to μ/U.

Assuming the first and last steps are deterministic:
    P(a|z,s,θ) = P(a|z) * P(z|ϕ) * P(ϕ|s,θ)
               = 1 * P(z|ϕ) * 1
               = P(z|ϕ)
So:
    ∇log P(a|z,s,θ) = ∇log P(z|ϕ)
=#

"""
Maps state to value: V(s)

We update the learnable parameters using gradient vector δ
"""
mutable struct ValueCritic{T,TRANS<:Learnable} <: Learnable
    trans::TRANS  # nS --> 1  transformation which outputs a value V(s) for state s
    γ::T            # discount
    lastv::T       # V(s)
    δ::T          # TD(0) delta: δ = r + γV(s′) - V(s)
    # lastδ::T
end

function ValueCritic(::Type{T}, trans::Learnable, γ::T) where T
    ValueCritic{T,typeof(trans)}(trans, γ, zero(T), zero(T))
end

# function transform!(critic::ValueCritic, s::AbstractArray)
#     # critic.lastv = output_value(critic.trans)[1]
#     transform!(critic.trans, s)
# end

# give reward r, compute output grad: δ = r + γV(s′) - V(s)
# then backprop to get ∇θ
function grad!(critic::ValueCritic{T}, r::Number) where T
    Vs′ = output_value(critic.trans)[1]
    Vs = critic.lastv
    # critic.lastδ = critic.δ

    # the loss function is L2 loss with "truth" (r+λVₛ′) and "estimate" Vₛ
    # the output gradient is:
    #   ∂(δ²/2)/∂Vₛ′

    # this is the discounted return δ:
    critic.δ = r + critic.γ * Vs′ - Vs
    output_grad(critic.trans)[1] = -critic.δ
    # output_grad(critic.trans)[1] = -critic.γ * critic.δ

    # # this tries to solve for the average future reward
    # critic.δ = critic.γ * r + (one(T) - critic.γ) * Vs′ - Vs
    # output_grad(critic.trans)[1] = -critic.γ * critic.δ

    # critic.lastv = Vs′
    grad!(critic.trans)
end

"""
Online Generalized Advantage Estimation for Actor-Critic Reinforcement Learning

Transforms states (input) to actions (output) using learnable parameters θ

We assume the general form of mapping states (s) to actions (a):
    s --> ϕ(s,θ) --> D(ϕ) --> a

ϕ is:
    - a learnable transformation with learnable parameters θ
    - the sufficient statistics of distribution D
    - some concatenated combination of μ/U/σ for multivariate normals
"""
mutable struct OnlineGAE{T      <: Number,
                         ASET   <: AbstractSet,
                         PHI    <: Learnable,
                         DIST   <: MvNormalTransformation,
                         CRITIC <: ValueCritic,
                       # P      <: Params,
                         PEN    <: Penalty,
                         AL     <: LearningStrategy,
                         CL     <: LearningStrategy
                        } <: Learnable
    A::ASET           # the action space
    ϕ::PHI            # learnable transformation to output sufficient statistics of D
    D::DIST           # generative transformation for sampling inputs to a
    critic::CRITIC    # the critic... wraps a learnable estimating value function V(s)
    γ::T              # the discount for the critic
    λ::T              # the extra discount for the actor
    ϵ::Vector{T}      # eligibility traces for the learnable params θ in transformation ϕ
    # a::Vector{T}      # the action
    # t::Int            # current timestep
    # ∇logP::Vector{T}  # policy gradient: ∇log P(a | s) == ∇log P(z | ϕ)
    # lastr::T          # most recent return
    # params::P         # the combined parameters from the actor transformation ϕ and the critic transformation
    penalty::PEN      # a penalty to add to param gradients
    actor_learner::AL
    critic_learner::CL
end

function OnlineGAE(A::AbstractSet,
                   ϕ::Learnable,
                   D::MvNormalTransformation,
                   critic_trans::Learnable,
                   γ::T,
                   λ::T,
                   actor_learner::LearningStrategy,
                   critic_learner::LearningStrategy;
                   penalty::Penalty = NoPenalty()) where T
    # connect transformations, init the critic
    link_nodes!(ϕ, D)
    critic = ValueCritic(T, critic_trans, γ)
    np = params_length(ϕ)
    ϵ = zeros(T, np)
    # ∇logP = zeros(T, np)
    # params = consolidate_params(T, ϕ, critic_trans)
    setup!(actor_learner, ϕ)
    setup!(critic_learner, critic.trans)
    OnlineGAE(A, ϕ, D, critic, γ, λ, ϵ, penalty, actor_learner, critic_learner)
end

# don't do anything here... we'll update later
LearnBase.update!(π::OnlineGAE, ::Void) = return

function Reinforce.reset!(π::OnlineGAE{T}) where T
    fill!(π.ϵ, zero(T))
    # π.critic.lastv = 0
    # setup!(π.actor_learner, π.ϕ)
    # setup!(π.critic_learner, π.critic.trans)
end

function Reinforce.action(π::OnlineGAE, r, s′, A′)
    # sample z ~ N(μ,Σ) which is determined by ϕ
    transform!(π.ϕ, s′)
    z = transform!(π.D)

    # project our squashed sample onto into the action space to get our actions
    # a = (â --> [lo,hi])
    a = A′.lo .+ logistic.(z) .* (A′.hi .- A′.lo)
    if !(a in A′)
        warn("a=$a not in A=$(A′)")
        a = rand(A′)
    end
    a
end

function learn!(π::OnlineGAE, s, a, r, s′)
    # update the critic.  we use the current model to get the lastv == Vₛ
    # as well as the current == Vₛ′
    π.critic.lastv = transform!(π.critic.trans, s)[1]
    transform!(π.critic.trans, s′)
    grad!(π.critic, r)
    # π.lastr = r
    t = π.critic.trans
    addgrad!(grad(t), π.penalty, params(t))

    learn!(t, π.critic_learner, nothing)

    #=
    update the actor using the OnlineGAE formulas:
        ϵₜ = (γλ)ϵₜ₋₁ + ∇
        ĝₜ = δₜϵₜ

    note: we use the grad-log-prob from the last timestep, since we
    can't update until we compute the critic's δ, which depends
    on the next timestep
    =#

    # update the grad-log-prob of distribution D, and store that for the next timestep
    # NOTE: grad(π.ϕ) now contains the grad-log-prob of this timestep... but we don't use this
    #   until the next timestep
    grad!(π.D)
    grad!(π.ϕ)

    # we use last timestep's ∇logP to update the eligibility trace of the last timestep ϵ
    γλ = π.γ * π.λ
    ϵ = π.ϵ
    ∇ = grad(π.ϕ)
    addgrad!(∇, π.penalty, params(π.ϕ))
    for i=1:length(ϵ)
        ϵ[i] = γλ * ϵ[i] + ∇[i]
    end

    # copy!(π.∇logP, grad(π.ϕ))

    # overwrite the gradient estimate: ĝ = δϵ
    δ = π.critic.δ
    for i=1:length(ϵ)
        ∇[i] = -δ * ϵ[i]
    end

    # # add the penalty to the gradient
    # addgrad!(∇, π.penalty, params(π.ϕ))

    # learn the actor
    learn!(π.ϕ, π.actor_learner, nothing)
end


text:
using GaussDCA

#module ReWeight # from GaussDCA

const PACKBITS = 64
const _u = 0x0084210842108421
const _alt = 0x007c1f07c1f07c1f
const packfactor = div(PACKBITS, 5)
const packrest = PACKBITS % 5
const _msk = ~(UInt64(0));

clength(l::Int) = div(l-1, packfactor) + 1
crest(l::Int)   = (l-1) % packfactor + 1
nnz_aux(x::UInt64) = ((x | (x >>> 1) | (x >>> 2) | (x >>> 3) | (x >>> 4)) & _u)
nz_aux(nx::UInt64) = (nx & (nx >>> 1) & (nx >>> 2) & (nx >>> 3) & (nx >>> 4) & _u)
nz_aux2(nx::UInt64, s) = nz_aux(nx) & (_msk >>> s)

function compress_Z(Z::Matrix{Int8})
    N, M = size(Z)
    ZZ = Vector{Int8}[Z[:,i] for i = 1:M]
    cl = clength(N)
    cZ = [zeros(UInt64, cl) for i=1:M]
    @inbounds for i = 1:M
        cZi = cZ[i]
        ZZi = ZZ[i]
        for k = 1:N
            k0 = div(k-1, packfactor)+1
            k1 = (k-1) % packfactor
            cZi[k0] |= (UInt64(ZZi[k]) << (5*k1))
        end
    end
    return cZ
end

function compute_weights(Z::Matrix{Int8}, q::Integer, theta::Real)
    @assert q ≤ 32
    N, M = size(Z)
    cZ = compress_Z(Z)
    return compute_weights(cZ, theta, N, M)
end

function compute_weights(cZ::Vector{Vector{UInt64}}, theta::Real, N::Int, M::Int)
    theta = Float64(theta)
    cl = clength(N)
    kmax = div(cl - 1, 31)
    rmax = (cl - 1) % 31 + 1
    Meff = 0.0
    W = ones(M)
    thresh = floor(theta * N)
    println("theta = $theta threshold = $thresh")
    if theta == 0
        println("M = $M N = $N Meff = $M")
        return W, Float64(M)
    end
    for i = 1:M-1
        cZi = unsafe(cZ[i])
        for j = i+1:M
            cZj = unsafe(cZ[j])
            czi = start(cZi)
            czj = start(cZj)
            dist::UInt64 = 0
            z::UInt64 = 0
            for k = 1:kmax
                z = 0
                for r = 1:31
                    zi, czi = next(cZi, czi)
                    zj, czj = next(cZj, czj)
                    y = zi ⊻ zj
                    z += nnz_aux(y)
                end
                t = @collapse(z)
                dist += t
                dist >= thresh && break
            end
            if dist < thresh
                z = 0
                for r = 1:rmax
                    zi, czi = next(cZi, czi)
                    zj, czj = next(cZj, czj)

                    y = zi ⊻ zj
                    z += nnz_aux(y)
                end
                t = @collapse(z)
                dist += t
            end
            if dist < thresh
                W[i] += 1
                W[j] += 1
            end
        end
    end
    for i = 1:M
        W[i] = 1 / W[i]
    end
    Meff = sum(W)
    println("M = $M N = $N Meff = $Meff")
    return W, Meff
end

function ReadFasta(filename::AbstractString,
                   max_gap_fraction::Real = 0.9,
                   theta::Any = :auto,
                   remove_dups::Bool=true)

    Z = GaussDCA.read_fasta_alignment(filename, max_gap_fraction)
    if remove_dups
        Z, _ = GaussDCA.remove_duplicate_seqs(Z)
    end

    N, M = size(Z)
    q = round(Int,maximum(Z))

    q > 32 && error("parameter q=$q is too big (max 31 is allowed)")
    W, Meff = GaussDCA.compute_weights(Z, q, theta)

    rmul!(W, 1.0/Meff)
    Zint = round.(Int,Z)
    return W, Zint,N,M,q
end

filename = "/home/cossio/Desktop/julia/test.fasta"
max_gap_fraction = 0.9
theta = :auto
remove_dups = true

Z = GaussDCA.read_fasta_alignment(filename, max_gap_fraction)
if remove_dups
    Z, _ = GaussDCA.remove_duplicate_seqs(Z)
end

N, M = size(Z)
q = round(Int,maximum(Z))

q > 32 && error("parameter q=$q is too big (max 31 is allowed)")
W, Meff = GaussDCA.compute_weights(Z, q, 0.1)

W ./= Meff

#rmul!(W, 1.0/Meff)
Zint=round.(Int,Z)

cZ = compress_Z(Z)

x = GaussDCA.ReadFastaAlignment.letter2num('A')

ReadFasta("/home/cossio/Desktop/julia/test.fasta")


text:
using DoubleFloats
using Random

using SemiLagrangian:
    Advection,
    AdvectionData,
    getstcoef,
    advection!,
    sizeall,
    getext,
    getdata,
    getcur_t,
    getindsplit,
    _getcurrentindice,
    getprecal,
    getitr,
    gett_split,
    nextstate!,
    UniformMesh,
    totuple,
    tovector,
    Lagrange,
    getpoissonvar,
    compute_ke,
    getinterp,
    points,
    dotprod

function initmesh(t_deb, t_end, t_size)
    t_step = (t_end - t_deb) ./ t_size
    return totuple(UniformMesh.(t_deb, t_end, t_size)), t_step
end

function test_adv(T::DataType)

    println("trace1")
    t_debsp = T.([-1, -10, -3])
    t_endsp = T.([3, 6, 5])
    t_szsp = (16, 8, 32)
    t_debv = T.([-3, -9, 1])
    t_endv = T.([1, 7, 1])
    t_szv = (4, 8, 4)
    base_dt = one(T) / 80

    Nsum = 6

    t_meshsp, t_stepsp = initmesh(t_debsp, t_endsp, t_szsp)
    t_meshv, t_stepv = initmesh(t_debv, t_endv, t_szv)

    interp = Lagrange(3, T)
    adv = Advection(
        (t_meshsp...,t_meshv...),
        map(x -> Lagrange(3, T), 1:6),
        base_dt,
        [
            ([1, 2, 3, 6, 5, 4], 1, 1, true),
            ([2, 1, 3, 4, 6, 5], 1, 1, true),
            ([3, 2, 1, 4, 5, 6], 1, 1, true), # space states
            ([4, 5, 6, 1, 2, 3], 1, 2, true),
            ([5, 4, 6, 1, 2, 3], 1, 2, true),
            ([6, 5, 4, 1, 2, 3], 1, 2, true), # velocity states
        ]
    )
    println("trace2")

    sref = (t_szsp..., t_szv...)
    @test sref == sizeall(adv)

    tab = rand(T, sizeall(adv))

    advd = AdvectionData(adv, tab, getpoissonvar(adv))

    @test compute_ke(t_meshsp, t_meshv, tab) == compute_ke(advd)

    @test advd.state_gen == 1

    @test advd.parext == getext(advd)
    @test advd.data == getdata(advd)
    @test adv.tab_coef[1] == getcur_t(adv, 1)
    @test adv.tab_coef[3] == getcur_t(adv, 1)
    @test adv.tab_coef[1] == getcur_t(adv, 2)
    @test adv.tab_coef[1] == getcur_t(adv, 3)
    @test adv.tab_coef[2] == getcur_t(adv, 4)
    @test adv.tab_coef[2] == getcur_t(adv, 5)
    @test adv.tab_coef[2] == getcur_t(adv, 6)

    @test getcur_t(advd) == advd.adv.tab_coef[1]

    println("trace3")

    t_coef = [1, 1, 1, 2, 2, 2, 3, 3, 3, 1]
#    t_dim = [1, 2, 3, 1, 2, 3, 1, 2, 3, 1]
    t_indice = [1, 2, 3, 4, 5, 6, 1, 2, 3, 1]
    t_v = [false, false, false, true, true, true, false, false, false, false]
    t_result = [true, true, true, true, true, true, true, true, false, true]
    # resfirst = [4, 4, 4, (4, 1, 1), (4, 1, 1), (4, 1, 1), 4, 4, 4, 4]
    # ressecond = [
    #     (:, 3, 1, 4, 1, 1),
    #     (3, :, 1, 1, 4, 1),
    #     (3, 1, :, 1, 1, 4),
    #     (4, 1, 1, :, 3, 1),
    #     (4, 1, 1, 3, :, 1),
    #     (4, 1, 1, 3, 1, :),
    #     (:, 3, 1, 4, 1, 1),
    #     (3, :, 1, 1, 4, 1),
    #     (3, 1, :, 1, 1, 4),
    #     (:, 3, 1, 4, 1, 1),
    # ]

    for i = 1:length(t_coef)
        @test getstcoef(advd) == t_coef[i]
        @test advd.state_gen == (i-1)%9 + 1
        @test t_indice[i] == _getcurrentindice(advd)

        @test getcur_t(advd) == adv.tab_coef[t_coef[i]]
        t = adv.t_interp
        @test t[t_indice[i]] == getinterp(advd)[1]

        # itrfirst = getitrfirst(advd)
        # (res, _) = Iterators.peel(Iterators.drop(itrfirst,3))

        # @test resfirst[i] == res

        # itrsecond = getitrsecond(advd,res)

        # (res2, _) = Iterators.peel(Iterators.drop(itrsecond,2))

        # @test ressecond[i] == res2

#        x = t_indice[i]
        #    @time @test addcolon.(x, Iterators.product(refitr[vcat(1:(x-1),(x+1):Nsum)]...)) == getitr(advd)

        ret = nextstate!(advd)
        @test ret == t_result[i]
    end
end

# @testset "test Advection1d Float" begin

#     test_adv(Float64)

# end
# @testset "test Advection1d BigFloat" begin

#     test_adv(BigFloat)

# end
# @testset "test Advection1d Double64" begin

#     test_adv(Double64)

# end
function test_ke(T::DataType)
    t_debsp = T.([-1 // 1, -10 // 1, -3 // 1])
    t_endsp = T.([3 // 1, 6 // 1, 5 // 1])
    t_szsp = [2, 4, 8]
    t_stepsp = (t_endsp - t_debsp) ./ t_szsp
    tt_meshsp = UniformMesh.(t_debsp, t_endsp, t_szsp)
    t_meshsp = totuple(tt_meshsp)
    szsp = totuple(t_szsp)

    t_debv = T.([-3 // 1, -9 // 1, 1 // 1, -1 // 1])
    t_endv = T.([1 // 1, 7 // 1, 5 // 1, 3 // 1])
    t_szv = [4, 8, 4, 2]
    t_stepv = (t_endv - t_debv) ./ t_szv
    tt_meshv = UniformMesh.(t_debv, t_endv, t_szv)
    t_meshv = totuple(tt_meshv)
    szv = totuple(t_szv)

    fxv = if T <: Rational
        rationalize.(BigInt, rand(Float64, (szsp..., szv...)))
    else
        rand(T, (szsp..., szv...))
    end

    Nsp = length(szsp)
    Nv = length(szv)
    Nsum = Nsp + Nv
    dx = prod(step, t_meshsp)
    dv = prod(step, t_meshv)
    sum_sp = Array{T,Nv}(undef, szv)
    sum_sp .= reshape(sum(fxv, dims = ntuple(x -> x, Nsp)), szv)
    refres = (dx * dv) * sum(dotprod(points.(t_meshv)) .^ 2 .* sum_sp)

    @test refres == compute_ke(t_meshsp, t_meshv, fxv)
end

# function test_itr(T::DataType)
#     t_debsp = T.([-1, -10, -3])
#     t_endsp = T.([3, 6, 5])
#     t_szsp = (16, 8, 32)
#     t_debv = T.([-3, -9, 1])
#     t_endv = T.([1, 7, 1])
#     t_szv = (4, 8, 4)
#     base_dt = one(T) / 80

#     Nsum = 6

#     t_meshsp, t_stepsp = initmesh(t_debsp, t_endsp, t_szsp)
#     t_meshv, t_stepv = initmesh(t_debv, t_endv, t_szv)

#     interp = Lagrange(3, T)
#     adv = Advection1d(
#         t_meshsp,
#         t_meshv,
#         ntuple(x -> Lagrange(3, T), 3),
#         ntuple(x -> Lagrange(3, T), 3),
#         base_dt,
#     )

#     sref = (t_szsp..., t_szv...)

#     refitr = ntuple(x -> 1:sref[x], size(sref, 1))

#     tab = rand(T, sizeall(adv))

#     @show size(tab)

#     advd = Advection1dData(adv, tab, getpoissonvar(adv))

#     #    @show advd.t_itrfirst

#     resfirst = [
#         CartesianIndex(2, 17, 1, 2, 1),
#         CartesianIndex(2, 9, 3, 1, 1),
#         CartesianIndex(2, 1, 2, 3, 1),
#         CartesianIndex(2, 1, 5, 3, 1),
#         CartesianIndex(2, 1, 9, 5, 1),
#         CartesianIndex(2, 1, 5, 3, 1),
#         CartesianIndex(2, 17, 1, 2, 1),
#         CartesianIndex(2, 9, 3, 1, 1),
#         CartesianIndex(2, 1, 2, 3, 1),
#         CartesianIndex(2, 17, 1, 2, 1),
#     ]
#     ressecond = [
#         CartesianIndex(6, 21, 1, 3, 1),
#         CartesianIndex(6, 11, 1, 2, 1),
#         CartesianIndex(6, 5, 2, 5, 1),
#         CartesianIndex(6, 1, 6, 5, 1),
#         CartesianIndex(2, 2, 11, 1, 2),
#         CartesianIndex(6, 1, 6, 5, 1),
#         CartesianIndex(6, 21, 1, 3, 1),
#         CartesianIndex(6, 11, 1, 2, 1),
#         CartesianIndex(6, 5, 2, 5, 1),
#         CartesianIndex(6, 21, 1, 3, 1),
#     ]

#     for i = 1:length(resfirst)
#         itr = getitr(advd)

#         (res, _) = Iterators.peel(Iterators.drop(itr, 1153))

#         #        @show res

#         @test resfirst[i] == res

#         (res2, _) = Iterators.peel(Iterators.drop(itr, 2213))

#         #        @show res2

#         @test ressecond[i] == res2

#         nextstate!(advd)
#     end

# end

@testset "test compute_ke" begin
    test_ke(Rational{BigInt})
    test_ke(Float64)
    test_ke(BigFloat)
end

@testset "test advection" begin
    test_adv(Float64)
end
# @testset "test itr" begin
#     test_itr(Float64)
#     test_itr(BigFloat)
# end


text:
using Plots, Plots.Measures

"""
    plot_lmiss_n_hwe(list)
---
Plot the loci missing results to determine thresholds to elimate SNP.
Note:
- I found it is better to plot accumulate distribution of allele missing rate.
  - Hence I totally rewrote the function
- This must be done after plink stats missing.
- Removing low quality SNP seems should be the first step
  - as there are many SNP, e.g., thousands.
- Low quality ID can be done after this procedure.
"""
function plot_lmiss_n_hwe(list)
    cd(work_dir)
    empty_dir("tmp")

    sdir = "data/genotypes/step-2.plk"
    tdir = "notebooks/fig"
    isdir(tdir) || mkdir(tdir)

    title("Plotting l-missing and Hwe stats")
    for pf in list
        item("Platform $pf")
        n = "t"                 # dummy, to read N_id. i don't parse it into Int here.
        y = Float64[]
        message("Missing data stats")
        miss_allele_stats("$sdir/$pf", "tmp/$pf")

        open("tmp/$pf.lmiss", "r") do io # plot lmiss
            _ = readline(io)
            for line in eachline(io)
                n, m = [split(line)[i] for i in [4, 5]]
                push!(y, parse(Float64, m))
            end
            sort!(y)
        end

        z = Float64[]
        message("HWE stats")
        hwe_stats("$sdir/$pf", "tmp/$pf")
        open("tmp/$pf.hwe", "r") do io # plot hwe
            _ = readline(io)
            for line in eachline(io)
                push!(z, parse(Float64, split(line)[9]))
            end
            sort!(z, rev = true)
        end

        message("Plot the figures")
        nlc = length(y)
        x = range(0, 100., length = nlc)
        h = 0.1                 # the threshold
        p1 = plot(x, y,
                  label = "$pf l-miss",
                  ylabel = "Missing freq., $n ID",
                  xlabel = "Accumulate SNP%",
                  width = 2,
                  bottom_margin=2mm)
        plot!(p1, [h], seriestype="hline", width=.5, label="Missing frq. = $h")
        v = round(count(x->x<h, y) / nlc * 100, digits = 2)
        plot!(p1, [v], seriestype="vline", width=.5, label="Resv. rate = $v%")

        t = -log10.(z)
        h = 0.0001              # the threshold
        p2 = plot(x, t,
                  label = "$pf HWE",
                  ylabel = "-log10 HWE P-value",
                  xlabel = "Accumulate SNP%, $nlc loci",
                  width = 2)
        plot!(p2, [-log10(h)], seriestype="hline", width=.5, label="P(hwe) = $h")
        v = round(count(x->x>h, z) /nlc *100, digits = 2)
        plot!(p2, [v], seriestype="vline", width=.5, label="Resv. rate = $v%")

        plot(p1, p2, size=(800,300), dpi=300)

        savefig("$tdir/$pf-ms-hwe.png")
        done()
    end
end

"""
    filter_lowQ_snp(list; geno::Float64=0.1, maf::Float64=0.05, hwe::Float64=0.0001)
---
Remove low quality SNP in each platform by country.

SNP of
- missing ratio > 0.1
- MAF < 0.01
- P-value of HWE statistics <0.0001

are removed.

Note this is a keyword argument example. Call example:

`filter_lowQ_snp(maf = 0.02`
"""
function filter_lowQ_snp(list;
                         geno::Float64 = 0.1,
                         maf::Float64 = 0.01,
                         hwe::Float64 = 0.0001)
    cd(work_dir)
    title("Remove low quality SNP")
    sdir = "data/genotypes/step-2.plk"
    tdir = "data/genotypes/step-3.plk"
    isdir(tdir) || mkdir(tdir)

    for pf in list
        item("Filtering platform $pf")
        plink_filter_snp("$sdir/$pf", geno, maf, hwe, "$tdir/$pf")
        done()
    end
end

"""
    plot_imiss(list, threshold::Float64 = 0.1)
---
Plot missing data per ID to decide a threshold to filter away low quality ID.
"""
function plot_imiss(list,
                    threshold::Float64 = 0.1)
    title("Plot allele missing per ID")
    sdir = "data/genotypes/step-3.plk"
    tdir = "notebooks/fig"
    isdir(tdir) || mkdir(tdir)
    empty_dir("tmp")

    for pf in list
        item("Dealing with platform $pf")
        miss_allele_stats("$sdir/$pf", "tmp/out")
        y = Float64[]
        open("tmp/out.imiss", "r") do io
            _ = readline(io)
            for line in eachline(io)
                r = parse(Float64, split(line)[6])
                push!(y, r)
            end
        end
        sort!(y)
        tid = length(y)
        nid = count(x->x<threshold, y)

        plot(1:tid, y, label=pf, width=2, dpi=300)
        plot!([nid], seriestype="vline", width=.5, label="$tid reduced to $nid")
        plot!([threshold], seriestype="hline", width=.5, label="Threshold = $threshold")
        savefig("$tdir/$pf-imiss.png")
        done()
    end
end

"""
    filter_id(list, mind::Float64 = 0.05)
---
Remove ID with missing alleles over `mind`
"""
function filter_id(list, mind::Float64 = 0.05)
    cd(work_dir)
    title("Remove ID with missing genotypes over $mind")
    sdir = "data/genotypes/step-3.plk"
    tdir = "data/genotypes/step-4.plk"
    isdir(tdir) || mkdir(tdir)
    for pf in list
        item("Dealing with platform $pf")
        plink_filter_id("$sdir/$pf", mind, "$tdir/$pf")
        done()
    end
end

"""
    remove_duplicates(list)
---
Remove duplicates SNP, i.e., different names, same chromosome and bp, according
to `data/maps/duplicated.snp`.
- Merge genotypes with the non-missing values.
- Use SNP names in the 2nd column

# Note
ID `XXXYYYF100000007336`, which was labeled so, was removed at this stage.  According to my judgement,
its DNA was sampled from `XXXYYYF100000006921`.
"""
function remove_duplicates(list)
    cd(work_dir)
    sdir = joinpath(work_dir, "data/genotypes/step-4.plk")
    tdir = joinpath(work_dir, "data/genotypes/step-5.plk")
    empty_dir(tdir)
    excl = joinpath(work_dir, "data/maps/excl.snp")
    aref = joinpath(work_dir, "data/maps/ref.allele")
    isdir(tdir) || mkdir(tdir)

    title("Remove duplicated SNP")
    for pf in list
        item("Platform $pf")
        _ = read(`plink --cow
                        --bfile $sdir/$pf
                        --exclude $excl
                        --reference-allele $aref
                        --make-bed
                        --out $tdir/$pf`,
                 String)
        done()
    end
    cd(tdir)
    _ = read(`plink
		--cow
		--bfile dutch-d1
		--remove ../../pedigree/dd1-excl.id
		--make-bed
		--out tmp`,
             String)
    for ext in ["bed", "bim", "fam"]
        mv("tmp.$ext", "dutch-d1.$ext", force=true)
    end
end

#=
"""
    unify_ref_allele(list)
---
With the list of `bim` files given, unify `bed`s with the reference alleles appeared last.
The benefit of this is the homozygote genotypes won't flip later for `G` matrix calculation.
"""
function unify_ref_allele(list)
    title("Unify reference allele across platforms")
    fra = joinpath(work_dir, "data/genotypes/step-5.plk")
    til = joinpath(work_dir, "data/genotypes/uniref.plk")
    isdir(til) || mkdir(til)
    empty_dir(til)
    ref = Dict()
    for l in list
        for line in eachline(joinpath(fra, "$l.bim"))
            snp, r = [split(line)[i] for i in [2, 5]]
            push!(ref, snp=>r)
        end
    end
    open(joinpath(til, "ref.allele"), "w") do io
        for (snp, allele) in ref
            write(io, "$snp $allele\n")
        end
    end
    for l in list
        print("$l ")
        force_reference_allele(joinpath(fra, l), joinpath(til, "ref.allele"), joinpath(til, l))
    end
    done()
end
=#