text:
# functions for populations
"""
   pop = generate_population(par::PTParameters, nm::Int, nf::Int)
   pop = generate_population(par::PTParameters, qmsimfile, nm::Int, nf::Int; gfile)

Generates a new population with `nm` males and `nf` females.
For a historical population, you should provide the number of males and females.
For a new population, you do not need the numbers because the animals will be
generated in a course of simulation.

When using genotypes, you have to provide external genotypes, generated
by QMSim. The data should be converted to an HDF5 file by the QMSimData
package. The QMSim file must be needed for any populations because the file
defines the genetic architecture. Also, you should provide `gfile` to store
gennotypes that going to be generated in simulation.
```
# hostorical population with 5 males and 10 females
hp = generate_population(par, nm=5, nf=10)

# hostorical population with genotypes
hp = generate_population(par, "qmsimdata.h5", nm=5, nf=10)

# new population with no animals
pop = generate_population(par)

# new population with genotypes
pop = generate_population(par, "qmsimdata.h5", nm=0, nf=0, gfile="pop.h5")
```
"""
function generate_population(par::PTParameters; nm::Int=0, nf::Int=0, map::Union{Nothing,QMSimMap}=nothing, gfile="")
   # check parameters
   check_parameters(par)

   # initialization
   maxAnimal = nm+nf
   hp = ifelse(maxAnimal>0, true, false)
   maxGroup = 0
   maxCG = 0

   # generate polygenic breeding values
   pbv = zeros(Float64,maxAnimal)
   if par.h2_poly > 0.0
      var_poly = par.var_p * par.h2_poly
      pbv .= randn(maxAnimal)*sqrt(var_poly)
   end

   # breeding value due to QTL
   qbv = zeros(Float64,maxAnimal)

   # total breeding value
   tbv = zeros(Float64,maxAnimal)
   tbv .= pbv .+ qbv

   # genotyping status
   if isnothing(map)
      genotyped = fill(false, maxAnimal)
   else
      genotyped = fill(true, maxAnimal)
   end

   df = DataFrame(
      id = collect(1:maxAnimal),
      male = vcat([true for i=1:nm],[false for i=1:nf]),
      year = zeros(Int, maxAnimal),
      alive = fill(true, maxAnimal),
      pregnant = fill(false, maxAnimal),
      genotyped = genotyped,
      candidate = fill(false, maxAnimal),
      status = zeros(Int, maxAnimal),
      sire = zeros(Int, maxAnimal),
      dam = zeros(Int, maxAnimal),
      siregroup = zeros(Int, maxAnimal),
      damgroup = zeros(Int, maxAnimal),
      nprog = zeros(Int, maxAnimal),
      nrecprog = zeros(Int, maxAnimal),
      inb = zeros(Float64,maxAnimal),
      pbv = pbv,
      qbv = qbv,
      tbv = tbv,
      ebv = zeros(Float64,maxAnimal),
      gebv = Vector{Union{Missing,Float64}}(missing,maxAnimal),
      rel = Vector{Union{Missing,Float64}}(missing,maxAnimal),
      nrec = zeros(Int, maxAnimal),
      firsty = Vector{Union{Missing,Float64}}(missing,maxAnimal),
      lasty = Vector{Union{Missing,Float64}}(missing,maxAnimal),
      avgy = Vector{Union{Missing,Float64}}(missing,maxAnimal))

   # generate empty entry
   animal = Vector{PTAnimal}()
   @inbounds for id in 1:maxAnimal
      cg = zeros(Int,0)
      y = zeros(Float64,0)
      pe = zeros(Float64,0)
      e = zeros(Float64,0)
      push!(animal, PTAnimal(cg,y,pe,e))
   end

   return PTPopulation(par,hp,maxAnimal,maxGroup,maxCG,df,animal,map,gfile)
end

function generate_population(par::PTParameters, qmsimfile::String; nm::Int=0, nf::Int=0, map=nothing, gfile="")
   map = read_qmsim_map_hdf5(qmsimfile)
   if nm>0 || nf>0
      if gfile==qmsimfile || gfile==""
         @info "hostorical population: qmsimfile is the same as gfile, which will be read-only."
         pop = generate_population(par, nm=nm, nf=nf, map=map, gfile=qmsimfile)
      else
         throw(ArgumentError("For a historical population, qmsimfile must be the same as gfile."))
      end
   else
      if gfile==qmsimfile
         throw(ArgumentError("gfile must be provided as a different file as qmsimfile."))
      end
      # genomic data file created by QMSim/QMSimFiles
      if gfile==""
         h5file = "snp_"*randstring(4)*".h5"
      else
         h5file = gfile
      end
      pop = generate_population(par, nm=nm, nf=nf, map=map, gfile=h5file)
      create_qmsim_map_hdf5(map, pop.gfile)
   end
   return pop
end

"""
    add_new_animal!(pop::PTPopulation)

Add an animal to population `pop`.
The added animals lose the pedigree information.
"""
function add_new_animal!(pop::PTPopulation, male::Bool, animal::PTAnimal;
                     genotyped=false, candidate=false, year=0, alive=true, pregnant=false, status=0, sire=0, dam=0, siregroup=0, damgroup=0, inb=0.0,
                     pbv=0.0, qbv=0.0, tbv=0.0, ebv=0.0, rel=missing, gebv=missing)
   pop.maxAnimal = pop.maxAnimal + 1
   (nrec,firsty,lasty,avgy) = get_phenotypic_information(animal.y)
   newdata = (pop.maxAnimal, male, year, alive, pregnant, genotyped, candidate, status, sire, dam, siregroup, damgroup, 0, 0, inb, pbv, qbv, tbv, ebv, gebv, rel, nrec, firsty, lasty, avgy)
   push!(pop.df, newdata)
   push!(pop.animal, animal)
   return nothing
end

function add_new_animal!(pop::PTPopulation, row::DataFrameRow{DataFrames.DataFrame,DataFrames.Index}, animal::PTAnimal)
   pop.maxAnimal = pop.maxAnimal + 1
   (nrec,firsty,lasty,avgy) = get_phenotypic_information(animal.y)
   newdata = (id=pop.maxAnimal, male=row.male, year=row.year, alive=row.alive, pregnant=row.pregnant, genotyped=row.genotyped,
              candidate=row.candidate, status=row.status, sire=row.sire, dam=row.dam, siregroup=row.siregroup, damgroup=row.damgroup, nprog=0, nrecprog=0,
              inb=row.inb, pbv=row.pbv, qbv=row.qbv, tbv=row.tbv, ebv=row.ebv, gebv=row.gebv, rel=row.rel,
              nrec=nrec, firsty=firsty, lasty=lasty, avgy=avgy)
   push!(pop.df, newdata)
   push!(pop.animal, animal)
   return nothing
end

function get_phenotypic_information(y::Vector{Float64})
   nrec = length(y)
   if nrec>0
      firsty = y[1]
      lasty = y[end]
      avgy = mean(y)
   else
      firsty = missing
      lasty = missing
      avgy = missing
   end
   return (nrec,firsty,lasty,avgy)
end

"""
    newidlist = migrate_from_hp!(hp::PTPopulation, pop::PTPopulation, idlist::Vector{Int}; year=0, halt=true)

Migrate some animals from a historical population `hp` to another population `pop`.
The migrated animals have new IDs in the new population.
The animals will have "dead" status in the historical population.
It returns a new ID list `newidlist` in the new population.
With `halt=true` (default), this function throws when there are no enough individuals in the historical population.
"""
function migrate_from_hp!(hp::PTPopulation, pop::PTPopulation, idlist::Vector{Int}; year::Union{Int,Vector{Int}}=0, halt::Bool=true)
   unique_idlist = unique(idlist)
   n = length(unique_idlist)
   if n==0 && halt
      error("empty ID list; maybe there is no more individual in the historical population.")
   end
   newidlist = Vector{Int}()
   @inbounds for id in unique_idlist
      # check id in the historiacl population
      if id > hp.maxAnimal
         error("No such ID $(id) in the first population")
      end
      if hp.df[id,:alive]
         # new animal for the new population
         add_new_animal!(pop, hp.df[id,:], hp.animal[id])
         pop.df[pop.maxAnimal,:id] = pop.maxAnimal
         # moved
         hp.df[id,:alive] = false
         # assigned new id
         push!(newidlist, pop.maxAnimal)
      else
         if halt
            error("id $(id) in hp is not alive; maybe there is no more individual in the historical population.")
         else
            @warn "id $(id) in hp is not alive, and has not been added to the population."
         end
      end

      # (hp) id to (pop) pop.maxAnimal
      # assigned new genotypes
      if hp.df[id,:genotyped]
         # copy genotypes from the database
         g = read_qmsim_individual_hdf5(hp.map,hp.gfile,id)
         add_qmsim_individual_hdf5(pop.map,pop.gfile,g)
         pop.df[pop.maxAnimal,:genotyped] = true
         pop.df[pop.maxAnimal,:qbv] = g.tbv
         pop.df[pop.maxAnimal,:tbv] = pop.df[pop.maxAnimal,:pbv] + pop.df[pop.maxAnimal,:qbv]
      end
   end
   # assign year
   assign_year!(pop, newidlist, year)
   return newidlist
end

"""
    assign_year!(pop::PTPopulation, idlist::Vector{Int}, year::Vector{Int})
    assign_year!(pop::PTPopulation, idlist::Vector{Int}, year::Int)

Assign birth year to animals defied by `idlist`.
By default, the "years" will be evenly assigned to the animals.
"""
function assign_year!(pop::PTPopulation, idlist::Vector{Int}, year::Vector{Int})
   m = length(year)
   j = 1
   @inbounds for i in idlist
      if 0<i && i<=pop.maxAnimal
         pop.df[i,:year] = year[j]
         j = j + 1
         if j>m
            j = 1
         end
      end
   end
   return nothing
end
function assign_year!(pop::PTPopulation, idlist::Vector{Int}, year::Int)
   assign_year!(pop, idlist, [year])
   return nothing
end

# new group with founder males and females
function generate_group(pop::PTPopulation; sires::Vector{Int}=Int[], dams::Vector{Int}=Int[], aliveonly::Bool=true, empty::Bool=false)
   if empty
      pop.maxGroup = pop.maxGroup + 1
      return PTGroup(pop, pop.maxGroup, 0, 0, 0, Int[], Int[], Int[], Int[])
   else
      return _generate_group(pop, sires=sires, dams=dams, aliveonly=aliveonly)
   end
end

function _generate_group(pop::PTPopulation; sires::Vector{Int}=Int[], dams::Vector{Int}=Int[], aliveonly::Bool=true)
   check_idlist(sires)
   check_idlist(dams)

   if aliveonly
      sirelist = sires[map(x->pop.df[x,:alive],sires)]
      damlist = dams[map(x->pop.df[x,:alive],dams)]
   else
      sirelist = sires
      damlist = dams
   end

   nsires = length(sirelist)
   ndams = length(damlist)
   hassire = ifelse(nsires>0, true, false)
   hasdam = ifelse(ndams>0, true, false)
   if !hassire && !hasdam
      throw(ArgumentError("no sire and dam lists; use generate_group(pop,empty=true) for an empty group"))
   end

   n = nsires + ndams
   pop.maxGroup = pop.maxGroup + 1
   idlist = [sirelist; damlist]
   sort!(idlist)
   generation = zeros(Int,n)
   return PTGroup(pop, pop.maxGroup, n, nsires, ndams, sirelist, damlist, idlist, generation)
end

function check_idlist(idlist::Vector{Int})
   n = length(idlist)
   if n>0 && n!=length(unique(idlist))
      throw(ArgumentError("possible duplicated animals in the list"))
   end
   if sum(idlist .< 1)>0
      throw(ArgumentError("zero or negative ID(s) in the list"))
   end
   return nothing
end

"""
    group = copy_group(orig::PTGroup)

Make a deep copy of the griginal group `orig`.
A new group code (gid) will be assigned.
"""
function copy_group(orig::PTGroup)
   pop = orig.pop
   pop.maxGroup = pop.maxGroup + 1
   n = length(orig.sires) + length(orig.dams)
   sires = copy(orig.sires)
   dams = copy(orig.dams)
   id = copy(orig.id)
   generation = copy(orig.generation)
   return PTGroup(pop, pop.maxGroup, n, orig.maxSire, orig.maxDam, sires, dams, id, generation)
end

# a general function to select IDs for a filter function
"""
    id = selectid(fun, pop, ...)
    id = selectid(fun, group, ...)
    id = selectid(fun, groups, ...)

Returns an array of IDs that meet the condition defined by a function `f` in a population, a group, or an array of groups.
The query function `f` should be reasonable for the dataframe stored in `pop`.
Some options are available.

- `aliveonly`: select only individuals alive (default=`true`).
- `limit`: maximum length of output (default=`nothing` i.e., all IDs)
- `idlist`: select IDs only in this list (default=`Int[]` i.e., no constrant).
- `allowempty`: allow an empty ID list (default=`true`); for `false`, the function throws an error if `length(id)==0`.
- `sortby`: sort IDs by items in the dataframe or particular condition
    - `:id`: no sorting
    - `:random`: randomized array
    - any comparable quantities (e.g., `:year`, `:tbv`, `:pbv`, `:qbv`, `:ebv`, `:gebv`, etc.)
- `rev`: revserse sort (default=`true`, i.e., descending "high to low" order)

When an array of groups is supplied, `aliveonly`, `idlist`, and `allowempty` will be aplied to each group; `limit`, `sortby` and `rev` will be applied only to the final output.

```julia
# males in the population
selectid(:male => x -> x==true, pop)

# culled males with unknown sire (coded by 0)
selectid([:male,:sire] => (x,y) -> x==true && y==0, group, aliveonly=false)

# top males with highest EBV in the group
selectid([:male] => x -> x==true, group, sortby=:ebv)
```
"""
function selectid(fun, pop::PTPopulation; limit::Union{Nothing,Int}=nothing, idlist::Vector{Int}=Int[], aliveonly::Bool=true,
   sortby::Symbol=:id, rev::Bool=true, allowempty::Bool=true)

   # symbols
   symb = get_symbol_array(fun)
   if !in(:alive,symb); push!(symb,:alive); end
   if !in(:id,symb);    push!(symb,:id);    end

   selected_id = _selectid(fun, symb, pop.df, idlist=idlist, aliveonly=aliveonly, sortby=sortby, rev=rev, allowempty=allowempty)

   leng = valid_length(selected_id, limit)
   if leng==length(selected_id)
      return selected_id
   else
      return selected_id[1:leng]
   end
end

# for groups
function selectid(fun, group::PTGroup; limit::Union{Nothing,Int}=nothing, idlist::Vector{Int}=Int[], aliveonly::Bool=true,
   sortby::Symbol=:id, rev::Bool=true, allowempty::Bool=true)

   # symbols
   symb = get_symbol_array(fun)
   if !in(:alive,symb); push!(symb,:alive); end
   if !in(:id,symb);    push!(symb,:id);    end

   if isempty(idlist)
      selected_id = _selectid(fun, symb, group.pop.df, idlist=group.id, aliveonly=aliveonly, sortby=sortby, rev=rev, allowempty=allowempty)
   else
      unique_id = sort(unique(intersect(idlist,group.id)))
      selected_id = _selectid(fun, symb, group.pop.df, idlist=unique_id, aliveonly=aliveonly, sortby=sortby, rev=rev, allowempty=allowempty)
   end

   leng = valid_length(selected_id, limit)
   if leng==length(selected_id)
      return selected_id
   else
      return selected_id[1:leng]
   end
end

function selectid(fun, groups::Vector{PTGroup}; limit::Union{Nothing,Int}=nothing, idlist::Vector{Int}=Int[], aliveonly::Bool=true,
   sortby::Symbol=:id, rev::Bool=true, allowempty::Bool=true)

   if length(groups)==0
      return Int[]
   end

   # symbols
   symb = get_symbol_array(fun)
   if !in(:alive,symb); push!(symb,:alive); end
   if !in(:id,symb);    push!(symb,:id);    end

   selected_id = Int[]
   for group in groups
      if isempty(idlist)
         selected_id = [selected_id; _selectid(fun, symb, group.pop.df, idlist=group.id, aliveonly=aliveonly, allowempty=allowempty)]
      else
         unique_id = sort(unique(intersect(idlist,group.id)))
         selected_id = [selected_id; _selectid(fun, symb, group.pop.df, idlist=unique_id, aliveonly=aliveonly, allowempty=allowempty)]
      end
   end

   selected_id = unique(selected_id)

   # sort
   if sortby==:id
      ret_id = selected_id
   elseif sortby==:random
      ret_id = shuffle(selected_id)
   else
      df = groups[1].pop.df
      perm = sortperm(df[selected_id,sortby],rev=rev)
      ret_id = selected_id[perm]
   end
   # check
   if !allowempty && isempty(ret_id)
      error("empty ID list generated by selectid")
   end
   # length check
   leng = valid_length(ret_id, limit)
   if leng==length(ret_id)
      return ret_id
   else
      return ret_id[1:leng]
   end
end

function valid_length(idlist::Vector{Int},limit::Union{Nothing,Int})
   listlen = length(idlist)
   if isnothing(limit)
      return listlen
   else
      if limit<0
         throw(ArgumentError("invalid limit < 0"))
      else
         return min(limit,listlen)
      end
   end
end

function selectid_general(fun, group::PTGroup; idlist::Vector{Int}=Int[], aliveonly::Bool=true, sortby::Symbol=:id, rev::Bool=true, allowempty::Bool=true)
   # symbols
   symb = get_symbol_array(fun)
   if !in(:alive,symb); push!(symb,:alive); end
   if !in(:id,symb);    push!(symb,:id);    end

   # including :generation
   if in(:generation,symb)
      if length(idlist)>0
         unique_id = sort(unique(intersect(idlist,group.id)))
         df = pop.df[unique_id,symb]
         df[:,:generation] = group.generation[unique_id]
      else
         df = pop.df[group.id,symb]
         df[:,:generation] = group.generation
      end
   else
      if length(idlist)>0
         unique_id = sort(unique(intersect(idlist,group.id)))
         df = group.pop.df[unique_id,symb]
      else
         df = group.pop.df[group.id,symb]
      end
   end

   return _selectid(fun, symb, df, aliveonly=aliveonly, sortby=sortby, rev=rev, allowempty=allowempty)
end

function _selectid(fun, symb::Vector{Symbol}, df::DataFrame; idlist::Vector{Int}=Int[], aliveonly::Bool=true, sortby::Symbol=:id, rev::Bool=true, allowempty::Bool=true)
   # filtering
   if aliveonly
      if length(idlist)>0
         alive_id = unique(idlist[findall(df[idlist,:alive])])
      else
         alive_id = df[df[!,:alive],:id]
      end
      selected_id = collect( filter(fun, df[alive_id,symb],view=true)[!,:id] )
   else
      if length(idlist)>0
         unique_id = unique(idlist)
         selected_id = collect( filter(fun, df[unique_id,symb],view=true)[!,:id] )
      else
         selected_id = collect( filter(fun, df[:,symb],view=true)[!,:id] )
      end
   end
   # sort
   if sortby==:id
      ret_id = selected_id
   elseif sortby==:random
      ret_id = shuffle(selected_id)
   else
      perm = sortperm(df[selected_id,sortby],rev=rev)
      ret_id = selected_id[perm]
   end
   # check
   if !allowempty && isempty(ret_id)
      error("empty ID list generated by selectid")
   end
   return ret_id
end

function get_symbol_array(fun)
   if typeof(fun[1])==Symbol
      symb = [fun[1]]
   else
      symb = copy(fun[1])
   end
   return symb
end

# random sampling
function random_sampling(pop::PTPopulation, n::Int; male::Bool=false, female::Bool=false, aliveonly::Bool=true, allowempty::Bool=true)
   if male && female
      idlist = selectid([:id] => x->x>0, pop, aliveonly=aliveonly, sortby=:random, allowempty=allowempty)
   elseif male
      idlist = selectid([:male] => x->x==true, pop, aliveonly=aliveonly, sortby=:random, allowempty=allowempty)
   elseif female
      idlist = selectid([:male] => x->x==!true, pop, aliveonly=aliveonly, sortby=:random, allowempty=allowempty)
   else
      throw(ArgumentError("needs male=true and/or female=true"))
   end
   if length(idlist)>=n
      return idlist[1:n]
   elseif allowempty
      return idlist
   else
      error("too few sampled individuals")
   end
end

function random_sampling(pop::PTPopulation; nm::Int=0, nf::Int=0, aliveonly::Bool=true, allowempty::Bool=true)
   if nm>0
      mlist = selectid([:male] => x->x==true, pop, aliveonly=aliveonly, sortby=:random, allowempty=allowempty)
      if length(mlist)>=nm
         return mlist[1:nm]
      elseif allowempty
         return mlist
      else
         error("too few sampled individuals")
      end
   else
      mlist = Int[]
   end
   if nf>0
      flist = selectid([:male] => x->x==!true, pop, aliveonly=aliveonly, sortby=:random, allowempty=allowempty)
      if length(flist)>=nf
         return flist[1:nf]
      elseif allowempty
         return flist
      else
         error("too few sampled individuals")
      end
   else
      flist = Int[]
   end
   idlist = [mlist; flist]
   return idlist
end

"""
    n = add_sires!(group::PTGroup, idlist::Vector{Int})

Add animals to a sire list in `group`.
The function reads the candidate list from the top, then adds the animal if it is not in the list.
If the size of list reaches the maximum (defined by `maxSire` in the group structure), this function stops adding any more.
It returns how many sires have been put into the sire list in the group.
"""
function add_sires!(group::PTGroup, idlist::Vector{Int}; verbose::Bool=false)
   n = length(group.sires)
   ntaken = 0
   @inbounds for i=1:length(idlist)
      candidate = idlist[i]
      if n>=group.maxSire
         break
      end
      if candidate>0
         if !group.pop.df[candidate,:alive]
            if verbose
               @warn "culled sire $(candidate); skip."
            end
            continue
         end
         if !in(candidate,group.sires) && group.pop.df[candidate,:male]
            push!(group.sires,candidate)
            n = n + 1
            ntaken = ntaken + 1
            if !in(candidate,group.id)
               push!(group.id,candidate)
               push!(group.generation,0)
            end
         end
      end
   end
   group.n = length(group.id)
   return ntaken
end

function add_sires!(group::PTGroup, id::Int; verbose::Bool=false)
   return add_sires!(group, [id], verbose=verbose)
end

"""
    n = add_dams!(group::PTGroup, idlist::Vector{Int})

Add animals to a sire list in `group`.
The function reads the candidate list from the top, then adds the animal if it is not in the list.
If the size of list reaches the maximum (defined by `maxDam` in the group structure), this function stops adding any more.
It returns how many dams have been taken.
"""
function add_dams!(group::PTGroup, idlist::Vector{Int}; verbose::Bool=false)
   n = length(group.dams)
   ntaken = 0
   @inbounds for i=1:length(idlist)
      candidate = idlist[i]
      if n>=group.maxDam
         break
      end
      if candidate>0
         if !group.pop.df[candidate,:alive]
            if verbose
               @warn "culled dam $(candidate); skip."
            end
            continue
         end
         if !in(candidate,group.dams) && !group.pop.df[candidate,:male]
            push!(group.dams,candidate)
            n = n + 1
            ntaken = ntaken + 1
            if !in(candidate,group.id)
               push!(group.id,candidate)
               push!(group.generation,0)
            end
         end
      end
   end
   group.n = length(group.id)
   return ntaken
end

function add_dams!(group::PTGroup, id::Int; verbose::Bool=false)
   return add_dams!(group, [id], verbose=verbose)
end

"""
    n = mark_candidate!(pop::PTPopulation, idlist::Vector{Int})

Mark individuals in `idlist` as candidates for selection.
It just puts a flag in the dataframe, and does nothing more.
"""
function mark_candidate!(pop::PTPopulation,idlist::Vector{Int}; check::Bool=true)
   for id in idlist
      if check
         if id<=0 || id>pop.maxAnimal || !pop.df[id,:alive]
            error("id $(id) out of range, dead, or inappropriate")
         end
      end
      pop.df[id,:candidate] = true
   end
   return nothing
end

"""
    n = change_status!(pop::PTPopulation, idlist::Vector{Int}, status=0)

Change status for individuals in `idlist`.
It just puts a number in the dataframe, and does nothing more.
The status code is arbitrary; the user can define it (the default code is 0).
This function does not check if the individual is alive or not.
"""
function change_status!(pop::PTPopulation,idlist::Vector{Int},status::Int=0; check::Bool=true)
   for id in idlist
      if check
         if id<=0 || id>pop.maxAnimal
            error("id $(id) out of range")
         end
      end
      pop.df[id,:status] = status
   end
   return nothing
end

"""
    n = vacancy_for_sires(group)

The number of vacant sires to be filled.
"""
function vacancy_for_sires(group)
   return group.maxSire - length(group.sires)
end

"""
    n = vacancy_for_dams(group)

The number of vacant dams to be filled.
"""
function vacancy_for_dams(group)
   return group.maxDam - length(group.dams)
end

"""
    n = cull!(pop::PTPopulation, idlist::Vector{Int})
    n = cull!(group::PTGroup, idlist::Vector{Int})

Remove the animal IDs in `idlist` from population `pop` or group `group`.
Returns the number of culled animals, `n`.
With `group`, the animals will also be removed from the population.
The status of the animals becomes "dead", i.e., `alive=false`.
"""
function cull!(pop::PTPopulation,idlist::Vector{Int})
   nculled = 0
   @inbounds for i in idlist
      if 0<i && i<=pop.maxAnimal && pop.df[i,:alive]
         pop.df[i,:alive] = false
         nculled = nculled + 1
      end
   end
   return nculled
end

function cull!(group::PTGroup,idlist::Vector{Int})
   pop = group.pop
   nculled = cull!(pop,idlist)
   setdiff!(group.sires,idlist)
   setdiff!(group.dams,idlist)
   return nculled
end

function cull!(groups::Vector{PTGroup},idlist::Vector{Int})
   nculled = 0
   for group in groups
      nculled = nculled + cull!(group,idlist)
   end
   return nculled
end


text:
using MinimaxCalibratedEBayes
using QuadGK
using OffsetArrays

const MCEB = MinimaxCalibratedEBayes

using ECOS

@testset "Check Sobolev to coefficients expansion" begin

target = PriorDensityTarget(0.0)
marginal_grid = -3:0.02:3
Zs_discr = DiscretizedStandardNormalSamples(marginal_grid)

bc = MCEB.ButuceaComteEstimator(target;  n=1000)
bc_marginal_affine = DiscretizedAffineEstimator(Zs_discr.mhist, bc)

herm_class = HermitePriorClass(qmax=30, sobolev_order=1.0, sobolev_bound=0.2,
                               solver=ECOS.Optimizer)

herm_worst_case_bias = worst_case_bias(bc_marginal_affine,
                                       Zs_discr, herm_class, target,
						               boundary_mass=0.1)

g_min = herm_worst_case_bias.min_g

@test quadgk(x->pdf(g_min, x), -Inf, Inf)[1] ≈ 1
@test cf(g_min, 0.0) ≈ 1 + 0im
@test cf(g_min, 1.0) ≈ quadgk(x->exp(im*x)*pdf(g_min,x), -Inf, Inf)[1]
@test cf(g_min, -1.0) ≈ quadgk(x->exp(-im*x)*pdf(g_min,x), -Inf, Inf)[1]
@test herm_class.sobolev_bound ≈ quadgk(x -> abs2(cf(g_min, x))*(x^2+1), -Inf, Inf)[1]/(2π)

#function _coefs_to_sobolev_norm(αs)
#	q = length(αs)
#	αstilde = zeros(q + 1)
#	αs_o = OffsetArray(αs, 0:(q-1))
#	αstilde_o = OffsetArray(αstilde, 0:q)
#	αstilde_o[0] = αs_o[1]/sqrt(2)
#	for j=1:q
#	end
#	αstilde_o
#end
function _coefs_to_sobolev_norm(αs)
	q = length(αs)
	αs_o = OffsetArray([αs;0], 0:q)
	quad_term = abs2(αs_o[1])/2 + sum(abs2, αs_o)
	sob_vec = [(sqrt(j)*αs_o[j-1]-sqrt(j+1)*αs_o[j+1]) for j in 1:(q-1)]
	sob_term = sum(abs2, sob_vec)/2
	sob_norm = quad_term + sob_term

	herm_vec = [(sqrt(j/2)*αs_o[j-1]+sqrt((j+1)/2)*αs_o[j+1]) for j in 1:(q-1)]
	hermite_norm = abs2(αs_o[1])/2 + sum(abs2, herm_vec)
	joint_norm_bound = sum([ (j*abs2(αs_o[j-1]) + (j+1)*abs2(αs_o[j+1])) for j in 1:(q-1)])
	joint_norm_bound_diff =  abs2(αs_o[1]) + sum(abs2, αs_o)
	(sob_norm = sob_norm, hermite_norm=hermite_norm,
	 joint_norm_bound = joint_norm_bound,
	 joint_norm_bound_diff = joint_norm_bound_diff)
end
sob_norms = _coefs_to_sobolev_norm(g_min.coefs)

@test sob_norms[:sob_norm] ≈ herm_class.sobolev_bound atol=0.00001
@test sob_norms[:hermite_norm] ≈ quadgk(x->abs2(pdf(g_min,x)*x), -Inf,+Inf)[1] atol=0.00001
@test sob_norms[:joint_norm_bound] < sob_norms[:hermite_norm] + sob_norms[:sob_norm]

@test sob_norms[:hermite_norm] + sob_norms[:sob_norm] - sob_norms[:joint_norm_bound] ≈ sob_norms[:joint_norm_bound_diff]

end

ext: jl


text:
export Frobenius

#Frobenius matrices or Gaussian elimination matrices
#
#Matrices of the form
#[ 1 0 ...     0 ]
#[ 0 1 ...     0 ]
#[ .........     ]
#[ ... 1 ...     ]
#[ ... c1 1 ...  ]
#[ ... c2 0 1 ...]
#[ ............. ]
#[ ... ck ...   1]
#
#i.e. an identity matrix with nonzero subdiagonal elements along a single
#column.
#
#In this implementation, the subdiagonal of the nonzero column is stored as a
#dense vector, so that the size can be inferred automatically as j+k where j is
#the index of the column and k is the number of subdiagonal elements.

immutable Frobenius{T} <: AbstractArray{T, 2}
    colidx :: Int
    c :: Vector{T}
end

#Basic property computations
size(F::Frobenius, r::Int) = (r==1 || r==2) ? F.colidx + length(F.c) :
    throw(ArgumentError("Frobenius matrix is of rank 2"))

function size(F::Frobenius)
    n = F.colidx + length(F.c)
    n, n
end

#XXX Inefficient but works
getindex(F::Frobenius, i, j) = getindex(full(F), i, j)
isassigned(F::Frobenius, i, j) = isassigned(full(F), i, j)

function full{T}(F::Frobenius{T})
    M = eye(T, size(F, 1))
    M[F.colidx+1:end, F.colidx] = F.c
    M
end

#Linear algebra stuff
function A_mul_B!{T}(F::Frobenius{T}, b::Vector{T})
    (n = size(F, 2)) == length(b) || throw(DimensionMismatch("$n $(length(b))"))
    for i=F.colidx+1:n
        b[i] += F.c[i-F.colidx] * b[F.colidx]
    end
    b
end
*{T}(F::Frobenius{T}, b::Vector{T}) = A_mul_B!(F, copy(b))

function *{T}(F::Frobenius{T}, G::Frobenius{T})
    (n = size(F, 2)) == size(G, 2) || throw(DimensionMismatch(""))
    if F.colidx == G.colidx #Answer is still expressible as a Frobenius
        return Frobenius(F.colidx, F.c + G.c)
    else
        M = full(F)
        M[G.colidx+1:end, G.colidx] = F.colidx < G.colidx ? G.c :
            Frobenius(F.colidx-G.colidx, F.c) * G.c
        return M
    end
end

inv(F::Frobenius) = Frobenius(F.colidx, -F.c)


text:
import FastTransforms: paduapoints, plan_chebyshevtransform, plan_ichebyshevtransform
import ClassicalOrthogonalPolynomials

const COP = ClassicalOrthogonalPolynomials

struct PoincareInvariant2ndCOP{DT,TT,ET,ΩT,CT} <: AbstractPoincareInvariant2nd{DT}
    equ::ET
    ω::ΩT
    Δt::TT
    nx::Int
    ny::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::OffsetArray{Double64,1,Vector{Double64}}
    J::OffsetArray{Double64,1,Vector{Double64}}
    ΔI::OffsetArray{Double64,1,Vector{Double64}}
    ΔJ::OffsetArray{Double64,1,Vector{Double64}}
    y::Array{DT,3}
    z̃::Array{DT,3}
    z::Array{DT,3}
    Ω::Matrix{DT}
    cheb::CT
end

function PoincareInvariant2ndCOP(f_equ::Function, f_surface::Function, ω::ΩT, Δt::TT, d::Int, nx::Int, ny::Int, ntime::Int, nsave::Int=1, DT=Float64) where {TT,ΩT}

    nx = 2nx
    ny = 2ny

    if get_config(:verbosity) > 1
        println()
        println("Second Euler-Poincaré Integral Invariant (OrthogonalPolynomials)")
        println("================================================================")
        println()
        println(" nx    = ", nx)
        println(" ny    = ", ny)
        println(" ntime = ", ntime)
        println(" nsave = ", nsave)
        println(" Δt    = ", Δt)
        println()
    end

    # compute initial conditions
    xGrid = COP.grid(COP.ChebyshevT()[:,1:nx])
    yGrid = COP.grid(COP.ChebyshevT()[:,1:ny])

    q₀ = reshape([f_surface((xGrid[i] + 1) / 2, (yGrid[j] + 1) / 2) for i in 1:nx, j in 1:ny], nx*ny)

    equ = f_equ(q₀)

    # compute initial conditions
    # _padua_length(N) = Int(cld(-3+sqrt(1+8N),2))

    # c = paduapoints(DT, _padua_length(nx*ny))

    # println(length(c))

    # q₀ = zeros(DT, (d, length(c)))

    # for i in axes(c,1)
    #     q₀[:,i] .= f_surface(c[i][1], c[i][2])
    # end

    # equ = f_equ(q₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I  = OffsetArray(zeros(Double64, nt+1), 0:nt)
    J  = OffsetArray(zeros(Double64, nt+1), 0:nt)
    ΔI = OffsetArray(zeros(Double64, nt+1), 0:nt)
    ΔJ = OffsetArray(zeros(Double64, nt+1), 0:nt)

    y  = zeros(DT, d, nx,   ny  )
    z̃  = zeros(DT, d, nx,   ny-1)
    z  = zeros(DT, d, nx-1, ny-1)
    Ω  = zeros(DT, d, d)

    Di = zeros(DT, d, nx-1, ny  )
    Dj = zeros(DT, d, nx,   ny-1)
    D1 = zeros(DT, d, nx-1, ny-1)
    D2 = zeros(DT, d, nx-1, ny-1)

    Tx = COP.ChebyshevT()[:,1:nx]
    Ty = COP.ChebyshevT()[:,1:ny]
    Ux = COP.ChebyshevU()[:,1:nx-1]
    Uy = COP.ChebyshevU()[:,1:ny-1]

    Dx = COP.Derivative(axes(Tx,1)) * Tx
    Dy = COP.Derivative(axes(Ty,1)) * Ty

    PTx = plan_chebyshevtransform(zeros(nx), Val(1))
    PTy = plan_chebyshevtransform(zeros(ny), Val(1))

    iPUx = plan_ichebyshevtransform(zeros(nx-1), Val(2))
    iPUy = plan_ichebyshevtransform(zeros(ny-1), Val(2))

    PTUx = Ux \ Tx
    PTUy = Uy \ Ty

    cheb = (Di = Di, Dj = Dj, D1 = D1, D2 = D2,
            Tx = Tx, Ty = Ty, Ux = Ux, Uy = Uy,
            Dx = Dx, Dy = Dy,
            PTx  = PTx,  PTy  = PTy,
            iPUx = iPUx, iPUy = iPUy,
            PTUx = PTUx, PTUy = PTUy)

    PoincareInvariant2ndCOP{DT,TT,typeof(equ),ΩT,typeof(cheb)}(equ, ω, Δt, nx, ny, ntime, nsave, nt, I, J, ΔI, ΔJ, y, z̃, z, Ω, cheb)
end

# function compute_position(x::AbstractArray{DT}, nx, ny, Ux, Uy, Px, Py, PUx, PUy, UxGrid, UyGrid, PUxi, PUyi) where {DT}
#     local y = copy(reshape(x, (size(x,1), nx, ny)))
#     local z̃ = zeros(DT, (size(x,1), nx,   ny-1))
#     local z = zeros(DT, (size(x,1), nx-1, ny-1))

#     for d in axes(y,1)
#         for i in axes(y,2)
#             @views y[d,i,:] .= Py \ y[d,i,:]
#             # @views y[d,i,:] .= Py * y[d,i,:]
#         end
#         for j in axes(y,3)
#             @views y[d,:,j] .= Px \ y[d,:,j]
#             # @views y[d,:,j] .= Px * y[d,:,j]
#         end
#     end

#     for d in axes(z̃,1)
#         for i in axes(z̃,2)
#             @views z̃[d,i,:] .= PUy * y[d,i,:]
#         end
#     end
#     for d in axes(z,1)
#         for j in axes(z,3)
#             @views z[d,:,j] .= PUx * z̃[d,:,j]
#         end
#     end

#     for d in axes(z,1)
#         for i in axes(z,2)
#             @views z[d,i,:] .= (Uy * z[d,i,:])[UyGrid]
#             # @views z[d,i,:] .= PUyi * z[d,i,:]
#         end
#         for j in axes(z,3)
#             @views z[d,:,j] .= (Ux * z[d,:,j])[UxGrid]
#             # @views z[d,:,j] .= PUxi * z[d,:,j]
#         end
#     end

#     return y, z
# end

# function compute_derivative_i!(Di, D1, y, Dx, Ux, Uy, PUy, UxGrid, UyGrid, PUxi, PUyi)
#     for d in axes(Di,1)
#         for j in axes(Di,3)
#             # @views Di[d,:,j] .= Ux \ ( Dx * y[d,:,j] )
#             @views Di[d,:,j] .= ( Dx * y[d,:,j] ).args[2].args[1]

#         end
#     end

#     for d in axes(D1,1)
#         for i in axes(D1,2)
#             @views D1[d,i,:] .= PUy * Di[d,i,:]
#         end
#     end

#     for d in axes(D1,1)
#         for i in axes(D1,2)
#             @views D1[d,i,:] .= (Uy * D1[d,i,:])[UyGrid]
#             # @views D1[d,i,:] .= PUyi * D1[d,i,:]
#         end
#         for j in axes(D1,3)
#             @views D1[d,:,j] .= (Ux * D1[d,:,j])[UxGrid]
#             # @views D1[d,:,j] .= PUxi * D1[d,:,j]
#         end
#     end
# end

# function compute_derivative_j!(Dj, D2, y, Dy, Ux, Uy, PUx, UxGrid, UyGrid, PUxi, PUyi)
#     for d in axes(Dj,1)
#         for i in axes(Dj,2)
# #            @views Dj[d,i,:] .= Uy \ ( Dy * y[d,i,:] )
#             @views Dj[d,i,:] .= ( Dy * y[d,i,:] ).args[2].args[1]
#         end
#     end

#     for d in axes(D2,1)
#         for j in axes(D2,3)
#             @views D2[d,:,j] .= PUx * Dj[d,:,j]
#         end
#     end

#     for d in axes(D2,1)
#         for i in axes(D2,2)
#             @views D2[d,i,:] .= (Uy * D2[d,i,:])[UyGrid]
#             # @views D2[d,i,:] .= PUyi * D2[d,i,:]
#         end
#         for j in axes(D2,3)
#             @views D2[d,:,j] .= (Ux * D2[d,:,j])[UxGrid]
#             # @views D2[d,:,j] .= PUxi * D2[d,:,j]
#         end
#     end
# end

# function compute_position(x::AbstractArray{DT}, nx, ny, PTx, PTy, PTUx, PTUy, iPUx, iPUy) where {DT}
#     local y = copy(reshape(x, (size(x,1), nx, ny)))
#     local z̃ = zeros(DT, (size(x,1), nx,   ny-1))
#     local z = zeros(DT, (size(x,1), nx-1, ny-1))

#     # Py = iPUy * PTUy * PTy
#     # Px = iPUx * PTUx * PTx

#     for d in axes(z̃,1)
#         for i in axes(z̃,2)
#             # @views z̃[d,i,:] .= Py * y[d,i,:]
#             @views z̃[d,i,:] .= iPUy * ( PTUy * ( PTy * y[d,i,:] ) )
#             # @views z̃[d,i,:] .= iPUy * ( PTUy * ( PTy \ y[d,i,:] ) )
#         end
#     end

#     for d in axes(z,1)
#         for j in axes(z,3)
#             # @views z[d,:,j] .= Px * z̃[d,:,j]
#             @views z[d,:,j] .= iPUx * ( PTUx * ( PTx * z̃[d,:,j] ) )
#             # @views z[d,:,j] .= iPUx * ( PTUx * ( PTx \ z̃[d,:,j] ) )
#         end
#     end

#     return z
# end

# function compute_derivative_i!(Di, D1, x, nx, ny, Dx, PTx, PTy, PTUy, iPUx, iPUy)
#     local y = copy(reshape(x, (size(x,1), nx, ny)))

#     # Py = iPUy * PTUy * PTy

#     for d in axes(Di,1)
#         for j in axes(Di,3)
#             @views Di[d,:,j] .= iPUx * ( Dx * (PTx * y[d,:,j]) ).args[2].args[1]
#             # @views Di[d,:,j] .= iPUx * ( Dx * (PTx \ y[d,:,j]) ).args[2].args[1]
#         end
#     end

#     for d in axes(D1,1)
#         for i in axes(D1,2)
#             # @views D1[d,i,:] .= Py * Di[d,i,:]
#             @views D1[d,i,:] .= iPUy * ( PTUy * ( PTy * Di[d,i,:] ) )
#             # @views D1[d,i,:] .= iPUy * ( PTUy * ( PTy \ Di[d,i,:] ) )
#         end
#     end
# end

# function compute_derivative_j!(Dj, D2, x, nx, ny, Dy, PTx, PTy, PTUx, iPUx, iPUy)
#     local y = copy(reshape(x, (size(x,1), nx, ny)))

#     # Px = iPUx * PUx * PTx

#     for d in axes(Dj,1)
#         for i in axes(Dj,2)
#             @views Dj[d,i,:] .= iPUy * ( Dy * (PTy * y[d,i,:]) ).args[2].args[1]
#             # @views Dj[d,i,:] .= iPUy * ( Dy * (PTy \ y[d,i,:]) ).args[2].args[1]
#         end
#     end

#     for d in axes(D2,1)
#         for j in axes(D2,3)
#             # @views D2[d,:,j] .= Px * Dj[d,:,j]
#             @views D2[d,:,j] .= iPUx * ( PTUx * ( PTx * Dj[d,:,j] ) )
#             # @views D2[d,:,j] .= iPUx * ( PTUx * ( PTx \ Dj[d,:,j] ) )
#         end
#     end
# end

# function _integrate_cop(t, ω, z, nx, ny, D1, D2, Ω::Matrix)
#     local result::Double64 = 0

#     for j in axes(z,3)
#         for i in axes(z,2)
#             ω(t, z[:,i,j], Ω)
#             @views result += D2[:,i,j]' * Ω * D1[:,i,j]
#         end
#     end

#     return 4 * result / nx / ny
# end

function _integrate_cop(t, x, pinv)
    Di, Dj, D1, D2, Tx, Ty, Ux, Uy, Dx, Dy, PTx, PTy, iPUx, iPUy, PTUx, PTUy = pinv.cheb

    pinv.y .= reshape(x, (size(x,1), pinv.nx, pinv.ny))

    for d in axes(pinv.y,1)
        for i in axes(pinv.y,2)
            @views pinv.y[d,i,:] .= PTy * pinv.y[d,i,:]
        end
        for j in axes(pinv.y,3)
            @views pinv.y[d,:,j] .= PTx * pinv.y[d,:,j]
        end
    end

    for d in axes(pinv.z̃,1)
        for i in axes(pinv.z̃,2)
            @views pinv.z̃[d,i,:] .= PTUy * pinv.y[d,i,:]
        end
    end
    for d in axes(pinv.z,1)
        for j in axes(pinv.z,3)
            @views pinv.z[d,:,j] .= PTUx * pinv.z̃[d,:,j]
        end
    end

    for d in axes(pinv.z,1)
        for i in axes(pinv.z,2)
            @views pinv.z[d,i,:] .= iPUy * pinv.z[d,i,:]
        end
        for j in axes(pinv.z,3)
            @views pinv.z[d,:,j] .= iPUx * pinv.z[d,:,j]
        end
    end

    for d in axes(Di,1)
        for j in axes(Di,3)
            @views Di[d,:,j] .= ( Dx * pinv.y[d,:,j] ).args[2].args[1]
        end
    end
    for d in axes(D1,1)
        for i in axes(D1,2)
            @views D1[d,i,:] .= iPUy * (PTUy * Di[d,i,:])
        end
        for j in axes(D1,3)
            @views D1[d,:,j] .= iPUx * D1[d,:,j]
        end
    end

    for d in axes(Dj,1)
        for i in axes(Dj,2)
            @views Dj[d,i,:] .= ( Dy * pinv.y[d,i,:] ).args[2].args[1]
        end
    end
    for d in axes(D2,1)
        for j in axes(D2,3)
            @views D2[d,:,j] .= iPUx * (PTUx * Dj[d,:,j])
        end
        for i in axes(D2,2)
            @views D2[d,i,:] .= iPUy * D2[d,i,:]
        end
    end

    local result::Double64 = 0

    for j in axes(pinv.z,3)
        for i in axes(pinv.z,2)
            @views pinv.ω(t, pinv.z[:,i,j], pinv.Ω)
            @views result += D2[:,i,j]' * pinv.Ω * D1[:,i,j]
        end
    end

    return 4 * result / pinv.nx / pinv.ny
end

function surface_integral(pinv, t, x::AbstractMatrix{DT}) where {DT}
    # local B = zeros(DT, size(x,1), size(x,1))
    # local I = zero(Double64)

    # y, z = compute_position(x, pinv.nx, pinv.ny, Ux, Uy, Px, Py, PUx, PUy, UxGrid, UyGrid)
    # y, z = compute_position(x, pinv.nx, pinv.ny, Ux, Uy, Px, Py, PUx, PUy, UxGrid, UyGrid, UxiPlan, UyiPlan)

    # z = compute_position(x, pinv.nx, pinv.ny, PTx, PTy, PTUx, PTUy, iPUx, iPUy)

    # compute_derivative_i!(Di, D1, x, pinv.nx, pinv.ny, Dx, PTx, PTy, PTUy, iPUx, iPUy)
    # compute_derivative_j!(Dj, D2, x, pinv.nx, pinv.ny, Dy, PTx, PTy, PTUx, iPUx, iPUy)

    # return _integrate_cop(t, pinv.ω, z, pinv.nx, pinv.ny, D1, D2, B)

    return _integrate_cop(t, x, pinv)
end

# function integrate_canonical(γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ, b::Vector{TT}, c::Vector{TT}, vᵢ::Vector{DT}, vⱼ::Vector{DT}) where {DT,TT}
#     @assert length(b) == length(c)

#     local result = zero(Double64)

#     for i in eachindex(b)
#         for j in eachindex(b)
#             Θ̇ᵢ(c[i], c[j], vᵢ)
#             γ̇ⱼ(c[i], c[j], vⱼ)
#             result += b[i] * b[j] * dot(vᵢ,vⱼ)

#             γ̇ᵢ(c[i], c[j], vᵢ)
#             Θ̇ⱼ(c[i], c[j], vⱼ)
#             result -= b[i] * b[j] * dot(vᵢ,vⱼ)
#         end
#     end

#     return result
# end

# function surface_integral_canonical(q::AbstractMatrix{DT}, p::AbstractMatrix{DT}, nx, ny) where {DT}
#     local b = [0.5, 0.5]
#     local c = [0.0, 1.0]

#     local vᵢ = zeros(DT, size(q,1))
#     local vⱼ = zeros(DT, size(q,1))
#     local I  = zero(Double64)

#     integrate_trapezoidal = (γ, γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ) -> integrate_canonical(γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ, b, c, vᵢ, vⱼ)

#     for j in 1:ny-1
#         for i in 1:nx-1
#             γ̇ᵢ = (λ, μ, y) -> interpolate_derivative_i(q, i, j, λ, μ, y, nx, ny)
#             γ̇ⱼ = (λ, μ, y) -> interpolate_derivative_j(q, i, j, λ, μ, y, nx, ny)
#             Θ̇ᵢ = (λ, μ, y) -> interpolate_derivative_i(p, i, j, λ, μ, y, nx, ny)
#             Θ̇ⱼ = (λ, μ, y) -> interpolate_derivative_j(p, i, j, λ, μ, y, nx, ny)
#             I += integrate_trapezoidal(γ, γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ)
#         end
#     end

#     return I
# end

function evaluate_poincare_invariant(pinv::PoincareInvariant2ndCOP{DT}, sol::Solution) where {DT}
    local verbosity = get_config(:verbosity)

    for i in axes(sol.q,1)
        verbosity > 1 ? println("      it = ", i) : nothing
        pinv.I[i]  = surface_integral(pinv, sol.t[i], hcat(sol.q[i,:]...))
        pinv.ΔI[i] = abs(pinv.I[0]) < sqrt(eps()) ? pinv.I[i] : (pinv.I[i] .- pinv.I[0]) ./ pinv.I[0]
        verbosity > 1 ? println("           I_q = ", pinv.I[i], ",   ε_q = ", pinv.ΔI[i]) : nothing

        # if hasproperty(sol, :p)
        #     pinv.J[i] = surface_integral_canonical(pinv, hcat(sol.q[i,:]...), hcat(sol.p[i,:]...))
        #     pinv.ΔJ[i] = abs(pinv.J[0]) < sqrt(eps()) ? pinv.J[i] : (pinv.J[i] .- pinv.J[0]) ./ pinv.J[0]
        #     verbosity > 1 ? println("           I_p = ", pinv.J[i], ",   ε_p = ", pinv.ΔJ[i]) : nothing
        # end
    end

    return (DT.(pinv.I), DT.(pinv.J), DT.(pinv.ΔI), DT.(pinv.ΔJ))
end

function write_to_hdf5(pinv::PoincareInvariant2ndCOP, sol::Solution, output_file::String)
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5

        write(h5, "t", sol.t)
        write(h5, "I", pinv.I)

        isdefined(sol, :p) ? write(h5, "J", pinv.J) : nothing

    end
end


text:
using LazySets

function projectile_specification()
    # initial set
    X0 = Singleton([0.0, 5.0, 100.0, 0.0])

    # time horizon: 20 time units
    time_horizon = 20.0

    # specification
    O = Dict{Symbol, Any}(:T => time_horizon)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
using YaoExperiment
using YaoExperiment: WEBCONFIG

WEBCONFIG[:maxn_dense] = 10
WEBCONFIG[:maxn_sparse] = 15
WEBCONFIG[:PORT] = 8080
SERVER = run_server(YaoExperiment.naive_handler)


text:
module ITensorMakie

using Colors
using Graphs
using NetworkLayout
using Reexport
using GraphMakie
@reexport using ITensorVisualizationBase

using GraphMakie.Makie:
  Makie,
  Figure,
  contents,
  hidedecorations!,
  hidespines!,
  deregister_interaction!,
  register_interaction!

using ITensorVisualizationBase:
  @Backend_str,
  default_vertex_labels,
  default_vertex_labels_prefix,
  default_vertex_size,
  default_vertex_textsize,
  default_edge_textsize,
  default_edge_widths,
  default_edge_labels,
  default_arrow_show,
  default_arrow_size,
  default_siteinds_direction,
  is_self_loop,
  _ndims

import ITensorVisualizationBase: visualize, visualize!, _graphplot

function __init__()
  return ITensorVisualizationBase.set_backend!(Backend"Makie"())
end

fill_number(a::AbstractVector, n::Integer) = a
fill_number(x::Number, n::Integer) = fill(x, n)

function visualize(b::Backend"Makie", g::AbstractGraph; kwargs...)
  f = Figure()
  visualize!(b, f[1, 1], g; kwargs...)
  return f
end

function visualize!(b::Backend"Makie", f::Figure, g::AbstractGraph; kwargs...)
  visualize!(b, f[1, 1], g; kwargs...)
  return f
end

function visualize!(
  b::Backend"Makie",
  f::Makie.GridPosition,
  g::AbstractGraph;
  interactive=true,
  ndims=2,
  layout=Spring(; dim=ndims),

  # vertex
  vertex_labels_prefix=default_vertex_labels_prefix(b, g),
  vertex_labels=default_vertex_labels(b, g, vertex_labels_prefix),
  vertex_size=default_vertex_size(b, g),
  vertex_textsize=default_vertex_textsize(b, g),

  # edge
  edge_textsize=default_edge_textsize(b),
  edge_widths=default_edge_widths(b, g),
  edge_labels=default_edge_labels(b, g),

  # arrow
  arrow_show=default_arrow_show(b, g),
  arrow_size=default_arrow_size(b, g),
  siteinds_direction=default_siteinds_direction(b, g),
)
  if ismissing(Makie.current_backend[])
    error(
      """
  You have not loaded a backend.  Please load one (`using GLMakie` or `using CairoMakie`)
  before trying to visualize a graph.
"""
    )
  end

  edge_labels = ITensorVisualizationBase.edge_labels(b, edge_labels, g)

  if length(vertex_labels) ≠ nv(g)
    throw(
      DimensionMismatch(
        "$(length(vertex_labels)) vertex labels $(vertex_labels) were specified but there are $(nv(g)) tensors in the diagram, please specify the correct number of labels.",
      ),
    )
  end

  graphplot_kwargs = (;
    layout=layout,

    # vertex
    node_size=fill_number(vertex_size, nv(g)),
    node_color=colorant"lightblue1", # TODO: store in vertex, make a default
    node_marker='●', # TODO: allow other options, like '◼'
    node_attr=(; strokecolor=:black, strokewidth=3),

    # vertex labels
    nlabels=vertex_labels,
    nlabels_textsize=vertex_textsize,
    nlabels_color=colorant"black",
    nlabels_align=(:center, :center),

    # edge
    edge_width=edge_widths,
    edge_color=colorant"black",

    # edge labels
    elabels=edge_labels,
    elabels_textsize=edge_textsize,
    elabels_color=colorant"red",

    # self-edge
    selfedge_width=1e-5, # Small enough so you can't see the loop, big enough for site label to show up
    selfedge_direction=siteinds_direction,
    selfedge_size=3,

    # arrow
    arrow_show=arrow_show,
    arrow_size=arrow_size,
    arrow_shift=0.49,
  )

  overwrite_axis = false
  if isempty(contents(f))
    axis_plot = graphplot(f, g; graphplot_kwargs...)
  else
    @warn "Visualizing a graph in the same axis as an existing graph. This feature is experimental and some features like interactivity might now work"
    overwrite_axis = true
    graphplot!(f, g; graphplot_kwargs...)
  end

  # Set rotation of edge labels to 0.
  axis_plot.plot.elabels_rotation[] = zeros(ne(g))

  if !overwrite_axis && (_ndims(layout) == 2)
    hidedecorations!(axis_plot.axis)
    # This would hide the box around the plot
    # TODO: make this optional
    #hidespines!(axis_plot.axis)
    if interactive
      deregister_interaction!(axis_plot.axis, :rectanglezoom)
      register_interaction!(axis_plot.axis, :nhover, NodeHoverHighlight(axis_plot.plot))
      register_interaction!(axis_plot.axis, :ehover, EdgeHoverHighlight(axis_plot.plot))
      register_interaction!(axis_plot.axis, :ndrag, NodeDrag(axis_plot.plot))
      register_interaction!(axis_plot.axis, :edrag, EdgeDrag(axis_plot.plot))
    end
  end
  return f
end

# For use in sequence visualization.
# TODO: Make this more generalizable to other backends.
function _graphplot(::Backend"Makie", graph; all_labels)
  fig, ax, plt = graphplot(
    reverse(graph); arrow_show=false, nlabels=all_labels, layout=Buchheim()
  )
  hidedecorations!(ax)
  #hidespines!(ax)
  return fig
end

end


text:
### Julia module for Network Flows algorithms
__precompile__() # Precompile option for the module. Is recompiled after changes

module NetworkFlows

# Items to be exported
export Network # constructors
export zero_to_one! # index the nodes from 1 to n instead of 0 to (n-1) for the vectors in Julia
export print # redefinition of the print function for Network
export edmondsKarp # max-flow
export mincut # min-cut
export connectivity # connectivity (max-flow on the unweighted graph)
export kishimoto # multiroute flow
export breakingPoints # of the extended multiroute flow function
export mixedMLA, successMLA # Solution to Network Interdiction/ MLA problems

include("links.jl")
include("network.jl")
include("io.jl")
include("search.jl")
include("flow.jl")
include("cut.jl")
include("extmrflow.jl")
include("mla.jl")

end # module


text:
using ModelingToolkit
using DifferentialEquations
using Plots
using BenchmarkTools
using Coexist

# initial state
nAge=9
nHS=8
nIso=4
nI=4
nTest=4
param=(nAge=nAge, nHS=nHS, nIso=nIso, nTest=nTest)
ssize=(nAge, nHS, nIso, nTest)

function dydt_Complete(stateTensor_flattened,p,t)
	T = eltype(stateTensor_flattened)
	stateTensor = reshape(stateTensor_flattened, (nTest, nIso, nHS, nAge))
	dydt = zeros(T,size(stateTensor)...)
	trTensor_complete = zeros(T,(nTest, nIso, nHS, nTest, nIso, nHS, nAge))
	trTensor_diseaseProgression = Coexist.trFunc_diseaseProgression(;param...)
	for k in 1:4
		shape = size(trTensor_diseaseProgression[:,k,:,:])
		expand_dims = reshape(trTensor_diseaseProgression[:,k,:,:],
					(shape[1:end-2]..., 1, shape[end-1:end]...)) # equal to np.exapand_dims
		# slice = trTensor_complete[:,k,:,:,k,:,:]
		# @einsum view[m,l,j,i] := slice[l,m,l,j,i]
		view = Coexist.einsum("ijlml->ijlm", trTensor_complete[:,k,:,:,k,:,:],T)
		view .+= expand_dims
		trTensor_complete[:,k,:,:,k,:,:] = Coexist._einsum12(trTensor_complete[:,k,:,:,k,:,:], view, T)
		# @einsum slice[l,m,l,j,i] = view[m,l,j,i]
		# trTensor_complete[:,k,:,:,k,:,:] = slice
	end

	# @einsum to_be_modified[l,k,j,i] := trTensor_complete[l,k,j,l,k,j,i]
	to_be_modified = Coexist.einsum("ijkljkl->ijkl", trTensor_complete, T)
	to_be_modified .-= Coexist.einsum("...jkl->...", trTensor_complete, T)
	# @einsum trTensor_complete[l,k,j,l,k,j,i] = view[l,k,j,i]
	trTensor_complete = Coexist._einsum11(trTensor_complete, to_be_modified, T)

	dydt = Coexist.einsum("ijkl,ijklmnp->imnp", stateTensor, trTensor_complete, T)
	return vec(dydt)
end

state = 50*ones(9*8*4*4)
dstate = similar(state)

using ModelingToolkit
@variables t u[1:1152](t)
@derivatives D'~t

dydt_symbolic = dydt_Complete(u,0,nothing)
dydt_symbolic = simplify.(dydt_symbolic)
sys = ODESystem(D.(u) .~ dydt_symbolic,t,u,[])
computed_f = generate_function(sys)[2]

open(joinpath(@__DIR__,"..","src","generated_dydt.jl"), "w") do io
   write(io, "const dydt = $computed_f")
end

f = eval(computed_f)

function solveSystem_unoptimized(
	state,
	timeSpan = (0.0, 80.0),
	p=nothing
	)
	prob = ODEProblem(dydt_Complete,state,(0.0,80.0),p=nothing)
	sol = solve(prob,Tsit5(),reltol=1e-3,abstol=1e-3)
	sol = convert(Array, sol)
	return reshape(sol, (4,4,8,9, length(sol) ÷ (4*4*8*9)))
end

function solveSystem(
	state,
	timeSpan = (0.0, 80.0),
	p=nothing
	)
	prob = ODEProblem(f,state,(0.0,80.0),p=nothing)
	sol = solve(prob,Tsit5(),reltol=1e-3,abstol=1e-3)
	sol = convert(Array, sol)
	reshape(sol, (4,4,8,9, length(sol) ÷ (4*4*8*9)))
end

using BenchmarkTools
dydt_benchmark = @benchmark dydt_Complete(state,0,nothing)
f_benchmark = @benchmark f(dstate,state,0,nothing)
solveSystem_unoptimized_benchmark = @benchmark solveSystem_unoptimized(state)
solveSystem_benchmark = @benchmark solveSystem(state)

println("BENCHMARK OF f")
display(dydt_benchmark)
display(f_benchmark)
println()
println("BENCHMARK OF solveSystem")
display(solveSystem_unoptimized_benchmark)
display(solveSystem_benchmark)
println()

max_stars_repo_name: vollmersj/coexist-julia

max_forks_repo_name: vollmersj/coexist-julia


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "moments" begin
    N = 3
    @testset "moments constructor" begin
        g = CircuitGate((2,), XGate())
        m = Moment(g)
        @test m[1] === g

        h = CircuitGate((3,), YGate())
        m = Moment([g,h])
        @test m[1] === g
        @test m[2] === h
    end

    @testset "moments constructor exceptions" begin
        g = CircuitGate((2,), XGate())
        h = CircuitGate((2,), YGate())
        @test_throws ErrorException("Only gates on different wires are allowed in a Moment") cgc_ref =
            Moment([g, h])
    end

    @testset "moments adjoint" begin
        g = CircuitGate((2,), ZGate())
        h = CircuitGate((1,), SGate())
        m = Moment([g, h])
        madj = adjoint(m)
        @test madj[1] ≈ adjoint(h)
        @test madj[2] ≈ adjoint(g)

        @test sparse_matrix(adjoint(m)) ≈ adjoint(sparse_matrix(m))
    end

    @testset "moments approx" begin
        g = CircuitGate((2,), ZGate())
        h = CircuitGate((1,), SGate())
        m = Moment([g, h])
        m1 = Moment([g, h])
        m2 = Moment([h, g])

        @test !(m1 ≈ m2)
        @test m1 ≈ m
    end

    @testset "moments array operations" begin
        m = Moment([CircuitGate((2,), ZGate()), CircuitGate((1,), SGate())])

        @test all(m[:] .≈ m.gates)
        @test all(m[1:2] .≈ m.gates)
        @test isempty(m) == false
        @test CircuitGate((1,), SGate()) ≈ pop!(m)
        @test isempty(m) == false
        @test CircuitGate((2,), ZGate()) ≈ pop!(m)
        @test isempty(m) == true
        @test_throws ArgumentError("array must be non-empty") pop!(m)
    end

    @testset "moments sparse_matrix" begin
        g = CircuitGate((2,), ZGate())
        h = CircuitGate((1,), SGate())
        m = Moment([g, h])

        @test sparse_matrix(m) ≈ sparse_matrix(g, 2) * sparse_matrix(h, 2)
    end

    @testset "moments reverse" begin
        g = CircuitGate((2,), ZGate())
        h = CircuitGate((1,), SGate())
        m = Moment([g, h])
        m_reverse = reverse(m)
        @test m_reverse[2] ≈ g
        @test m_reverse[1] ≈ h
        @test m[1] ≈ g
        @test m[2] ≈ h

        @test all(reverse!(m) .≈ m_reverse)
        @test m[1] ≈ h
        @test m[2] ≈ g
    end
end


text:
__precompile__()
module Astellarn
    const ASTELLARN_VERSION = "v0.2.9"

    using Crayons
    using StaticArrays
    using Printf
    #using SIMD

    import Base.&, Base.|, Base.~, Base.<<, Base.>>, Base.⊻, Base.!, Base.bswap
    import Base.isempty, Base.isone, Base.isequal
    import Base.getindex, Base.setindex!, Base.push!
    import Base.iterate, Base.length, Base.eltype, Base.size, Base.IndexStyle
    import Base.copy!
    import Base.show

    include("../deps/config.jl")
    include("bitboard.jl")
    include("pieces.jl")
    include("zobrist.jl")
    include("timeman.jl")
    include("utils.jl")
    include("parameters.jl")
    include("board.jl")
    include("types.jl")
    include("fen.jl")
    include("magic.jl")
    include("attacks.jl")
    include("move.jl")
    include("movegen.jl")
    include("movecount.jl")
    include("perft.jl")
    include("masks.jl")
    include("evaluate.jl")
    include("syzygy.jl")
    include("transposition.jl")
    include("thread.jl")
    include("moveorder.jl")
    include("history.jl")
    include("search.jl")
    include("repl.jl")
    include("uci.jl")

    export Bitboard, Board, Piece, PieceType, Color, Magic, Move, Undo, MoveStack, UndoStack
    export @newgame, @move, @random, @engine, @importfen, @perft

    export importfen, exportfen
    export pawns, kings, bishops, knights, rooks, queens, enemy, friendly, occupied, empty, rooklike, bishoplike
    export checkers, pinned, cancastlekingside, cancastlequeenside
    export ischeck, islegal, ischeckmate, isstalemate, isdrawbymaterial
    export monkey!, perft, engine!

    export static_exchange_evaluator

    export WHITE, BLACK
    export PAWN, KNIGHT, BISHOP, ROOK, KING, QUEEN
    export WHITEPAWN, WHITEKNIGHT, WHITEBISHOP, WHITEROOK, WHITEQUEEN, WHITEKING
    export BLACKPAWN, BLACKKNIGHT, BLACKBISHOP, BLACKROOK, BLACKQUEEN, BLACKKING

    export __KNIGHT_PROMO, __BISHOP_PROMO, __ROOK_PROMO, __QUEEN_PROMO, SEE_VALUES, PawnKingTable, evaluate, initEvalInfo

    export evaluate_knights, evaluate_bishops, evaluate_rooks, evaluate_queens, evaluate_kings, evaluate_pins, evaluate_threats, evaluate_initiative, evaluate_space, evaluate_passed

    export uci_main

end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
obs(env::CartPoleEnv, ::Nothing) = obs(env, zeros(4))
obs(env::CartPoleEnv, (x, x̄, θ, θ̄)) =
    Dict("x" => x, "theta"=>θ)

function Ctx(env::CartPoleEnv)
    path = (p) -> normpath("$(@__DIR__)/$p")
    s = Scope(imports=path.(["../../assets/cartpole/js/Board.js", "../../assets/cartpole/css/cartpole.css"]))
    config = Dict(
        "cart_height"=> env.length/60,
        "cart_length"=> 2*env.length,
        "pole_length"=> env.polemass_length,
        "pole_diameter"=> env.polemass_length/5,
        "x_threshold"=> env.x_threshold)

    o = Observable(s, "obs", obs(env, Flux.data(env.state)))
    onimport(s, @js () -> begin
        window.pick = (e) -> document.querySelector(e)
        window.container = window.pick(".wio-scope") || window.pick(".webio-scope")
        window.board = @new Board(window.container, $(config))
        board.render($o[])
    end)

    onjs(o, @js function (x)
        board.render(x)
    end)

    Ctx(s, o)
end

function render(env::CartPoleEnv, ctx::Ctx)
    ctx.o[] = obs(env, Flux.data(env.state))
end


text:
using JSON

struct XJld2 end
struct XJson end
struct XJld end

_decode_dict(D::Dict) = convert(Val(Symbol(D["__id__"])), D)

# ------------ writing and reading NBodyIP ------------

NBodyIP(D::Dict) = NBodyIP(_decode_dict.(D["components"]))

Dict(IP::NBodyIP) = Dict("__id__" => "NBodyIP",
                         "components" => Dict.(IP.components))

Base.convert(::Val{:NBodyIP}, D::Dict) = NBodyIP(D)

# ------------ writing and reading nothing ------------

Dict(::Nothing) = Dict("__id__" => "Void")

Base.convert(::Val{:Nothing}, D::Dict) = nothing

# ------------ FileIO -------------------

function _checkextension(fname)
   if fname[end-3:end] == "jld2"
      return XJld2()
   elseif fname[end-2:end] == "jld"
      return XJld()
   elseif fname[end-3:end] == "json"
      return XJson()
   end
   error("unknown extension for `$fname`")
end

save_ip_args(IP) = IP, Dict{String, Any}()
save_ip_args(IP, info) = IP, info

"""
`save_ip(fname, IP, info)`

where `IP` must be an `NBodyIP` and `info` a `Dict`
"""
save_ip(fname::AbstractString, args...) =
   save_ip(_checkextension(fname), fname, save_ip_args(args...)...)

"""
`load_ip(fname) -> IP, info`

where `IP`is an `NBodyIP` and `info` a `Dict`
"""
load_ip(fname::AbstractString) =
   load_ip(_checkextension(fname), fname)

save_ip(::Union{XJld2, XJld}, args...) = err_jld()
load_ip(::Union{XJld2, XJld}, args...) = err_jld()

err_jld() =
   @warn("""`load_ip` and `save_ip` do not directly support jld and jld2. In order
           to load an `NBodyIP` stored in one of those formats, please use `FileIO`,
           load the IP as a `ipdict::Dict` and then decode it using
           `NBodyIP(ipdict)`. To save an `ip::NBodyIP` as a jld or jld2 file,
           convert it to a `Dict` using `Dict(ip)` and then save it using
           `FileIO.save`.""")

function save_ip(::XJson, fname, IP, info)
   f = open(fname, "w")
   print(f, JSON.json(Dict("IP" => Dict(IP), "info" => info)))
   close(f)
end

function load_ip(::XJson, fname)
   IPj = JSON.parsefile(fname)
   if haskey(IPj, "__id__")
      @info("`load_ip`: Old IP filetype")
      @assert IPj["__id__"] == "NBodyIP"
      return NBodyIP(IPj), IPj["info"], Dict{String, Any}()
   end

   # new filetype
   @assert IPj["IP"]["__id__"] == "NBodyIP"
   IP = NBodyIP(IPj["IP"])
   if haskey(IPj, "info")
      info = IPj["info"]
   else
      info = Dict{String,Any}()
   end
   return IP, info
end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# 01-speech-blstm.jl
#
# See Graves & Schmidhuber ([Graves, A., &
# Schmidhuber, J. (2005). Framewise phoneme classification with
# bidirectional LSTM and other neural network architectures. Neural
# Networks, 18(5-6), 602-610.]).

using Flux
using Flux: crossentropy, softmax, flip, sigmoid, LSTM, @epochs
using BSON
using Random

# Paths to the training and test data directories
const TRAINDIR = "train"
const TESTDIR = "test"
const EPOCHS = 20

# Component layers of the bidirectional LSTM layer
forward = LSTM(26, 93)
backward = LSTM(26, 93)
output = Dense(186, 61)

"""
  BLSTM(x)

BLSTM layer using above LSTM layers

# Parameters
* **x** A 2-tuple containing the forward and backward time samples;
the first is from processing the sequence forward, and the second
is from processing it backward

# Returns
* The concatenation of the forward and backward LSTM predictions
"""
BLSTM(x) = vcat.(forward.(x), flip(backward, x))

"""
  model(x)

The chain of functions representing the trained model.

# Parameters
* **x** The utterance that the model should process

# Returns
* The model's predictions for each time step in `x`
"""
model(x) = softmax.(output.(BLSTM(x)))

"""
   loss(x, y)

Calculates the categorical cross-entropy loss for an utterance

# Parameters
* **x** Iterable containing the frames to classify
* **y** Iterable containing the labels corresponding to the frames
in `x`

# Returns
* The calculated loss value

# Side-effects
* Resets the state in the BLSTM layer
"""
function loss(x, y)
  l = sum(crossentropy.(model(x), y))
  Flux.reset!((forward, backward))
  return l
end

"""
  readData(dataDir)

Reads in the data contained in a specified directory

# Parameters
* **dataDir** String of the path to the directory containing the data

# Return
* **Xs** Vector where each element is a vector of the frames for
one utterance
* **Ys** A vector where each element is a vector of the labels for
the frames for one utterance
"""
function readData(dataDir)
  fnames = readdir(dataDir)

  Xs = Vector()
  Ys = Vector()

  for (i, fname) in enumerate(fnames)
    print(string(i) * "/" * string(length(fnames)) * "\r")
    BSON.@load joinpath(dataDir, fname) x y
    x = [x[i,:] for i in 1:size(x,1)]
    y = [y[:,i] for i in 1:size(y,2)]
    push!(Xs, x)
    push!(Ys, y)
  end

  return (Xs, Ys)
end

"""
  evaluateAccuracy(data)

Evaluates the accuracy of the model on a set of data; can be used
either for validation or test accuracy

# Parameters
* **data** An iterable of paired values where the first element is
all the frames for a single utterance, and the second is the
associated frame labels to compare the model's predictions against

# Returns
* The predicted accuracy value as a proportion of the number of
correct predictions over the total number of predictions made
"""
function evaluateAccuracy(data)
  correct = Vector()
  for (x, y) in data
    y = argmax.(y)
    ŷ = argmax.(model(x))
    Flux.reset!((forward, backward))
    append!(correct, [ŷ_n == y_n for (ŷ_n, y_n) in zip(ŷ, y)])
  end
  sum(correct) / length(correct)
end

function main()

  println("Loading files")
  Xs, Ys = readData(TRAINDIR)
  data = collect(zip(Xs, Ys))

  valData = data[1:184]
  data = data[185:end]

  # Begin training
  println("Beginning training")

  opt = Momentum(params((forward, backward, output)), 10.0^-5; ρ=0.9)

  i = 0

  @epochs EPOCHS begin

    i += 1

    shuffle!(data)
    valData = valData[shuffle(1:length(valData))]

    Flux.train!(loss, data, opt)

    BSON.@save "model_epoch$(i).bson" forward backward output

    print("Validating\r")
    val_acc = evaluateAccuracy(valData)
    println("Val acc. " * string(val_acc))
    println()
  end

  # Clean up some memory
  valData = nothing
  data = nothing
  Xs = nothing
  Ys = nothing
  GC.gc()

  # Test model
  print("Testing\r")
  Xs_test, Ys_test = readData(TESTDIR)
  test_data = collect(zip(Xs_test, Ys_test))
  test_acc = evaluateAccuracy(test_data)
  println("Test acc. " * string(test_acc))
  println()
end

main()

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
using Test, SegmentIntersections

@testset "Test event creation" begin
    e1 = Event(1,2)
    e2 = Event(0,5)
    e3 = Event(1,2)
    @test e1.point == Point(1,2)
    @test e2 < e1
    @test e1 == e3
end


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 2, Zadanie 1.
#

# assert ARGS
err_mess = "provide two arguments: [64, 32] [up, down, max_to_min, min_to_max]"
if length(ARGS) != 2
    println(err_mess)
    exit(1)
end

# name the arguments
fl_str = ARGS[1]
mode = ARGS[2]

# parse the arguments
if !in(fl_str, ["64", "32"]) || !in(mode, ["up", "down", "max_to_min", "min_to_max"])
    println(err_mess)
    exit(1)
end

fl = Float64

# apply the first argument
if fl_str== "32"
    global fl = Float32
end
# end assert

# x₄ and x₅ were altered
x = [2.718281828, -3.141592654, 1.414213562, 0.577215664, 0.301029995]
y = [1486.2497, 878366.9879, -22.37492, 4773714.647, 0.000185049]

function calc_scalar_product(x, y)
    S = fl(0)
    # (a), (b)
    if in(mode, ["up", "down"])
        # change the direction of the range depending on the provided program argument
        for i in (mode == "down" ? (5:-1:1) : (1:5))
            S = fl(S + fl(x[i] * y[i]))
        end
    # (c), (d)
    elseif in(mode, ["max_to_min", "min_to_max"])
        # save partial products
        mid_positive::Array{fl} = []
        mid_negative::Array{fl} = []
        # iterate over the vectors
        for i in (1:5)
            out = fl(x[i] * y[i])
            if out < fl(0)
                push!(mid_negative, out)
            else
                push!(mid_positive, out)
            end
        end

        # sort the partial products
        sort!(mid_negative, rev=(mode == "max_to_min"))
        sort!(mid_positive, rev=(mode == "max_to_min"))

        # add the partial products
        S_positive = fl(0)
        S_negative = fl(0)
        for i in (1:length(mid_negative))
            S_negative = fl(S_negative + mid_negative[i])
        end
        for i in (1:length(mid_positive))
            S_positive = fl(S_positive + mid_positive[i])
        end
        # finally, add the partial sums
        S = fl(S_positive + S_negative)
    end
    return S
end

# print out the results
println(calc_scalar_product(x,y))

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
"""
    Aerosol activation scheme, which includes:

  - mean hygroscopicity for each mode of the aerosol size distribution
  - critical supersaturation for each mode of the aerosol size distribution
  - maximum supersaturation
  - total number of particles actived
  - total mass of particles actived
"""
module AerosolActivation

import SpecialFunctions

import Thermodynamics
import CloudMicrophysics
import CLIMAParameters

const SF = SpecialFunctions
const TD = Thermodynamics
const AM = CloudMicrophysics.AerosolModel
const CO = CloudMicrophysics.Common
const CP = CLIMAParameters
const CP_planet = CLIMAParameters.Planet
const APS = CP.AbstractParameterSet

export mean_hygroscopicity_parameter
export max_supersaturation

export N_activated_per_mode
export M_activated_per_mode

export total_N_activated
export total_M_activated

"""
    coeff_of_curvature(param_set, T)

  - `param_set` - abstract set with Earth's parameters
  - `T` - air temperature

Returns a curvature coefficient.
"""
function coeff_of_curvature(param_set::APS, T::FT) where {FT <: Real}

    _molmass_water::FT = CP_planet.molmass_water(param_set)
    _gas_constant::FT = CP.gas_constant()
    _ρ_cloud_liq::FT = CP_planet.ρ_cloud_liq(param_set)
    _surface_tension::FT = CP_planet.surface_tension_coeff(param_set)

    return 2 * _surface_tension * _molmass_water / _ρ_cloud_liq /
           _gas_constant / T
end

"""
    mean_hygroscopicity_parameter(param_set, ad)

  - `param_set` - abstract set with Earth's parameters
  - `ad` - aerosol distribution struct

Returns a tuple of hygroscopicity parameters
(one tuple element for each aerosol size distribution mode).
The tuple is computed either as mass-weighted B parameters
(Abdul-Razzak and Ghan 2000)
or volume weighted kappa parameters (Petters and Kreidenweis 2007).
Implemented via a dispatch based on aerosol distribution mode type.
"""
function mean_hygroscopicity_parameter(
    param_set::APS,
    ad::AM.AerosolDistribution{NTuple{N, T}},
) where {N, T <: AM.Mode_B}

    _molmass_water = CP_planet.molmass_water(param_set)
    _ρ_cloud_liq = CP_planet.ρ_cloud_liq(param_set)

    return ntuple(length(ad.Modes)) do i

        mode_i = ad.Modes[i]

        nom = sum(1:(mode_i.n_components)) do j
            mode_i.mass_mix_ratio[j] *
            mode_i.dissoc[j] *
            mode_i.osmotic_coeff[j] *
            mode_i.soluble_mass_frac[j] / mode_i.molar_mass[j]
        end

        den = sum(1:(mode_i.n_components)) do j
            mode_i.mass_mix_ratio[j] / mode_i.aerosol_density[j]
        end

        nom / den * _molmass_water / _ρ_cloud_liq
    end
end
function mean_hygroscopicity_parameter(
    param_set::APS,
    ad::AM.AerosolDistribution{NTuple{N, T}},
) where {N, T <: AM.Mode_κ}

    return ntuple(length(ad.Modes)) do i

        modei = ad.Modes[i]
        sum(1:(modei.n_components)) do j
            modei.vol_mix_ratio[j] * modei.kappa[j]
        end
    end
end

"""
    critical_supersaturation(param_set, ad, T)

  - `param_set` - abstract set with Earth's parameters
  - `ad` - aerosol distribution struct
  - `T` - air temperature

Returns a tuple of critical supersaturations
(one tuple element for each aerosol size distribution mode).
"""
function critical_supersaturation(
    param_set::APS,
    ad::AM.AbstractAerosolDistribution,
    T::FT,
) where {FT <: Real}

    A::FT = coeff_of_curvature(param_set, T)
    hygro = mean_hygroscopicity_parameter(param_set, ad)

    return ntuple(length(ad.Modes)) do i
        2 / sqrt(hygro[i]) * (A / 3 / ad.Modes[i].r_dry)^(3 / 2)
    end
end

"""
    max_supersaturation(param_set, ad, T, p, w, q)

  - `param_set` - abstract set with Earth's parameters
  - `ad` - aerosol distribution struct
  - `T` - air temperature
  - `p` - air pressure
  - `w` - vertical velocity
  - `q` - phase partition

Returns the maximum supersaturation.
"""
function max_supersaturation(
    param_set::APS,
    ad::AM.AbstractAerosolDistribution,
    T::FT,
    p::FT,
    w::FT,
    q::TD.PhasePartition{FT},
) where {FT <: Real}

    _grav::FT = CP_planet.grav(param_set)
    _ρ_cloud_liq::FT = CP_planet.ρ_cloud_liq(param_set)

    ϵ::FT = 1 / CP_planet.molmass_ratio(param_set)
    R_m::FT = TD.gas_constant_air(param_set, q)
    cp_m::FT = TD.cp_m(param_set, q)

    L::FT = TD.latent_heat_vapor(param_set, T)
    p_vs::FT = TD.saturation_vapor_pressure(param_set, T, TD.Liquid())
    G::FT = CO.G_func(param_set, T, TD.Liquid()) / _ρ_cloud_liq

    # eq 11, 12 in Razzak et al 1998
    # but following eq 10 from Rogers 1975
    α::FT = L * _grav * ϵ / R_m / cp_m / T^2 - _grav / R_m / T
    γ::FT = R_m * T / ϵ / p_vs + ϵ * L^2 / cp_m / T / p

    A::FT = coeff_of_curvature(param_set, T)
    ζ::FT = 2 * A / 3 * sqrt(α * w / G)

    Sm = critical_supersaturation(param_set, ad, T)

    tmp::FT = sum(1:length(ad.Modes)) do i

        mode_i = ad.Modes[i]

        f::FT = 0.5 * exp(2.5 * (log(mode_i.stdev))^2)
        g::FT = 1 + 0.25 * log(mode_i.stdev)
        η::FT = (α * w / G)^(3 / 2) / (2 * pi * _ρ_cloud_liq * γ * mode_i.N)

        1 / (Sm[i])^2 *
        (f * (ζ / η)^(3 / 2) + g * (Sm[i]^2 / (η + 3 * ζ))^(3 / 4))
    end

    return FT(1) / sqrt(tmp)
end

"""
    N_activated_per_mode(param_set, ad, T, p, w, q)

  - `param_set` - abstract set with Earth's parameters
  - `ad` - aerosol distribution struct
  - `T` - air temperature
  - `p` - air pressure
  - `w` - vertical velocity
  - `q` - phase partition

Returns the number of activated aerosol particles
in each aerosol size distribution mode.
"""
function N_activated_per_mode(
    param_set::APS,
    ad::AM.AbstractAerosolDistribution,
    T::FT,
    p::FT,
    w::FT,
    q::TD.PhasePartition{FT},
) where {FT <: Real}

    smax::FT = max_supersaturation(param_set, ad, T, p, w, q)
    sm = critical_supersaturation(param_set, ad, T)

    return ntuple(length(ad.Modes)) do i

        mode_i = ad.Modes[i]
        u_i::FT = 2 * log(sm[i] / smax) / 3 / sqrt(2) / log(mode_i.stdev)

        mode_i.N * (1 / 2) * (1 - SF.erf(u_i))
    end
end

"""
    M_activated_per_mode(param_set, ad, T, p, w, q)

  - `param_set` - abstract set with Earth's parameters
  - `ad` - aerosol distribution struct
  - `T` - air temperature
  - `p` - air pressure
  - `w` - vertical velocity
  - `q` - phase partition

Returns the mass of activated aerosol particles
per mode of the aerosol size distribution.
"""
function M_activated_per_mode(
    param_set::APS,
    ad::AM.AbstractAerosolDistribution,
    T::FT,
    p::FT,
    w::FT,
    q::TD.PhasePartition{FT},
) where {FT <: Real}

    smax = max_supersaturation(param_set, ad, T, p, w, q)
    sm = critical_supersaturation(param_set, ad, T)

    return ntuple(length(ad.Modes)) do i

        mode_i = ad.Modes[i]

        avg_molar_mass_i = sum(1:(mode_i.n_components)) do j
            mode_i.molar_mass[j] * mode_i.mass_mix_ratio[j]
        end

        u_i = 2 * log(sm[i] / smax) / 3 / sqrt(2) / log(mode_i.stdev)

        avg_molar_mass_i * 1 / 2 *
        (1 - SF.erf(u_i - 3 * sqrt(2) / 2 * log(mode_i.stdev)))
    end
end

"""
    total_N_activated(param_set, ad, T, p, w, q)

  - `param_set` - abstract set with Earth's parameters
  - `ad` - aerosol distribution struct
  - `T` - air temperature
  - `p` - air pressure
  - `w` - vertical velocity
  - `q` - phase partition

Returns the total number of activated aerosol particles.
"""
function total_N_activated(
    param_set::APS,
    ad::AM.AbstractAerosolDistribution,
    T::FT,
    p::FT,
    w::FT,
    q::TD.PhasePartition{FT},
) where {FT <: Real}

    return sum(N_activated_per_mode(param_set, ad, T, p, w, q))

end

"""
    total_M_activated(param_set, ad, T, p, w, q)

  - `param_set` - abstract set with Earth's parameters
  - `ad` - aerosol distribution struct
  - `T` - air temperature
  - `p` - air pressure
  - `w` - vertical velocity
  - `q` - phase partition

Returns the total mass of activated aerosol particles.
"""
function total_M_activated(
    param_set::APS,
    ad::AM.AbstractAerosolDistribution,
    T::FT,
    p::FT,
    w::FT,
    q::TD.PhasePartition{FT},
) where {FT <: Real}

    return sum(M_activated_per_mode(param_set, ad, T, p, w, q))

end

end # module AerosolActivation.jl


text:
# ======================================
#    Row of data
# ======================================

"""
A `Row` stores multiple elements of data referenced by their column names.
"""
immutable Row{Names, Types <: Tuple}
    data::Types

    function Row{T <: Tuple}(data_in::T)
        check_row(Val{Names}, Types)
        new(convert(Types, data_in))
    end
end

@generated function check_row{Names, Types}(::Type{Val{Names}}, ::Type{Types})
    if !isa(Names, Tuple) || eltype(Names) != Symbol || length(Names) != length(unique(Names))
        str = "Row parameter 1 (Names) is expected to be a tuple of unique symbols, got $Names" # TODO: reinsert colons in symbols?
        return :(error($str))
    end
    if :Row ∈ Names
        return :( error("Field name cannot be :Row") )
    end
    N = length(Names)
    if length(Types.parameters) != N || reduce((a,b) -> a | !isa(b, DataType), false, Types.parameters)
        str = "Row parameter 2 (Types) is expected to be a Tuple{} of $N DataTypes, got $Types"
        return :(error($str))
    end

    return nothing
end

@compat @generated function (::Type{Row{Names}}){Names}(data::Tuple)
    if !isa(Names, Tuple) || eltype(Names) != Symbol || length(Names) != length(unique(Names))
        str = "Row parameter 1 (Names) is expected to be a tuple of unique symbols, got $Names" # TODO: reinsert colons in symbols?
        return :(error($str))
    end

    if length(Names) == length(data.parameters)
        return quote
            $(Expr(:meta,:inline))
            Row{Names, $data}(data)
        end
    else
        return :(error("Can't construct Row with $(length(Names)) columns with input $data"))
    end
end

@compat (::Type{Row{Names,Types_new}}){Names, Types, Types_new}(r::Row{Names,Types}) = convert(Row{Names,Types_new}, r)

@generated function Base.convert{Names, Names_new, Types, Types_new <: Tuple}(::Type{Row{Names_new,Types_new}}, r::Row{Names,Types})
    if !isa(Names_new, Tuple) || eltype(Names_new) != Symbol || length(Names_new) != length(Names) || length(Names_new) != length(unique(Names_new))
        str = "Cannot convert $(length(Names)) columns to new names $(Names_new)."
        return :(error($str))
    end
    if length(Types_new.parameters) != length(Types.parameters)
        str = "Cannot convert $(length(Types.pareters)) columns to $(length(Types_new.parameters)) new types."
        return :(error($str))
    end
    exprs = [:(convert(Types_new.parameters[$j], r.data[$j])) for j = 1:length(Names)]
    return Expr(:call, Row{Names_new,Types_new}, Expr(:tuple, exprs...))
end

@compat Base.:(==){Names, Types1, Types2}(row1::Row{Names, Types1}, row2::Row{Names, Types2}) = (row1.data == row2.data)
@compat @generated function Base.:(==){Names1, Types1, Names2, Types2}(row1::Row{Names1,Types1}, row2::Row{Names2,Types2})
    try
        order = permutator(Names1, Names2)
        expr = :( row1.data[$(order[1])] == row2.data[1] )
        for j = 2:length(Names1)
            expr = Expr(:call, :(&), expr, :( row1.data[$(order[j])] == row2.data[$j] ))
        end
        return expr
    catch
        return false
    end
end

rename{Names, NewNames}(row::Row{Names}, ::Type{Val{NewNames}}) = Row{NewNames}(row.data)

@generated function rename{Names, OldName, NewName}(row::Row{Names}, ::Type{Val{OldName}}, ::Type{Val{NewName}})
    j = columnindex(Names, OldName)

    NewNames = [Names...]
    NewNames[j] = NewName
    NewNames = (NewNames...)

    return :(Row{$NewNames}(row.data))
end

columnindex{N}(names::NTuple{N,Symbol}, name) = error("Can't search for columns $name")
function columnindex{N}(names::NTuple{N,Symbol}, name::Symbol)
    for i = 1:length(names)
        if names[i] == name
            return i
        end
    end
    error("Can't find column with name :$name")
end
function columnindex{N, M}(names::NTuple{N,Symbol}, searchnames::NTuple{M,Symbol})
    return ntuple(i->columnindex(names, searchnames[i]), M)
end

@inline Base.names{Names}(::Row{Names}) = Names
@inline Base.names{Names, Types <: Tuple}(::Type{Row{Names,Types}}) = Names
@inline Base.names{Names}(::Type{Row{Names}}) = Names
@inline eltypes{Names, Types <: Tuple}(::Row{Names,Types}) = Types
@inline eltypes{Names, Types <: Tuple}(::Type{Row{Names,Types}}) = Types

@inline nrow(::Row) = 1
@inline ncol{Names}(::Row{Names}) = length(Names)
@inline ncol{Names,Types}(::Type{Row{Names,Types}}) = length(Names)
@inline ncol{Names}(::Type{Row{Names}}) = length(Names)

@generated function Base.getindex{Names, GetName}(row::Row{Names}, ::Type{Val{GetName}})
    if isa(GetName, Symbol)
        j = columnindex(Names, GetName)

        return quote
            $(Expr(:meta, :inline))
            row.data[$j]
        end
    elseif isa(GetName, Tuple)
        inds = columnindex(Names, GetName)
        exprs = [:(row.data[$(inds[j])]) for j = 1:length(inds)]
        expr = Expr(:call, Row{GetName}, Expr(:tuple, exprs...))

        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    else
        str = "Can't get column(s) named $Name"
        return :(error($str))
    end
end

@inline Base.length(::Row) = 1
@inline Base.endof(::Row) = 1

Base.start(r::Row) = false
Base.next(r::Row, state) = (r, true)
Base.done(r::Row, state) = state

Base.getindex(r::Row) = r
Base.getindex(r::Row, i::Integer) = i == 1 ? r : error("Cannot index Row at $i")
Base.getindex(r::Row, ::Colon) = r

# reordering
@generated function permutecols{Names1,Names2,Types}(r::Row{Names1,Types}, ::Type{Val{Names2}})
    if Names1 == Names2
        return :(r)
    else
        if !(isa(Names2, Tuple)) || eltype(Names2) != Symbol || length(Names2) != length(Names1) || length(Names2) != length(unique(Names2))
            str = "New column names $Names2 do not match existing names $Names1"
            return :(error($str))
        end

        order = permutator(Names1, Names2)

        exprs = [:(r.data[$(order[j])]) for j = 1:length(Names1)]
        return Expr(:call, Row{Names2}, Expr(:tuple, exprs...))
    end
end

function permutator{N}(names1::NTuple{N,Symbol}, names2::NTuple{N,Symbol})
    order = zeros(Int, N)
    for i = 1:N
        isfound = false
        for j = 1:N
            if names1[i] == names2[j]
                isfound = true
                order[j] = i
                break
            end
        end

        if !isfound
            str = "New column names $names2 do not match existing names $names1"
            return :(error($str))
        end
    end

    return order
end

# Horizontally concatenate cells and rows into rows
@generated Base.hcat{Name}(c::Cell{Name}) = :( Row{$((Name,))}((c.data,)) )
Base.hcat(r::Row) = r

@generated function Base.hcat(r1::Union{Cell,Row}, r2::Union{Cell,Row})
    names1 = (r1 <: Cell ? (name(r1),) : names(r1))
    names2 = (r2 <: Cell ? (name(r2),) : names(r2))

    if length(intersect(names1, names2)) != 0
        str = "Column names are not distinct. Got $names1 and $names2"
        return :(error($str))
    end

    newnames = (names1..., names2...)
    exprs = vcat([:(r1.data[$j]) for j = 1:length(names1)], [:(r2.data[$j]) for j = 1:length(names2)])

    return Expr(:call, Row{newnames}, Expr(:tuple, exprs...))
end

Base.hcat(r1::Union{Cell,Row}, r2::Union{Cell,Row}, rs::Union{Cell,Row}...) = hcat(hcat(r1, r2), rs...)

# copy
@generated function Base.copy{Names}(r::Row{Names})
    exprs = [:(copy(r.data[$j])) for j = 1:length(Names)]
    return Expr(:call, Row{Names}, Expr(:tuple, exprs...))
end

macro Row(exprs...)
    N = length(exprs)
    names = Vector{Any}(N)
    values = Vector{Any}(N)
    for i = 1:N
        expr = exprs[i]
        if expr.head != :(=) && expr.head != :(kw) # strange Julia bug, see issue 7669
            error("A Expecting expression like @Row(name1::Type1 = value1, name2::Type2 = value2) or @Row(name1 = value1, name2 = value2)")
        end
        if isa(expr.args[1],Symbol)
            names[i] = (expr.args[1])
            values[i] = esc(expr.args[2])
        elseif isa(expr.args[1],Expr)
            if expr.args[1].head != :(::) || length(expr.args[1].args) != 2
                error("A Expecting expression like @Row(name1::Type1 = value1, name2::Type2 = value2) or @Row(name1 = value1, name2 = value2)")
            end
            names[i] = (expr.args[1].args[1])
            values[i] = esc(Expr(:call, :convert, expr.args[1].args[2], expr.args[2]))
        else
            error("A Expecting expression like @Row(name1::Type1 = value1, name2::Type2 = value2) or @Row(name1 = value1, name2 = value2)")
        end
    end

    rowtype = TypedTables.Row{(names...)}
    return Expr(:call, rowtype, Expr(:tuple, values...))
end

#=
# Some interrogation
@inline Base.names{Index,DataTypes}(row::Union{Row{Index,DataTypes},Type{Row{Index,DataTypes}}}) = names(Index)
@inline eltypes{Index,DataTypes}(row::Union{Row{Index,DataTypes},Type{Row{Index,DataTypes}}}) = DataTypes
@inline index{Index,DataTypes}(row::Union{Row{Index,DataTypes},Type{Row{Index,DataTypes}}}) = Index
@generated Base.length{Index,DataTypes}(row::Union{Row{Index,DataTypes},Type{Row{Index,DataTypes}}}) = :($(length(Index)))
Base.endof(row::Row) = length(row)
@generated ncol{Index,DataTypes}(row::Union{Row{Index,DataTypes},Type{Row{Index,DataTypes}}}) = :($(length(Index)))
nrow{Index,DataTypes}(row::Union{Row{Index,DataTypes},Type{Row{Index,DataTypes}}}) = 1

samefields{Index1,Index2}(::Row{Index1},::Row{Index2}) = samefields(Index1,Index2)
samefields{Index1<:FieldIndex,Index2}(::Index1,::Row{Index2}) = samefields(Index1(),Index2)
samefields{Index1,Index2<:FieldIndex}(::Row{Index1},::Index2) = samefields(Index1,Index2())

=={Index,ElTypes}(row1::Row{Index,ElTypes},row2::Row{Index,ElTypes}) = (row1.data == row2.data)
function =={Index1,ElTypes1,Index2,ElTypes2}(row1::Row{Index1,ElTypes1},row2::Row{Index2,ElTypes2})
    if samefields(Index1,Index2)
        idx = Index2[Index1]
        for i = 1:length(Index1)
            if row1.data[i] != row2.data[idx[i]]
                return false
            end
        end
        return true
    else
        return false
    end
end

rename{Index,DataTypes}(row::Row{Index,DataTypes}, new_names::FieldIndex) = Row(rename(Index,new_names),row.data)
rename{Index,DataTypes}(row::Row{Index,DataTypes}, old_names::Union{FieldIndex,Field}, new_names::Union{FieldIndex,Field}) = Row(rename(Index,old_names,new_names),row.data)

#function Base.show{Index,DataTypes}(io::IO,row::Row{Index,DataTypes})
#    print(io,"Row(")
#    for i = 1:length(Index)
#        print(io,"$(name(Index[i]))=$(compactstring(row.data[i]))")
#        if i < length(Index)
#            print(io,", ")
#        end
#    end
#    print(io,")")
#end

# Can index with integers or fields
# A single field returns a scalar value, an integer a cell, a FieldIndex or other collection returns a smaller Row.
@inline Base.getindex{Index,DataTypes}(row::Row{Index,DataTypes},i::Field) = row.data[Index[i]] #::DataTypes.parameters[i] # Is this considered "type safe"??
@generated function Base.getindex{R<:Row,I<:FieldIndex}(row::R,idx::I)
    tmp = ntuple(i->:(row.data[$(index(R)[I()[i]])]), length(I))
    tmp2 = Expr(:tuple,tmp...)
    return :(idx($tmp2))
end

# Getting indices (immutable, so no setindex)
@generated function Base.getindex{Index,DataTypes,I}(row::Row{Index,DataTypes},::Type{Val{I}})
    if isa(I, Int)
        return :( $(Cell{Index[I],DataTypes.parameters[I]})(row.data[$I]) )
    elseif isa(I, Symbol)
        return :( row.data[$(Index[Val{I}])] )
    elseif isa(I, Tuple)
        l_I = length(I)
        if isa(I, NTuple{l_I, Int})
            expr = Expr(:tuple, ntuple(i-> :(row.data[$(I[i])]), l_I)...)
            return :( Row($(Index[Val{I}]), $expr) )
        elseif isa(I, NTuple{l_I, Symbol})
            expr = Expr(:tuple, ntuple(i -> :(row.data[$(Index[Val{I[i]}])]), l_I)...)
            return :( Row($(Index[Val{Index[Val{I}]}]), $expr) )
        else
            str = "Can't index Row with fields $Fields with a Val{$I}"
            return :(error($str))
        end
    else # e.g. UnitRange{Int} and other methods of indexing a Tuple
        str = "Can't index Row with fields $Fields with a Val{$I}"
        return :(error($str))
    end
end

#@inline Base.getindex{Index,DataTypes}(row::Row{Index,DataTypes},i) = Index[i](row.data[i])
@inline Base.getindex{Index,DataTypes}(row::Row{Index,DataTypes},::Colon) = row

# TODO setindex broadcasts to `[] =` (assume its like a zero-dimensional sub-array) or use SubRow class

# copies
Base.copy{Index,DataTypes}(row::Row{Index,DataTypes}) = Row{Index,DataTypes}(copy(row.data))
Base.deepcopy{Index,DataTypes}(row::Row{Index,DataTypes}) = Row{Index,DataTypes}(deepcopy(row.data))

# Concatenate cells and rows into rows
# Generated functions appear to be needed for speed...
Base.hcat{F,ElType}(cell::Cell{F,ElType}) = Row(F,cell.data)
Base.hcat{Index,ElTypes}(row::Row{Index,ElTypes}) = Row{Index,ElTypes}(row.data)

@generated Base.hcat{F1,ElType1,F2,ElType2}(cell1::Cell{F1,ElType1},cell2::Cell{F2,ElType2}) = :( Row{$(FieldIndex{(F1,F2)}()),$(Tuple{ElType1,ElType2})}((cell1.data,cell2.data)) )
@generated function Base.hcat{F1,ElType1,Index2,ElTypes2}(cell1::Cell{F1,ElType1},row2::Row{Index2,ElTypes2})
    Index = F1 + Index2
    ElTypes = eltypes(Index)
    :(Row{$Index,$ElTypes}((cell1.data, row2.data...)) )
end
@generated function Base.hcat{Index1,ElTypes1,F2,ElType2}(row1::Row{Index1,ElTypes1},cell2::Cell{F2,ElType2})
    Index = Index1 + F2
    ElTypes = eltypes(Index)
    :(Row{$Index,$ElTypes}((row1.data..., cell2.data)) )
end
@generated function Base.hcat{Index1,ElTypes1,Index2,ElTypes2}(row1::Row{Index1,ElTypes1},row2::Row{Index2,ElTypes2})
    Index = Index1 + Index2
    ElTypes = eltypes(Index)
    :(Row{$Index,$ElTypes}((row1.data..., row2.data...)) )
end

@inline Base.hcat(c1::Union{Cell,Row},c2::Union{Cell,Row},cs::Union{Cell,Row}...) = hcat(hcat(c1,c2),cs...) # Splatting is expensive but it works.

macro row(exprs...)
    N = length(exprs)
    field = Vector{Any}(N)
    value = Vector{Any}(N)
    for i = 1:N
        expr = exprs[i]
        if expr.head != :(=) && expr.head != :(kw) # strange Julia bug, see issue 7669
            error("A Expecting expression like @row(name1::Type1 = value1, name2::Type2 = value2) or @row(field1 = value1, field2 = value2)")
        end
        if isa(expr.args[1],Symbol)
            field[i] = expr.args[1]
        elseif isa(expr.args[1],Expr)
            if expr.args[1].head != :(::) || length(expr.args[1].args) != 2
                error("A Expecting expression like @row(name1::Type1 = value1, name2::Type2 = value2) or @row(field1 = value1, field2 = value2)")
            end
            field[i] = :(TypedTables.Field{$(Expr(:quote,expr.args[1].args[1])),$(expr.args[1].args[2])}())
        else
            error("A Expecting expression like @row(name1::Type1 = value1, name2::Type2 = value2) or @row(field1 = value1, field2 = value2)")
        end
        value[i] = expr.args[2]
    end

    fields = Expr(:tuple,field...)
    values = Expr(:tuple,value...)

    return :(TypedTables.Row(TypedTables.FieldIndex{$(esc(fields))}(),$(esc(values))))
end
=#


text:
const SOL_SOCKET = 1
const SO_PROTOCOL = 38
const SO_DETACH_BPF = 27
const SO_ATTACH_BPF = 50
function getsockopt(sock, level, optname, T)
    optval_ref = Ref{T}(zero(T))
    optlen_ref = Ref{Cint}(sizeof(T))
    GC.@preserve optval_ref begin
    ret = ccall(:getsockopt, Cint,
                             (Cint, Cint, Cint, Ptr{Cvoid}, Ptr{Cint}),
                             Base._fd(sock), level, optname, optval_ref, optlen_ref)
    end
    ret != 0 && Base.systemerror(ret)
    optval_ref[]
end
function setsockopt(sock, level, optname, optval)
    optval_ref = Ref{typeof(optval)}(optval)
    GC.@preserve optval_ref begin
    ret = ccall(:setsockopt, Cint,
                             (Cint, Cint, Cint, Ptr{Cvoid}, Cint),
                             Base._fd(sock), level, optname, optval_ref, sizeof(optval))
    end
    ret != 0 && Base.systemerror(ret)
    ret
end