text:
# This file currently contains both the objective function and preprocessing
# steps, plus helper functions SSE and norm

# some global parameters -
BIG = 1e10
num_objs = 11

# How we calculate error
function SSE(model_values,data_values)
  y = model_values
  m = data_values
  # (max(y)-10)^2
  E = sum((m-y).^2)
end

# How we normalize the data
function norm_dat(data)
  min = minimum(data)
  max = maximum(data)
  normalized_data = (data-min)/(max-min)
end

#     #    #   # Prep objective data
# Load and Normalize each blot individually
fpath="/home/dbassen/Dropbox/server_swap_space/gen_2_model/poets/medici_blot_combined.jl"
include(fpath)
Fig1B = norm_dat(Fig1B)
Fig1F = norm_dat(Fig1F)
Fig2A = norm_dat(Fig2A)
Fig2E = norm_dat(Fig2E)
Fig3D = norm_dat(Fig3D)
Fig3F = norm_dat(Fig3F)
Fig4D = norm_dat(Fig4D)
Fig4E = norm_dat(Fig4E)
Fig4G = norm_dat(Fig4G)
Fig5A = norm_dat(Fig5A)
Fig5C = norm_dat(Fig5C)

#     #    #   # Initialize and  Run model to steady state
 # EXECUTE RunAllSimulations UP TO LOOP

 # Load data dictionary off disk

 # Initialize simulation functions
 fpath="/home/dbassen/Dropbox/server_swap_space/gen_2_model/src/new_SimulationFunctions.jl" # Loads up DataFile , SolveBalances , Balances, Kinetics, Control
 include(fpath)

 # Establish simulation parameters
 TSTART = 0.0;
 Ts = 0.1;
 TSTOP = 100.0;

 # Initialize intial data dictionary
 data_dictionary = DataFile(TSTART,TSTOP,Ts)

 # # # #

# update control parameter array here for now
GAIN_MASTER = 0.3
num_regs = 17;
control_parameter_array=data_dictionary["CONTROL_PARAMETER_ARRAY"]
for i in 1:num_regs
  control_parameter_array[i,1]=GAIN_MASTER # Gain for transcriptional activity
end
# TGFB3_inhibition: TGFB3_Ab inhibition; TGFB3_binding
control_parameter_array[1,1]=30.0 # GAIN # _tgfb3_Ab
# GSK3_P_inhibition_Binding_APC_AXIN_forw: GSK3_P inhibition; Binding_APC_AXIN_forw
control_parameter_array[2,1]=3.0 # GAIN # GSK3_P_inhibition_Binding_APC_AXIN_forw
control_parameter_array[2,2]=1.0 # ORDER# GSK3_P_inhibition_Binding_APC_AXIN_forw
# Active_TCF4_inhibition_Ecad: Ecadherin_Bcatenin inhibition; TCF4_Bcatenin_complex_formation
control_parameter_array[3,1]=0.5 # GAIN # Active_TCF4_inhibition_Ecad
# Active_TCF4_inhibition_APC: APC_AXIN_Bcatenin inhibition; TCF4_Bcatenin_complex_formation
control_parameter_array[4,1]=0.5 # GAIN # Active_TCF4_inhibition_APC
# activation_Ecadherin_by_Active_NFATc: Active_NFATc induction induction_gene_Ecadherin
control_parameter_array[5,1]=0.3 # GAIN # activation_Ecadherin_by_Active_NFATc
# repression_Ecadherin_by_pSmad: pSmad repression induction_gene_Ecadherin
control_parameter_array[6,1]=1.0 # GAIN # REP_ECAD_PSMAD
# repression_Ecadherin_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_Ecadherin
control_parameter_array[7,1]=2.0 # GAIN_SNAIL_SLUG_ECAD
# repression_Ecadherin_by_Active_LEF1: Active_LEF1 repression induction_gene_Ecadherin
control_parameter_array[8,1]=12.0# GAIN_ACTIVE_LEF1_ECAD
#  control_parameter_array[6,1]=0.3 # GAIN_BCAT_TGFB3
# activation_Ecadherin_by_ecad_virus: ecad_virus induction induction_gene_Ecadherin
control_parameter_array[9,1]=0.5 # ecad_virus ***
# activation_TGFB3_by_TCF4: Active_TCF4 induction induction_gene_TGFB3
control_parameter_array[10,1]=0.5 # GAIN_TCF4_TGFB3
# repression_TGFB3_by_YREG1: YREG1 repression induction_gene_TGFB3
control_parameter_array[11,1]=8.0 # YREG_TGFB3
# repression_YREG1_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_YREG1
control_parameter_array[12,1]=6.0 # SNAIL_YREG1
# activation_LEF1_by_pSmad: pSmad induction induction_gene_LEF1
control_parameter_array[13,1]=0.8# GAIN # _PSMAD_LEF1 # 0.5
control_parameter_array[13,2]=1.0# ORDER# _PSMAD_LEF1
# repression_LEF1_by_YREG1: YREG1 repression induction_gene_LEF1
control_parameter_array[14,1]=5.0 # YREG1_LEF1
# activation_vimentin_by_Active_LEF1: Active_LEF1 induction induction_gene_vimentin
control_parameter_array[15,1]=0.5# GAIN_LEF1_VIM
control_parameter_array[15,2]=2.0# ORDER_LEF1_VIM
# activation_SNAIL_SLUG_by_ss_virus: ss_virus induction induction_gene_SNAIL_SLUG
control_parameter_array[17,1]=1.0# SS_VIRUS ***

# # # #

 data_dictionary["CONTROL_PARAMETER_ARRAY"]=control_parameter_array

 # Also use edit_file_global here
 fpath = "/home/dbassen/Dropbox/server_swap_space/gen_2_model/src/edit_file_global.jl"
 include(fpath)
 UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
 UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)

  # False, this will break the simulations
#  const pa103=parameter_array[103]
#  const pa104=parameter_array[104]
#  const pa105=parameter_array[105]
#  const pa106=parameter_array[106]
#  const pa107=parameter_array[107]
#  const pa108=parameter_array[108]
#  const pa109=parameter_array[109]
#  const pa110=parameter_array[110]
#  const pa111=parameter_array[111]
#  const pa112=parameter_array[112]
#  const pa113=parameter_array[113]
#  const pa114=parameter_array[114]
#  const pa115=parameter_array[115]
#  const pa116=parameter_array[116]
#  const pa117=parameter_array[117]
#  const pa118=parameter_array[118]

 println("Running to steady state ... ")

 # # # STEADY STATE
 println(now())
 XSS = EstimateSteadyState(data_dictionary)
 println(now())
 initial_condition_array = vec(XSS);
 data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
 xss_data_dict = deepcopy(data_dictionary) # may not be necessary

#     #    #   # Prep data dictionary for rapid deconstruction / reassembly
len_cpa = length(data_dictionary["CONTROL_PARAMETER_ARRAY"])
len_rca = length(data_dictionary["RATE_CONSTANT_ARRAY"])
len_sca = length(data_dictionary["SATURATION_CONSTANT_ARRAY"])

# We cannot simply dump each array into the parameter array, since each
# sub array as a different shape (especially the SCA)
 # get lengths of data dictoinary components to be included in simulation

 # SCA is sparse, so let's gather linear indices for its non-zero terms
 i = 1 #index to start at
 sca_indices,sca_values = [],[]
 for i in 1:len_sca
   val = data_dictionary["SATURATION_CONSTANT_ARRAY"][i]
 if val !=0
   push!(sca_indices,i)
   push!(sca_values,val)
 end
 end
 # package for quick reconstruction
 # this needs to be reassembled at runtime anyway :/
 sca_pack = zip(sca_indices,sca_values)

  # CPA, let's give access to gain AND power
  cpa_gain = data_dictionary["CONTROL_PARAMETER_ARRAY"][:,1] # Gain
  cpa_order= data_dictionary["CONTROL_PARAMETER_ARRAY"][:,2] # Power

# Initialize parameter_array
 parameter_array = []
 append!(parameter_array,data_dictionary["RATE_CONSTANT_ARRAY"])
 append!(parameter_array,cpa_gain)
 append!(parameter_array,cpa_order)
 append!(parameter_array,sca_values)

 # get lengths of new arrays
 len_cpa_gain = length(cpa_gain)
 len_cpa_order = length(cpa_order)
 len_sca_values = length(sca_values)
 # array order RCA,CPAGAIN,CPAORDER,SCA
   # need indices for reconstruction
len_sca_prime = length(sca_values)
  A = len_rca + len_cpa_gain
  B = A + len_cpa_order
  C = B + len_sca_prime

# parameter_array gets sent to estimate estimate_ensemble as initial_state

#Evaluates the objective function values -
function objective_function(parameter_array)

 # Alias the species -
 # in order of appearance below!

# Descriptions include the CURRENT indexed value,
# NOT necessarily the correct one
# Indices are from the DataFile and are NOT automously updated!!!
 # LEF 1 protein (technically this not not total, due to active amt)
 index_lef = 58
 # LEF 1 active
 index_lef_act = 59
 # Snail_slug total protein
 index_ss = 62
 # ecadherin protein (unbound)
 index_ecad = 51
 # vimentin protein
 index_vim = 63
 # tgfb3 protein
 index_tgfb3 = 91
 # This is fit to ACTIVE_TCF4, data is IP from bcat and tcf-4
 index_bcat_tcf = 57 # i.e. this index is for Active_TCF4
 # Should probably be total protein, so sum index_bcat with index_tcf
 index_tcf = 56
 #
 #
 #

 # Alias species as they will appear in the output array

# indices for simulation data array sim_data
# 1 "edit_file_none.jl"
# 2 "edit_file_tgfb.jl"
# 3 "edit_file_mapkin.jl"
# 4 "edit_file_ss.jl"
# 5 "edit_file_ssecad.jl"
# 6 "edit_file_ss_tgfb3ab.jl"
# 7 "edit_file_ss_dnlef1.jl"
# 8 "edit_file_tgfb3.jl"
# 9 "edit_file_dnsmad.jl"
# 10 "edit_file_tgfb3_dnlef1.jl"
# 11 "edit_file_tgfb_tgfb3ab.jl"

#  MOD,BLOT = norm_dat(MOD),norm_dat(BLOT)

 # Execute simulations

 # Calculate the objective function array -
 obj_array = BIG*ones(num_objs,1)

 # Refactor data dictionary for model execution
 # need to specify implicit strucutre of parameter_array
 # structure: rca;cpa(1);cpa(2);sca
 data_dictionary["RATE_CONSTANT_ARRAY"] = parameter_array[1:len_rca]
 data_dictionary["CONTROL_PARAMETER_ARRAY"][:,1] = parameter_array[len_rca+1:A]
 data_dictionary["CONTROL_PARAMETER_ARRAY"][:,2] = parameter_array[A+1:B]
# Rebuild SCA
 sca_pack = zip(sca_indices,parameter_array[B+1:C]) # new values
 for (i,val) in sca_pack
 data_dictionary["SATURATION_CONSTANT_ARRAY"][i]=val
 end

### RUN TO NEW STEADY STATE!
# Make new xss_data_dict

# # # STEADY STATE
# WAIT, what initial condiiton array do we use???
# use last steady state, we should still be stored in the data_dictionary object
# as long as it is global
# # # # THIS DATA DICTIONARY CONTAINS THE MASTER STATE
# # # # FOR A GIVEN SET OF PARAMETERS
XSS = EstimateSteadyState(data_dictionary) # BE CAREFUL WITH DATA DICTIONARY NAME
initial_condition_array = vec(XSS);
data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
xss_data_dict = deepcopy(data_dictionary) # BE CAREFUL WITH DATA DICTIONARY NAME

### EXECUTE RunAllSimulations loop (without file writes)

# ??? This should be read off of ram SolveModelGeneric needs to be updated
simulations = ["edit_file_none.jl","edit_file_tgfb.jl","edit_file_mapkin.jl","edit_file_ss.jl","edit_file_ssecad.jl","edit_file_ss_tgfb3ab.jl","edit_file_ss_dnlef1.jl","edit_file_tgfb3.jl","edit_file_dnsmad.jl","edit_file_tgfb3_dnlef1.jl","edit_file_tgfb_tgfb3ab.jl"]
sim_data = [] # Keep data in RAM
for sim in simulations
  ##### What data_dictionary should be loaded here??? clearly the active one
  data_dictionary = deepcopy(xss_data_dict) # RENAME THIS TO RETAIN ORIGINAL call it data_dictionary_active
  (TSIM,X) =  SolveModelGenericPOETS(TSTART,TSTOP,Ts,data_dictionary,string("/home/dbassen/Dropbox/server_swap_space/gen_2_model/src/",sim))
  push!(sim_data,X)
  # now reset the data dictionary
end

###
t48 = 528 # 48 hours(24hrs = 264)
 #  sim_data[1][526,55] -> [t48,index_species]
 # specific objectives
 # EXP [conditions] RO [readouts]
 #1 EXP SS RO SS LEF

 MOD  = [sim_data[1][t48,index_ss],sim_data[1][t48,index_lef],
         sim_data[4][t48,index_ss],sim_data[4][t48,index_lef]]
 BLOT = [Fig1B[1],Fig1B[9],
         Fig1B[2],Fig1B[10]]
 MOD = norm_dat(MOD)
 obj_array[1]  = SSE(MOD,BLOT)

 #2 EXP SS DNLEF RO ECAD VIM
 MOD  = [sim_data[1][t48,index_ecad],sim_data[1][t48,index_vim],
         sim_data[4][t48,index_ecad],sim_data[4][t48,index_vim],
         sim_data[7][t48,index_ecad],sim_data[7][t48,index_vim]]
 BLOT = [Fig1F[1],Fig1F[5],
         Fig1F[2],Fig1F[6],
         Fig1F[4],Fig1F[8]]
 MOD = norm_dat(MOD)

 obj_array[2]  = SSE(MOD,BLOT)

 #3 EXP SS RO TGFB3
 MOD  = [sim_data[1][t48,index_tgfb3],
         sim_data[4][t48,index_tgfb3]]
 BLOT = [Fig2A[1],
         Fig2A[2]]
 MOD = norm_dat(MOD)
 obj_array[3]  = SSE(MOD,BLOT)

 #4 EXP SS TGFBab RO LEF ECAD VIM
 MOD  = [sim_data[1][t48,index_ecad],sim_data[1][t48,index_vim],sim_data[1][t48,index_lef],
         sim_data[4][t48,index_ecad],sim_data[4][t48,index_vim],sim_data[4][t48,index_lef],
         sim_data[6][t48,index_ecad],sim_data[6][t48,index_vim],sim_data[6][t48,index_lef]]
 BLOT = [Fig2E[6],Fig2E[11],Fig2E[1],
         Fig2E[7],Fig2E[12],Fig2E[2],
         Fig2E[10],Fig2E[15],Fig2E[5]]
 MOD = norm_dat(MOD)
 obj_array[4] = SSE(MOD,BLOT)

 #5 EXP SS ECAD RO BCAT TCF4
 # CAREFUL in setting index_bcat & index_tcf, total BCAT & TCF?
 # consistent with total TCF4 and free/nuclear bcat, but need to verify
 MOD  = [sim_data[1][t48,index_bcat_tcf],sim_data[1][t48,index_tcf]+sim_data[1][t48,index_bcat_tcf],
         sim_data[4][t48,index_bcat_tcf],sim_data[4][t48,index_tcf]+sim_data[4][t48,index_bcat_tcf],
         sim_data[5][t48,index_bcat_tcf],sim_data[5][t48,index_tcf]+sim_data[5][t48,index_bcat_tcf]]
 BLOT = [Fig3D[1],Fig3D[6],
         Fig3D[2],Fig3D[7],
         Fig3D[5],Fig3D[10]]
 MOD = norm_dat(MOD)
 obj_array[5]  = SSE(MOD,BLOT)

 #6 EXP SS ECAD RO TGFB3 ECAD VIM
 MOD  = [sim_data[1][t48,index_tgfb3],sim_data[1][t48,index_vim],sim_data[1][t48,index_ecad],
         sim_data[4][t48,index_tgfb3],sim_data[4][t48,index_vim],sim_data[4][t48,index_ecad],
         sim_data[5][t48,index_tgfb3],sim_data[5][t48,index_vim],sim_data[5][t48,index_ecad]]
 BLOT = [Fig3F[1],Fig3F[11],Fig3F[6],
         Fig3F[2],Fig3F[12],Fig3F[7],
         Fig3F[5],Fig3F[15],Fig3F[10]]
 MOD = norm_dat(MOD)
 obj_array[6]  = SSE(MOD,BLOT)

 #7 EXP TGFB3 RO SS
 MOD  = [sim_data[1][t48,index_tgfb3],
         sim_data[8][t48,index_tgfb3]]
 BLOT = [Fig4D[1],
         Fig4D[2]]
 MOD = norm_dat(MOD)
 obj_array[7]  = SSE(MOD,BLOT)

 #8 EXP TGFB3 DNSMAD RO LEF
 MOD  = [sim_data[1][t48,index_lef],
         sim_data[8][t48,index_lef],
         sim_data[9][t48,index_lef]]
 BLOT = [Fig4E[1],
         Fig4E[2],
         Fig4E[3]]
 MOD = norm_dat(MOD)
 obj_array[8]  = SSE(MOD,BLOT)

 #9 EXP TGFB3 DNSMAD DNLEF RO ECAD VIM
 # RESTART HERE ..
 MOD  = [sim_data[1][t48,index_ecad],sim_data[1][t48,index_vim],
         sim_data[8][t48,index_ecad],sim_data[8][t48,index_vim],
         sim_data[9][t48,index_ecad],sim_data[9][t48,index_vim],
         sim_data[10][t48,index_ecad],sim_data[10][t48,index_vim]]
 BLOT = [Fig4G[1],Fig4G[5],
         Fig4G[2],Fig4G[6],
         Fig4G[3],Fig4G[7],
         Fig4G[4],Fig4G[8]]
 MOD = norm_dat(MOD)
 obj_array[9]  = SSE(MOD,BLOT)

 #10 EXP TGFB2 DNMEK RO SS
 MOD  = [sim_data[1][t48,index_ss],
         sim_data[2][t48,index_ss],
         sim_data[3][t48,index_ss]]
 BLOT = [Fig5A[1],
         Fig5A[2],
         Fig5A[3]]
 MOD = norm_dat(MOD)
 obj_array[10] = SSE(MOD,BLOT)

 #11 EXP TGFB TGFB3_Abn RO LEF1
  # generate new experiment EXP TGFB + TGFB3_Ab RO: LEF1
 MOD  = [sim_data[1][t48,index_lef],
         sim_data[2][t48,index_lef],
         sim_data[11][t48,index_lef]] # patched
 BLOT = [Fig5C[1],
         Fig5C[2],
         Fig5C[3]]
 MOD = norm_dat(MOD)
 obj_array[11] = SSE(MOD,BLOT)

 # Constraints are implemented as a penaltly on obj value
 lambda_value = 100.0

 # How much do we violate the constraints?
 #violation_constraint_1 = 25 - (x-5.0)^2 - y^2
 #violation_constraint_2 = (x-8.0)^2 + (y-3.0)^2 - 7.7
 penaltly_array = zeros(num_objs)
 #penaltly_array[1] = lambda_value*(min(0,violation_constraint_1))^2
 #penaltly_array[2] = lambda_value*(min(0,violation_constraint_2))^2

 # return the obj_array -
 return obj_array+penaltly_array
end

#     #    #   # Provide other model specific functions for POETS
# FROM JV ...
# Generates new parameter array, given current array -
function neighbor_function(parameter_array)

  SIGMA = 0.05
  number_of_parameters = length(parameter_array)

  # calculate new parameters -
  new_parameter_array = parameter_array.*(1+SIGMA*randn(number_of_parameters))

  # Check the bound constraints -
  LOWER_BOUND = zeros(parameter_array)
  UPPER_BOUND = LOWER_BOUND + 100

  # return the corrected parameter arrays -
  return parameter_bounds_function(new_parameter_array,LOWER_BOUND,UPPER_BOUND)
end

function acceptance_probability_function(rank_array,temperature)

  return (exp(-rank_array[end]/temperature))

end

function cooling_function(temperature)

  # define my new temperature -
  alpha = 0.8 #0.9
  return alpha*temperature
end

# Helper functions -
function parameter_bounds_function(parameter_array,lower_bound_array,upper_bound_array)

  # reflection_factor -
  epsilon = 0.01

  # iterate through and fix the parameters -
  new_parameter_array = copy(parameter_array)
  for (index,value) in enumerate(parameter_array)

    lower_bound = lower_bound_array[index]
    upper_bound = upper_bound_array[index]
    if (value<lower_bound)
      new_parameter_array[index] = lower_bound+epsilon*upper_bound
    elseif (value>upper_bound)
      new_parameter_array[index] = upper_bound - epsilon*lower_bound
    end
  end

  # fix parameters
  new_parameter_array[len_rca+1]=30.0# GAIN_tgfb3_Ab

  # addition terms
#  parameter_array[103]=pa103
#  parameter_array[104]=pa104
#  parameter_array[105]=pa105
#  parameter_array[106]=pa106
#  parameter_array[107]=pa107
#  parameter_array[108]=pa108
#  parameter_array[109]=pa109
#  parameter_array[110]=pa110
#  parameter_array[111]=pa111
#  parameter_array[112]=pa112
#  parameter_array[113]=pa113
#  parameter_array[114]=pa114
#  parameter_array[115]=pa115
#  parameter_array[116]=pa116
#  parameter_array[117]=pa117
#  parameter_array[118]=pa118

  return new_parameter_array
end

function n_alpha(alpha,Tf=0.0001;Ti=1.0)
Ti = 1.0
(log(Tf)-log(Ti))/log(alpha)
end


text:
using StarAlgebras

@testset "affordable real degrees/dot" begin
    G = SmallPermGroups[10][2] # C₂⊕C₅
    tbl = SymbolicWedderburn.CharacterTable(Rational{Int}, G)
    chars = SymbolicWedderburn.irreducible_characters(tbl)

    @test all(isone ∘ SymbolicWedderburn.degree, chars)
    @test all(χ -> isone(dot(χ, χ)), chars)
    chars_fl = SymbolicWedderburn.Character{ComplexF64}.(chars)
    @test all(isone ∘ SymbolicWedderburn.degree, chars_fl)
    @test all(χ -> isapprox(dot(χ, χ), 1), chars_fl)

    charsR, _ = SymbolicWedderburn.affordable_real(chars)
    @test SymbolicWedderburn.degree.(charsR) ==
        [dot(χ, χ) for χ in charsR] ==
        [1,1,2,2,2,2]
    charsR_fl = SymbolicWedderburn.Character{Float64}.(charsR)
    @test all([1,1,2,2,2,2] .==
        SymbolicWedderburn.degree.(charsR_fl) .≈
        [dot(χ, χ) for χ in charsR_fl]
        )
end

@testset "Symmetry adapted basis" begin
    @testset "step by step" begin
        G = PermGroup(perm"(1,2)", perm"(1,2,3)")
        irr = SymbolicWedderburn.irreducible_characters(G)
        @test irr isa AbstractVector{<:SymbolicWedderburn.Character{<:Cyclotomic}}
        @test SymbolicWedderburn.degree.(irr) == [2,1,1]

        RG = let G = G
            b = StarAlgebras.Basis{UInt16}(collect(G))
            StarAlgebra(G, b, (length(b), length(b)))
        end

        µ, s = SymbolicWedderburn.minimal_rank_projection(irr[1], RG)
        @test µ isa AlgebraElement
        @test s

        @test SymbolicWedderburn.Character{Rational{Int}}(irr[1]) isa
        SymbolicWedderburn.Character{Rational{Int}}
        @test collect(values(SymbolicWedderburn.Character{Float64}(irr[1]))) isa Vector{Float64}

        mps, simple = SymbolicWedderburn.minimal_projection_system(irr, RG)
        @test all(simple)

        @test rank(float.(SymbolicWedderburn.matrix_projection(irr[1]))) == 2
        @test rank(float.(SymbolicWedderburn.matrix_projection(mps[1]))) == 1

        sa_basis_ssimple = symmetry_adapted_basis(Rational{Int}, G, Rational{Int}, semisimple=true)

        @test issimple.(sa_basis_ssimple) == [false, true]
        @test rank.(convert.(Matrix, sa_basis_ssimple)) == [2, 1]
        @test dot(
            multiplicity.(sa_basis_ssimple),
            SymbolicWedderburn.degree.(sa_basis_ssimple)
        ) == PermutationGroups.degree(G)

        sa_basis = symmetry_adapted_basis(Rational{Int}, G, Rational{Int}, semisimple=false)
        @test all(issimple.(sa_basis))
        @test rank.(convert.(Matrix{Float64}, sa_basis)) == [1,1]
    end

    C₄ = PermGroup([perm"(1,2,3,4)"])
    A₄ = PermGroup([perm"(1,2,3)", perm"(2,3,4)"])
    S₄ = PermGroup([perm"(1,2,3,4)", perm"(1,2)"])

    sa_basis = symmetry_adapted_basis(C₄, semisimple=true)
    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(C₄)
    sa_basis = symmetry_adapted_basis(A₄, semisimple=true)
    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(A₄)
    sa_basis = symmetry_adapted_basis(S₄, semisimple=true)
    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(S₄)

    sa_basis = symmetry_adapted_basis(ComplexF64, C₄, semisimple=true)
    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(C₄)
    sa_basis = symmetry_adapted_basis(ComplexF64, A₄, semisimple=true)
    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(A₄)
    sa_basis = symmetry_adapted_basis(ComplexF64, S₄, semisimple=true)
    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(S₄)

    sa_basis = symmetry_adapted_basis(Float64, C₄, semisimple=true)
    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(C₄)
    sa_basis = symmetry_adapted_basis(Float64, A₄, semisimple=true)
    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(A₄)
    sa_basis = symmetry_adapted_basis(Float64, S₄, semisimple=true)
    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(S₄)
end

@testset "Symmetry adapted basis: small groups" begin
    @time @testset "Rational" begin
        for ord = 2:30
            @testset "SmallGroup($ord, $n)" for (n, G) in enumerate(SmallPermGroups[ord])
                @testset "complex" begin
                    S = Rational{Int}

                    sa_basis = symmetry_adapted_basis(G, S, semisimple=true)

                    @test dot(SymbolicWedderburn.degree.(sa_basis), multiplicity.(sa_basis)) ==
                        PermutationGroups.degree(G)

                    @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(G)

                    S = if (ord, n) in ((21,1),)
                        Rational{BigInt}
                    else
                        Rational{Int}
                    end
                    sa_basis = symmetry_adapted_basis(G, S, semisimple=false)
                    for b in sa_basis
                        if issimple(b)
                            @test multiplicity(b) == size(b, 1)
                        else
                            @test multiplicity(b)*SymbolicWedderburn.degree(b) == size(b, 1)
                        end
                    end
                end

                @testset "real" begin
                    S = Rational{Int}

                    sa_basisR = symmetry_adapted_basis(Float64, G, S, semisimple=true)
                    @test dot(SymbolicWedderburn.degree.(sa_basisR), multiplicity.(sa_basisR)) ==
                        PermutationGroups.degree(G)
                    @test sum(first ∘ size, sa_basisR) == PermutationGroups.degree(G)

                    sa_basisR = symmetry_adapted_basis(Float64, G, S, semisimple=false)
                    for b in sa_basisR
                        if issimple(b)
                            @test multiplicity(b) == size(b, 1) || 2*multiplicity(b) == size(b, 1)
                            # the first condiditon doesn't hold for realified characters;
                        else
                            @test multiplicity(b)*SymbolicWedderburn.degree(b) == size(b, 1)
                        end
                    end
                end
            end
        end
    end

    @time @testset "Float64/Complex64" begin
        for ord = 2:30
            @testset "SmallGroup($ord, $n)" for (n, G) in enumerate(SmallPermGroups[ord])

                T = Float64
                sa_basis = symmetry_adapted_basis(T, G, semisimple=true)
                @test sa_basis isa Vector{<:AbstractMatrix{T}}

                @test all(sa_basis) do b
                    rank(b) == size(b, 1)
                end

                @test sum(first ∘ size, sa_basis) == PermutationGroups.degree(G)

                basisC = symmetry_adapted_basis(Complex{T}, G, semisimple=true)
                @test sum(first ∘ size, basisC) == PermutationGroups.degree(G)
            end
        end
    end
end

ext: jl


text:
"""
    start_fabin3(model)

Return a vector of FABIN 3 starting values (see Hägglund 1982).
Not available for ensemble models.
"""
function start_fabin3 end

# splice model and loss functions
function start_fabin3(
        model::AbstractSemSingle;
        kwargs...)
    return start_fabin3(
        model.observed,
        model.imply,
        model.optimizer,
        model.loss.functions...,
        kwargs...)
end

function start_fabin3(
        observed,
        imply,
        optimizer,
        args...;
        kwargs...)
    return start_fabin3(
        ram_matrices(imply),
        obs_cov(observed),
        obs_mean(observed))
end

# SemObservedMissing
function start_fabin3(
        observed::SemObservedMissing,
        imply,
        optimizer,
        args...;
        kwargs...)

    if !observed.em_model.fitted
        em_mvn(observed; kwargs...)
    end

    return start_fabin3(
        ram_matrices(imply),
        observed.em_model.Σ,
        observed.em_model.μ)
end

function start_fabin3(ram_matrices::RAMMatrices, Σ, μ)

    A_ind, S_ind, F_ind, M_ind, parameters =
        ram_matrices.A_ind,
        ram_matrices.S_ind,
        ram_matrices.F_ind,
        ram_matrices.M_ind,
        ram_matrices.parameters

    n_par = length(parameters)
    start_val = zeros(n_par)
    n_var, n_nod = ram_matrices.size_F
    n_latent = n_nod - n_var

    C_indices = CartesianIndices((n_nod, n_nod))

    # check in which matrix each parameter appears

    indices = Vector{CartesianIndex}(undef, n_par)

    start_val = zeros(n_par)

#=     in_S = length.(S_ind) .!= 0
    in_A = length.(A_ind) .!= 0
    A_ind_c = [linear2cartesian(ind, (n_nod, n_nod)) for ind in A_ind]
    in_Λ = [any(ind[2] .∈ F_ind) for ind in A_ind_c]

    if !isnothing(M)
        in_M = length.(M_ind) .!= 0
        in_any = in_A .| in_S .| in_M
    else
        in_any = in_A .| in_S
    end

    if !all(in_any)
        @warn "Could not determine fabin3 starting values for some parameters, default to 0."
    end =#

    # set undirected parameters in S
    for (i, S_ind) in enumerate(S_ind)
        if length(S_ind) != 0
            c_ind = C_indices[S_ind][1]
            if c_ind[1] == c_ind[2]
                if c_ind[1] ∈ F_ind
                    index_position = findall(c_ind[1] .== F_ind)
                    start_val[i] = Σ[index_position[1], index_position[1]]/2
                else
                    start_val[i] = 0.05
                end
            end # covariances stay 0
        end
    end

    # set loadings
    constants = ram_matrices.constants
    A_ind_c = [linear2cartesian(ind, (n_nod, n_nod)) for ind in A_ind]
    # ind_Λ = findall([is_in_Λ(ind_vec, F_ind) for ind_vec in A_ind_c])

    for i ∈ findall(.!(1:n_nod .∈ [F_ind]))
        reference = Int64[]
        indicators = Int64[]
        loadings = Symbol[]

        for (j, ind_c_vec) in enumerate(A_ind_c)
            for ind_c in ind_c_vec
                if (ind_c[2] == i) & (ind_c[1] ∈ F_ind)
                    push!(indicators, ind_c[1])
                    push!(loadings, parameters[j])
                end
            end
        end

        for ram_constant in constants
            if (ram_constant.matrix == :A) & (ram_constant.index[2] == i) & (ram_constant.index[1] ∈ F_ind)
                push!(loadings, Symbol(""))
                if isone(ram_constant.value)
                    push!(reference, ram_constant.index[1])
                else
                    push!(indicators, ram_constant.index[1])
                end
            end
        end

        reference = [findfirst(x -> x == ref, F_ind) for ref in reference]
        indicators = [findfirst(x -> x == ind, F_ind) for ind in indicators]

        # is there at least one reference indicator?
        if size(reference, 1) > 0
            if size(reference, 1) > 1
                @warn "You have more than 1 scaling indicator"
                reference = reference[1]
            else
                reference = reference[1]
            end

            for (j, indicator) in enumerate(indicators)
                if indicator != reference
                    instruments = indicators[.!(indicators .∈ [[reference; indicator]])]

                    s32 = Σ[instruments, indicator]
                    s13 = Σ[reference, instruments]
                    S33 = Σ[instruments, instruments]

                    if size(instruments, 1) == 1
                        temp = S33[1]/s13[1]
                        λ = s32[1]*temp/(s13[1]*temp)
                        start_val[isequal.(parameters, loadings[j])] .= λ
                    else
                        temp = S33\s13
                        λ = s32'*temp/(s13'*temp)
                        start_val[isequal.(parameters, loadings[j])] .= λ
                    end

                end
            end
        # no reference indicator:
        else
            reference = indicators[1]
            λ = zeros(size(indicators, 1)); λ[1] = 1.0
            for (j, indicator) in enumerate(indicators)
                if indicator != reference
                    instruments = indicators[.!(indicators .∈ [[reference; indicator]])]

                    s32 = Σ[instruments, indicator]
                    s13 = Σ[reference, instruments]
                    S33 = Σ[instruments, instruments]

                    if size(instruments, 1) == 1
                        temp = S33[1]/s13[1]
                        λ[j] = s32[1]*temp/(s13[1]*temp)
                    else
                        temp = S33\s13
                        λ[j] = s32'*temp/(s13'*temp)
                    end

                end
            end

            Σ_λ = Σ[indicators, indicators]
            D = λ*λ' ./ sum(λ.^2)
            θ = (I - D.^2)\(diag(Σ_λ - D*Σ_λ*D))

            # 3. psi
            Σ₁ = Σ_λ - Diagonal(θ)
            l₂ = sum(λ.^2)
            Ψ = sum( sum(λ.*Σ₁, dims = 1).*λ') ./ l₂^2

            λ = λ .* sign(Ψ) .* sqrt(abs(Ψ))

            for (j, indicator) ∈ enumerate(indicators)
                start_val[isequal.(parameters, loadings[j])] .= λ[j]
            end
        end
    end

    # set means
    if !isnothing(M_ind)
        for (i, M_ind) in enumerate(M_ind)
            if length(M_ind) != 0
                ind = M_ind[1]
                if ind[1] ∈ F_ind
                    index_position = findfirst(ind[1] .== F_ind)
                    start_val[i] = μ[index_position]
                end # latent means stay 0
            end
        end
    end

    return start_val
end

function is_in_Λ(ind_vec, F_ind)
    res = [!(ind[2] ∈ F_ind) & (ind[1] ∈ F_ind) for ind in ind_vec]
    return any(res)
end


text:
# This file stores data structures used in imputation and phasing.
# For JLSO compression, the relevant data structures are in compress.jl

"""
Data structure for recording haplotype mosaic of one strand:
`start[i]` to `start[i+1]` has haplotype `haplotypelabel[i]`
in `window[i]`. The haplotype label is the column index of
`CW_typed[window[i]].uniqueH`. `start[end]` to`length` has
haplotype `haplotypelabel[end]` in `window[end]`.
"""
struct HaplotypeMosaic
    length::Int64
    start::Vector{Int64}
    window::Vector{Int32}
    haplotypelabel::Vector{Int32}
end
HaplotypeMosaic(len) = HaplotypeMosaic(len, Int64[], Int32[], Int32[])

function push_Mosaic!(x::HaplotypeMosaic, y::Tuple{Int64, T}) where T <: Integer
    newstart, newlabel = y[1], y[2]
    xlen = length(x.start)
    if xlen != 0
        # check if start occurs before previous start position, since searching
        # breakpoints for 2 consecutive windows can cause "overlaps"
        xstart = last(x.start)
        if newstart ≤ xstart
            deleteat!(x.start, xlen)
            deleteat!(x.haplotypelabel, xlen)
        end
    end
    push!(x.start, newstart)
    push!(x.haplotypelabel, newlabel)
end

function push_Mosaic!(x::HaplotypeMosaic, y::Tuple{Int64, T, T}) where T <: Integer
    newstart, newlabel, newwindow = y[1], y[2], y[3]
    xlen = length(x.start)
    if xlen != 0
        # check if start occurs before previous start position, since searching
        # breakpoints for 2 consecutive windows can cause "overlaps"
        xstart = last(x.start)
        if newstart ≤ xstart
            deleteat!(x.start, xlen)
            deleteat!(x.haplotypelabel, xlen)
            deleteat!(x.window, xlen)
        end
    end
    push!(x.start, newstart)
    push!(x.haplotypelabel, newlabel)
    push!(x.window, newwindow)
end

"""
data structure for recording haplotype mosaic of two strands
"""
struct HaplotypeMosaicPair
    strand1::HaplotypeMosaic
    strand2::HaplotypeMosaic
end
HaplotypeMosaicPair(len) = HaplotypeMosaicPair(HaplotypeMosaic(len),
    HaplotypeMosaic(len))

# reset function for dynamic programming phasing
function initialize!(x::Vector{Vector{Vector{Tuple{Int32, Int32}}}})
    n = length(x)
    win = length(x[1])
    @inbounds for i in 1:n, w in 1:win
        empty!(x[i][w])
    end
end

# resize function for dynamic programming phasing
function resize!(x::Vector{Vector{Vector{Tuple{Int32, Int32}}}}, windows::Int)
    n = length(x)
    @inbounds for i in 1:n
        Base.resize!(x[i], windows)
        sizehint!(x[i], windows)
    end
end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
module Classes

using DataStructures
using MacroTools
using MacroTools:combinedef, combinestructdef
using InteractiveUtils: subtypes

export @class, Class, AbstractClass, isclass, classof, superclass, superclasses, issubclass, subclasses, absclass

#
# Functional interface to MacroTools' dict-based expression creation functions
# TBD: suggest adding these (without requiring dict) to MacroTools
#
function emit_struct(name::Symbol, supertype::Symbol, mutable::Bool, params::Vector, fields::Vector, ctors::Vector)
    fieldtups = [(tup[1], tup[2]) for tup in map(splitarg, fields)]
    d = Dict(:name=>name, :supertype=>supertype, :mutable=>mutable, :params=>params, :fields=>fieldtups, :constructors=>ctors)
    combinestructdef(d)
end

function emit_function(name::Symbol, body; args::Vector=[], kwargs::Vector=[], rtype=nothing, params::Vector=[], wparams::Vector=[])
    d = Dict(:name=>name, :args=>args, :kwargs=>kwargs, :body=>body, :params=>params, :whereparams=>wparams)
    if rtype !== nothing
        d[:rtype] = rtype
    end
    combinedef(d)
end

abstract type AbstractClass end            # supertype of all shadow class types
abstract type Class <: AbstractClass end   # superclass of all concrete classes

abs_symbol(cls::Symbol) = Symbol("Abstract", cls)

# Since nameof() doesn't cover all the cases we need, we define our own
typename(t) = t
typename(t::TypeVar) = t.name

# fieldnames(DataType)
# (:name, :super, :parameters, :types, :names, :instance, :layout, :size, :ninitialized, :uid, :abstract, :mutable, :hasfreetypevars,
# :isconcretetype, :isdispatchtuple, :isbitstype, :zeroinit, :isinlinealloc, Symbol("llvm::StructType"), Symbol("llvm::DIType"))
#
# if dtype.hasfreetypevars, dtype.types is like svec(XYZ<:ABC,...)

function _translate_ivar(d::Dict, ivar)
    if ! @capture(ivar, vname_::vtype_ | vname_)
        error("Expected field definition, got $ivar")
    end

    if vtype === nothing
        return ivar    # no type, nothing to translate
    end

    vtype = get(d, vtype, vtype)    # translate parameterized types
    return :($vname::$vtype)
end

function _translate_where(d::Dict, wparam::TypeVar)
    # supname = :Any
    supname = wparam.ub          # TBD: not sure this suffices
    name = wparam.name
    name = get(d, name, name)    # translate, if a type parameter, else pass through

    return :($name <: $supname)
end

# If a symbol is already a gensym, extract the symbol and re-gensym with it
regensym(s) = MacroTools.isgensym(s) ? gensym(Symbol(MacroTools.gensymname(s))) : gensym(s)

_cache = nothing

# Return info about a class in a named tuple
function _class_info(::Type{T}) where {T <: AbstractClass}
    global _cache
    _cache === nothing && (_cache = Dict())
    haskey(_cache, T) && return _cache[T]

    # @info "_class_info($T)"

    typ = (typeof(T) === UnionAll ? Base.unwrap_unionall(T) : T)

    # note: must extract symbol from type to create required expression
    ivars = (isabstracttype(typ) ? Expr[] : [:($vname::$(typename(vtype))) for (vname, vtype) in zip(fieldnames(typ), typ.types)])
    wheres = typ.parameters

    d = Dict(t.name=>regensym(t.name) for t in wheres)    # create mapping of type params to gen'd symbols
    ivars = [_translate_ivar(d, iv) for iv in ivars]      # translate types to use gensyms
    wheres = [_translate_where(d, w) for w in wheres]

    result = (wheres=wheres, ivars=ivars, super=superclass(typ))
    _cache[T] = result
    return result
end

"""
    superclass(t::Type{Class})

Returns the type of the concrete superclass of the given class, or `nothing`
for `Class`, which is the root of the class hierarchy.
"""
superclass(::Type{Class}) = nothing

"""
    superclasses(::Type{T}) where {T <: AbstractClass}

Returns a vector of superclasses from the superclass of the current class
to `Class`, in order.
"""
superclasses(::Type{Class}) = []

function superclasses(::Type{T}) where {T <: AbstractClass}
    super = superclass(T)
    [super, superclasses(super)...]
end

"""
    isclass(X)

Return `true` if `X` is a concrete subclass of `AbstractClass`, or is `Class`, which is abstract.
"""
isclass(any) = false

# Note that !isabstracttype(T) != isconcretetype(T): parameterized types return false for both
isclass(::Type{T}) where {T <: AbstractClass} = (T === Class || !isabstracttype(T))

"""
    issubclass(t1::DataType, t2::DataType)

Returns `true` if `t1` is a subclass of `t2`, else false.
"""
# identity
issubclass(::Type{T}, ::Type{T}) where {T <: AbstractClass} = true

issubclass(::Type{T1}, ::Type{T2}) where {T1 <: AbstractClass, T2 <: AbstractClass} = T1 in Set(superclasses(T2))

"""
    classof(::Type{T}) where {T <: AbstractClass}

Compute the concrete class associated with abstract class `T`, which must
be a subclass of `AbstractClass`.
"""
function classof(::Type{T}) where {T <: AbstractClass}
    if isclass(T)
        return T
    end

    # Abstract types should have only one concrete subtype
    concrete = filter(isconcretetype, subtypes(T))

    if length(concrete) == 1
        return concrete[1]
    end

    # Should never happen unless user manually creates errant subtypes
    error("Abstract class supertype $T has multiple concrete subtypes: $concrete")
end

"""
    subclasses(::Type{T}) where {T <: AbstractClass}

Compute the vector of subclasses for a given class.
"""
function subclasses(::Type{T}) where {T <: AbstractClass}
    # immediate supertype is "our" entry in the type hierarchy
    super = supertype(T)

    # collect immediate subclasses
    subs = [classof(t) for t in subtypes(super) if isabstracttype(t)]

    # recurse on subclasses
    return [subs; [subclasses(t) for t in subs]...]
end

"""
    absclass(::Type{T}) where {T <: AbstractClass}

Returns the abstract type associated with the concrete class `T`.
"""
function absclass(::Type{T}) where {T <: AbstractClass}
    isclass(T) ? supertype(T) : error("absclass(T) must be called on concrete classes; $T is abstract.")
end

function _argnames(fields)
    return [sym for (sym, arg_type, slurp, default) in map(splitarg, fields)]
end

# We generate two initializer functions: one takes all fields, cumulative through superclasses,
# and another initializes only locally-defined fields. This function produces either, depending
# on the fields passed by _constructors().
function _initializer(class, fields, wheres)
    args = _argnames(fields)
    assigns = [:(_self.$arg = $arg) for arg in args]
    T = gensym(class)

    funcdef = :(
        function $class(_self::$T, $(fields...)) where {$T <: $(abs_symbol(class)), $(wheres...)}
            $(assigns...)
            _self
        end
    )

    return funcdef
end

function _constructors(clsname, super, super_info, local_fields, all_fields, wheres)
    all_wheres = [super_info.wheres; wheres]
    init_all = _initializer(clsname, all_fields, all_wheres)
    inits = [init_all]

    # If clsname is a direct subclasses of Classes.Class, it has no fields
    # other than those defined locally, so the two methods would be identical.
    # In this case, we emit only one of them.
    if all_fields != local_fields
        init_local = _initializer(clsname, local_fields, wheres)
        push!(inits, init_local)
    end

    # extract parameter names from where clauses
    params = [(clause isa Expr ? clause.args[1] : clause) for clause in all_wheres]
    has_params = length(params) != 0

    args = _argnames(all_fields)
    body = has_params ? :(new{$(params...)}($(args...))) : :(new($(args...)))
    dflt = emit_function(clsname, body, args=all_fields, params=params, wparams=all_wheres, rtype=clsname)

    methods = [dflt]

    # Primarily for immutable classes, we emit a constructor that takes an instance
    # of the direct superclass and copies values when creating a new object.
    super_fields = super_info.ivars
    if length(super_fields) != 0
        super_args = [:(_super.$arg) for arg in _argnames(super_fields)]
        local_args = _argnames(local_fields)
        all_args = [super_args; local_args]

        body = has_params ? :(new{$(params...)}($(all_args...))) : :(new($(all_args...)))
        args = [:(_super::$super); local_fields]
        immut_init = emit_function(clsname, body; args=args, params=params, wparams=all_wheres, rtype=clsname)
        push!(methods, immut_init)
    end

    return methods, inits
end

function _defclass(clsname, supercls, mutable, wheres, exprs)
    wheres   = (wheres === nothing ? [] : wheres)
    # @info "clsname:$clsname supercls:$supercls mutable:$mutable wheres:$wheres exprs:$exprs"

    # partition expressions into constructors and field defs
    ctors  = Vector{Expr}()
    fields = Vector{Expr}()
    for ex in exprs
        try
            splitdef(ex)        # throws AssertionError if not a func def
            push!(ctors, ex)
        catch
            push!(fields, ex)
        end
    end

    super_info = _class_info(supercls)
    all_fields = [super_info.ivars; fields]
    all_wheres = [super_info.wheres; wheres]

    # add default constructors
    inner, outer = _constructors(clsname, supercls, super_info, fields, all_fields, wheres)
    append!(ctors, inner)

    abs_class = abs_symbol(clsname)
    abs_super = nameof(absclass(supercls))

    struct_def = emit_struct(clsname, abs_class, mutable, all_wheres, all_fields, ctors)

    # set mutability flag
    struct_def.args[1] = mutable

    result = quote
        abstract type $abs_class <: $abs_super end
        $struct_def
        $(outer...)
        Classes.superclass(::Type{$clsname}) = $supercls
        $clsname    # return the struct type
    end

    return result
end

macro class(elements...)
    if (mutable = (elements[1] == :mutable))
        elements = elements[2:end]
    end

    if (len = length(elements)) == 1                       # no fields defined
        name_expr = elements[1]
        definition = quote end
    elseif len == 2
        (name_expr, definition) = elements
    else
        error("Unrecognized form for @class definition: $elements")
    end

    # @info "name_expr: $name_expr, definition: $definition"

    # initialize the "captured" vars to avoid "unknown var" warnings
    cls = clsname = exprs = wheres = nothing

    @capture(definition, begin exprs__ end)

    # allow for optional type params and supertype
    if ! (@capture(name_expr, ((cls_{wheres__} | cls_) <: supername_) | (cls_{wheres__} | cls_)) && cls isa Symbol)
        error("Unrecognized class name expression: `$name_expr`")
    end

    supername = (supername === nothing ? :Class : supername)

    # __module__ is a "hidden" arg passed to macros with the caller's Module
    expr = _defclass(cls, getproperty(__module__, supername), mutable, wheres, exprs)
    return esc(expr)
end

end # module


text:
using RecipesBase: apply_recipe
using Statistics

APERTURES = [
    CircularAperture(3, 3, 3),
    CircularAnnulus(3, 3, 2, 4),
    EllipticalAperture(3, 3, 4, 2, 45),
    EllipticalAnnulus(3, 3, 3, 4, 2, -26),
    RectangularAperture(3, 3, 3, 4, 15),
    RectangularAnnulus(3, 3, 3, 4, 2, -5)
]

@testset "Aperture Plots - $(typeof(ap))" for ap in APERTURES
    rec = apply_recipe(Dict{Symbol,Any}(), ap)
    for i in length(rec)
        @test getfield(rec[i], 1) == Dict{Symbol,Any}(
            :seriescolor => :red,
            :label => "",
            :seriestype => :path,
            :aspect_ratio => :equal)

        # test to make sure our position is correct (should be +0.5 the given (x,y))
        x, y = rec[i].args
        @test mean(x[1:end - 1]) ≈ 3.5
        @test mean(y[1:end - 1]) ≈ 3.5
    end
end


text:
type Packet
    ptr::Ptr{Void}

    function Packet(ptr::Ptr{Void})
        p = new(ptr)
        finalizer(p, destroy)
        p
    end
end

function Packet()
    Packet(ccall((:sfPacket_create, libcsfml_network), Ptr{Void}, ()))
end

function copy(packet::Packet)
    return Packet(ccall((:sfPacket_copy, libcsfml_network), Ptr{Void}, (Ptr{Void},), packet.ptr))
end

function destroy(packet::Packet)
    ccall((:sfPacket_destroy, libcsfml_network), Void, (Ptr{Void},), packet.ptr)
end

function clear(packet::Packet)
    ccall((:sfPacket_clear, libcsfml_network), Void, (Ptr{Void},), packet.ptr)
end

function get_data_size(packet::Packet)
    return ccall((:sfPacket_getDataSize, libcsfml_network), Int64, (Ptr{Void},), packet.ptr)
end

function read_bool(packet::Packet)
    return Bool(ccall((:sfPacket_readBool, libcsfml_network), UInt8, (Ptr{Void},), packet.ptr))
end

function read_string(packet::Packet)
    str = ""
    str = bytestring(ccall((:sjPacket_readString, "libjuliasfml"), Ptr{Cchar}, (Ptr{Void}, Ptr{Cchar},), packet.ptr, str))
    return str
end

function read_int(packet::Packet)
    return Int(ccall((:sfPacket_readInt32, libcsfml_network), Int32, (Ptr{Void},), packet.ptr))
end

function read_float(packet::Packet)
    return ccall((:sfPacket_readFloat, libcsfml_network), Cfloat, (Ptr{Void},), packet.ptr)
end

function read_double(packet::Packet)
    return ccall((:sfPacket_readDouble, libcsfml_network), Cdouble, (Ptr{Void},), packet.ptr)
end

function write(packet::Packet, value::Bool)
    ccall((:sfPacket_writeBool, libcsfml_network), Void, (Ptr{Void}, Int32,), packet.ptr, value)
end

function write(packet::Packet, val::Integer)
    ccall((:sfPacket_writeInt32, libcsfml_network), Void, (Ptr{Void}, Int32,), packet.ptr, val)
end

function write(packet::Packet, val::Cfloat)
    ccall((:sfPacket_writeFloat, libcsfml_network), Void, (Ptr{Void}, Cfloat,), packet.ptr, val)
end

function write(packet::Packet, val::Cdouble)
    ccall((:sfPacket_writeDouble, libcsfml_network), Void, (Ptr{Void}, Cdouble,), packet.ptr, val)
end

function write(packet::Packet, string::AbstractString)
    ccall((:sfPacket_writeString, libcsfml_network), Void, (Ptr{Void}, Ptr{Cchar},), packet.ptr, string)
end


text:
function parametrized_strogatz(du,u,p,t)
    #NumberOfOscillators,CouplingConstant, and NaturalFreqs
    N,K,w = p
    ##Centroid
    du[1] = 0
    for i in 1:N
        du[1] += du[i+1]
    end
    du[1] = du[1]/N
    #Oscillators
    for i in 1:N
        du[i+1] = (1.0-abs2(u[i+1])+w[i]*im)*u[i+1] + K*(u[1]-u[i+1])
    end
end

function getparameters_strogatz(N,K,γ,distribution)
    distribution == "Uniform" ? nothing : error("$(distribution) is not yet supported")
    j = collect(1:N)
    w = broadcast(j->((2j-N-1)/(N-1))*((pi*γ)/2),j)
    p = (N=N, K=K, w=w)
    return p
end

function getstartpoints_strogatz(N)
    u0 = ComplexF64[]
    d = Uniform(-1,1)
    for i in 1:N
        push!(u0,rand(d)+rand(d)im)
    end
    pushfirst!(u0,mean(u0))
    return u0
end

function strogatz(N,K,T,fs,γ,distribution)
    tspan = (0.0,T)
    u0 = getstartpoints_strogatz(N)
    p = getparameters_strogatz(N,K,γ,distribution)
    prob = ODEProblem(parametrized_strogatz,u0,tspan,p)
    sol = solve(prob,RK4(),saveat=collect(0:1/fs:T),progress=true)
end


text:
# 2D diffusion problem

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test,OrdinaryDiffEq
using ModelingToolkit: Differential

# Tests
@testset "Test 00: Dt(u(t,x)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dt = Differential(t)
    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 2.
    y_min = 0.
    y_max = 2.

    # 3D PDE
    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))

    analytic_sol_func(t,x,y) = exp(x+y)*cos(x+y+4t)
    # Initial and boundary conditions
    bcs = [u(t_min,x,y) ~ analytic_sol_func(t_min,x,y),
        u(t,x_min,y) ~ analytic_sol_func(t,x_min,y),
        u(t,x_max,y) ~ analytic_sol_func(t,x_max,y),
        u(t,x,y_min) ~ analytic_sol_func(t,x,y_min),
        u(t,x,y_max) ~ analytic_sol_func(t,x,y_max)]

    # Space and time domains
    domains = [t ∈ IntervalDomain(t_min,t_max),
            x ∈ IntervalDomain(x_min,x_max),
            y ∈ IntervalDomain(y_min,y_max)]
    pdesys = PDESystem([eq],bcs,domains,[t,x,y],[u(t,x,y)])

    # Method of lines discretization
    dx = 0.1; dy = 0.2
    for order in [2,4,6]
        discretization = MOLFiniteDifference([x=>dx,y=>dy],t;centered_order=order)
        prob = ModelingToolkit.discretize(pdesys,discretization)
        sol = solve(prob,Tsit5())

        # Test against exact solution
        # TODO: do this properly when sol[u] with reshape etc works
        @test sol.u[1][1] ≈ analytic_sol_func(sol.t[1],0.1,0.2)
        @test sol.u[1][2] ≈ analytic_sol_func(sol.t[1],0.2,0.2)

    end
end


text:
VERSION < v"0.7.0-beta2.199" && __precompile__(true)
module Nullables

if !isdefined(Base, :NullSafeTypes)
    include("nullable.jl")
else
    using Base: NullSafeTypes
end

export Nullable, NullException, isnull, unsafe_get

end # module


text:
########################################################################
# Codes for the paper:
# An Anthropocene Map of Genetic Diversity
#
# Andreia Miraldo, Sen Li, Michael K. Borregaard, Alexander Floréz-Rodriguéz, Shyam Gopalakrishnan, Mirneza Risvanovic, Zhiheng Wang, Carsten Rahbek, Katharine A. Marske & David Nogués-Bravo
#
# Submitted to Science, 2016
# Code in this file by Sen Li and Michael K. Borregaard
#########################################################################

# load the necessary library
using DataFrames

"""
A function to compute the GD (Genetic Diversity) value from a set of sequences. The algorithm compares all pairwise combinations of sequences and calculates the proportion of loci that differ between the pair. Returns a DataFrame with the identity of sequences, the lengths of each sequence, the overlap, and the computed pairwise divergence values.

**Parameters**
* 'species_seqs': A matrix where the rows are aligned genetic sequences, and columns are loci. Basepairs must be coded as 1, 2, 3 or 4, or with a 0 signifying that the locus is absent from the alignment.
"""
function compare_pairwise(species_seqs::Matrix{Int})
    num_seqs = size(species_seqs,1)  # count the number of sequences
    if(num_seqs < 2) # return an empty state if there is only one sequence
        return(DataFrame(seq1 = -1, seq2 = -1, length_seq1 = -1, length_seq2 = -1, overlap = -1., commons = -1, num_per_bp = -1.))
    end

    #Preinitialize containers
    seq1 = Int[]
    seq2 = Int[]
    length_seq1 = Int[]
    length_seq2 = Int[]
    overlap = Float64[]
    commons = Int[]
    num_per_bp = Float64[]
    for iter = 1:(num_seqs-1)
        for iter_sub = (iter + 1):num_seqs # For each pairwise combination of sequences
            pair_seqs = species_seqs[[iter, iter_sub], :]
            common_sites = minimum(pair_seqs, 1) .> 0 #Ensure that only sites existing in both sequences are compared (i.e. ignore loci missing from the alignment)
            if sum(common_sites) == 0.
                push!(num_per_bp, 0.)
            else
                num_segregating = sum(diff(pair_seqs[:, common_sites], 1) .!= 0)
                push!(num_per_bp, num_segregating/sum(common_sites))
            end
            push!(seq1, iter)
            push!(seq2, iter_sub)
            ls1 = sum(species_seqs[iter, :] .!= 0)
            ls2 = sum(species_seqs[iter_sub, :] .!= 0)
            push!(length_seq1, ls1)
            push!(length_seq2, ls2)
            push!(overlap, sum(common_sites)/max(ls1, ls2))  #Compute the overlap between the sequences
            push!(commons, sum(common_sites))
        end
    end
    DataFrame(seq1 = seq1, seq2 = seq2, length_seq1 = length_seq1, length_seq2 = length_seq2, overlap = overlap, commons = commons, num_per_bp = num_per_bp)
end

max_stars_repo_path: Miraldo et al. 2016/Compare_Pairwise_Function.jl

max_stars_repo_name: mkborregaard/codes-for-published-papers

max_issues_repo_name: DalavanCloud/codes-for-published-papers

max_forks_repo_name: DalavanCloud/codes-for-published-papers


text:
#!/usr/bin/env julia
include("src/QHull.jl")
using QHull

cd("test")
include("test/runtests.jl")


text:
import CUDA.CUDNN: cudnnDropoutForwardAD
using CUDA.CUDNN: cudnnDropoutBackward, handle
using AutoGrad: AutoGrad, @primitive1, value

@primitive1((cudnnDropoutForwardAD(x; xDesc, y, yDesc, dropoutDesc, reserveSpace),
             _dy, _y),
            ((x,y,dy,dx) = (value(x),value(_y),value(_dy),similar(x));
             cudnnDropoutBackward(handle(), dropoutDesc, yDesc, dy, xDesc, dx, reserveSpace, sizeof(reserveSpace));
             dx))


text:
using Flux
using Flux.Data: DataLoader
using BSON: @save, @load # for data saving
using ProgressMeter
using LinearAlgebra
using Dates

pdename = ARGS[1]
@info "Learning Green function of $pdename equation"

cp = "checkpoints_$(pdename)_$(now())"
mkdir(cp)
@info "$cp folder created"

# Load mat file
using MAT
file = matopen("greenlearning/examples/datasets/$(pdename).mat")
varnames = collect(names(file))
for varname in varnames
    val = read(file, varname)
    @eval $(Symbol(varname)) = $val
end
close(file)
X = reshape(X, 100)
Y = reshape(Y, 200)
xmin, xmax = X[begin], X[end]
ymin, ymax = Y[begin], Y[end]
dx = X[2]-X[1]
dy = Y[2]-Y[1]
Nu, Nf = length(X), length(Y)
@info "Train data loaded"

hidden_size = 50
G_network = Chain(
    Dense(2, hidden_size, relu),
    Dense(hidden_size, hidden_size, relu),
    Dense(hidden_size, hidden_size, relu),
    Dense(hidden_size, hidden_size, relu),
    Dense(hidden_size, 1)
)
U_hom_network = Chain(
    Dense(1, hidden_size, relu),
    Dense(hidden_size, hidden_size, relu),
    Dense(hidden_size, hidden_size, relu),
    Dense(hidden_size, hidden_size, relu),
    Dense(hidden_size, 1)
)
parameters = Flux.params(G_network, U_hom_network)
@info "Model" G_network U_hom_network

# integrator weights
trapezoidal_x = ones(Nu); trapezoidal_x[[begin,end]] .= 0.5
trapezoidal_y = ones(Nf); trapezoidal_y[[begin,end]] .= 0.5

function loss(U_train, F_train)
    l = 0.0
    U_hom_predict = [U_hom_network([x])[1] for x in X]
    # G_predict = [G_network([x,y])[1] for x in X, y in Y]
    G_predict = [[G_network([x,y])[1] for x in X] for y in Y]
    G_predict = hcat(G_predict...)
    num = size(U_train, 2)
    for n in 1:num
        @views u = U_train[:,n] # 100×1
        @views f = F_train[:,n] # 200×1
        Gf = dy*G_predict*(trapezoidal_y.*f)
        integ = (u-U_hom_predict-Gf).^2
        loss_u = dx*dot(integ, trapezoidal_x)
        u_L2 = dx*dot(u.^2, trapezoidal_x)
        l += loss_u/u_L2
    end
    return l/num
end

opt = ADAM()
@info "Optimiser" opt

batch_size = 100
train_loader = DataLoader((U, F), batchsize=batch_size, shuffle=true)
data_num = size(U, 2)
epochs = parse(Int, ARGS[2])
digit_epochs = length(string(epochs))
iteration_num = Int(epochs*batch_size/data_num)
@info "Training start" data_num epochs batch_size iteration_num

start = now()

learning = Dict([
    ("Equation name", pdename),
    ("Total epochs", epochs),
    ("Batch size", batch_size),
    ("Number of data", data_num),
    ("Current epoch number", 0),
    ("Current loss value", []),
    ("Start time", start),
    ("Time", 0)
])

skip = 5
counter = 0
@showprogress for epoch in 1:epochs
    for (U_train, F_train) in train_loader
        gs = gradient(()->loss(U_train, F_train), parameters)
        Flux.Optimise.update!(opt, parameters, gs)
    end
    global counter += 1
    if counter % skip == 0
        d = now()
        learning["Current epoch number"] = epoch
        learning["Time"] = d
        current_loss = loss(U, F)
        push!(learning["Current loss value"], current_loss)
        fname1 = "$cp/$(pdename)_description_$(lpad(epoch, digit_epochs, "0")).bson"
        fname2 = "$cp/$(pdename)_model_$(lpad(epoch, digit_epochs, "0")).bson"
        @save fname1 learning
        @save fname2 G_network U_hom_network opt
    end
end

d = now()
learning["Current epoch number"] = epochs
learning["Time"] = d
fname1 = "$cp/$(pdename)_description_$epochs.bson"
fname2 = "$cp/$(pdename)_model_$epochs.bson"
@save fname1 learning
@save fname2 G_network U_hom_network opt
@info "Training ended\nModel saved to $fname2"


text:
@safetestset "SVG API support in Renderer.Html" begin
  @safetestset "SVG API is available" begin
    using Genie
    using Genie.Renderer.Html

    @test svg() == "<svg></svg>"

    @test_throws UndefVarError clippath()
    @test clipPath() == "<clipPath></clipPath>"
  end;

end;


text:
#!/usr/bin/env julia

using Luxor, Test

function drawcircles(Δα)
    @layer begin
        setline(8)
        for r = 0:20:1500
            α = r / 100  + Δα
            circle(O + sqrt(r ) .* (cos(α), sin(α)), r, :stroke)
            settext("<span font = '16' foreground = 'darkblue' > $(string(r)) </span>",
                 O + (0, 8 - r ); markup = true)
        end
    end
end

fna(i, suf) = "test_snap_$i.$suf"

function testsnapshot(fname; cb = missing, scalefactor = missing, kwargs...)
    if ismissing(cb)
        sn = if ismissing(scalefactor)
            snapshot(;fname = fname, kwargs...)
        else
            snapshot(;fname = fname, scalefactor = scalefactor, kwargs...)
        end
        @test typeof(sn) == Drawing
    else
        if ismissing(scalefactor)
            sn = snapshot(;fname = fname, cb = cb, kwargs...)
            @test sn.width == boxwidth(cb) && sn.height == boxheight(cb)
        else
            sn = snapshot(;fname = fname, cb = cb, scalefactor = scalefactor, kwargs...)
            @test sn.width == round(boxwidth(cb) * scalefactor)
            @test sn.height == round(boxheight(cb) * scalefactor)
        end
    end
    println("...finished test: output in $fname")
    sn
end
function circles_then_snapshot(fno, suf; kwargs...)
    drawcircles(0.5 * fno)
    testsnapshot(fna(fno, suf); kwargs...)
end

function frame_then_snapshot(fno, suf, cb; kwargs...)
    @layer begin
        setline(2)
        sethue("red")
        box(cb, :stroke)
    end
    testsnapshot(fna(fno, suf); cb = cb, kwargs...)
end

function ripples_then_snapshot(fno, suf; kwargs...)
    ctm = cairotojuliamatrix(getmatrix())
    # User space coordinates for device top
    midtopx, midtopy, _ = ctm^-1 * [-ctm[1, 3] + currentdrawing().width / 2, -ctm[2, 3], 0.0]
    colos = Luxor.Colors.colormap("blues")
    for ulr in range(0, 1, length = 3 * 577)
        r = 2 * 577 * ulr^1.5
        # phase
        ϕ = 20 * sin(40ulr)
        sethue(colos[80 + Int(round(ϕ))])
        p = Point(midtopx, midtopy + 2ϕ)
        circle(p, r, :stroke)
        # Add some rays for showing shear transformation
        sethue("yellow")
        [circle(p + (r * cos(α), -r * sin(α)), 2, :stroke) for α in [0.25, 0.75, 1.25, 1.75]π]
    end
    testsnapshot(fna(fno, suf); kwargs...)
end

# 1-4: Full drawing snapshot from bounded recording surface
Drawing(640, 480, :rec)
origin(400, 280)
sethue("gold4")
circles_then_snapshot(1, "png")
sethue("gold3")
circles_then_snapshot(2, "pdf")
sethue("gold2")
circles_then_snapshot(3, "svg")
sethue("gold1")
circles_then_snapshot(4, "png", scalefactor = 0.5)

# 5-6: Crop snapshots from bounded recording surface
ptul = -Point(getmatrix()[5:6]...) # origin is not at center
ptbr = ptul + (currentdrawing().width, currentdrawing().height)
cb = BoundingBox(ptul, ptbr) - 50
frame_then_snapshot(5, "png", cb)

cb = BoundingBox(box(Point(10, -115), Point(42, -91), :stroke))
frame_then_snapshot(6, "png", cb)

# 7: Crop and scale up
frame_then_snapshot(7, "png", cb, scalefactor = 20)

# 8: Shrink full image
origin()
cb = BoundingBox()
frame_then_snapshot(8, "png", cb, scalefactor = 0.5)

# 9-10: Crop snapshot from free extents recording surface

# It seems that the svg device called by Cairo fails to find
# a size reference from unbounded recording surfaces. Botched up svgs
# include "-1.#QNAN", or very long strings. This occurs with and
# without 'Pro api' text on the drawing.
# However, it seems that if an svg has previously been rendered from a surface
# with bounds, the svg renderer works normally with unbounded surfaces.

Drawing(NaN, NaN, :rec)
@test snapshot() == false
origin(400, 280)
sethue("gold4")
cb = BoundingBox(Point(-128, -96), Point(128, 96))
circles_then_snapshot(9, :png; cb = cb)

# Large image, more effective storage as vector type svg
cb = BoundingBox(Point(-1280, -960), Point(1280, 960))
circles_then_snapshot(10, "svg"; cb = cb)

# 11-13: Full drawing snapshot from transformed recording surface
Drawing(640, 480, :rec)
origin()
scale(0.5)
ripples_then_snapshot(11, "png")
scale(2)
ripples_then_snapshot(12, "png")
origin()
scale(1, 0.5)
ctm = getmatrix()
ctm[3] = 0.1     # shear y
setmatrix(ctm)
rotate(π / 15)
ripples_then_snapshot(13, "png")

# 14-16: Crop and scale snapshot from transformed recording surface
pul = Point(-240, -380)
fontsize(36)
text("ul", pul, valign = :top)
pur = Point(80, -380)
text("ur", pur, valign = :top)
pbr = Point(80, -140)
text("br", pbr, valign = :top)
pbl = Point(-240, -140)
text("Bottom left", pbl, valign = :bottom)
cb = BoundingBox(pul, pbr)
frame_then_snapshot(14, "png", cb)
frame_then_snapshot(15, "png", cb, scalefactor = 2)
testsnapshot(fna(16, "png"))
finish()


text:
# This functionality ought to be moved to another package.
# Perhaps ImageMorphology.jl ?

"""
    y = complement(x)

Take the complement `1-x` of `x`.  If `x` is a color with an alpha channel,
the alpha channel is left untouched. Don't forget to add a dot when `x` is
an array: `complement.(x)`
"""
complement(x::Union{Number,Colorant}) = oneunit(x)-x
complement(x::TransparentColor) = typeof(x)(complement(color(x)), alpha(x))
@deprecate complement(x::AbstractArray) complement.(x)


text:
using ZStackAnalysis
using Test

@testset "ZStackAnalysis.jl" begin
    # Write your tests here.
end


text:
@testset "DataCollection" begin
    mutable struct Nested
        data::Vector{Float64}
    end

    function initialize()
        model = ABM(
            Agent3,
            GridSpace((10, 10));
            properties = Dict(
                :year => 0,
                :tick => 0,
                :flag => false,
                :container => Float64[],
                :deep => Nested([20.0, 52.1]),
            ),
        )
        add_agent!((4, 3), model, 0.1)
        add_agent!((7, 5), model, 0.35)
        add_agent!((2, 9), model, 0.67)
        return model
    end

    function agent_step!(agent, model)
        agent.weight += 0.05
        if model.tick % 365 == 0
            agent.weight *= 2
        end
    end
    function model_step!(model)
        model.tick += 1
        model.flag = !model.flag
        if model.tick % 365 == 0
            model.year += 1
            model.deep.data[1] += 0.5
        end
    end

    x_position(agent) = first(agent.pos)
    model = initialize()

    @testset "DataFrame init" begin
        @test init_agent_dataframe(model, nothing) == DataFrame()
        @test collect_agent_data!(DataFrame(), model, nothing, 1) == DataFrame()

        @test init_model_dataframe(model, nothing) == DataFrame()
        @test collect_model_data!(DataFrame(), model, nothing, 1) == DataFrame()

        props = [:weight]
        @test sprint(
            show,
            "text/csv",
            describe(init_agent_dataframe(model, props), :eltype),
        ) ==
              "\"variable\",\"eltype\"\n\"step\",\"Int64\"\n\"id\",\"Int64\"\n\"weight\",\"Float64\"\n"
        props = [:year]
        @test sprint(
            show,
            "text/csv",
            describe(init_model_dataframe(model, props), :eltype),
        ) == "\"variable\",\"eltype\"\n\"step\",\"Int64\"\n\"year\",\"Int64\"\n"

        @test_throws ErrorException init_agent_dataframe(model, [:UNKNOWN])
    end

    @testset "dataname" begin
        adata = [(:weight, mean), (x_position, length)]
        @test dataname(:weight) == "weight"
        @test dataname((:weight, mean)) == "mean_weight"
        @test dataname(adata[2]) == "length_x_position"
        @test dataname((x_position, length)) == "length_x_position"
        ypos(a) = a.pos[2] > 5
        @test dataname((x_position, length, ypos)) == "length_x_position_ypos"

        funcs = Vector{Function}(undef, 2)
        for i in 1:length(funcs)
            inline_func(x) = i * x
            funcs[i] = inline_func
        end
        @test dataname(funcs[2]) == "inline_func_i=2"
    end

    @testset "Aggregate Collections" begin
        props = [:weight]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        # Expecting weight values of all three agents. ID and step included.
        @test size(df) == (3, 3)
        @test propertynames(df) == [:step, :id, :weight]
        @test mean(df[!, :weight]) ≈ 0.37333333333

        props = [(:weight, mean)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        # Activate aggregation. Weight column is expected to be one value for this step,
        # renamed mean(weight). ID is meaningless and will therefore be dropped.
        @test size(df) == (1, 2)
        @test propertynames(df) == [:step, :mean_weight]
        @test df[1, dataname((:weight, mean))] ≈ 0.37333333333

        # Add a function as a property
        props = [:weight, x_position]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (3, 4)
        @test propertynames(df) == [:step, :id, :weight, :x_position]
        @test mean(df[!, :x_position]) ≈ 4.3333333

        props = [(:weight, mean), (x_position, mean)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (1, 3)
        @test propertynames(df) == [:step, :mean_weight, :mean_x_position]
        @test df[1, dataname(props[2])] ≈ 4.3333333

        xtest(agent) = agent.pos[1] > 5
        ytest(agent) = agent.pos[2] > 5

        props = [(:weight, mean, ytest)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (1, 2)
        @test propertynames(df) == [:step, :mean_weight_ytest]
        @test df[1, dataname((:weight, mean, ytest))] ≈ 0.67

        props = [(:weight, mean), (:weight, mean, ytest)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (1, 3)
        @test propertynames(df) == [:step, :mean_weight, :mean_weight_ytest]
        @test df[1, dataname(props[1])] ≈ 0.37333333333
        @test df[1, dataname(props[2])] ≈ 0.67

        props = [(:weight, mean, xtest), (:weight, mean, ytest)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props, 1)
        @test size(df) == (1, 3)
        @test propertynames(df) == [:step, :mean_weight_xtest, :mean_weight_ytest]
        @test df[1, dataname(props[1])] ≈ 0.35
        @test df[1, dataname(props[2])] ≈ 0.67
    end

    @testset "High-level API for Collections" begin
        # Extract data from the model every year for five years,
        # with the average `weight` of all agents every six months.
        each_year(model, step) = step % 365 == 0
        six_months(model, step) = step % 182 == 0
        agent_data, model_data = run!(
            model,
            agent_step!,
            model_step!,
            365 * 5;
            when_model = each_year,
            when = six_months,
            mdata = [:flag, :year],
            adata = [(:weight, mean)],
        )

        @test size(agent_data) == (11, 2)
        @test propertynames(agent_data) == [:step, :mean_weight]
        @test maximum(agent_data[!, :step]) == 1820

        @test size(model_data) == (6, 3)
        @test propertynames(model_data) == [:step, :flag, :year]
        @test maximum(model_data[!, :step]) == 1825

        agent_data, model_data = run!(
            model,
            agent_step!,
            model_step!,
            365 * 5;
            when_model = [1, 365 * 5],
            when = false,
            mdata = [:flag, :year],
            adata = [(:weight, mean)],
        )
        @test size(agent_data) == (0, 2)
        @test size(model_data) == (2, 3)

        _, model_data = run!(
            model,
            agent_step!,
            model_step!,
            365 * 5;
            when_model = [1, 365 * 5],
            when = false,
            mdata = [:deep],
            obtainer = deepcopy,
        )
        @test model_data[1, :deep].data[1] < model_data[end, :deep].data[1]

        _, model_data = run!(
            model,
            agent_step!,
            model_step!,
            365 * 5;
            when_model = [365 * 5],
            when = false,
            mdata = [(m) -> (m.deep.data[i]) for i in 1:length(model.deep.data)],
        )
        @test Array{Float64,1}(model_data[1, 2:end]) == model.deep.data
    end

    @testset "Low-level API for Collections" begin
        # Generate three separate dataframes using the low level API.
        # User controls the evolution, as well as the identity of each `step`.
        # Note in this example daily data are labelled with a daily `step`,
        # and yearly data with a yearly `step`.
        model = initialize()
        model_props = [:flag, :year]
        agent_agg = [(:weight, mean)]
        agent_props = [:weight]
        daily_model_data = init_model_dataframe(model, model_props)
        daily_agent_aggregate = init_agent_dataframe(model, agent_agg)
        yearly_agent_data = init_agent_dataframe(model, agent_props)

        for year in 1:5
            for day in 1:365
                step!(model, agent_step!, model_step!, 1)
                collect_model_data!(daily_model_data, model, model_props, day * year)
                collect_agent_data!(daily_agent_aggregate, model, agent_agg, day * year)
            end
            collect_agent_data!(yearly_agent_data, model, agent_props, year)
        end

        @test size(daily_model_data) == (1825, 3)
        @test propertynames(daily_model_data) == [:step, :flag, :year]
        @test maximum(daily_model_data[!, :step]) == 1825

        @test size(daily_agent_aggregate) == (1825, 2)
        @test propertynames(daily_agent_aggregate) == [:step, :mean_weight]
        @test maximum(daily_agent_aggregate[!, :step]) == 1825

        @test size(yearly_agent_data) == (15, 3)
        @test propertynames(yearly_agent_data) == [:step, :id, :weight]
        @test maximum(yearly_agent_data[!, :step]) == 5

        @test dummystep(model) === nothing
        @test dummystep(model[1], model) === nothing
        tick = model.tick
        step!(model, agent_step!, 1)
        @test tick == model.tick
        stop(m, s) = m.year == 6
        step!(model, agent_step!, model_step!, stop)
        @test model.tick == 365 * 6
        step!(model, agent_step!, model_step!, 1, false)
    end

    @testset "Mixed-ABM collections" begin
        model = ABM(Union{Agent3,Agent4}, GridSpace((10, 10)); warn = false)
        add_agent_pos!(Agent3(1, (6, 8), 54.65), model)
        add_agent_pos!(Agent4(2, (10, 7), 5), model)

        # Expect position type (both agents have it)
        props = [:pos]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (2, 4)
        @test propertynames(df) == [:step, :id, :agent_type, :pos]
        @test df[1, :pos] == model[1].pos
        @test df[2, :pos] == model[2].pos

        # Expect weight for Agent3, but missing for Agent4
        props = [:weight]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (2, 4)
        @test df[1, :weight] == model[1].weight
        @test ismissing(df[2, :weight])

        # Expect similar output, but using functions
        wpos(a) = a.pos[1] + a.weight
        props = [wpos]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (2, 4)
        @test df[1, :wpos] == model[1].pos[1] + model[1].weight
        @test ismissing(df[2, :wpos])

        # Expect similar output, but using anonymous accessor functions
        props = [(a) -> (a.pos[i] + a.weight) for i in 1:2]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (2, 5)
        @test df[1, dataname(props[1])] == model[1].pos[1] + model[1].weight
        @test ismissing(df[2, dataname(props[2])])

        add_agent_pos!(Agent3(3, (2, 4), 19.81), model)
        add_agent_pos!(Agent4(4, (4, 1), 3), model)

        props = [:pos, :weight, :p, wpos]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (4, 7)
        @test typeof(df.pos) <: Vector{Tuple{Int,Int}}
        @test typeof(df.weight) <: Vector{Union{Missing,Float64}}
        @test typeof(df.p) <: Vector{Union{Missing,Int}}
        @test typeof(df.wpos) <: Vector{Union{Missing,Float64}}

        # Expect something completely unknown to fail
        props = [:UNKNOWN]
        @test_throws ErrorException init_agent_dataframe(model, props)

        # Aggregates should behave in a similiar fashion
        pos1(a) = a.pos[1]
        props = [(:pos, length), (pos1, sum)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (1, 3)
        @test typeof(df.length_pos) <: Vector{Int}
        @test df[1, :length_pos] == 4
        @test typeof(df.sum_pos1) <: Vector{Int}
        @test df[1, :sum_pos1] == 22

        # Expect aggregate to filter out Agent4's
        a3(a) = a isa Agent3
        props = [(pos1, sum, a3)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (1, 2)
        @test typeof(df.sum_pos1_a3) <: Vector{Int}
        @test df[1, :sum_pos1_a3] == 8

        # When aggregating, missing data must be handled explicitly.
        props = [(:weight, sum)]
        df = init_agent_dataframe(model, props)
        @test_throws ErrorException collect_agent_data!(df, model, props)

        # Filtering out missings makes this work.
        props = [(:weight, sum, a3)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (1, 2)
        @test typeof(df.sum_weight_a3) <: Vector{Float64}
        @test df[1, :sum_weight_a3] ≈ 74.46

        # Handle missmatches
        # In this example, weight exists in both agents, but they have different types
        mutable struct Agent3Int <: AbstractAgent
            id::Int
            pos::Dims{2}
            weight::Int
        end
        model = ABM(Union{Agent3,Agent3Int}, GridSpace((10, 10)); warn = false)
        add_agent_pos!(Agent3(1, (6, 8), 54.65), model)
        add_agent_pos!(Agent3Int(2, (10, 7), 5), model)
        add_agent_pos!(Agent3(3, (2, 4), 19.81), model)
        add_agent_pos!(Agent3Int(4, (4, 1), 3), model)

        props = [:weight]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (4, 4)
        @test typeof(df.weight) <: Vector{Union{Float64,Int}}

        # Expect a1.weight <: Float64, a2.weight <: Int64 to fail in aggregate
        props = [(:weight, sum)]
        @test_throws ErrorException init_agent_dataframe(model, props)

        # Promotion is the fix in this case
        fweight(a) = Float64(a.weight)
        props = [(fweight, sum)]
        df = init_agent_dataframe(model, props)
        collect_agent_data!(df, model, props)
        @test size(df) == (1, 2)
        @test df[1, :sum_fweight] ≈ 82.46
    end

    @testset "Observers" begin
        model = initialize()
        model_props = [:container]
        model_data = init_model_dataframe(model, model_props)
        push!(model.container, 50.0)
        collect_model_data!(model_data, model, model_props, 0; obtainer = copy)
        push!(model.container, 37.2)
        collect_model_data!(model_data, model, model_props, 1; obtainer = copy)
        model.container[1] += 21.9
        collect_model_data!(model_data, model, model_props, 2; obtainer = copy)
        @test model_data.container[1][1] ≈ 50.0
        @test model_data.container[3][1] ≈ 71.9
        @test length.(model_data.container) == [1, 2, 2]

        model = initialize()
        model_props = [:deep]
        model_data = init_model_dataframe(model, model_props)
        push!(model.deep.data, 17.5)
        collect_model_data!(model_data, model, model_props, 0; obtainer = deepcopy)
        push!(model.deep.data, 1.2)
        collect_model_data!(model_data, model, model_props, 1; obtainer = deepcopy)
        model.deep.data[1] += 0.9
        collect_model_data!(model_data, model, model_props, 2; obtainer = deepcopy)
        @test model_data[1, :deep].data[1] ≈ 20.0
        @test model_data[3, :deep].data[1] ≈ 20.9
        @test [length(d.data) for d in model_data[!, :deep]] == [3, 4, 4]

        model = initialize()
        agent_data, model_data = run!(
            model,
            agent_step!,
            model_step!,
            365 * 5;
            when_model = [1, 365 * 5],
            when = false,
            mdata = [:flag, :year, :container, :deep],
            adata = [(:weight, mean)],
            obtainer = deepcopy,
        )
        @test size(agent_data) == (0, 2)
        @test size(model_data) == (2, 5)
    end
end

@testset "Parameter scan" begin
    n = 10
    parameters = Dict(:density => [0.6, 0.7, 0.8], :griddims => (20, 20))

    forest, agent_step!, forest_step! = Models.forest_fire()
    forest_initiation(; kwargs...) = Models.forest_fire(; kwargs...)[1]

    burnt(f) = count(t == 3 for t in f.trees)
    unburnt(f) = count(t == 1 for t in f.trees)
    @testset "Standard Scan" begin
        mdata = [unburnt, burnt]
        _, data = paramscan(
            parameters,
            forest_initiation;
            n = n,
            agent_step! = agent_step!,
            model_step! = forest_step!,
            mdata,
        )
        # 3 is the number of combinations of changing params
        @test size(data) == ((n + 1) * 3, 4)
        _, data = paramscan(
            parameters,
            forest_initiation;
            n = n,
            agent_step! = agent_step!,
            model_step! = forest_step!,
            include_constants = true,
            mdata,
        )
        # 3 is the number of combinations of changing params,
        # 5 is 3+2, where 2 is the number of constant parameters
        @test size(data) == ((n + 1) * 3, 5)

        mdata = [burnt]
        _, data = paramscan(
            parameters,
            forest_initiation;
            n = n,
            agent_step! = agent_step!,
            model_step! = forest_step!,
            mdata,
        )
        @test unique(data.step) == 0:10
        @test unique(data.density) == [0.6, 0.7, 0.8]
    end
end

@testset "Issue #179 fix" begin
    # only ids sorted, not properties
    model = ABM(Agent2)
    for i in 1:5
        add_agent!(model, i * 0.2)
    end
    data, _ = run!(model, dummystep, 2; adata = [:weight])
    @test data[1, :id] == 1 && data[1, :weight] ≈ 0.2
    @test data[3, :id] == 3 && data[3, :weight] ≈ 0.6
    @test data[6, :id] == 1 && data[6, :weight] ≈ 0.2
end

@testset "ensemblerun! and different seeds" begin
    _, as!, ms! = Models.daisyworld(; griddims = (4, 4), init_black = 0.5, init_white = 0.5)
    generator(seed) =
        Models.daisyworld(; griddims = (4, 4), init_black = 0.5, init_white = 0.5, seed)[1]
    seeds = [1234, 563, 211]
    daisy(a) = a isa Models.Daisy
    adata = [(:age, sum, daisy)]
    adf, _ = ensemblerun!(generator, as!, ms!, 2; adata, seeds)
    @test adf[!, :sum_age_daisy] == unique(adf[!, :sum_age_daisy])
    @test sort!(adf[:, :ensemble]) == [1, 1, 1, 2, 2, 2, 3, 3, 3]
end


text:
module Spaces

using ..PCG.Types
using ..PCG.Universes
using ..PCG.Topologies
using ..PCG.Topologies.SquareGreedTopologies
using ..PCG.Topologies.HexGreedTopologies
using ..PCG.Storages
using ..PCG.Storages.LinearStorages

export initialize

# TODO: add dispatch by storage type

# TODO: does name correct
# TODO: does topology attribute requred (SquareGreedIndex should contain all required information)
function Universes.storage_index(storage::LinearStorage, topology::SquareGreedTopology, i::SquareGreedIndex)
    return LinearStorageIndex((i.y - 1) * topology.height + i.x)
end

function Universes.storage_size(::Type{LinearStorage}, topology::SquareGreedTopology)
    return topology.height * topology.width
end

# TODO: does name correct
# TODO: does topology attribute requred (HexGreedIndex should contain all required information)
function Universes.storage_index(storage::LinearStorage, topology::HexGreedTopology, i::HexGreedIndex)
    # see hex_greed_topologies.cells_hexagon
    height = topology.radius * 2 # q
    width = topology.radius  * 4 # r

    q = i.q + topology.radius
    r = i.r + topology.radius * 2

    return LinearStorageIndex((q - 1) * height + r + 1)
end

function Universes.storage_size(::Type{LinearStorage}, topology::HexGreedTopology)
    # see hex_greed_topologies.cells_hexagon
    height = topology.radius * 2 # q
    width = topology.radius  * 4 # r
    return height * width
end

function initialize(topology::Topology, base_property, recorders::Recorders)
    storage = LinearStorage(base_property, storage_size(LinearStorage, topology))

    universe = Universe(storage, topology, recorders)

    universe.cache = AreaCache{typeof(universe),
                               index_type(topology),
                               LinearStorageIndex,
                               typeof(base_property)}()

    return universe
end

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
println(ENV["JULIA_REVISE"])

using Revise

using ReviseLoop

using Test

# These do not influence the outcome, I just tried adding them on a lark
pre_sleep = get(ENV, "pre_sleep", nothing) !== nothing
post_sleep = get(ENV, "post_sleep", nothing) !== nothing

function revise_now()
    pre_sleep && sleep(5)
    println("revising")
    Revise.revise()
    println("revised")
    post_sleep && sleep(5)
end

function more_than_revisable(x)
    revise_now()
    # n.b. new version of `more_than` doesn't get called until second revise call
    more_than(x)
end

function revisable(f)
    return function(args...; kwargs...)
        revise_now()
        # n.b. new version of `f` doesn't get called until second revise call here, either
        f(args...; kwargs...)
    end
end

atexit(() -> edit_pos_num(1))

n = more_than(1)

@test n - 1 != 2
edit_pos_num(2)

@test more_than(1) == n
revise_now()
@test more_than(1) == 3 # Issue '1': fails when run as script, running rest in repl

rmt = revisable(more_than)

@test rmt(1) == 3
edit_pos_num(3)

# Issue '2a'
@test rmt(1) == 3 # stale
@test rmt(1) == 4 # delayed update

@test more_than_revisable(1) == 4
edit_pos_num(4)

# Issue '2b'
@test more_than_revisable(1) == 4 # stale
@test more_than_revisable(1) == 5 # delayed update


text:
abstract type Langmodel end
abstract type gammamodel <: Langmodel end #BaseNgram with Add-one smoothing algo
abstract type InterpolatedLanguageModel <: Langmodel end #Interpolated language model with smoothing

#DataType MLE
#Type for providing MLE ngram model scores.
#Implementation of Base Ngram Model.

struct MLE <: Langmodel
    vocab::Vocabulary
end

"""
    MLE(word::Vector{T}, unk_cutoff=1, unk_label="<unk>") where {T <: AbstractString}

Initiate Type for providing MLE ngram model scores.

Implementation of Base Ngram Model.

"""
function MLE(word::Vector{T}, unk_cutoff=1, unk_label="<unk>") where {T <: AbstractString}
    MLE(Vocabulary(word, unk_cutoff, unk_label))
end

function (lm::MLE)(text::Vector{T}, min::Integer, max::Integer) where {T <: AbstractString}
    text = lookup(lm.vocab, text)
    text=convert(Array{String}, text)
    return counter2(text, min, max)
end

struct Lidstone <: gammamodel
    vocab::Vocabulary
    gamma::Float64
end

"""
    Lidstone(word::Vector{T}, gamma:: Float64, unk_cutoff=1, unk_label="<unk>") where {T <: AbstractString}

Function to initiate Type(Lidstone) for providing Lidstone-smoothed scores.

In addition to initialization arguments from BaseNgramModel also requires
a number by which to increase the counts, gamma.
"""
function Lidstone(word::Vector{T}, gamma = 1.0, unk_cutoff=1, unk_label="<unk>") where {T <: AbstractString}
    Lidstone(Vocabulary(word, unk_cutoff, unk_label), gamma)
end

function (lm::Lidstone)(text::Vector{T}, min::Integer, max::Integer) where {T <: AbstractString}
    text = lookup(lm.vocab, text)
    text=convert(Array{String}, text)
    return counter2(text, min, max)
end

"""
    Laplace(word::Vector{T}, unk_cutoff=1, unk_label="<unk>") where {T <: AbstractString}
Function to initiate Type(Laplace) for providing Laplace-smoothed scores.

In addition to initialization arguments from BaseNgramModel also requires
a number by which to increase the counts, gamma = 1.
"""
struct Laplace <: gammamodel
    vocab::Vocabulary
    gamma::Float64
end

function Laplace(word::Vector{T}, unk_cutoff=1, unk_label="<unk>") where {T <: AbstractString}
    Laplace(Vocabulary(word, unk_cutoff, unk_label), 1.0)
end

function (lm::Laplace)(text, min::Integer, max::Integer)
    text = lookup(lm.vocab, text)
    text = convert(Array{String}, text)
    return counter2(text, min, max)
end

"""
	score(m::gammamodel, temp_lm::DefaultDict, word::AbstractString, context::AbstractString)

score is used to output probablity of word given that context

Add-one smoothing to Lidstone or Laplace(gammamodel) models

"""
function score(m::gammamodel, temp_lm::DefaultDict, word, context) #score for gammamodel output probabl
    accum = temp_lm[context]
    #print(accum)
    s = float(sum(accum)+(m.gamma)*length(m.vocab.vocab))
    for (text, count) in accum
        if text == word
            return(float(count+m.gamma)/s)
        end
    end
    return(float(m.gamma)/s)
end

"""
To get probability of word given that context

In otherwords, for given context calculate frequency distribution of word

"""
function prob(m::Langmodel, templ_lm::DefaultDict, word, context=nothing)
    if context == nothing || context == ""
        return(1/float(length(templ_lm))) #provide distribution
    else
        accum = templ_lm[context]
    end
    s = float(sum(accum))
    for (text, count) in accum
        if text == word
            return(float(count) / s)
        end
    end
    if context in keys(m.vocab.vocab)
        return(0)
    end
    return(Inf)
end

"""
	score(m::MLE, temp_lm::DefaultDict, word::AbstractString, context::AbstractString)

score is used to output probablity of word given that context in MLE

"""
function score(m::MLE, temp_lm::DefaultDict, word, context=nothing)
    prob(m, temp_lm, word, context)
end

struct WittenBellInterpolated <: InterpolatedLanguageModel
    vocab ::Vocabulary
end

"""
    WittenBellInterpolated(word::Vector{T}, unk_cutoff=1, unk_label="<unk>") where { T <: AbstractString}

Initiate Type for providing Interpolated version of Witten-Bell smoothing.

The idea to abstract this comes from Chen & Goodman 1995.

"""
function WittenBellInterpolated(word::Vector{T}, unk_cutoff=1, unk_label="<unk>") where {T <: AbstractString}
    WittenBellInterpolated(Vocabulary(word, unk_cutoff, unk_label))
end

function (lm::WittenBellInterpolated)(text::Vector{T}, min::Integer, max::Integer) where {T <: AbstractString}
    text = lookup(lm.vocab, text)
    text=convert(Array{String}, text)
    return counter2(text, min, max)
end
# alpha_gamma function for KneserNeyInterpolated
function alpha_gammma(m::WittenBellInterpolated, templ_lm::DefaultDict, word, context)
    local alpha
    local gam
    accum = templ_lm[context]
    s = float(sum(accum))
    for (text,count) in accum
        if text == word
            alpha=(float(count) / s)
            break
        else
            alpha = 1/s
        end
    end

    gam = gamma(accum)
    return alpha*(1- gam), gam
end

function count_non_zero_vals(accum::Accumulator{})
    return(length(accum))
end

function gamma(accum)
    nplus=count_non_zero_vals(accum)
    return(nplus/(nplus+float(sum(accum))))
end

"""
	score(m::InterpolatedLanguageModel, temp_lm::DefaultDict, word::AbstractString, context::AbstractString)

score is used to output probablity of word given that context in InterpolatedLanguageModel

Apply Kneserney and WittenBell smoothing
depending upon the sub-Type

"""
function score(m::InterpolatedLanguageModel, temp_lm::DefaultDict, word, context=nothing)
    if context == nothing || context == ""
        return prob(m, temp_lm, word, context)
    end
    if context in keys(temp_lm)
        alpha,gamma = alpha_gammma(m, temp_lm, word, context)
        return (alpha + gamma*score(m, temp_lm, word, context_reduce(context)))
    else
        return score(m, temp_lm, word, context_reduce(context))
    end
end

function context_reduce(context)
    context = split(context)
    join(context[2:end], " ")
end

struct KneserNeyInterpolated <: InterpolatedLanguageModel
    vocab::Vocabulary
    discount::Float64
end

"""
    KneserNeyInterpolated(word::Vector{T}, discount:: Float64,unk_cutoff=1, unk_label="<unk>") where {T <: AbstractString}

Initiate Type for providing KneserNey Interpolated language model.

The idea to abstract this comes from Chen & Goodman 1995.

"""
function KneserNeyInterpolated(word::Vector{T}, disc = 0.1, unk_cutoff=1, unk_label="<unk>") where {T <: AbstractString}
    KneserNeyInterpolated(Vocabulary(word, unk_cutoff, unk_label) ,disc)
end

function (lm::KneserNeyInterpolated)(text::Vector{T}, min::Integer, max::Integer) where {T <: AbstractString}
    text = lookup(lm.vocab, text)
    text=convert(Array{String}, text)
    return counter2(text, min, max)
end
# alpha_gamma function for KneserNeyInterpolated
function alpha_gammma(m::KneserNeyInterpolated, templ_lm::DefaultDict, word, context)
    local alpha
    local gamma
    accum = templ_lm[context]
    s = float(sum(accum))
    for (text, count) in accum
        if text == word
            alpha=(max(float(count)-m.discount, 0.0) / s)
            break
        else
            alpha = 1/length(m.vocab.vocab)
        end
    end
    gamma = (m.discount * count_non_zero_vals(accum) /s)
    return alpha, gamma
end


text:
using Revise
using ZonalFlow

dn      = "validate/stochastic/linear/"

tspan = (0.0,1000.0);
tsargs  = ( dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=100000,
            save_everystep=false,
            save_noise=false,
            saveat=5
           );

eqs     = [GQL(Λ=0),CE2(),CE2()];

"""
   Single mode excitation
"""

coeffs  = Coefficients(Ω=5.0,θ=0.0,μ=0.01,ν=0.0,ν₄=0.0,linear=true);
forcing = Stochastic(kf=3,dk=0,ε=0.01); #note dk = 0 => only kf is excited

labels  = ["ql_6x6_m3","ce2_6x6_m3","ce2_6x6_qlic_m3"];

domain  = Domain(extent=(2π,2π),res=(6,6));
prob    = BetaPlane(domain,coeffs,forcing);

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

"""
   Three-mode excitations
"""

# 6x6

domain  = Domain(extent=(2π,2π),res=(6,6));
prob    = BetaPlane(domain,coeffs,forcing);

coeffs  = Coefficients(Ω=5.0,θ=0.0,μ=0.01,ν=0.0,ν₄=0.0,linear=true);
forcing = Stochastic(kf=3,dk=1,ε=0.01);

ob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

# 6x6 with dt = 0.005
domain  = Domain(extent=(2π,2π),res=(6,6));
prob    = BetaPlane(domain,coeffs,forcing);

tsargs  = ( dt=0.0005,
            adaptive=false,
            progress=true,
            progress_steps=100000,
            save_everystep=false,
            saveat=5
           );

labels  = ["ql_6x6_dt0005","ce2_6x6_dt0005","ce2_6x6_qlic_dt0005"];

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

# 8x8
domain  = Domain(extent=(2π,2π),res=(8,8));
prob    = BetaPlane(domain,coeffs,forcing);

tsargs  = ( dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=100000,
            save_everystep=false,
            saveat=5
           );

labels  = ["ql_8x8","ce2_8x8","ce2_8x8_qlic"];

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

# 10x10

domain  = Domain(extent=(2π,2π),res=(10,10));
prob    = BetaPlane(domain,coeffs,forcing);

labels  = ["ql_10x10","ce2_10x10","ce2_10x10_qlic"];

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

# 12x12

domain  = Domain(extent=(2π,2π),res=(12,12));
prob    = BetaPlane(domain,coeffs,forcing);

labels  = ["ql_12x12","ce2_12x12","ce2_12x12_qlic"];

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

# GCE2
# eqs     = [NL(),GQL(Λ=prob.d.nx-1),GCE2(Λ=prob.d.nx-1)];
# labels  = ["nl_linear","gql_5_linear","gce2_5_linear"];
# sols    = integrate(prob,eqs,tspan;tsargs...);
# write(prob,eqs,sols,dn="validate/stochastic/6x6/",labels=labels)


text:
############################################################################################
## This file is part of RLinearAlgebra.jl
##
## Overview: abstractions and methods for performing various deploying Random Linear System
## solvers
##
## Contents
## - `RLSSolver` struct
## - `rsolve` functions
## - Randomized linear Solvers
## - Export Statements
##
############################################################################################

# Dependencies:

#############################################
# `RLSSolver` stucture
#############################################

"""
    RLSSolver{S<:LinSysSampler, R<:LinSysSolveRoutine,
        L<:LinSysSolverLog, C<:LinSysStopCriterion}

An encapsulation for a (R)andomized (L)inear (S)ysterm Solver as a mutable structure.

# Fields
- `sampler::S`, a procedure for sampling, sketching or selecting with respect to a linear
    system.
- `routine::R`, a routine for randomly generating the subproblem generated by the sampling,
    sketching or selecting operation
- `log::L`, a logger for the progress and behavior of the randomized linear system solver.
- `stop::C`, a stopping criterion for the randomized linear system solver.
- `x::Union{Vector{Float64},Nothing}`, the result produced by the randomized linear system
    solver.

# Constructions
- `RLSSolver(iter::Int64)` specifies a solver with: random cyclic row sampling;  partially
    orthogonalized row projections (memory of five) subproblem solver; a full residual
    logger; and a maximum iteration stopping criterion specified by `iter`.
"""
mutable struct RLSSolver{S<:LinSysSampler, R<:LinSysSolveRoutine,
    L<:LinSysSolverLog, C<:LinSysStopCriterion}

    sampler::S
    routine::R
    log::L
    stop::C
    x::Union{Vector{Float64},Nothing}
end

RLSSolver(iter::Int64) = RLSSolver(
    LinSysVecRowRandCyclic(),   # Random Cyclic Sampling
    LinSysVecRowProjPO(),       # Partially Orthogonalized Row Projection, 5 vector memory
    LSLogFull(),                # Full Logger: maintains residual history
    LSStopMaxIterations(iter),  # Maximum iterations stopping criterion
    nothing                     # System solution
)

#############################################
# `rsolve` functions
#############################################
"""
    rsolve([solver::RLSSolver,] A, b)

Deploys an iterative randomized solver for the linear system whose coefficient matrix is
    encapsulated in `A` and whose constant vector is encapsulated in `b`. If the `solver`
    is not specified, then the default `RLSSolver` is constructed and used with a maximum
    number of iterations of `10 * length(b)`. See `RLSSolver` for more details on the
    default solver.

The function returns the solution.
"""
function rsolve(A, b)
    solver = RLSSolver(10 * length(b))
    x = zeros(size(A, 2))
    rsolve!(solver, A, b, x)

    return x
end

function rsolve(solver::RLSSolver, A, b)
    x = zeros(size(A, 2))
    rsolve!(solver, A, b, x)

    return x
end

"""
    rsolve!([solver::RLSSolver,] A, b, x::AbstractVector)

Identical to `rsolve` except an initial iterate is supplied by the argument `x`. The
    function overwrites the solution into the vector `x`.
"""
function rsolve!(A, b, x::AbstractVector)
    solver = RLSSolver(10 * length(b))
    rsolve!(solver, A, b, x)

    return nothing
end

function rsolve!(solver::RLSSolver, A, b, x::AbstractVector)

    iter = 0
    log_update!(solver.log, solver.sampler, x, (), iter, A, b)

    while check_stop_criterion(solver.log, solver.stop) == false
        iter += 1
        samp = sample(solver.sampler, A, b, x, iter)
        rsubsolve!(solver.routine, x, samp, iter)
        log_update!(solver.log, solver.sampler, x, samp, iter, A, b)
    end

    solver.x = x

    return nothing
end

#############################################
# Randomized Linear Solvers
#############################################

#############################################
# Export Statements
#############################################
# See RLinearAlgebra.jl
#export RLSSolver, rsolve, rsolve!