text:
################################################
#
# Calibration exercise to create balanced SAM
#
################################################

using SLiDE
using CSV
using JuMP
using DataFrames
using Ipopt

##########
# FUNCTIONS
##########

#replace here with "collect"
  function key_to_vec(d::Dict,index_num::Int64)
    return [k[index_num] for k in keys(d)]
  end

  function df_to_dict(df::DataFrame,remove_column::Symbol,value_column::Symbol,year::Int64,sub::Bool)
    colnames = setdiff(names(df),[value_column,remove_column])
    #subset on year
    df = df[df[!,:yr].==year,:]

    #filter out negative values
    if sub
      df[df[!,value_column].<0,value_column] .= 0
    end

    #only want non-zero values
    df = df[df[!,value_column].!=0.0,:]

    #df[!,value_column] .= round.(df[!,value_column],digits=3)

    if length(colnames) == 1
        return Dict((row[colnames]...)=>row[value_column] for row in eachrow(df))
    end

    if length(colnames) > 1
        return Dict(tuple(row[colnames]...)=>row[value_column] for row in eachrow(df))
    end

  end

##################
# -- Load Data --
##################

mod_year = 2016;

symbols_calibration = [:y0,:ys0,:fs0,:id0,:fd0,:va0,:m0,:x0,:ms0,:md0,:a0,:ta0,:tm0];

iot = copy(io)
io[:x0] = io[:x0][[:yr,:i,:value]]
io[:m0] = io[:m0][[:yr,:i,:value]]
io[:fs0] = io[:fs0][[:yr,:i,:value]]
io[:ys0] = io[:ys0][[:yr,:j,:i,:value]]

#need to set all negative pce values in :fd0 to zero
# but others can remain
kk = copy(io[:fd0])
k1 = kk[(kk[:fd].=="pce"),:]
k1[(k1[!,:value].<0),:value] .= 0
k2 = kk[(kk[!,:fd].!="pce"),:]

kkt = [k1;k2]
io[:fd0] = kkt

cal = Dict();

skpneg = [:fd0,:ta0, :tm0]

for i in symbols_calibration
  if !(i in skpneg)
    push!(cal,i=>df_to_dict(io[i],:yr,:value,mod_year,true));
  end

  if i in skpneg
    push!(cal,i=>df_to_dict(io[i],:yr,:value,mod_year,false));
  end
end

i_set = unique(io[:y0][!,:i]);
j_set = copy(i_set);

fd_set = ["pce","equipment","intelprop","residential","changinv","structures","def_equipment","defense",
          "def_intelprop","def_structures","nondefense","fed_equipment","fed_intelprop","fed_structures",
          "state_consume","state_equipment","state_intelprop","state_invest"];

ts_set = ["taxes","subsidies"];

va_set = ["compen","surplus","othtax"];

m_set = ["trn","trd"];

output_use_set = ["use","oth"];

#fill_zero not working for single dimensions..
for i in i_set
  if haskey(cal[:ta0],i)==false
    push!(cal[:ta0],i=>0)
  end
  if haskey(cal[:tm0],i)==false
    push!(cal[:tm0],i=>0)
  end
end

####################
# -- Calibration --
####################

calib = Model(with_optimizer(Ipopt.Optimizer,nlp_scaling_method="gradient-based"))

@variable(calib,ys0_est[j in j_set,i in i_set]>=0,start=0);
@variable(calib,fs0_est[i in i_set]>=0,start=0);
@variable(calib,ms0_est[i in i_set,m in m_set]>=0,start=0);
@variable(calib,y0_est[i in i_set]>=0,start=0);
@variable(calib,id0_est[i in i_set,j in j_set]>=0,start=0);
@variable(calib,fd0_est[i in i_set,fd in fd_set]>=0,start=0);
@variable(calib,va0_est[va in va_set,j in j_set]>=0,start=0);
@variable(calib,a0_est[i in i_set]>=0,start=0);
@variable(calib,x0_est[i in i_set]>=0,start=0);
@variable(calib,m0_est[i in i_set]>=0,start=0);
@variable(calib,md0_est[m in m_set,i in i_set]>=0,start=0);

```
mkt_py(i)..	sum(j, ys0_(j,i)) +  fs0_(i) =E= sum(m, ms0_(i,m)) + y0_(i);
mkt_pa(i)..	a0_(i) =E= sum(j, id0_(i,j)) + sum(fd,fd0_(i,fd));
mkt_pm(m)..	sum(i,ms0_(i,m)) =E= sum(i, md0_(m,i));
prf_y(j)..	sum(i, ys0_(j,i)) =E= sum(i, id0_(i,j)) + sum(va,va0_(va,j));
prf_a(i)..	a0_(i)*(1-ta0(i)) + x0_(i) =E= y0_(i) + m0_(i)*(1+tm0(i)) + sum(m, md0_(m,i));
```

@constraint(calib,mkt_py[i in i_set],
  sum(ys0_est[j,i] for j in j_set) + fs0_est[i] == sum(ms0_est[i,m] for m in m_set) + y0_est[i]
);

@constraint(calib,mkt_pa[i in i_set],
  a0_est[i] == sum(id0_est[i,j] for j in j_set) + sum(fd0_est[i,fd] for fd in fd_set)
);

@constraint(calib,mkt_pm[m in m_set],
  sum(ms0_est[i,m] for i in i_set) == sum(md0_est[m,i] for i in i_set)
);

@constraint(calib,prf_y[j in j_set],
  sum(ys0_est[j,i] for i in i_set) == sum(id0_est[i,j] for i in i_set) + sum(va0_est[va,j] for va in va_set)
);

@constraint(calib,prf_a[i in i_set],
  a0_est[i] * (1-cal[:ta0][i]) + x0_est[i] == y0_est[i] + m0_est[i]*(1+cal[:tm0][i]) + sum(md0_est[m,i] for m in m_set)
);

penalty_nokey = 1e4

@objective(calib,Min,
  + sum(abs(cal[:ys0][j,i]) * (ys0_est[j,i] / cal[:ys0][j,i] - 1)^2 for i in i_set for j in j_set if haskey(cal[:ys0],(j,i))  )
  + sum(abs(cal[:id0][i,j]) * (id0_est[i,j] / cal[:id0][i,j] - 1)^2 for i in i_set for j in j_set if haskey(cal[:id0],(i,j)) )
  + sum(abs(cal[:fs0][i]) * (fs0_est[i] / cal[:fs0][i] - 1)^2 for i in i_set if haskey(cal[:fs0],i) )
  + sum(abs(cal[:ms0][i,m]) * (ms0_est[i,m] / cal[:ms0][i,m] - 1)^2 for i in i_set for m in m_set if haskey(cal[:ms0],(i,m)) )
  + sum(abs(cal[:y0][i]) * (y0_est[i] / cal[:y0][i] - 1)^2 for i in i_set if haskey(cal[:y0],i) )
  + sum(abs(cal[:fd0][i,fd]) * (fd0_est[i,fd] / cal[:fd0][i,fd] - 1)^2 for i in i_set for fd in fd_set if haskey(cal[:fd0],(i,fd)) )
  + sum(abs(cal[:va0][va,j]) * (va0_est[va,j] / cal[:va0][va,j] - 1)^2 for va in va_set for j in j_set if haskey(cal[:va0],(va,j)) )
  + sum(abs(cal[:a0][i]) * (a0_est[i] / cal[:a0][i] - 1)^2  for i in i_set if haskey(cal[:a0],i) )
  + sum(abs(cal[:x0][i]) * (x0_est[i] / cal[:x0][i] - 1)^2  for i in i_set if haskey(cal[:x0],i) )
  + sum(abs(cal[:m0][i]) * (m0_est[i] / cal[:m0][i] - 1)^2  for i in i_set if haskey(cal[:m0],i) )
  + sum(abs(cal[:md0][m,i]) * (md0_est[m,i] / cal[:md0][m,i] - 1)^2 for m in m_set for i in i_set if haskey(cal[:md0],(m,i)) )

+  penalty_nokey * (
  + sum(ys0_est[j,i]^2 for i in i_set for j in j_set if !haskey(cal[:ys0],(j,i)) )
  + sum(id0_est[i,j]^2  for i in i_set for j in j_set if !haskey(cal[:id0],(i,j)) )
  + sum(fs0_est[i]^2  for i in i_set if !haskey(cal[:fs0],i) )
  + sum(ms0_est[i,m]^2  for i in i_set for m in m_set if !haskey(cal[:ms0],(i,m)) )
  + sum(y0_est[i]^2 for i in i_set if !haskey(cal[:y0],i) )
  + sum(fd0_est[i,fd]^2 for i in i_set for fd in fd_set if !haskey(cal[:fd0],(i,fd)) )
  + sum(va0_est[va,j]^2 for va in va_set for j in j_set if !haskey(cal[:va0],(va,j)) )
  + sum(a0_est[i]^2 for i in i_set if !haskey(cal[:a0],i) )
  + sum(x0_est[i]^2 for i in i_set if !haskey(cal[:x0],i) )
  + sum(m0_est[i]^2 for i in i_set if !haskey(cal[:m0],i) )
  + sum(md0_est[m,i]^2  for m in m_set for i in i_set if !haskey(cal[:md0],(m,i)) )
  )
);

[set_start_value(ys0_est[j,i],cal[:ys0][j,i]) for i in i_set for j in j_set if haskey(cal[:ys0],(j,i)) ] ;
[set_start_value(id0_est[i,j],cal[:id0][i,j]) for i in i_set for j in j_set if haskey(cal[:id0],(i,j)) ] ;
[set_start_value(fs0_est[i],cal[:fs0][i]) for i in i_set if haskey(cal[:fs0],i) ];
[set_start_value(ms0_est[i,m],cal[:ms0][i,m]) for i in i_set for m in m_set if haskey(cal[:ms0],(i,m)) ];
[set_start_value(y0_est[i],cal[:y0][i]) for i in i_set if haskey(cal[:y0],i) ];
[set_start_value(fd0_est[i,fd],cal[:fd0][i,fd]) for i in i_set for fd in fd_set if haskey(cal[:fd0],(i,fd)) ];
[set_start_value(va0_est[va,j],cal[:va0][va,j]) for va in va_set for j in j_set if haskey(cal[:va0],(va,j)) ];
[set_start_value(a0_est[i],cal[:a0][i]) for i in i_set if haskey(cal[:a0],(i)) ];
[set_start_value(x0_est[i],cal[:x0][i]) for i in i_set if haskey(cal[:x0],i) ];
[set_start_value(m0_est[i],cal[:m0][i]) for i in i_set if haskey(cal[:m0],i) ];
[set_start_value(md0_est[m,i],cal[:md0][m,i]) for m in m_set for i in i_set if haskey(cal[:md0],(m,i)) ] ;

# multipliers for lower and upper bound relative
# to each respective variables reference parameter
lb = 0.1
ub = 5

[set_lower_bound(ys0_est[j,i],max(0,lb * cal[:ys0][j,i] )) for i in i_set for j in j_set if haskey(cal[:ys0],(j,i)) ] ;
[set_lower_bound(id0_est[i,j],max(0,lb * cal[:id0][i,j] )) for i in i_set for j in j_set if haskey(cal[:id0],(i,j)) ] ;
[set_lower_bound(fs0_est[i],max(0,lb * cal[:fs0][i])) for i in i_set if haskey(cal[:fs0],i) ];
[set_lower_bound(ms0_est[i,m],max(0,lb * cal[:ms0][i,m])) for i in i_set for m in m_set if haskey(cal[:ms0],(i,m)) ];
[set_lower_bound(y0_est[i],max(0,lb * cal[:y0][i])) for i in i_set if haskey(cal[:y0],i) ];
[set_lower_bound(fd0_est[i,fd],max(0,lb * cal[:fd0][i,fd])) for i in i_set for fd in fd_set if haskey(cal[:fd0],(i,fd)) ];
[set_lower_bound(va0_est[va,j],max(0,lb * cal[:va0][va,j])) for va in va_set for j in j_set if haskey(cal[:va0],(va,j)) ];
[set_lower_bound(a0_est[i],max(0,lb * cal[:a0][i])) for i in i_set if haskey(cal[:a0],(i)) ];
[set_lower_bound(x0_est[i],max(0,lb * cal[:x0][i])) for i in i_set if haskey(cal[:x0],i) ];
[set_lower_bound(m0_est[i],max(0,lb * cal[:m0][i])) for i in i_set if haskey(cal[:m0],i) ];
[set_lower_bound(md0_est[m,i],max(0,lb * cal[:md0][m,i])) for m in m_set for i in i_set if haskey(cal[:md0],(m,i)) ] ;

[set_upper_bound(ys0_est[j,i],abs(ub * cal[:ys0][j,i] )) for i in i_set for j in j_set if haskey(cal[:ys0],(j,i)) ] ;
[set_upper_bound(id0_est[i,j],abs(ub * cal[:id0][i,j] )) for i in i_set for j in j_set if haskey(cal[:id0],(i,j)) ] ;
[set_upper_bound(fs0_est[i],abs(ub * cal[:fs0][i])) for i in i_set if haskey(cal[:fs0],i) ];
[set_upper_bound(ms0_est[i,m],abs(ub * cal[:ms0][i,m])) for i in i_set for m in m_set if haskey(cal[:ms0],(i,m)) ];
[set_upper_bound(y0_est[i],abs(ub * cal[:y0][i])) for i in i_set if haskey(cal[:y0],i) ];
[set_upper_bound(fd0_est[i,fd],abs(ub * cal[:fd0][i,fd])) for i in i_set for fd in fd_set if haskey(cal[:fd0],(i,fd)) ];
[set_upper_bound(va0_est[va,j],abs(ub * cal[:va0][va,j])) for va in va_set for j in j_set if haskey(cal[:va0],(va,j)) ];
[set_upper_bound(a0_est[i],abs(ub * cal[:a0][i])) for i in i_set if haskey(cal[:a0],(i)) ];
[set_upper_bound(x0_est[i],abs(ub * cal[:x0][i])) for i in i_set if haskey(cal[:x0],i) ];
[set_upper_bound(m0_est[i],abs(ub * cal[:m0][i])) for i in i_set if haskey(cal[:m0],i) ];
[set_upper_bound(md0_est[m,i],abs(ub * cal[:md0][m,i])) for m in m_set for i in i_set if haskey(cal[:md0],(m,i)) ] ;

[fix(fs0_est[i],cal[:fs0][i],force=true) for i in i_set if haskey(cal[:fs0],i)];
[fix(fs0_est[i],0,force=true) for i in i_set if !haskey(cal[:fs0],i)];

[fix(va0_est[va,j],cal[:va0][va,j],force=true) for va in va_set for j in j_set if haskey(cal[:va0],(va,j))];
[fix(va0_est[va,j],0,force=true) for va in va_set for j in j_set if !haskey(cal[:va0],(va,j))];

[fix(m0_est[i],cal[:m0][i],force=true) for i in i_set if haskey(cal[:m0],i)];
[fix(m0_est[i],0,force=true) for i in i_set if !haskey(cal[:m0],i)];

#fd_temp = filter!(x->x≠"pce",fd_set)
#[fix(fd0_est[i,fd],cal[:fd0][i,fd],force=true) for i in i_set for fd in fd_temp if haskey(cal[:fd0],(i,fd))]
#[fix(fd0_est[i,fd],0,force=true) for i in i_set for fd in fd_temp if !haskey(cal[:fd0],(i,fd))]

[fix(ys0_est[j,i],0,force=true) for j in output_use_set for i in i_set]

JuMP.optimize!(calib)

##################
# -- Reporting --
##################

Obj_Diagnostics = Dict();

Obj_Diagnostics[:Z] = - objective_value(calib);

Obj_Diagnostics[:ys0_withkey] = sum(abs(cal[:ys0][j,i]) * (JuMP.value(ys0_est[j,i]) / cal[:ys0][j,i] - 1)^2 for i in i_set for j in j_set if haskey(cal[:ys0],(j,i))  ) ;
Obj_Diagnostics[:id0_withkey] = sum(abs(cal[:id0][i,j]) * (JuMP.value(id0_est[i,j]) / cal[:id0][i,j] - 1)^2 for i in i_set for j in j_set if haskey(cal[:id0],(i,j)) ) ;
Obj_Diagnostics[:fs0_withkey] = sum(abs(cal[:fs0][i]) * (JuMP.value(fs0_est[i]) / cal[:fs0][i] - 1)^2 for i in i_set if haskey(cal[:fs0],i) ) ;
Obj_Diagnostics[:ms0_withkey] = sum(abs(cal[:ms0][i,m]) * (JuMP.value(ms0_est[i,m]) / cal[:ms0][i,m] - 1)^2 for i in i_set for m in m_set if haskey(cal[:ms0],(i,m)) ) ;
Obj_Diagnostics[:y0_withkey] = sum(abs(cal[:y0][i]) * (JuMP.value(y0_est[i]) / cal[:y0][i] - 1)^2 for i in i_set if haskey(cal[:y0],i) ) ;
Obj_Diagnostics[:fd0_withkey] = sum(abs(cal[:fd0][i,fd]) * (JuMP.value(fd0_est[i,fd]) / cal[:fd0][i,fd] - 1)^2 for i in i_set for fd in fd_set if haskey(cal[:fd0],(i,fd)) ) ;
Obj_Diagnostics[:va0_withkey] = sum(abs(cal[:va0][va,j]) * (JuMP.value(va0_est[va,j]) / cal[:va0][va,j] - 1)^2 for va in va_set for j in j_set if haskey(cal[:va0],(va,j)) ) ;
Obj_Diagnostics[:a0_withkey] = sum(abs(cal[:a0][i]) * (JuMP.value(a0_est[i]) / cal[:a0][i] - 1)^2  for i in i_set if haskey(cal[:a0],i) ) ;
Obj_Diagnostics[:x0_withkey] = sum(abs(cal[:x0][i]) * (JuMP.value(x0_est[i]) / cal[:x0][i] - 1)^2  for i in i_set if haskey(cal[:x0],i) ) ;
Obj_Diagnostics[:m0_withkey] = sum(abs(cal[:m0][i]) * (JuMP.value(m0_est[i]) / cal[:m0][i] - 1)^2  for i in i_set if haskey(cal[:m0],i) ) ;
Obj_Diagnostics[:md0_withkey] = sum(abs(cal[:md0][m,i]) * (JuMP.value(md0_est[m,i]) / cal[:md0][m,i] - 1)^2 for m in m_set for i in i_set if haskey(cal[:md0],(m,i)) )  ;

Obj_Diagnostics[:ys0_nokey] = penalty_nokey * sum(JuMP.value(ys0_est[j,i]) for i in i_set for j in j_set if !haskey(cal[:ys0],(j,i)) ) ;
Obj_Diagnostics[:id0_nokey] = penalty_nokey * sum(JuMP.value(id0_est[i,j])  for i in i_set for j in j_set if !haskey(cal[:id0],(i,j)) ) ;
Obj_Diagnostics[:fs0_nokey] = penalty_nokey * sum(JuMP.value(fs0_est[i])  for i in i_set if !haskey(cal[:fs0],i) ) ;
Obj_Diagnostics[:ms0_nokey] = penalty_nokey * sum(JuMP.value(ms0_est[i,m])  for i in i_set for m in m_set if !haskey(cal[:ms0],(i,m)) ) ;
Obj_Diagnostics[:y0_nokey] = penalty_nokey * sum(JuMP.value(y0_est[i]) for i in i_set if !haskey(cal[:y0],i) ) ;
Obj_Diagnostics[:fd0_nokey] = penalty_nokey * sum(JuMP.value(fd0_est[i,fd]) for i in i_set for fd in fd_set if !haskey(cal[:fd0],(i,fd)) ) ;
Obj_Diagnostics[:va0_nokey] = penalty_nokey * sum(JuMP.value(va0_est[va,j]) for va in va_set for j in j_set if !haskey(cal[:va0],(va,j)) ) ;
Obj_Diagnostics[:a0_nokey] = penalty_nokey * sum(JuMP.value(a0_est[i]) for i in i_set if !haskey(cal[:a0],i) ) ;
Obj_Diagnostics[:x0_nokey] = penalty_nokey * sum(JuMP.value(x0_est[i]) for i in i_set if !haskey(cal[:x0],i) ) ;
Obj_Diagnostics[:m0_nokey] = penalty_nokey * sum(JuMP.value(m0_est[i]) for i in i_set if !haskey(cal[:m0],i) ) ;
Obj_Diagnostics[:md0_nokey] = penalty_nokey * sum(JuMP.value(md0_est[m,i])  for m in m_set for i in i_set if !haskey(cal[:md0],(m,i)) )  ;

Obj_Diagnostics[:difference] = sum(values(Obj_Diagnostics)) ;

#############################
# -- Variable diagnostics --
#############################

variable_diagnostics = Dict()

variable_diagnostics[:ys0_min] = minimum(JuMP.value.(ys0_est))
variable_diagnostics[:id0_min] = minimum(JuMP.value.(id0_est))
variable_diagnostics[:fs0_min] = minimum(JuMP.value.(fs0_est))
variable_diagnostics[:ms0_min] = minimum(JuMP.value.(ms0_est))
variable_diagnostics[:y0_min] = minimum(JuMP.value.(y0_est))
variable_diagnostics[:fd0_min] = minimum(JuMP.value.(fd0_est))
variable_diagnostics[:va0_min] = minimum(JuMP.value.(va0_est))
variable_diagnostics[:a0_min] = minimum(JuMP.value.(a0_est))
variable_diagnostics[:x0_min] = minimum(JuMP.value.(x0_est))
variable_diagnostics[:m0_min] = minimum(JuMP.value.(m0_est))
variable_diagnostics[:md0_min] = minimum(JuMP.value.(md0_est))

variable_diagnostics[:ys0_max] = maximum(JuMP.value.(ys0_est))
variable_diagnostics[:id0_max] = maximum(JuMP.value.(id0_est))
variable_diagnostics[:fs0_max] = maximum(JuMP.value.(fs0_est))
variable_diagnostics[:ms0_max] = maximum(JuMP.value.(ms0_est))
variable_diagnostics[:y0_max] = maximum(JuMP.value.(y0_est))
variable_diagnostics[:fd0_max] = maximum(JuMP.value.(fd0_est))
variable_diagnostics[:va0_max] = maximum(JuMP.value.(va0_est))
variable_diagnostics[:a0_max] = maximum(JuMP.value.(a0_est))
variable_diagnostics[:x0_max] = maximum(JuMP.value.(x0_est))
variable_diagnostics[:m0_max] = maximum(JuMP.value.(m0_est))
variable_diagnostics[:md0_max] = maximum(JuMP.value.(md0_est))

variable_diagnostics[:ys_rat] = Dict()
variable_diagnostics[:id_rat] = Dict()
variable_diagnostics[:fs_rat] = Dict()
variable_diagnostics[:ms_rat] = Dict()
variable_diagnostics[:y0_rat] = Dict()
variable_diagnostics[:fd_rat] = Dict()
variable_diagnostics[:va_rat] = Dict()
variable_diagnostics[:a0_rat] = Dict()
variable_diagnostics[:x0_rat] = Dict()
variable_diagnostics[:m0_rat] = Dict()
variable_diagnostics[:md_rat] = Dict()

for i in i_set
  for j in j_set
      if haskey(cal[:ys0],(j,i))
          val = (JuMP.value(ys0_est[j,i]) / cal[:ys0][j,i] )
          push!(variable_diagnostics[:ys_rat],[j,i]=>val)
      end
  end
end

for i in i_set
  for j in j_set
      if haskey(cal[:id0],(i,j))
          val = (JuMP.value(id0_est[i,j]) / cal[:id0][i,j] )
          push!(variable_diagnostics[:id_rat],[i,j]=>val)
      end
  end
end

for i in i_set
  if haskey(cal[:fs0],i)
      val = (JuMP.value(fs0_est[i]) / cal[:fs0][i])
      push!(variable_diagnostics[:fs_rat],[i]=>val)
  end
end

for i in i_set
  for m in m_set
      if haskey(cal[:ms0],(i,m))
      val = (JuMP.value(ms0_est[i,m]) / cal[:ms0][i,m])
      push!(variable_diagnostics[:ms_rat],[i,m]=>val)
      end
  end
end

for i in i_set
  if haskey(cal[:y0],i)
      val = (JuMP.value(y0_est[i]) / cal[:y0][i])
      push!(variable_diagnostics[:y0_rat],[i]=>val)
  end
end

for i in i_set
  for fd in fd_set
      if haskey(cal[:fd0],(i,fd))
          val = (JuMP.value(fd0_est[i,fd]) / cal[:fd0][i,fd])
          push!(variable_diagnostics[:fd_rat],[i,fd]=>val)
      end
  end
end

for va in va_set
  for j in j_set
      if haskey(cal[:va0],(va,j))
          val = (JuMP.value(va0_est[va,j]) / cal[:va0][va,j])
          push!(variable_diagnostics[:va_rat],[va,j]=>val)
      end
  end
end

for i in i_set
  if haskey(cal[:a0],(i))
      val = (JuMP.value(a0_est[i]) / cal[:a0][i])
      push!(variable_diagnostics[:a0_rat],[i]=>val)
  end
end

for i in i_set
  if haskey(cal[:x0],i)
      val = (JuMP.value(x0_est[i]) / cal[:x0][i])
      push!(variable_diagnostics[:x0_rat],[i]=>val)
  end
end

for i in i_set
  if haskey(cal[:m0],i)
      val = (JuMP.value(m0_est[i]) / cal[:m0][i])
      push!(variable_diagnostics[:m0_rat],[i]=>val)
  end
end

for m in m_set
  for i in i_set
      if haskey(cal[:md0],(m,i))
          val = (JuMP.value(md0_est[m,i]) / cal[:md0][m,i])
          push!(variable_diagnostics[:md_rat],[m,i]=>val)
      end
  end
end

variable_diagnostics[:ys_rat_max] = findmax(variable_diagnostics[:ys_rat])
variable_diagnostics[:id_rat_max] = findmax(variable_diagnostics[:id_rat])
variable_diagnostics[:fs_rat_max] = findmax(variable_diagnostics[:fs_rat])
variable_diagnostics[:ms_rat_max] = findmax(variable_diagnostics[:ms_rat])
variable_diagnostics[:y0_rat_max] = findmax(variable_diagnostics[:y0_rat])
variable_diagnostics[:fd_rat_max] = findmax(variable_diagnostics[:fd_rat])
variable_diagnostics[:va_rat_max] = findmax(variable_diagnostics[:va_rat])
variable_diagnostics[:a0_rat_max] = findmax(variable_diagnostics[:a0_rat])
variable_diagnostics[:x0_rat_max] = findmax(variable_diagnostics[:x0_rat])
variable_diagnostics[:m0_rat_max] = findmax(variable_diagnostics[:m0_rat])
variable_diagnostics[:md_rat_max] = findmax(variable_diagnostics[:md_rat])

variable_diagnostics[:ys_rat_min] = findmin(variable_diagnostics[:ys_rat])
variable_diagnostics[:id_rat_min] = findmin(variable_diagnostics[:id_rat])
variable_diagnostics[:fs_rat_min] = findmin(variable_diagnostics[:fs_rat])
variable_diagnostics[:ms_rat_min] = findmin(variable_diagnostics[:ms_rat])
variable_diagnostics[:y0_rat_min] = findmin(variable_diagnostics[:y0_rat])
variable_diagnostics[:fd_rat_min] = findmin(variable_diagnostics[:fd_rat])
variable_diagnostics[:va_rat_min] = findmin(variable_diagnostics[:va_rat])
variable_diagnostics[:a0_rat_min] = findmin(variable_diagnostics[:a0_rat])
variable_diagnostics[:x0_rat_min] = findmin(variable_diagnostics[:x0_rat])
variable_diagnostics[:m0_rat_min] = findmin(variable_diagnostics[:m0_rat])
variable_diagnostics[:md_rat_min] = findmin(variable_diagnostics[:md_rat])

################################
# -- Constraint diagnostics --
################################

constraint_diagnostic = Dict()

constraint_diagnostic[:mkt_py] = Dict()
constraint_diagnostic[:mkt_pa] = Dict()
constraint_diagnostic[:mkt_pm] = Dict()
constraint_diagnostic[:prf_y] = Dict()
constraint_diagnostic[:prf_a] = Dict()

for i in i_set
  push!(constraint_diagnostic[:mkt_py],i=>sum(JuMP.value(ys0_est[j,i]) for j in j_set) + JuMP.value(fs0_est[i])  - (sum(JuMP.value(ms0_est[i,m]) for m in m_set) + JuMP.value(y0_est[i])))
  push!(constraint_diagnostic[:mkt_pa],i=>JuMP.value(a0_est[i])  - (sum(JuMP.value(id0_est[i,j]) for j in j_set) + sum(JuMP.value(fd0_est[i,fd]) for fd in fd_set)))
  push!(constraint_diagnostic[:prf_a],i=>JuMP.value(a0_est[i]) * (1-cal[:ta0][i]) + JuMP.value(x0_est[i]) - ( JuMP.value(y0_est[i]) + JuMP.value(m0_est[i])*(1+cal[:tm0][i]) + sum(JuMP.value(md0_est[m,i]) for m in m_set) )  )
end

for m in m_set
  push!(constraint_diagnostic[:mkt_pm],m=>sum(JuMP.value(ms0_est[i,m]) for i in i_set)  - (sum(JuMP.value(md0_est[m,i]) for i in i_set)))
end

for j in j_set
  push!(constraint_diagnostic[:prf_y],j=>sum(JuMP.value(ys0_est[j,i]) for i in i_set)  - (sum(JuMP.value(id0_est[i,j]) for i in i_set) + sum(JuMP.value(va0_est[va,j]) for va in va_set)))
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
module TestInterface

using Phylo
using DataFrames
using Compat.Test

@testset "Build and tear down trees" begin
    @testset "For $TreeType" for TreeType in
        [NamedTree, NamedBinaryTree,
         BinaryTree{ManyRoots, DataFrame, Vector{Float64}}]

        species = ["Dog", "Cat", "Human", "Potato", "Apple"]
        tree = TreeType(species)
        othernodes = ["Some 1", "Some 2"]
        @test othernodes == createnodes!(tree, othernodes)
        extra = createnodes!(tree, 1)
        @test isa(extra[1], String)
        append!(othernodes, extra)
        also = createnode!(tree)
        @test also ∉ othernodes
        @test_throws ErrorException createnode!(tree, also)
        @test_throws ErrorException createnodes!(tree, [also])
        @test_throws ErrorException deletenode!(tree, "really not")
        @test also == deletenode!(tree, also)
        also = createnode!(tree)
        push!(othernodes, also)
        @test Set(othernodes) ⊆ Set(getnodenames(tree))
        innodes = append!(copy(species), othernodes)
        allnodes = reverse(innodes)
        pop!(innodes)
        branches =
            map(innodes) do node
                itr = Iterators.filter(name ->
                                       hasoutboundspace(tree, name) &&
                                       name != node, allnodes)
                createbranch!(tree, first(itr), node)
            end
        @test Set(branches) == Set(getbranchnames(tree))
        @test validate(tree)
        @test_throws ErrorException createbranch!(tree, allnodes[1], allnodes[2])
        who = getparent(tree, species[1])
        b = getinbound(tree, species[1])
        node1 = getnode(tree, species[1])
        @test !isunattached(tree, node1) &&
            hasoutboundspace(tree, node1) && !hasinboundspace(tree, node1) &&
            outdegree(tree, node1) == 0 && indegree(tree, node1) == 1
        @test hasbranch(tree, b)
        @test dst(getbranch(tree, b)) == species[1]
        @test src(getbranch(tree, b)) == who
        @test b == deletebranch!(tree, b)
        @test !hasbranch(tree, b)
        @test isunattached(tree, species[1]) &&
            hasoutboundspace(tree, species[1]) &&
            hasinboundspace(tree, species[1]) &&
            outdegree(tree, species[1]) == 0 &&
            indegree(tree, species[1]) == 0
        @test_throws ErrorException getbranch(tree, b)
        branches = [name for name in branches if name != b]
        @test Set(branches) == Set(getbranches(tree))
        b3 = getinbound(tree, species[2])
        source = src(tree, b3)
        destination = dst(tree, b3)
        @test b3 == deletebranch!(tree, b3)
        b3 = createbranch!(tree, who, species[1])
        @test who == src(tree, b3)
        @test species[1] == dst(tree, b3)
        b2 = createbranch!(tree, source, destination)
        push!(branches, b2)
        @test Set(branches) == Set(getbranchnames(tree))
        @test species[1] == deletenode!(tree, species[1])
        @test !hasnode(tree, species[1])
        @test species[1] == createnode!(tree, species[1])
        @test createbranch!(tree, who, species[1]) isa Int
        @test hasnode(tree, species[1])
        @test validate(tree)
        @test all(isleaf(tree, node) for node in species)
        @test all((!isroot(tree, node) & !isunattached(tree, node) &
                   !isinternal(tree, node)) for node in species)
        @test Set(src(tree, branch) for branch in getbranches(tree)) ∪
            Set(dst(tree, branch) for branch in getbranches(tree)) ==
            Set(allnodes)
    end
end

end

ext: jl


text:
# ------------ includes ------------ #

function primitive end # Declare for documentation.

include("macros/blackbox.jl")
include("macros/jaynes.jl")

# ------------ Documentation -------------- #

@doc(
"""
```julia
@primitive function logpdf(fn::typeof(foo), args, foo_ret)
    ...
end
```

`@primitive` is a convenience metaprogramming construct which derives contextual dispatch definitions for functions which you'd like the tracer to "summarize" to a single site. This is one mechanism which gives the user more control over the structure of the choice map structure in their programs.

Example:

```julia
geo(p::Float64) = rand(:flip, Bernoulli(p)) ? 1 : 1 + rand(:geo, geo, p)

# Define as primitive.
@primitive function logpdf(fn::typeof(geo), p, count)
    return Distributions.logpdf(Geometric(p), count)
end

function foo()
    ret = rand(:geo, geo, 0.3)
    ret
end
```

The above program would summarize the trace into a single choice site for `:geo`, as if `geo` was a primitive distribution.

```julia
  __________________________________

               Addresses

 geo : 42
  __________________________________
```
""", primitive)


text:
export MeanFieldGame
"""
    description of mean field game

    min_θ ∫∫ α1⋅L(x,-∇Φ) dx + α2⋅F(ρ(⋅,t)) dt + α3⋅G(ρ(⋅,1))
             + α4⋅∫∫ |-∂t ϕ(z(x,t),t) - δF(z(x,t),t) + ∇ϕ(z(x,t),t)| dx dt
			 + α5⋅∫ |ϕ(z(x,1),1) -  G(z(x,t),ρ(z(x,t),1))| dx

    s.t. ∂_t u = odefun(u,θ), u(0) = [x;0;0;0;0]

"""
mutable struct MeanFieldGame{R}
    F # function or array for interaction term
    G # function or array for terminal cost
    X0::AbstractArray{R} # training points
    rho0 # function to compute initial density rho0
    w::AbstractVector{R} # quadrature weights for X0
    rho0x::AbstractVector{R} # rho0(X0), stored for efficiency
    Φ # machine learning model for the potential
    α::Vector{R}  # vector containing penalties for objective functions [L, HJ]
    stepper # time integrator
    tspan::Vector{R} # time interval
    nt # number of time steps in ODE solve
    UN  # [X0,ldet,costL,costF,costHJ]
    cs  # [costL, costF, costG, costHJ,costHJfinal]
end

MeanFieldGame(F,G,X0::AbstractArray{R},rho0,w;
                Φ=PotentialNN(),α=ones(R,5),
                rho0x=rho0(X0),stepper=RK1Step(),tspan=R.([0.; 1.]),nt=2) where R <: Real =
            MeanFieldGame(F,G,X0,rho0,w,rho0x,Φ,α,stepper,tspan,nt,[],zeros(R,5))

"""
Pad X0 with zeros to initialize running costs
"""
function initUN(J,X0::AbstractArray{R}) where R <: Real
    J.UN = [X0; typeof(X0)(zeros(4,size(X0,2)))]
    return J.UN
end

function (J::MeanFieldGame{R})(Θ) where R <: Real
    (d,nex) = size(J.X0)

    # costF       = zero(R)
    h     = (J.tspan[2]-J.tspan[1])/J.nt

    UN = initUN(J,J.X0)

    tk    = J.tspan[1]
    for k=1:J.nt
        UN     = step(J.stepper,odefun,J,UN,Θ,tk,tk+h)
        tk    +=h
    end
    J.UN = UN

    # compute running costs
    costL = dot(vec(UN[end-2,:]),J.w)
    costF = dot(vec(UN[end-1,:]),J.w)
    # compute final costs
    costG = dot(J.G(UN),J.w)

	# compute HJB penalty
    costHJ = dot(vec(UN[end,:]),J.w)
    phi1 = vec(J.Φ([UN[1:d,:]; fill(R(1.0),1,size(J.X0,2))],Θ))
    costHJf = dot(abs.(phi1 - J.α[3].*vec(getDeltaG(J.G,UN))),J.w)

    cs = [costL, costF, costG, costHJ, costHJf]
    Jc = dot(J.α,cs)

    # store intermediate results for plotting and printing
    J.cs = cs .* J.α
    return Jc
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
using Bonsai, Test
using Bonsai: ws_upgrade
using HTTP

# @testset "start and stop" begin
# 	router = Router()
# 	start(router, port=9999, verbose=true)
# 	stop(router)
# 	@test true
# end

# @testset "web socket" begin
# 	function f(stream)
# 		ws = ws_upgrade(stream)
# 		try
# 			for i in ["hello", "mate"]
# 				write(ws, i)
# 			end
# 		finally
# 			close(ws)
# 		end
# 	end

# 	router = Router()
# 	get!(router, "/", f)
# 	start(router, port=9999, verbose=true)
# 	l = []

# 	HTTP.WebSockets.open("ws://127.0.0.1:9999") do ws
# 		while !eof(ws) && isopen(ws)
# 			x = readavailable(ws)
# 			if isempty(x)
# 				continue
# 			end
# 			push!(l, String(x))
# 		end
# 	end

# 	stop(router)
# 	@test l == ["hello", "mate"]
# end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
# All these functions were borrowed from the
# unregistered package LocalScopeBenchmarks.jl

function collect_symbols(expr)
    assignments = OrderedDict{Symbol, Expr}()
    prewalk(expr) do x
        if x isa Symbol
            assignments[x] = Expr(:$, x)
            return nothing
        elseif x isa Expr && x.head == :$
            # Don't recurse inside $() interpolations,
            # since those will already be interpolated
            return nothing
        else
            return x
        end
    end
    assignments
end

function parse_setup(setup::Expr)
    assignments = OrderedDict()
    postwalk(setup) do x
        if @capture(x, a_ = b_)
            assignments[a] = b
        end
        x
    end
    assignments
end

function lower_setup(assignments::AbstractDict)
    Expr(:block, [Expr(:(=), k, v) for (k, v) in assignments]...)
end

function parse_params(kwargs)
    params_dict = OrderedDict((@assert x.head == :kw; x.args[1] => x.args[2]) for x in kwargs)
end

function lower_params(params::AbstractDict)
    [Expr(:kw, k, v) for (k, v) in params]
end

function interpolate_locals_into_setup(args...)
    core, kwargs = BenchmarkTools.prunekwargs(args...)
    params = parse_params(kwargs)
    setup_assignments = parse_setup(get(() -> Expr(:block), params, :setup))
    local_assignments = collect_symbols(core)
    setup = merge(local_assignments, setup_assignments)
    params[:setup] = lower_setup(setup)
    core, lower_params(params)
end

macro localbenchmark(args...)
    core, params = interpolate_locals_into_setup(args...)
    quote
        BenchmarkTools.@benchmark(Ref($(core))[], $(params...))
    end
end


text:
using MyMCMCSampler
using SafeTestsets, Test


text:
function corner(r, side1, side2)
    path = []
    while (isborder(r,side1)&&isborder(r,side2))==false
        if isborder(r,side1)==false
            move!(r,side1)
            push!(path,rev(side1))
        end
        if isborder(r,side2)==false
            move!(r,side2)
            push!(path,rev(side2))
        end
    end
    return path
end

function rev(side)
    if side == Nord
        return Sud
    end
    if side == Sud
        return Nord
    end
    if side == West
        return Ost
    end
    if side == Ost
        return West
    end
end

function move_to_wall(r, side)
    while isborder(r,side)==false
        move!(r,side)
    end
end

function way_back(r,path)
    for i in 0:(length(path)-1)
        move!(r,path[length(path)-i])
    end
end

function marks_in_corners(r)
    path = corner(r,Sud,West)
    for i in [Nord,Ost,Sud,West]
        move_to_wall(r,i)
        putmarker!(r)
    end
    way_back(r,path)
end


text:
include("go.jl")
using Test

@testset "Illigal Moves" begin
    set_game(["b9","a8"], [], white)
    @test process_input("a9") == false

    state[27] = state[35] = state[45] = 2
    @test process_input("i6") == false

    state[57] = state[59] = state[49] = state[67] = 2
    @test process_input("d4") == false

    set_game(["c7", "d8", "e7", "e6", "d5", "c6"], ["d7"], white)
    @test process_input("d6") == false

    set_game(["c7", "d8", "e7", "e6", "d5", "c6"], ["d6"], white)
    @test process_input("d7") == false

    set_game(["c7", "d8", "e7", "f6", "e5", "d5", "c6"], ["d7", "e6"], white)
    @test process_input("d6") == false
end

@testset "Captures" begin
    set_game(["b9"], ["a9"], black)
    @test process_input("a8") == true
    @test stone_on("a9") == blank

    set_game(["d1", "e2"], ["e1"], black)
    @test process_input("f1") == true
    @test stone_on("e1") == blank

    set_game(["c7", "d8", "e7"], ["d7"], black)
    @test process_input("d6") == true
    @test stone_on("d7") == blank

    set_game(["c7", "d8", "e7", "f6", "e5", "d5"], ["d7", "e6", "d6"], black)
    @test process_input("c6") == true
    @test stone_on("d7") == blank
    @test stone_on("e6") == blank
    @test stone_on("d6") == blank
end

@testset "KO" begin
    set_game(["c7", "d8", "e7", "d6"], ["c6", "e6", "d5"], white)
    @test process_input("d7") == true
    @test stone_on("d6") == blank

    @test process_input("d6") == false
    @test stone_on("d6") == blank

    process_input("a9")
    process_input("b9")
    @test process_input("d6") == true
    @test stone_on("d7") == blank
end

max_stars_repo_name: aole/Games-in-Julia

max_issues_repo_name: aole/Games-in-Julia

max_forks_repo_name: aole/Games-in-Julia


text:
SLS2 = RheoModelClass(
        # Model name
        name="SLS2",
        # Model parameters,
        p = [:G₀, :G₁, :η₁, :G₂, :η₂],
        # Relaxation modulus
        G = quote
              G₀ + G₁*exp(-t*G₁/η₁) + G₂*exp(-t*G₂/η₂)
            end,
        # Creep modulus
        J = quote
              tau1 = G₁/η₁
              tau2 = G₂/η₂
              Jbar(s) = (1/s^2)*(G₀/s+G₁*1/(s+tau1)+G₂*1/(s+tau2))
              invLaplace(s -> Jbar(s), t)
            end,
        # Network
        info= "
                                       ___
                 ___╱╲  ╱╲  ╱╲  ________| |_____
                 |    ╲╱  ╲╱  ╲╱       _|_|     |
                 |      G₂              η₂      |
                 |                     ___      |
              ___|__╱╲  ╱╲  ╱╲  ________| |_____|___
                 |    ╲╱  ╲╱  ╲╱       _|_|     |
                 |      G₁              η₁      |
                 |                              |
                 |__________╱╲  ╱╲  ╱╲  ________|
                              ╲╱  ╲╱  ╲╱
                                G₀
               "
        )

#
#  Generate array of parameter symbols and expression of the relaxation modulus
#  for the generalised maxwell model
#

function G_GSLS(n)

    params=[:k0]
    G_terms=[:(k0)]

    for i=1:n
        k=Symbol("k_"*string(i))
        η=Symbol("η_"*string(i))
        e=:($k * exp.( (- $k/$η) * t))

        params=vcat(params, [k,η])
        G_terms=vcat(G_terms,e)
    end
    G_Expr=Expr(:call,:+,G_terms...)

    return (params, G_Expr)
end

# (p,e)=G_GSLS(3)

# cd("/home/alexandre/.julia/dev/RHEOS/src/"); include("RHEOS.jl"); using Main.RHEOS; m=RheoModel(SLS2,(G₀=1., G₁=2., η₁=3., G₂=4., η₂=5.))


text:
fit = perc_pos_to_num_tests_scale(model; lookback = 150)
θs = model |> flatten_chain

## Project PCR test results from the infection projections

sampled_PCRs = map((inf, θ) -> generate_PCR_predictions(inf, θ, model, fit; M_PCR = 150), infections, θs)
sampled_pos_tests = [sample[:, 1] for sample in sampled_PCRs]
sampled_prop_pos_tests = [[sample[i, 2] == 0 ? 0.0 : sample[i, 1] / sample[i, 2] for i = 1:size(sample, 1)] for sample in sampled_PCRs]

post_pred_pos_tests = get_credible_intervals(sampled_pos_tests)
post_pred_prop_pos_tests = get_credible_intervals(sampled_prop_pos_tests)

scatter(model.PCR_cases[model.model_start_time:end, 1])
plot!(post_pred_pos_tests.pred, ribbon = (post_pred_pos_tests.lb, post_pred_pos_tests.ub))

scatter(model.PCR_cases[model.model_start_time:end, 1] ./ model.PCR_cases[model.model_start_time:end, 2])
plot!(post_pred_prop_pos_tests.pred, ribbon = (post_pred_prop_pos_tests.lb, post_pred_prop_pos_tests.ub))

## Project serological test results from the infection projections

n = length(sero_prevalence_by_age_preds[1].pred)
uncertaintyfirstinfections_by_age = [[infs[:, a] for infs in firstinfections_by_age] for a = 1:size(firstinfections_by_age[1], 2)]
sero_prevalence_by_age = [map((ι_sero, θ) -> generate_sero_predictions(ι_sero, θ, age_grp, model), uncertaintyfirstinfections_by_age[age_grp], θs)
                          for age_grp = 1:6]

sero_data_by_age = [model.sero_cases[:, age_grp, :] for age_grp = 1:6]
dates = [Date(2020, 2, 19) + Day(k) for k = 1:size(model.sero_cases, 1)]
sero_prop_by_age = [monthly_proportions_with_errors(sero, dates, Date(2020, 2, 19)) for sero in sero_data_by_age]

##
k = 1
plot(model.model_start_time:(model.model_start_time-1+n), sero_prevalence_by_age_preds[k].pred,
    lab = "",
    ribbon = (sero_prevalence_by_age_preds[k].lb, sero_prevalence_by_age_preds[k].ub))

scatter!(sero_prop_by_age[k].ts, sero_prop_by_age[k].prop,
    yerrors = (sero_prop_by_age[k].lb, sero_prop_by_age[k].ub),
    lab = "")

##Project variant frequency

n = size(post_pred_pos_tests.pred, 1)
var_frqs = θs .|> θ -> get_variant_prediction(θ, model, n)

wt_frqs = [data.prop_inc_wt for data in var_frqs]
alpha_beta_frqs = [data.prop_inc_alpha_beta for data in var_frqs]
delta_frqs = [data.prop_inc_delta for data in var_frqs]

wt_frq_preds = get_credible_intervals(wt_frqs)
alpha_beta_frq_preds = get_credible_intervals(alpha_beta_frqs)
delta_frq_preds = get_credible_intervals(delta_frqs)

plot(wt_frq_preds.pred, ribbon = (wt_frq_preds.lb, wt_frq_preds.ub), lab = "WT")
plot!(alpha_beta_frq_preds.pred, ribbon = (alpha_beta_frq_preds.lb, alpha_beta_frq_preds.ub), lab = "Alpha/Beta")
plot!(delta_frq_preds.pred, ribbon = (delta_frq_preds.lb, delta_frq_preds.ub), lab = "Delta")

##

X = Matrix(RecursiveArrayTools.vecvec_to_mat(infections)')
pred = get_credible_intervals(X)

plot(pred.pred, ribbon = (pred.lb, pred.ub))

param_names = names(model.MCMC_results.chain)

get(model.MCMC_results.chain[1, :, :], param_names; flatten = true)

##
include("generation_methods.jl");

## Combine incidence data
fitfiles = readdir("modelfits", join = true)
filename = fitfiles[1]
fit_dict = load(filename)
model = fit_dict[first(keys(fit_dict))]
name = model.areaname
inc_by_type = get_daily_incidence_by_type(model, Date(2022, 6, 1),
    1.2,
    1.2,
    0.1,
    0.1,
    0.1)

kenya_inc_asymp = inc_by_type.asymp_infections
kenya_inc_M = inc_by_type.mild_infections
kenya_inc_severe = inc_by_type.severe_infections
kenya_inc_critical = inc_by_type.critical_infections

for filename in fitfiles[2:end]
    fit_dict = load(filename)
    model = fit_dict[first(keys(fit_dict))]
    name = model.areaname
    inc_by_type = get_daily_incidence_by_type(model, Date(2022, 6, 1), 1.2, 1.2, 0.1, 0.1, 0.1)
    kenya_inc_asymp .+= inc_by_type.asymp_infections
    kenya_inc_M .+= inc_by_type.mild_infections
    kenya_inc_severe .+= inc_by_type.severe_infections
    kenya_inc_critical .+= inc_by_type.critical_infections
end


text:
"""
Helps to maintain the pool of reusable arrays of different sizes
and reduce the burden on garbage collection.
"""
struct ArrayPool{T<:Number}
    length_pools::Dict{Int, Vector{Vector{T}}}

    ArrayPool{T}() where {T<:Number} = new{T}(Dict{Int, Vector{Vector{T}}}())
end

"""
Gets an array of specific size from the pool.
The returned array should be returned back to the pool using `release!()`.
"""
function acquire!(pool::ArrayPool{T}, size) where T
    len = prod(size)
    len_pool = haskey(pool.length_pools, len) ?
               pool.length_pools[len] :
               get!(() -> Vector{Vector{T}}(), pool.length_pools, len)
    return isempty(len_pool) ?
            Array{T}(undef, size) :
            reshape(pop!(len_pool), size)
end

"""
Releases an array returned by `acquire!()` back into the pool.
"""
function release!(pool::ArrayPool{T}, arr::Array{T}) where T
    len = length(arr)
    len_pool = get(pool.length_pools, len, nothing)
    if len_pool !== nothing
        #@info "release($(size(arr))) ($(length(len_pool)))"
        (length(len_pool) <= 100) || error("Overflow of $len-sized vectors pool")
        push!(len_pool, vec(arr))
    else
        throw(DimensionMismatch("No $len-element arrays were acquired before"))
    end
    return pool
end


text:
"""
A basis of the classical Jacobi polynomials on the interval `[-1,1]`.
These polynomials are orthogonal with respect to the weight function
```
w(x) = (1-x)^α (1+x)^β.
```
"""
struct Jacobi{T} <: OPS{T}
    n       ::  Int
    α       ::  T
    β       ::  T

    Jacobi{T}(n, α = 0, β = 0) where {T} = new{T}(n, α, β)
end

Jacobi(n::Int; α = 0, β = 0) = Jacobi(n, α, β)
Jacobi(n::Int, α, β) = Jacobi(n, promote(α, β)...)
Jacobi(n::Int, α::T, β::T) where {T <: AbstractFloat} = Jacobi{T}(n, α, β)
Jacobi(n::Int, α::S, β::S) where {S} = Jacobi(n, float(α), float(β))

similar(b::Jacobi, ::Type{T}, n::Int) where {T} = Jacobi{T}(n, b.α, b.β)

support(b::Jacobi{T}) where {T} = ChebyshevInterval{T}()

first_moment(b::Jacobi{T}) where {T} = (b.α+b.β+1≈0) ?
    T(2).^(b.α+b.β+1)*gamma(b.α+1)*gamma(b.β+1) :
    T(2).^(b.α+b.β+1)*gamma(b.α+1)*gamma(b.β+1)/(b.α+b.β+1)/gamma(b.α+b.β+1)
    # 2^(b.α+b.β) / (b.α+b.β+1) * gamma(b.α+1) * gamma(b.β+1) / gamma(b.α+b.β+1)

jacobi_α(b::Jacobi) = b.α
jacobi_β(b::Jacobi) = b.β

measure(b::Jacobi) = JacobiWeight(b.α, b.β)

iscompatible(d1::Jacobi, d2::Jacobi) = d1.α == d2.α && d1.β == d2.β

isorthogonal(dict::Jacobi, measure::JacobiWeight) =
	dict.α == measure.α && dict.β == measure.β

gauss_rule(dict::Jacobi) = GaussJacobi(length(dict), dict.α, dict.β)

interpolation_grid(dict::Jacobi) = JacobiNodes(length(dict), dict.α, dict.β)
iscompatible(dict::Jacobi, grid::JacobiNodes) =  length(dict) == length(grid) && dict.α ≈ grid.α && dict.β ≈ grid.β
isorthogonal(dict::Jacobi, measure::GaussJacobi) = jacobi_α(dict) ≈ jacobi_α(measure) && jacobi_β(dict) ≈ jacobi_β(measure) && opsorthogonal(dict, measure)

issymmetric(dict::Jacobi) = jacobi_α(dict)≈jacobi_β(dict)

# See DLMF (18.9.2)
# http://dlmf.nist.gov/18.9#i
function rec_An(b::Jacobi{T}, n::Int) where {T}
    if (n == 0) && (b.α + b.β+1 == 0)
        one(T)/2*(b.α+b.β)+1
    else
        T(2*n + b.α + b.β + 1) * (2*n + b.α + b.β + 2) / T(2 * (n+1) * (n + b.α + b.β + 1))
    end
end

function rec_Bn(b::Jacobi{T}, n::Int) where {T}
    if (n == 0) && ((b.α + b.β + 1 == 0) || (b.α+b.β == 0))
        one(T)/2*(b.α-b.β)
    else
        T(b.α^2 - b.β^2) * (2*n + b.α + b.β + 1) / T(2 * (n+1) * (n + b.α + b.β + 1) * (2*n + b.α + b.β))
    end
end

rec_Cn(b::Jacobi{T}, n::Int) where {T} =
    T(n + b.α) * (n + b.β) * (2*n + b.α + b.β + 2) / T((n+1) * (n + b.α + b.β + 1) * (2*n + b.α + b.β))

function innerproduct_native(d1::Jacobi, i::PolynomialDegree, d2::Jacobi, j::PolynomialDegree, measure::JacobiWeight; options...)
	T = coefficienttype(d1)
	if iscompatible(d1, d2) && isorthogonal(d1, measure)
		if i == j
			a = d1.α
			b = d1.β
			n = convert(T, value(i))
			2^(a+b+1)/(2n+a+b+1) * gamma(n+a+1)*gamma(n+b+1)/factorial(n)/gamma(n+a+b+1)
		else
			zero(T)
		end
	else
		innerproduct1(d1, i, d2, j, measure; options...)
	end
end

## Printing
function show(io::IO, b::Jacobi{Float64})
	if jacobi_α(b) == 0 && jacobi_β(b) == 0
		print(io, "Jacobi($(length(b)))")
	else
		print(io, "Jacobi($(length(b)); α = $(jacobi_α(b)), β = $(jacobi_β(b)))")
	end
end

function show(io::IO, b::Jacobi{T}) where T
	if jacobi_α(b) == 0 && jacobi_β(b) == 0
		print(io, "Jacobi{$(T)}($(length(b)))")
	else
		print(io, "Jacobi{$(T)}($(length(b)); α = $(jacobi_α(b)), β = $(jacobi_β(b)))")
	end
end

# TODO: move to its own file and make more complete
# Or better yet: implement in terms of Jacobi polynomials
struct UltrasphericalBasis{T} <: OPS{T}
	n		::	Int
	alpha	::	T
end

jacobi_α(b::UltrasphericalBasis) = b.α
jacobi_β(b::UltrasphericalBasis) = b.α

weightfun(b::UltrasphericalBasis, x) = (1-x)^(b.α) * (1+x)^(b.α)


text:
function Well44497a()
    return Well44497a(())
end

function Well44497a(arg0::Vector{jint})
    return Well44497a((Vector{jint},), arg0)
end

function Well44497a(arg0::jint)
    return Well44497a((jint,), arg0)
end

function Well44497a(arg0::jlong)
    return Well44497a((jlong,), arg0)
end

function next_int(obj::Well44497a)
    return jcall(obj, "nextInt", jint, ())
end


text:
#
# Extended Kalman Filter Types and Methods
#

"""
    ExtendedKalmanFilter(sys::NonlinearSystem, obs::NonlinearObserver,
                         estimate::AbstractState[, consider_states::Vector{UInt}])

Extended Kalman Filter type. The internal filter model takes on a specific form
depending on the type of the initial estimate. For a `DiscreteState`,

\$x_k = F(k, x_{k-1}) + w_k\$
\$y_k = H(k, x_k) + v_k\$

where \$w_k \\sim N(0, Q)\$ and \$v_k \\sim N(0, R)\$. For a `ContinuousState`,

\$\\dot{x}(t_k) = F(t_k, x(t_k)) + w(t_k)\$
\$y(t_k) = H(t_k, x(t_k)) + v(t_k)\$

Construction of an ExtendedKalmanFilter requires an initial estimate.
Internally, the initial estimate is an uncertain state type. The constructor
automatically converts absolute initial estimates to uncertain states with zero
covariance.

The `consider_states` input contains a list of indices of state elements to be
considered in the filtering process and not updated.
"""
immutable ExtendedKalmanFilter{T,S<:AbstractUncertainState{T}} <:
                                                       AbstractKalmanFilter{T,S}
    sys::NonlinearSystem{T}
    obs::NonlinearObserver{T}
    estimate::S
    consider_states::Vector{UInt16}

    function ExtendedKalmanFilter(sys, obs, estimate::S, consider_states::Vector
                                  ) where {T, S<:AbstractUncertainState{T}}
        if !allunique(consider_states)
            throw(ArgumentError("Consider state indices must be unique"))
        end
        for idx = 1:length(consider_states)
            if consider_states[idx] > length(estimate.x)
                throw(DimensionMismatch(
                    "Consider indices extend beyond length of initial state."))
            end
        end
        new{T,S}(NonlinearSystem(sys), NonlinearObserver(obs), estimate,
                 consider_states)
    end
end
ExtendedKalmanFilter(sys, obs, estimate) =
    ExtendedKalmanFilter(sys, obs, estimate, [])
ExtendedKalmanFilter(sys, obs, estimate::AbstractAbsoluteState, consider_states) =
    ExtendedKalmanFilter(sys, obs, make_uncertain(estimate), consider_states)

const DiscreteEKF{T} =
    ExtendedKalmanFilter{T,UncertainDiscreteState{T}} where T
const ContinuousEKF{T} =
    ExtendedKalmanFilter{T,UncertainContinuousState{T}} where T

"""
    kalman_predict(kf::KalmanFilter, t)
"""
function kalman_predict(ekf::ExtendedKalmanFilter, t)
    xk = predict(ekf.estimate, ekf.sys, t)
    yk = predict(xk, ekf.obs)
    Pxy = xk.P * ekf.obs.dH_dx(t, ekf.estimate.x)'
    return xk, yk, Pxy
end

"""
    extended_kalman_filter!(state, system, observer, measurement[, archive])

Extended Kalman Filter convenience function. Performs a single Kalman Filter
iteration, updating the provided state in-place. This function is not as
efficient as constructing creating a static `ExtendedKalmanFilter` and using the
`process!` method, but allows for increased flexibility.

Use `extended_kalman_filter` for not-in-place state updates.
"""
function extended_kalman_filter!(state, system, observer, measurement)
    ekf = ExtendedKalmanFilter(system, observer, state)
    process!(ekf, measurement)
end
function extended_kalman_filter!(state, system, observer, measurement, archive)
    ekf = ExtendedKalmanFilter(system, observer, state)
    process!(ekf, measurement, archive)
end

"""
extended_kalman_filter(state, system, observer, measurement[, archive])

Extended Kalman Filter convenience function. Performs a single Kalman Filter
iteration. This function is not as efficient as constructing creating a static
`ExtendedKalmanFilter` and using the `process!` method, but allows for increased
flexibility.

Use `extended_kalman_filter!` for in-place state updates.
"""
function extended_kalman_filter(state, system, observer, measurement)
    out_state = deepcopy(state)
    extended_kalman_filter!(out_state, system, observer, measurement)
    return out_state
end
function extended_kalman_filter(state, system, observer, measurement, archive)
    out_state = deepcopy(state)
    extended_kalman_filter!(out_state, system, observer, measurement, archive)
    return out_state
end


text:
using BinDeps, Compat
using BinDeps: libdir

modified_defaults = false
if !in(BinDeps.Binaries, BinDeps.defaults)
    pushfirst!(BinDeps.defaults, BinDeps.Binaries)
    modified_defaults = true
end

BinDeps.@setup
did_setup = true

const OSF_VERS = v"0.5.3"

openspecfun = library_dependency("libopenspecfun")

const URL = "https://github.com/ararslan/openspecfun-builder/releases/download/v$OSF_VERS" *
            "/libopenspecfun-$OSF_VERS"

const DOWNLOADS = Dict(
    "Linux-x86_64"   => ("$URL-linux-x86_64.tar.gz",
                         "d70a2a391915f64f44da21915bf93ce08d054127028088addca36e16ac53bcb1"),
    "Linux-i686"     => ("$URL-linux-i686.tar.gz",
                         "e5418b170b537af2f7f1f1d06eee9be01555404f5d22a47e18bc06a540321478"),
    "Darwin-x86_64"  => ("$URL-osx-x86_64.tar.gz",
                         "e57f5f84439757a2fd1d3821a6e19a3fa69b5b1e181cc40fec0d1652fbb9efdc"),
    "Windows-x86_64" => ("$URL-win-x86_64.zip",
                         "7a5f7be4ed46d7f9d6d18a599157075512c50a372da2b2908079a3dcab9a0f25"),
    "Windows-i686"   => ("$URL-win-i686.zip",
                         "2f63a08d80e67964e2c368367f4caef7039080828e217d288669416cd46f4584"),
)

const SYSTEM = string(BinDeps.OSNAME, '-', Sys.ARCH)

if haskey(DOWNLOADS, SYSTEM)
    url, sha = DOWNLOADS[SYSTEM]
    provides(Binaries, URI(url), openspecfun, SHA=sha, os=BinDeps.OSNAME,
             unpacked_dir=joinpath("usr", "lib"), installed_libpath=libdir(openspecfun))
else
    info("No precompiled binaries found for your system. Building from scratch...")
    include("scratch.jl")
end

BinDeps.@install Dict(:libopenspecfun => :openspecfun)

if modified_defaults
    popfirst!(BinDeps.defaults)
end


text:
using BlockArrays
using LinearAlgebra
using MixedModels
using PooledArrays
using Random
using SparseArrays
using Statistics
using StatsModels
using Tables
using Test

using MixedModels: dataset, likelihoodratiotest

const io = IOBuffer()

include("modelcache.jl")

@testset "Dyestuff" begin
    fm1 = only(models(:dyestuff))

    @test length(fm1.allterms) == 3
    @test size(fm1.reterms) == (1, )
    @test lowerbd(fm1) == zeros(1)
    @test fm1.lowerbd == zeros(1)
    @test fm1.optsum.initial == ones(1)
    fm1.θ = ones(1)
    @test fm1.θ == ones(1)

    @test_throws ArgumentError fit!(fm1)

    fm1.optsum.feval = -1
    @test_logs (:warn, "Model has not been fit") show(fm1)

    @test objective(updateL!(setθ!(fm1, [0.713]))) ≈ 327.34216280955366
    @test_deprecated MixedModels.describeblocks(IOBuffer(), fm1)

    show(io, BlockDescription(fm1))
    @test countlines(seekstart(io)) == 3
    output = String(take!(io))
    @test startswith(output, "rows:")

    refit!(fm1)

    @test :θ in propertynames(fm1)
    @test objective(fm1) ≈ 327.3270598811428 atol=0.001
    @test fm1.θ ≈ [0.752580] atol=1.e-5
    @test fm1.λ ≈ [LowerTriangular(reshape(fm1.θ, (1,1)))] atol=1.e-5
    @test deviance(fm1) ≈ 327.32705988 atol=0.001
    @test aic(fm1) ≈ 333.3270598811394 atol=0.001
    @test bic(fm1) ≈ 337.5306520261259 atol=0.001
    @test fixef(fm1) ≈ [1527.5]
    @test dispersion_parameter(fm1)
    @test first(first(fm1.σs)) ≈ 37.26034462135931 atol=0.0001
    @test fm1.β ≈ [1527.5]
    @test dof(fm1) == 3
    @test nobs(fm1) == 30
    @test MixedModels.fixef!(zeros(1),fm1) ≈ [1527.5]
    @test coef(fm1) ≈ [1527.5]
    fm1β = fm1.βs
    @test fm1β isa NamedTuple
    @test isone(length(fm1β))
    @test first(values(fm1β)) ≈ 1527.5
    fm1σρ = fm1.σρs
    @test fm1σρ isa NamedTuple
    @test isone(length(fm1σρ))
    @test isone(length(getproperty(first(fm1σρ), :σ)))
    @test isempty(getproperty(first(fm1σρ), :ρ))
    @test fm1.σ == sdest(fm1)
    @test fm1.b == ranef(fm1)
    @test fm1.u == ranef(fm1, uscale=true)
    @test fm1.stderror == stderror(fm1)
    @test isone(length(fm1.pvalues))
    @test fm1.objective == objective(fm1)
    @test fm1.σ ≈ 49.510099986291145 atol=1.e-5
    @test fm1.X == ones(30,1)
    ds = MixedModels.dataset(:dyestuff)
    @test fm1.y == ds[:yield]
    @test cond(fm1) == ones(1)
    @test first(leverage(fm1)) ≈ 0.15650534392640486 rtol=1.e-5
    @test sum(leverage(fm1)) ≈ 4.695160317792145 rtol=1.e-5
    cm = coeftable(fm1)
    @test length(cm.rownms) == 1
    @test length(cm.colnms) == 4
    @test fnames(fm1) == (:batch,)
    @test response(fm1) == ds[:yield]
    rfu = ranef(fm1, uscale = true)
    rfb = ranef(fm1)
    @test abs(sum(rfu[1])) < 1.e-5
    cv = condVar(fm1)
    @test length(cv) == 1
    @test size(first(cv)) == (1, 1, 6)
    show(IOBuffer(), fm1.optsum)

    @test logdet(fm1) ≈ 8.06014522999825 atol=0.001
    @test varest(fm1) ≈ 2451.2501089607676 atol=0.001
    @test pwrss(fm1) ≈ 73537.49947885796 atol=0.001
    @test stderror(fm1) ≈ [17.69455188898009] atol=0.0001

    vc = VarCorr(fm1)
    show(io, vc)
    str = String(take!(io))
    @test startswith(str, "Variance components:")
    @test vc.s == sdest(fm1)

    refit!(fm1, REML=true)
    @test objective(fm1) ≈ 319.65427684225216 atol=0.0001
    @test_throws ArgumentError loglikelihood(fm1)
    @test dof_residual(fm1) ≥ 0

    print(io, fm1)
    @test startswith(String(take!(io)), "Linear mixed model fit by REML")

    fm1.optsum.maxfeval = 5
    fm1.optsum.feval = -1
    @test_logs (:warn, "NLopt optimization failure: MAXEVAL_REACHED") fit!(fm1)
    fm1.optsum.maxfeval = -1

    vc = fm1.vcov
    @test isa(vc, Matrix{Float64})
    @test only(vc) ≈ 409.79495436473167 rtol=1.e-6
    # since we're caching the fits, we should get it back to being correctly fitted
    refit!(fm1)
end

@testset "Dyestuff2" begin
    ds2 = MixedModels.dataset(:dyestuff2)
    fm = only(models(:dyestuff2))
    @test lowerbd(fm) == zeros(1)
    show(IOBuffer(), fm)
    @test fm.θ ≈ zeros(1)
    @test objective(fm) ≈ 162.87303665382575
    @test abs(std(fm)[1][1]) < 1.0e-9
    @test std(fm)[2] ≈ [3.653231351374652]
    @test stderror(fm) ≈ [0.6669857396443261]
    @test coef(fm) ≈ [5.6656]
    @test logdet(fm) ≈ 0.0
    @test issingular(fm)
    #### modifies the model
    refit!(fm, float(MixedModels.dataset(:dyestuff)[:yield]))
    @test objective(fm) ≈ 327.3270598811428 atol=0.001
    refit!(fm, float(MixedModels.dataset(:dyestuff2)[:yield])) # restore the model in the cache
end

@testset "penicillin" begin
    fm = only(models(:penicillin))
    @test size(fm) == (144, 1, 30, 2)
    @test fm.optsum.initial == ones(2)
    @test lowerbd(fm) == zeros(2)

    @test objective(fm) ≈ 332.18834867227616 atol=0.001
    @test coef(fm) ≈ [22.97222222222222] atol=0.001
    @test fixef(fm) ≈ [22.97222222222222] atol=0.001
    @test coef(fm)[1] ≈ mean(dataset(:penicillin).diameter)
    @test stderror(fm) ≈ [0.7445960346851368] atol=0.0001
    @test fm.θ ≈ [1.5375772376554968, 3.219751321180035] atol=0.001
    @test first(std(fm)) ≈ [0.8455645948223015] atol=0.0001
    @test std(fm)[2] ≈ [1.770647779277388] atol=0.0001
    @test varest(fm) ≈ 0.3024263987592062 atol=0.0001
    @test logdet(fm) ≈ 95.74614821367786 atol=0.001

    @test_throws ArgumentError condVar(fm)

    rfu = ranef(fm, uscale=true)
    @test length(rfu) == 2
    @test first(first(rfu)) ≈ 0.523162392717432 rtol=1.e-4

    rfb = ranef(fm)
    @test length(rfb) == 2
    @test last(last(rfb)) ≈ -3.001823834230942 rtol=1.e-4

    show(io, BlockDescription(fm))
    @test countlines(seekstart(io)) == 4
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "Diag/Dense" in tokens
    @test "Diagonal" in tokens
end

@testset "pastes" begin
    fm = last(models(:pastes))
    @test size(fm) == (60, 1, 40, 2)
    @test fm.optsum.initial == ones(2)
    @test lowerbd(fm) == zeros(2)

    @test objective(fm) ≈ 247.99446586289676 atol=0.001
    @test coef(fm) ≈ [60.05333333333329] atol=0.001
    @test fixef(fm) ≈ [60.05333333333329] atol=0.001
    @test stderror(fm) ≈ [0.6421359883527029] atol=0.0001
    @test fm.θ ≈ [3.5268858714382905, 1.3299230213750168] atol=0.001
    @test first(std(fm)) ≈ [2.904069002535747] atol=0.001
    @test std(fm)[2] ≈ [1.095070371687089] atol=0.0001
    @test std(fm)[3] ≈ [0.8234088395243269] atol=0.0001
    @test varest(fm) ≈ 0.6780020742644107 atol=0.0001
    @test logdet(fm) ≈ 101.0381339953986 atol=0.001

    show(io, BlockDescription(fm))
    @test countlines(seekstart(io)) == 4
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "Sparse" in tokens
    @test "Diagonal" in tokens

    lrt = likelihoodratiotest(models(:pastes)...)
    @test length(lrt.deviance) == length(lrt.formulas) == length(lrt.models )== 2
    @test first(lrt.tests.pvalues) ≈ 0.5233767966395597 atol=0.0001
end

@testset "InstEval" begin
    fm1 = first(models(:insteval))
    @test size(fm1) == (73421, 2, 4114, 3)
    @test fm1.optsum.initial == ones(3)
    @test lowerbd(fm1) == zeros(3)

    @test objective(fm1) ≈ 237721.7687745563 atol=0.001
    ftd1 = fitted(fm1);
    @test size(ftd1) == (73421, )
    @test ftd1 == predict(fm1)
    @test first(ftd1) ≈ 3.17876 atol=0.0001
    resid1 = residuals(fm1);
    @test size(resid1) == (73421, )
    @test first(resid1) ≈ 1.82124 atol=0.00001

    @testset "PCA" begin
        @test length(fm1.rePCA) == 3
        pca = MixedModels.PCA(fm1)
        @test length(pca) == 3
        @test :covcor in propertynames(first(pca))
        str = String(take!(io))
        show(io, first(pca), stddevs=true, variances=true)
        str = String(take!(io))
        @test !isempty(findall("Standard deviations:", str))
        @test !isempty(findall("Variances:", str))
    end

    show(io, BlockDescription(fm1))
    @test countlines(seekstart(io)) == 5
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "Sparse" in tokens
    @test "Sparse/Dense" in tokens
    @test "Diag/Dense" in tokens

    fm2 = last(models(:insteval))
    @test objective(fm2) ≈ 237585.5534151694 atol=0.001
    @test size(fm2) == (73421, 28, 4100, 2)
end

@testset "sleep" begin
    fm = last(models(:sleepstudy))
    @test lowerbd(fm) == [0.0, -Inf, 0.0]
    A11 = getblock(fm.A, 1,1)
    @test isa(A11, UniformBlockDiagonal{Float64})
    @test isa(getblock(fm.L, 1, 1), UniformBlockDiagonal{Float64})
    @test size(A11) == (36, 36)
    a11 = view(A11.data, :, :, 1)
    @test a11 == [10. 45.; 45. 285.]
    @test size(A11.data, 3) == 18
    λ = first(fm.λ)
    b11 = LowerTriangular(view(getblock(fm.L, 1, 1).data, :, :, 1))
    @test b11 * b11' ≈ λ'a11*λ + I rtol=1e-5
    @test count(!iszero, Matrix(getblock(fm.L, 1, 1))) == 18 * 4
    @test rank(fm) == 2

    @test objective(fm) ≈ 1751.9393444647046
    @test fm.θ ≈ [0.929221307, 0.01816838, 0.22264487096] atol=1.e-6
    @test pwrss(fm) ≈ 117889.46144025437
    @test logdet(fm) ≈ 73.90322021999222 atol=0.001
    @test stderror(fm) ≈ [6.632257721914501, 1.5022354739749826] atol=0.0001
    @test coef(fm) ≈ [251.40510484848477,10.4672859595959]
    @test fixef(fm) ≈  [251.40510484848477,10.4672859595959]
    @test std(fm)[1] ≈ [23.780468100188497, 5.716827903196682] atol=0.01
    @test logdet(fm) ≈ 73.90337187545992 atol=0.001
    @test cond(fm) ≈ [4.175251] atol=0.0001
    @test loglikelihood(fm) ≈ -875.9696722323523
    @test sum(leverage(fm)) ≈ 28.611525700136877 rtol=1.e-5
    σs = fm.σs
    @test length(σs) == 1
    @test keys(σs) == (:subj,)
    @test length(σs.subj) == 2
    @test first(values(σs.subj)) ≈ 23.7804686 atol=0.0001
    @test last(values(first(σs))) ≈ 5.7168278 atol=0.0001
    @test fm.corr ≈ [1.0 -0.1375451787621904; -0.1375451787621904 1.0] atol=0.0001

    u3 = ranef(fm, uscale=true)
    @test length(u3) == 1
    @test size(first(u3)) == (2, 18)
    @test first(u3)[1, 1] ≈ 3.030300122575336 atol=0.001

    b3 = ranef(fm)
    @test length(b3) == 1
    @test size(first(b3)) == (2, 18)
    @test first(first(b3)) ≈ 2.815819441982976 atol=0.001

    b3tbl = raneftables(fm)
    @test length(b3tbl) == 1
    @test keys(b3tbl) == (:subj,)
    @test isa(b3tbl, NamedTuple)
    @test Tables.istable(only(b3tbl))

    simulate!(fm)  # to test one of the unscaledre methods
    # must restore state of fm as it is cached in the global fittedmodels
    slp = MixedModels.dataset(:sleepstudy)
    copyto!(fm.y, slp.reaction)
    updateL!(MixedModels.reevaluateAend!(fm))
    @test objective(fm) ≈ 1751.9393444647046 # check the model is properly restored

    fmnc = models(:sleepstudy)[2]
    @test size(fmnc) == (180,2,36,1)
    @test fmnc.optsum.initial == ones(2)
    @test lowerbd(fmnc) == zeros(2)

    @testset "zerocorr PCA" begin
        @test length(fmnc.rePCA) == 1
        @test fmnc.rePCA.subj ≈ [0.5, 1.0]
        @test any(Ref(fmnc.PCA.subj.loadings) .≈ (I(2), I(2)[:, [2,1]]))
        @test show(IOBuffer(), MixedModels.PCA(fmnc)) === nothing
    end

    @test deviance(fmnc) ≈ 1752.0032551398835 atol=0.001
    @test objective(fmnc) ≈ 1752.0032551398835 atol=0.001
    @test coef(fmnc) ≈ [251.40510484848585, 10.467285959595715]
    @test fixef(fmnc) ≈ [251.40510484848477, 10.467285959595715]
    @test stderror(fmnc) ≈ [6.707710260366577, 1.5193083237479683] atol=0.001
    @test fmnc.θ ≈ [0.9458106880922268, 0.22692826607677266] atol=0.0001
    @test first(std(fmnc)) ≈ [24.171449463289047, 5.799379721123582]
    @test last(std(fmnc)) ≈ [25.556130034081047]
    @test logdet(fmnc) ≈ 74.46952585564611 atol=0.001
    ρ = first(fmnc.σρs.subj.ρ)
    @test ρ === -0.0   # test that systematic zero correlations are returned as -0.0

    MixedModels.likelihoodratiotest(fm, fmnc)
    fmrs = fit(MixedModel, @formula(reaction ~ 1+days + (0+days|subj)), slp);
    @test objective(fmrs) ≈ 1774.080315280528 rtol=0.00001
    @test fmrs.θ ≈ [0.24353985679033105] rtol=0.00001

    fm_ind = models(:sleepstudy)[3]
    @test objective(fm_ind) ≈ objective(fmnc)
    @test coef(fm_ind) ≈ coef(fmnc)
    @test fixef(fm_ind) ≈ fixef(fmnc)
    @test stderror(fm_ind) ≈ stderror(fmnc)
    @test fm_ind.θ ≈ fmnc.θ
    @test std(fm_ind) ≈ std(fmnc)
    @test logdet(fm_ind) ≈ logdet(fmnc)

    # combining [ReMat{T,S1}, ReMat{T,S2}] for S1 ≠ S2
    slpcat = (subj = slp.subj, days = PooledArray(string.(slp.days)), reaction = slp.reaction)
    fm_cat = fit(MixedModel, @formula(reaction ~ 1+days+(1|subj)+(0+days|subj)),slpcat)
    @test fm_cat isa LinearMixedModel
    σρ = fm_cat.σρs
    @test σρ isa NamedTuple
    @test isone(length(σρ))
    @test first(keys(σρ)) == :subj
    @test keys(σρ.subj) == (:σ, :ρ)
    @test length(σρ.subj) == 2
    @test length(first(σρ.subj)) == 10
    @test length(σρ.subj.ρ) == 45
    # test that there's no correlation between the intercept and days columns
    ρs_intercept = σρ.subj.ρ[1 .+ cumsum(0:8)]
    @test all(iszero.(ρs_intercept))
    # amalgamate should set these to -0.0 to indicate structural zeros
    @test all(ρs_intercept .=== -0.0)

    # also works without explicitly dropped intercept
    fm_cat2 = fit(MixedModel, @formula(reaction ~ 1+days+(1|subj)+(days|subj)),slpcat)
    @test fm_cat2 isa LinearMixedModel
    σρ = fm_cat2.σρs
    @test σρ isa NamedTuple
    @test isone(length(σρ))
    @test first(keys(σρ)) == :subj
    @test keys(σρ.subj) == (:σ, :ρ)
    @test length(σρ.subj) == 2
    @test length(first(σρ.subj)) == 10
    @test length(σρ.subj.ρ) == 45
    # test that there's no correlation between the intercept and days columns
    ρs_intercept = σρ.subj.ρ[1 .+ cumsum(0:8)]
    @test all(iszero.(ρs_intercept))
    # amalgamate should set these to -0.0 to indicate structural zeros
    @test all(ρs_intercept .=== -0.0)

    show(io, BlockDescription(first(models(:sleepstudy))))
    @test countlines(seekstart(io)) == 3
    @test "Diagonal" in Set(split(String(take!(io)), r"\s+"))

    show(io, BlockDescription(last(models(:sleepstudy))))
    @test countlines(seekstart(io)) == 3
    @test "BlkDiag" in Set(split(String(take!(io)), r"\s+"))

end

@testset "d3" begin
    fm = only(models(:d3))
    @test pwrss(fm) ≈ 5.30480294295329e6 rtol=1.e-4
    @test objective(fm) ≈ 884957.5540213 rtol = 1e-4
    @test coef(fm) ≈ [0.4991229873, 0.31130780953] atol = 1.e-4
    @test length(ranef(fm)) == 3

    show(io, BlockDescription(fm))
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "BlkDiag/Dense" in tokens
end

@testset "kb07" begin
    global io
    pca = last(models(:kb07)).PCA
    @test keys(pca) == (:subj, :item)
    show(io, models(:kb07)[2])
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "Corr." in tokens
    @test "-0.89" in tokens
end

@testset "oxide" begin
    # this model has an interesting structure with two diagonal blocks
    m = first(models(:oxide))
    @test all(isapprox.(m.θ, [1.689182746, 2.98504262]; atol=1e-4))
    m = last(models(:oxide))
    # NB: this is a poorly defined fit
    # lme4 gives all sorts of convergence warnings for the different
    # optimizers and even quite different values
    # the overall estimates of the standard deviations are similar-ish
    # but the correlation structure seems particular unstable
    θneldermead = [1.6454, 8.6373e-02, 8.2128e-05, 8.9552e-01, 1.2014, 2.9286]
    # two different BOBYQA implementations
    θnlopt = [1.645, -0.221, 0.986, 0.895, 2.511, 1.169]
    θminqa = [1.6455, -0.2430, 1.0160, 0.8955, 2.7054, 0.0898]
    # very loose tolerance for unstable fit
    # but this is a convenient test of rankUpdate!(::UniformBlockDiagonal)
    @test all(isapprox.(m.θ, θnlopt; atol=1e-2))
end

@testset "Rank deficient" begin
    rng = MersenneTwister(0);
    x = rand(rng, 100);
    data = (x = x, x2 = 1.5 .* x, y = rand(rng, 100), z = repeat('A':'T', 5))
    model = fit(MixedModel, @formula(y ~ x + x2 + (1|z)), data)
    @test length(fixef(model)) == 2
    @test rank(model) == 2
    @test length(coef(model)) == 3
    ct = coeftable(model)
    @test ct.rownms ==  ["(Intercept)", "x", "x2"]
    @test length(fixefnames(model)) == 2
    @test coefnames(model) == ["(Intercept)", "x", "x2"]
    piv = first(model.feterms).piv
    r = first(model.feterms).rank
    @test coefnames(model)[piv][1:r] == fixefnames(model)
end

@testset "coeftable" begin
    ct = coeftable(only(models(:dyestuff)));
    @test [3,4] == [ct.teststatcol, ct.pvalcol]
end

@testset "wts" begin
    # example from https://github.com/JuliaStats/MixedModels.jl/issues/194
    data = (
        a = [1.55945122,0.004391538,0.005554163,-0.173029772,4.586284429,0.259493671,-0.091735715,5.546487603,0.457734831,-0.030169602],
        b = [0.24520519,0.080624178,0.228083467,0.2471453,0.398994279,0.037213859,0.102144973,0.241380251,0.206570975,0.15980803],
        c = PooledArray(["H","F","K","P","P","P","D","M","I","D"]),
        w1 = [20,40,35,12,29,25,65,105,30,75],
        w2 = [0.04587156,0.091743119,0.080275229,0.027522936,0.066513761,0.05733945,0.149082569,0.240825688,0.068807339,0.172018349],
    )

    #= no need to fit yet another model without weights, but here are the reference values from lme4
    m1 = fit(MixedModel, @formula(a ~ 1 + b + (1|c)), data)
    @test m1.θ ≈ [0.0]
    @test stderror(m1) ≈  [1.084912, 4.966336] atol = 1.e-4
    @test vcov(m1) ≈ [1.177035 -4.802598; -4.802598 24.664497] atol = 1.e-4
    =#

    m2 = fit(MixedModel, @formula(a ~ 1 + b + (1|c)), data, wts = data.w1)
    @test m2.θ ≈ [0.295181729258352]  atol = 1.e-4
    @test stderror(m2) ≈  [0.9640167, 3.6309696] atol = 1.e-4
    @test vcov(m2) ≈ [0.9293282 -2.557527; -2.5575267 13.183940] atol = 1.e-4
end

@testset "unifying ReMat eltypes" begin
    sleepstudy = dataset(:sleepstudy)

    re = LinearMixedModel(@formula(reaction ~ 1 + days + (1|subj) + (days|subj)), sleepstudy).reterms
    # make sure that the eltypes are still correct
    # otherwise this test isn't checking what it should be
    @test eltype(sleepstudy.days) == Int8
    @test eltype(sleepstudy.reaction) == Float64

    # use explicit typeof() and == is to remind us that things may break
    # if we change things and don't check their type implications now
    # that we're starting to support a non trivial type hierarchy
    @test typeof(re) == Vector{AbstractReMat{Float64}}
end


text:
struct VarType
  field :: Symbol
  type :: Type
  dim :: Int64
end


text:
module MapUtils
using Plots

export create_data_dict

"""
a Map is a representation of a 2D scatter plot with labels and sizes
"""
type ScatterMap
	xlims::Tuple{Float64, Float64}
	ylims::Tuple{Float64, Float64}
	xys::Array{Tuple{Float64, Float64}, 1}
	colors::Array{ColorTypes.RGBA{Float64}, 1}
	sizes::Array{Float64, 1}
	labels::Array{Tuple{Float64, Float64, String}, 1}
end
ScatterMap(xlims::Tuple{Float64, Float64}, ylims::Tuple{Float64, Float64}) = ScatterMap(xlims, ylims, [], [], [], [])
function ScatterMap(xys::Array{Tuple{Float64, Float64}, 1},
		 color_scheme::Symbol,
		 community_labels::Array{Int64, 1},
		 sizes::Array{Float64, 1},
		 labels::Array{Tuple{Float64, Float64, String}, 1})

	# defining min and max for the xy data
	xs = [xy[1] for xy in xys]
	ys = [xy[2] for xy in xys]
	xlims = (minimum(xs), maximum(xs))
	ylims = (minimum(ys), maximum(ys))

	# mapping communities to colors, creating the color scheme
	n_labels = maximum(community_labels)
	cg = ColorGradient(color_scheme)
	color_vec = [cg[i] for i in linspace(0, 1, n_labels)]
	color_data = [color_vec[i] for i in community_labels]

	ScatterMap(xlims, ylims, xys, color_data, sizes, labels)
end

"""
A TiledMap is a dictionary of the quadtree of nested maps indexed by
(zoom_level, x_index, y_index)
"""
type TiledMap
	map_dict::Dict{Tuple{Int64, Int64, Int64}, ScatterMap}
	depth::Int64
end
TiledMap(depth::Int64) = TiledMap(Dict{Tuple{Int64, Int64, Int64}, ScatterMap}(), depth)

"""
This constructor takes in a number of zoom levels, a set of xlims and ylims,
an array of tuples representing the data, an array of colors representing
some factor of the data, an array of sizes representing the data,
and a set of labels.
It returnes a TiledMap object, containing notably a dictionary indexed by
(zoom_level, xindex, yindex).
"""
function TiledMap(zoom_levels::Int64,
			  imap::ScatterMap,
			  markersize::Array{Float64, 1})
	# initializing an empty tilemap
	tm = TiledMap(zoom_levels)

	marker_size_limit = quantile(markersize, [max(0,1-10*4^(zoom_level+1)/1000000) for zoom_level in 0:zoom_levels])

	# we first need to initialize the dictionary so we can push to any tile.
	# we also calculate the boundaries for each tile in this
	# fist loop

	for zoom_level in 1:zoom_levels
		nb_subs = 2^(zoom_level - 1)
		delta_x = (imap.xlims[2] - imap.xlims[1]) / nb_subs
		delta_y = (imap.ylims[2] - imap.ylims[1]) / nb_subs
		for xindex in 1:nb_subs, yindex in 1:nb_subs
			local_xlims = (imap.xlims[1] + (xindex-1)*delta_x,  imap.xlims[1] + (xindex)*delta_x)
			local_ylims = (imap.ylims[2] - (yindex)*delta_y, imap.ylims[2] - (yindex-1)*delta_y)
			# this creates empty arrays for all the map properties
			tm.map_dict[(zoom_level, xindex, yindex)] = ScatterMap(local_xlims, local_ylims)
		end
	end

	for zoom_level in 1:zoom_levels
		nb_subs = 2^(zoom_level-1)
		delta_x = (imap.xlims[2] - imap.xlims[1]) / nb_subs
		delta_y = (imap.ylims[2] - imap.ylims[1]) / nb_subs
		for xindex in 1:nb_subs, yindex in 1:nb_subs
			for i in 1:size(imap.xys, 1)
				(xmin, xmax) = tm.map_dict[(zoom_level, xindex, yindex)].xlims
				(ymin, ymax) = tm.map_dict[(zoom_level, xindex, yindex)].ylims

				if xmin - delta_x/10 <= imap.xys[i][1] && imap.xys[i][1] <= xmax+delta_x/10 &&
					ymin -delta_x/10 <= imap.xys[i][2] && imap.xys[i][2] <= ymax+delta_x/10
					push!(tm.map_dict[(zoom_level, xindex, yindex)].xys, (imap.xys[i][1], imap.xys[i][2]))
					push!(tm.map_dict[(zoom_level, xindex, yindex)].colors, imap.colors[i])
					push!(tm.map_dict[(zoom_level, xindex, yindex)].sizes, imap.sizes[i])
				end

				# NB: for text, to prevent clipping, we also add all the labels
				# in the neighboring tiles!
				if markersize[i] > marker_size_limit[max(1, zoom_level-1)] +1 &&
					xmin - delta_x < imap.xys[i][1] &&
					imap.xys[i][1] <= xmax + delta_x &&
					ymin - delta_y <= imap.xys[i][2] &&
					imap.xys[i][2] < ymax + delta_y

					push!(tm.map_dict[(zoom_level, xindex, yindex)].labels,
	   					(imap.labels[i]))
				end
			end
		end
	end
	tm
end

function plot_sm(sm::ScatterMap, filename::String; kwargs...)
	scatter(sm.xys,
	 color = sm.colors,
	 markersize = sm.sizes,
	 xlims = sm.xlims,
	 ylims = sm.ylims,
	 grid = false,
	 axis = false,
	 legend = false,
	 margin = 0mm;
	 kwargs...);
	annotate!(sm.labels);
	png(filename)
end
"""
returns an empty png!
"""
function plot_sm(filename::String)
	scatter([(0,0)], grid = false, axis = false, legend = false, margin = 0mm);
	png(filename)
end

function plot_tm(tm::TiledMap, base_directory::String; kwargs...)
	file_number, directory_index = 0, 0

	# normalizing directory string and checking for existence
	endswith(base_directory, "/") && (base_directory = base_directory[1:end-1])
	isdir(base_directory) || error(base_directory * " is not a directory")

	isdir(base_directory * "/TileGroup" * string(directory_index)) || mkdir(base_directory * "/TileGroup" * string(directory_index))

	for d in 1:tm.depth
		nb_subs = 2^(d-1)
		for xi in 1:nb_subs, yi in 1:nb_subs
			filename = base_directory * "/TileGroup" * string(directory_index) * "/"  *
				string(d-1) * "-" * string(xi-1) * "-" * string(yi-1)
			if length(tm.map_dict[d, xi, yi].xys) > 0
				plot_sm(tm.map_dict[d, xi, yi], filename; kwargs...)
			else
				plot_sm(filename) # empty plot
			end
			file_number += 1
			if Int(floor(file_number/1024)) != directory_index
				directory_index += 1
				isdir(base_directory * "/TileGroup" * string(directory_index)) || mkdir(base_directory * "/TileGroup" * string(directory_index))
			end
		end
	end
end
end


text:
using Polyhedra
using SwitchOnSafety

@testset "Conitope" begin
    p = polyhedron(Conitope(2, [[1.0, 1.0]]))
    @test nallrays(p) == 0
    @test npoints(p) == 4
end


text:
"""
    min_occupancy(x, m=300)

Given an indicator signal this function uses run-length encoding to enforce minimum
lengths. If an indicator segment is less than the minimum it is added to the length
of the previous segment.
"""
function min_occupancy(x, m=300)
    # Use run-length encoding to find lengths and enforce minimum
    a, b = rle(x)
    a_new, b_new = Bool[], Int[]
    for i in eachindex(b)
        if b[i] > m || isempty(b_new)
            push!(a_new, a[i])
            push!(b_new, b[i])
        else
            b_new[end] += b[i]
        end
    end
    return BitVector(inverse_rle(a_new, b_new))
end

"""
    breath_availability(x)

Implements the breathing availbaility indicator from Long-Term Sleep Assessment by
Unobtrusive Pressure Sensor Arrays - 2018, Soleimani.

Expects an array with dimensions n_mats x 8 x samples. Returns the maximum distance
between mattress sensors and the mode of those distances.
"""
function breath_availability(x::AbstractArray{T, 3}) where T
    # Needs unfiltered data in a 9x8xN format in the proper order
    out = Vector{T}(undef, size(x, 3))
    for i in eachindex(out)
        active = findall(view(x, :, :, i) .> (2046 * 0.4))
        if !isempty(active)
            out[i] = pairwise(Euclidean(), Tuple.(active)) |> maximum
        else
            out[i] = Inf
        end
    end
    return mode(out), out
end

"""
    occupancy_detection(x::AbstractVector; kwargs)

Implements signle sensor simple occupancy detection from Long-Term Sleep Assessment by
Unobtrusive Pressure Sensor Arrays - 2018, Soleimani.

# Arguments
 - `β=350` : Average value of a sensor from an empty bed.
 - `m=300` : Minimum number of samples to be declared an occupancy.
 - `n=4` : Sensitivity parameter
"""
function occupancy_detection(x::AbstractVector; β=350, m=300, n=4)
    τ = cityblock(extrema(x)...) / n # Sareh uses n=4 in her code
    occupancy = (x .> (β + τ))
    return any(occupancy) ? min_occupancy(occupancy, m) : occupancy
end

"""
    occupancy_detection(x::AbstractMatrix; kwargs)

Implements simple occupancy detection from Long-Term Sleep Assessment by Unobtrusive
Pressure Sensor Arrays - 2018, Soleimani.

# Arguments
 - `max_dist=false` : Use breathing availability signal from the maximum distance method.
 - `β=25000` : Average value of the sum of all sensors from an empty bed.
 - `m=300` : Minimum number of samples to be declared an occupancy.
 - `n=4` : Sensitivity parameter
"""
function occupancy_detection(x::AbstractMatrix; max_dist=false, β=25000, m=300, n=4)
    Z = sum(x, dims=2) |> vec
    τ = cityblock(extrema(Z)...) / n # Sareh uses n=4 in her code

    if max_dist # Sareh's post-processing
        k, breath = breath_availability(reshape_psm(x))
        occupancy = (Z .> (β + τ)) .& (0.85*k .< breath .< 1.15*k)
    else
        occupancy = (Z .> (β + τ))
    end
    return any(occupancy) ? min_occupancy(occupancy, m) : occupancy
end


text:
# Julia wrapper for header: /opt/ffmpeg/include/libavformat/avformat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_packet,
    av_append_packet,
    av_stream_get_r_frame_rate,
    av_stream_set_r_frame_rate,
    av_stream_get_parser,
    av_stream_get_end_pts,
    av_format_get_probe_score,
    av_format_get_video_codec,
    av_format_set_video_codec,
    av_format_get_audio_codec,
    av_format_set_audio_codec,
    av_format_get_subtitle_codec,
    av_format_set_subtitle_codec,
    av_format_get_metadata_header_padding,
    av_format_set_metadata_header_padding,
    av_format_get_opaque,
    av_format_set_opaque,
    av_format_get_control_message_cb,
    av_format_set_control_message_cb,
    av_format_inject_global_side_data,
    av_fmt_ctx_get_duration_estimation_method,
    #avformat_version,
    avformat_configuration,
    avformat_license,
    av_register_all,
    av_register_input_format,
    av_register_output_format,
    avformat_network_init,
    avformat_network_deinit,
    av_iformat_next,
    av_oformat_next,
    avformat_alloc_context,
    avformat_free_context,
    avformat_get_class,
    avformat_new_stream,
    av_stream_get_side_data,
    av_new_program,
    avformat_alloc_output_context,
    avformat_alloc_output_context2,
    av_find_input_format,
    av_probe_input_format,
    av_probe_input_format2,
    av_probe_input_format3,
    av_probe_input_buffer2,
    av_probe_input_buffer,
    avformat_open_input,
    av_demuxer_open,
    av_find_stream_info,
    avformat_find_stream_info,
    av_find_program_from_stream,
    av_find_best_stream,
    av_read_packet,
    av_read_frame,
    av_seek_frame,
    avformat_seek_file,
    av_read_play,
    av_read_pause,
    av_close_input_file,
    avformat_close_input,
    av_new_stream,
    av_set_pts_info,
    avformat_write_header,
    av_write_frame,
    av_interleaved_write_frame,
    av_write_uncoded_frame,
    av_interleaved_write_uncoded_frame,
    av_write_uncoded_frame_query,
    av_write_trailer,
    av_guess_format,
    av_guess_codec,
    av_get_output_timestamp,
    av_hex_dump,
    av_hex_dump_log,
    av_pkt_dump2,
    av_pkt_dump_log2,
    av_codec_get_id,
    av_codec_get_tag,
    av_codec_get_tag2,
    av_find_default_stream_index,
    av_index_search_timestamp,
    av_add_index_entry,
    av_url_split,
    av_dump_format,
    av_get_frame_filename,
    av_filename_number_test,
    av_sdp_create,
    av_match_ext,
    avformat_query_codec,
    avformat_get_riff_video_tags,
    avformat_get_riff_audio_tags,
    avformat_get_mov_video_tags,
    avformat_get_mov_audio_tags,
    av_guess_sample_aspect_ratio,
    av_guess_frame_rate,
    avformat_match_stream_specifier,
    avformat_queue_attached_pictures

function av_get_packet(s,pkt,size::Integer)
    ccall((:av_get_packet,libavformat),Cint,(Ptr{AVIOContext},Ptr{AVPacket},Cint),s,pkt,size)
end

function av_append_packet(s,pkt,size::Integer)
    ccall((:av_append_packet,libavformat),Cint,(Ptr{AVIOContext},Ptr{AVPacket},Cint),s,pkt,size)
end

function av_stream_get_r_frame_rate(s)
    ccall((:av_stream_get_r_frame_rate,libavformat),AVRational,(Ptr{AVStream},),s)
end

function av_stream_set_r_frame_rate(s,r::AVRational)
    ccall((:av_stream_set_r_frame_rate,libavformat),Void,(Ptr{AVStream},AVRational),s,r)
end

function av_stream_get_parser(s)
    ccall((:av_stream_get_parser,libavformat),Ptr{AVCodecParserContext},(Ptr{AVStream},),s)
end

function av_stream_get_end_pts(st)
    ccall((:av_stream_get_end_pts,libavformat),Int64,(Ptr{AVStream},),st)
end

function av_format_get_probe_score(s)
    ccall((:av_format_get_probe_score,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_format_get_video_codec(s)
    ccall((:av_format_get_video_codec,libavformat),Ptr{AVCodec},(Ptr{AVFormatContext},),s)
end

function av_format_set_video_codec(s,c)
    ccall((:av_format_set_video_codec,libavformat),Void,(Ptr{AVFormatContext},Ptr{AVCodec}),s,c)
end

function av_format_get_audio_codec(s)
    ccall((:av_format_get_audio_codec,libavformat),Ptr{AVCodec},(Ptr{AVFormatContext},),s)
end

function av_format_set_audio_codec(s,c)
    ccall((:av_format_set_audio_codec,libavformat),Void,(Ptr{AVFormatContext},Ptr{AVCodec}),s,c)
end

function av_format_get_subtitle_codec(s)
    ccall((:av_format_get_subtitle_codec,libavformat),Ptr{AVCodec},(Ptr{AVFormatContext},),s)
end

function av_format_set_subtitle_codec(s,c)
    ccall((:av_format_set_subtitle_codec,libavformat),Void,(Ptr{AVFormatContext},Ptr{AVCodec}),s,c)
end

function av_format_get_metadata_header_padding(s)
    ccall((:av_format_get_metadata_header_padding,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_format_set_metadata_header_padding(s,c::Integer)
    ccall((:av_format_set_metadata_header_padding,libavformat),Void,(Ptr{AVFormatContext},Cint),s,c)
end

function av_format_get_opaque(s)
    ccall((:av_format_get_opaque,libavformat),Ptr{Void},(Ptr{AVFormatContext},),s)
end

function av_format_set_opaque(s,opaque)
    ccall((:av_format_set_opaque,libavformat),Void,(Ptr{AVFormatContext},Ptr{Void}),s,opaque)
end

function av_format_get_control_message_cb(s)
    ccall((:av_format_get_control_message_cb,libavformat),av_format_control_message,(Ptr{AVFormatContext},),s)
end

function av_format_set_control_message_cb(s,callback::av_format_control_message)
    ccall((:av_format_set_control_message_cb,libavformat),Void,(Ptr{AVFormatContext},av_format_control_message),s,callback)
end

function av_format_inject_global_side_data(s)
    ccall((:av_format_inject_global_side_data,libavformat),Void,(Ptr{AVFormatContext},),s)
end

function av_fmt_ctx_get_duration_estimation_method(ctx)
    ccall((:av_fmt_ctx_get_duration_estimation_method,libavformat),Cint,(Ptr{AVFormatContext},),ctx)
end

function avformat_version()
    ccall((:avformat_version,libavformat),UInt32,())
end

function avformat_configuration()
    ccall((:avformat_configuration,libavformat),Ptr{UInt8},())
end

function avformat_license()
    ccall((:avformat_license,libavformat),Ptr{UInt8},())
end

function av_register_all()
    ccall((:av_register_all,libavformat),Void,())
end

function av_register_input_format(format)
    ccall((:av_register_input_format,libavformat),Void,(Ptr{AVInputFormat},),format)
end

function av_register_output_format(format)
    ccall((:av_register_output_format,libavformat),Void,(Ptr{AVOutputFormat},),format)
end

function avformat_network_init()
    ccall((:avformat_network_init,libavformat),Cint,())
end

function avformat_network_deinit()
    ccall((:avformat_network_deinit,libavformat),Cint,())
end

function av_iformat_next(f)
    ccall((:av_iformat_next,libavformat),Ptr{AVInputFormat},(Ptr{AVInputFormat},),f)
end

function av_oformat_next(f)
    ccall((:av_oformat_next,libavformat),Ptr{AVOutputFormat},(Ptr{AVOutputFormat},),f)
end

function avformat_alloc_context()
    ccall((:avformat_alloc_context,libavformat),Ptr{AVFormatContext},())
end

function avformat_free_context(s)
    ccall((:avformat_free_context,libavformat),Void,(Ptr{AVFormatContext},),s)
end

function avformat_get_class()
    ccall((:avformat_get_class,libavformat),Ptr{AVClass},())
end

function avformat_new_stream(s,c)
    ccall((:avformat_new_stream,libavformat),Ptr{AVStream},(Ptr{AVFormatContext},Ptr{AVCodec}),s,c)
end

function av_stream_get_side_data(stream,_type::AVPacketSideDataType,size)
    ccall((:av_stream_get_side_data,libavformat),Ptr{UInt8},(Ptr{AVStream},AVPacketSideDataType,Ptr{Cint}),stream,_type,size)
end

function av_new_program(s,id::Integer)
    ccall((:av_new_program,libavformat),Ptr{AVProgram},(Ptr{AVFormatContext},Cint),s,id)
end

function avformat_alloc_output_context(format,oformat,filename)
    ccall((:avformat_alloc_output_context,libavformat),Ptr{AVFormatContext},(Ptr{UInt8},Ptr{AVOutputFormat},Ptr{UInt8}),format,oformat,filename)
end

function avformat_alloc_output_context2(ctx,oformat,format_name,filename)
    ccall((:avformat_alloc_output_context2,libavformat),Cint,(Ptr{Ptr{AVFormatContext}},Ptr{AVOutputFormat},Ptr{UInt8},Ptr{UInt8}),ctx,oformat,format_name,filename)
end

function av_find_input_format(short_name)
    ccall((:av_find_input_format,libavformat),Ptr{AVInputFormat},(Ptr{UInt8},),short_name)
end

function av_probe_input_format(pd,is_opened::Integer)
    ccall((:av_probe_input_format,libavformat),Ptr{AVInputFormat},(Ptr{AVProbeData},Cint),pd,is_opened)
end

function av_probe_input_format2(pd,is_opened::Integer,score_max)
    ccall((:av_probe_input_format2,libavformat),Ptr{AVInputFormat},(Ptr{AVProbeData},Cint,Ptr{Cint}),pd,is_opened,score_max)
end

function av_probe_input_format3(pd,is_opened::Integer,score_ret)
    ccall((:av_probe_input_format3,libavformat),Ptr{AVInputFormat},(Ptr{AVProbeData},Cint,Ptr{Cint}),pd,is_opened,score_ret)
end

function av_probe_input_buffer2(pb,fmt,filename,logctx,offset::Integer,max_probe_size::Integer)
    ccall((:av_probe_input_buffer2,libavformat),Cint,(Ptr{AVIOContext},Ptr{Ptr{AVInputFormat}},Ptr{UInt8},Ptr{Void},UInt32,UInt32),pb,fmt,filename,logctx,offset,max_probe_size)
end

function av_probe_input_buffer(pb,fmt,filename,logctx,offset::Integer,max_probe_size::Integer)
    ccall((:av_probe_input_buffer,libavformat),Cint,(Ptr{AVIOContext},Ptr{Ptr{AVInputFormat}},Ptr{UInt8},Ptr{Void},UInt32,UInt32),pb,fmt,filename,logctx,offset,max_probe_size)
end

function avformat_open_input(ps,filename,fmt,options)
    ccall((:avformat_open_input,libavformat),Cint,(Ptr{Ptr{AVFormatContext}},Ptr{UInt8},Ptr{AVInputFormat},Ptr{Ptr{AVDictionary}}),ps,filename,fmt,options)
    #ccall((:avformat_open_input,libavformat),Cint,(Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),ps,filename,fmt,options)
end

function av_demuxer_open(ic)
    ccall((:av_demuxer_open,libavformat),Cint,(Ptr{AVFormatContext},),ic)
end

function av_find_stream_info(ic)
    ccall((:av_find_stream_info,libavformat),Cint,(Ptr{AVFormatContext},),ic)
end

function avformat_find_stream_info(ic,options)
    ccall((:avformat_find_stream_info,libavformat),Cint,(Ptr{AVFormatContext},Ptr{Ptr{AVDictionary}}),ic,options)
end

function av_find_program_from_stream(ic,last,s::Integer)
    ccall((:av_find_program_from_stream,libavformat),Ptr{AVProgram},(Ptr{AVFormatContext},Ptr{AVProgram},Cint),ic,last,s)
end

function av_find_best_stream(ic,_type::AVMediaType,wanted_stream_nb::Integer,related_stream::Integer,decoder_ret,flags::Integer)
    ccall((:av_find_best_stream,libavformat),Cint,(Ptr{AVFormatContext},AVMediaType,Cint,Cint,Ptr{Ptr{AVCodec}},Cint),ic,_type,wanted_stream_nb,related_stream,decoder_ret,flags)
end

function av_read_packet(s,pkt)
    ccall((:av_read_packet,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVPacket}),s,pkt)
end

function av_read_frame(s,pkt)
    ccall((:av_read_frame,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVPacket}),s,pkt)
end

function av_seek_frame(s,stream_index::Integer,timestamp::Int64,flags::Integer)
    ccall((:av_seek_frame,libavformat),Cint,(Ptr{AVFormatContext},Cint,Int64,Cint),s,stream_index,timestamp,flags)
end

function avformat_seek_file(s,stream_index::Integer,min_ts::Int64,ts::Int64,max_ts::Int64,flags::Integer)
    ccall((:avformat_seek_file,libavformat),Cint,(Ptr{AVFormatContext},Cint,Int64,Int64,Int64,Cint),s,stream_index,min_ts,ts,max_ts,flags)
end

function av_read_play(s)
    ccall((:av_read_play,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_read_pause(s)
    ccall((:av_read_pause,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_close_input_file(s)
    ccall((:av_close_input_file,libavformat),Void,(Ptr{AVFormatContext},),s)
end

function avformat_close_input(s)
    ccall((:avformat_close_input,libavformat),Void,(Ptr{Ptr{AVFormatContext}},),s)
end

function av_new_stream(s,id::Integer)
    ccall((:av_new_stream,libavformat),Ptr{AVStream},(Ptr{AVFormatContext},Cint),s,id)
end

function av_set_pts_info(s,pts_wrap_bits::Integer,pts_num::Integer,pts_den::Integer)
    ccall((:av_set_pts_info,libavformat),Void,(Ptr{AVStream},Cint,UInt32,UInt32),s,pts_wrap_bits,pts_num,pts_den)
end

function avformat_write_header(s,options)
    ccall((:avformat_write_header,libavformat),Cint,(Ptr{AVFormatContext},Ptr{Ptr{AVDictionary}}),s,options)
end

function av_write_frame(s,pkt)
    ccall((:av_write_frame,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVPacket}),s,pkt)
end

function av_interleaved_write_frame(s,pkt)
    ccall((:av_interleaved_write_frame,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVPacket}),s,pkt)
end

function av_write_uncoded_frame(s,stream_index::Integer,frame)
    ccall((:av_write_uncoded_frame,libavformat),Cint,(Ptr{AVFormatContext},Cint,Ptr{AVFrame}),s,stream_index,frame)
end

function av_interleaved_write_uncoded_frame(s,stream_index::Integer,frame)
    ccall((:av_interleaved_write_uncoded_frame,libavformat),Cint,(Ptr{AVFormatContext},Cint,Ptr{AVFrame}),s,stream_index,frame)
end

function av_write_uncoded_frame_query(s,stream_index::Integer)
    ccall((:av_write_uncoded_frame_query,libavformat),Cint,(Ptr{AVFormatContext},Cint),s,stream_index)
end

function av_write_trailer(s)
    ccall((:av_write_trailer,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_guess_format(short_name,filename,mime_type)
    ccall((:av_guess_format,libavformat),Ptr{AVOutputFormat},(Ptr{UInt8},Ptr{UInt8},Ptr{UInt8}),short_name,filename,mime_type)
end

function av_guess_codec(fmt,short_name,filename,mime_type,_type::AVMediaType)
    ccall((:av_guess_codec,libavformat),Cint,(Ptr{AVOutputFormat},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},AVMediaType),fmt,short_name,filename,mime_type,_type)
end

function av_get_output_timestamp(s,stream::Integer,dts,wall)
    ccall((:av_get_output_timestamp,libavformat),Cint,(Ptr{AVFormatContext},Cint,Ptr{Int64},Ptr{Int64}),s,stream,dts,wall)
end

function av_hex_dump(f,buf,size::Integer)
    ccall((:av_hex_dump,libavformat),Void,(Ptr{FILE},Ptr{UInt8},Cint),f,buf,size)
end

function av_hex_dump_log(avcl,level::Integer,buf,size::Integer)
    ccall((:av_hex_dump_log,libavformat),Void,(Ptr{Void},Cint,Ptr{UInt8},Cint),avcl,level,buf,size)
end

function av_pkt_dump2(f,pkt,dump_payload::Integer,st)
    ccall((:av_pkt_dump2,libavformat),Void,(Ptr{FILE},Ptr{AVPacket},Cint,Ptr{AVStream}),f,pkt,dump_payload,st)
end

function av_pkt_dump_log2(avcl,level::Integer,pkt,dump_payload::Integer,st)
    ccall((:av_pkt_dump_log2,libavformat),Void,(Ptr{Void},Cint,Ptr{AVPacket},Cint,Ptr{AVStream}),avcl,level,pkt,dump_payload,st)
end

function av_codec_get_id(tags,tag::Integer)
    ccall((:av_codec_get_id,libavformat),Cint,(Ptr{Ptr{AVCodecTag}},UInt32),tags,tag)
end

function av_codec_get_tag(tags,id::AVCodecID)
    ccall((:av_codec_get_tag,libavformat),UInt32,(Ptr{Ptr{AVCodecTag}},AVCodecID),tags,id)
end

function av_codec_get_tag2(tags,id::AVCodecID,tag)
    ccall((:av_codec_get_tag2,libavformat),Cint,(Ptr{Ptr{AVCodecTag}},AVCodecID,Ptr{UInt32}),tags,id,tag)
end

function av_find_default_stream_index(s)
    ccall((:av_find_default_stream_index,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_index_search_timestamp(st,timestamp::Int64,flags::Integer)
    ccall((:av_index_search_timestamp,libavformat),Cint,(Ptr{AVStream},Int64,Cint),st,timestamp,flags)
end

function av_add_index_entry(st,pos::Int64,timestamp::Int64,size::Integer,distance::Integer,flags::Integer)
    ccall((:av_add_index_entry,libavformat),Cint,(Ptr{AVStream},Int64,Int64,Cint,Cint,Cint),st,pos,timestamp,size,distance,flags)
end

function av_url_split(proto,proto_size::Integer,authorization,authorization_size::Integer,hostname,hostname_size::Integer,port_ptr,path,path_size::Integer,url)
    ccall((:av_url_split,libavformat),Void,(Ptr{UInt8},Cint,Ptr{UInt8},Cint,Ptr{UInt8},Cint,Ptr{Cint},Ptr{UInt8},Cint,Ptr{UInt8}),proto,proto_size,authorization,authorization_size,hostname,hostname_size,port_ptr,path,path_size,url)
end

function av_dump_format(ic,index::Integer,url,is_output::Integer)
    ccall((:av_dump_format,libavformat),Void,(Ptr{AVFormatContext},Cint,Ptr{UInt8},Cint),ic,index,url,is_output)
end

function av_get_frame_filename(buf,buf_size::Integer,path,number::Integer)
    ccall((:av_get_frame_filename,libavformat),Cint,(Ptr{UInt8},Cint,Ptr{UInt8},Cint),buf,buf_size,path,number)
end

function av_filename_number_test(filename)
    ccall((:av_filename_number_test,libavformat),Cint,(Ptr{UInt8},),filename)
end

function av_sdp_create(ac,n_files::Integer,buf,size::Integer)
    ccall((:av_sdp_create,libavformat),Cint,(Ptr{Ptr{AVFormatContext}},Cint,Ptr{UInt8},Cint),ac,n_files,buf,size)
end

function av_match_ext(filename,extensions)
    ccall((:av_match_ext,libavformat),Cint,(Ptr{UInt8},Ptr{UInt8}),filename,extensions)
end

function avformat_query_codec(ofmt,codec_id::AVCodecID,std_compliance::Integer)
    ccall((:avformat_query_codec,libavformat),Cint,(Ptr{AVOutputFormat},AVCodecID,Cint),ofmt,codec_id,std_compliance)
end

function avformat_get_riff_video_tags()
    ccall((:avformat_get_riff_video_tags,libavformat),Ptr{AVCodecTag},())
end

function avformat_get_riff_audio_tags()
    ccall((:avformat_get_riff_audio_tags,libavformat),Ptr{AVCodecTag},())
end

function avformat_get_mov_video_tags()
    ccall((:avformat_get_mov_video_tags,libavformat),Ptr{AVCodecTag},())
end

function avformat_get_mov_audio_tags()
    ccall((:avformat_get_mov_audio_tags,libavformat),Ptr{AVCodecTag},())
end

function av_guess_sample_aspect_ratio(format,stream,frame)
    ccall((:av_guess_sample_aspect_ratio,libavformat),AVRational,(Ptr{AVFormatContext},Ptr{AVStream},Ptr{AVFrame}),format,stream,frame)
end

function av_guess_frame_rate(ctx,stream,frame)
    ccall((:av_guess_frame_rate,libavformat),AVRational,(Ptr{AVFormatContext},Ptr{AVStream},Ptr{AVFrame}),ctx,stream,frame)
end

function avformat_match_stream_specifier(s,st,spec)
    ccall((:avformat_match_stream_specifier,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVStream},Ptr{UInt8}),s,st,spec)
end

function avformat_queue_attached_pictures(s)
    ccall((:avformat_queue_attached_pictures,libavformat),Cint,(Ptr{AVFormatContext},),s)
end


text:
using Documenter, Contour

makedocs(
    modules = [Contour],
    format = Documenter.HTML(),
    sitename = "Contour.jl",
    pages = Any[
        "Introduction" => "index.md",
        "Tutorial" => "tutorial.md",
        "Reference" => "reference.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaGeometry/Contour.jl.git",
)


text:
@testset "Test renamefile" begin
    @test renamefile("textprocessing/df2latex.jl", "123.jl") == "textprocessing\\123.jl"
	@test renamefile(raw"D:\GoogleDrive\Sites\CGRG\doc-archive\content\en\docs\magtip", "new_folder") == "D:\\GoogleDrive\\Sites\\CGRG\\doc-archive\\content\\en\\docs\\new_folder"
end

max_stars_repo_name: CGRG-lab/HypertextTools

max_issues_repo_name: CGRG-lab/HypertextTools

max_forks_repo_name: CGRG-lab/HypertextTools