text:
using Schistoxpkg
using Test
using Distributions
using Random
using JLD

# filename to save population in =#
# filename = "equ_runs_1.jld"
N = 1000
max_age = 100
initial_worms = 10
time_step = 10
worm_stages = 1
female_factor = 1
male_factor = 1
contact_rate = 0.03
ages_per_index = 5

# if more than one community, then specify how many here
#N_communities = 3
N_communities = 1

# next parameter is the relative probabilities of being in each community
# if entries are all equal, then all communities are equally likely and will
# be roughly the same size
#community_probs = [1,1,1]
community_probs = 1

# community_contact_rate = [1,1,1]
community_contact_rate = 1

# parameter for proportion of people who are given mda who will take it
mda_adherence = .9
mda_access = .9

heavy_burden_threshold = 16
# number of days after which miracidia become cercariae
miracidia_maturity_time = 24 # for S. mansoni
# miracidia_maturity_time = 21 # for S. haemotobium

env_cercariae = 0
initial_miracidia = 50000*N/1000
init_env_cercariae = 50000*N/1000
initial_miracidia_days = trunc(Int,ceil(miracidia_maturity_time/time_step, digits = 0))

# how long to run simulation for
number_years_equ = 200

max_fecundity = 0.34  # for S. mansoni [Toor et al JID paper SI]
#max_fecundity = 0.3  # for S. haematobium [Toor et al JID paper SI]

density_dependent_fecundity = 0.0007 # for S. mansoni [Toor et al JID paper SI]
#density_dependent_fecundity = 0.0006 # for S. haematobium [Toor et al JID paper SI]

r = 0.03 # aggregation parameter for negative binomial for egg production
num_time_steps_equ = trunc(Int, 365*number_years_equ / time_step)

# human birth rate
birth_rate = 28*time_step/(1000*365)

average_worm_lifespan = 5.7 # years for S. mansoni [Toor et al JID paper SI]
#average_worm_lifespan = 4 # years for S. haematobium [Toor et al JID paper SI]

# this is the aggregation parameter for the predisposition
predis_aggregation = 0.24
predis_weight = 1

# what proportion of miracidias and cercariae survive each round
env_miracidia_survival_prop = 1/2
env_cercariae_survival_prop = 1/2
mda_coverage = 0.8 # proportion of target age group reached by mda
mda_round = 0

# proportion of cercariae which can infect humans
human_cercariae_prop = 1

# gamma distribution for Kato-Katz method
gamma_k = Gamma(0.87,1/0.87)

vaccine_effectiveness = 0.95
num_sims = 1

# record the state of the population this often in years
record_frequency = 1/24

#= this is the number of thousands of people in 5 year (0-4, 5-9,...) intervals in Kenya
and will be used to give a specified age structure when we run to equilibrium =#
spec_ages = 7639, 7082, 6524, 5674, 4725, 4147, 3928, 3362,
            2636, 1970, 1468, 1166, 943, 718, 455, 244

#= number of deaths per 1000 individuals by age
    first entry is for under 1's, then for 5 year intervals from then on =#

death_prob_by_age = [0.0656, 0.0093, 0.003, 0.0023, 0.0027, 0.0038, 0.0044, 0.0048, 0.0053,
                                                0.0065, 0.0088, 0.0106, 0.0144, 0.021, 0.0333, 0.0529, 0.0851, 0.1366, 0.2183, 0.2998 , 0.3698, 1]

ages_for_deaths = [1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60,
                                              65, 70, 75, 80, 85, 90, 95, 100, 110]

use_kato_katz = 0
kato_katz_par = 0.87
drug_efficacy = 0.863
scenario = "moderate "
@testset "update_contact_rate" begin
    @test isapprox(update_contact_rate([5,10,3], [0,0,0], [1,2,3,4,5,6,7,8,9,10,11,12]), [6,11,4])
end

@testset "miracidia_death" begin
    @test miracidia_death([1000,1421,43523,55],1/2) == [1000,1421,43523,28]
end

@testset "miracidia_death" begin
    @test miracidia_death([1000,1421,43523,55],1/3) == [1000,1421,43523,18]
end

@testset "cercariae_death" begin
    @test cercariae_death(1000,1/2,1) == 500
end

@testset "cercariae_death" begin
    @test cercariae_death(1000,1/3,1) == 333
end

@testset "worm_pairs" begin
    @test isapprox(calculate_worm_pairs([[1,3],[17,0],[400,12312]],[[100,3000],[1,6],[1,1]]),[4,7,2] )
end

@testset "total_worms" begin
    @test isapprox(calculate_total_worms([[1,3],[17,0],[400,1]],[[10,30],[1,6],[1,1]])[1], [4,17,401])
end

@testset "total_worms" begin
    @test isapprox(calculate_total_worms([[1,3],[17,0],[400,1]],[[10,30],[1,6],[1,1]])[2], [40,7,2])
end

time_step = 10

death_prob_by_age = [0.0656, 0.0093, 0.003, 0.0023, 0.0027, 0.0038, 0.0044, 0.0048, 0.0053,
                                                0.0065, 0.0088, 0.0106, 0.0144, 0.021, 0.0333, 0.0529, 0.0851, 0.1366, 0.2183, 0.2998 , 0.3698, 1]

ages_for_deaths = [1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60,
                                              65, 70, 75, 80, 85, 90, 95, 100, 110]

@testset "gen_ages_and_deaths" begin
    @test generate_ages_and_deaths(1, [1.01,2.3,3.1,4.2,5.3], [2.1,1.2,4.3,2.3,6.0], death_prob_by_age, ages_for_deaths, time_step)[1][1] == 1.01 + time_step/365
end

@testset "gen_ages_and_deaths" begin
    @test generate_ages_and_deaths(1, [1.01,2.3,3.1,4.2,5.3], [2.1,1.2,4.3,2.3,6.0], death_prob_by_age, ages_for_deaths, time_step)[1][end-1] == time_step/365
end

@test_throws ErrorException("max_age must be greater than 60") make_age_contact_rate_array(10, "high adult", [], [])

@testset "make_age_contact_rate_array(max_age,scenario)" begin
    @test contact_rates_by_age = make_age_contact_rate_array(max_age, scenario, [1,2,3], [3,5,6])[end] == 6
end

@testset "make_age_contact_rate_array(max_age,scenario)" begin
    @test make_age_contact_rate_array(100, "high adult", [], [])[1] == 0.01
end

@testset "make_age_contact_rate_array(max_age, scenario)" begin
    @test make_age_contact_rate_array(100, "high adult", [], [])[100] == 0.12
end

@testset "make_age_contact_rate_array(max_age, scenario)" begin
    @test make_age_contact_rate_array(100, "low adult", [], [])[100] == 0.02
end

@testset "make_age_contact_rate_array(max_age, scenario)" begin
    @test make_age_contact_rate_array(100, "moderate adult", [], [])[100] == 0.06
end
time_step = 10
human_cercariae = [[4],[5],[6],[7],[8]]
env_miracidia = [10, 10, 10, 10, 1]
env_cercariae = 100
contact_rate = 0.01
predisposition  = [0.1, 0.21, 1.4, 3.1, 100]
age_contact_rate = [0.0, 0.06, 0.12, 0.02, 0.06]
vac_status = [0,0,0,0,0]
vaccine_effectiveness = 0.95
community = 1,2,3,4,1
community_contact_rate = 1,1,1,0,1
female_worms= [[1,1], [2,1], [2,4],[3,1], [2,4]]
male_worms= [[1,1], [2,1], [2,5],[2,1], [2,4]]

@testset "cercariae_uptake(max_age, scenario)" begin
    @test isapprox(cercariae_uptake(copy( env_miracidia), copy(env_cercariae), copy(time_step), copy(contact_rate),
    community, community_contact_rate, copy(female_worms), copy(male_worms),
        copy(predisposition), copy(age_contact_rate), copy(vac_status), copy(vaccine_effectiveness), 1, 24)[2],  [10, 10, 10, 1])
end

time_step = 10
human_cercariae = [[4],[5],[6],[7],[8]]
env_miracidia = [10, 10, 10, 10, 1]
env_cercariae = 100
contact_rate = 0.01
predisposition  = [0.1, 0.21, 1.4, 3.1, 100]
age_contact_rate = [0.0, 0.06, 0.12, 0.02, 0.06]
vac_status = [0,0,0,0,0]
vaccine_effectiveness = 0.95

@testset "cercariae_uptake" begin
    @test isapprox(cercariae_uptake(copy( env_miracidia), copy(env_cercariae), copy(time_step), copy(contact_rate),
    community, community_contact_rate, copy(female_worms), copy(male_worms),
        copy(predisposition), copy(age_contact_rate), copy(vac_status), copy(vaccine_effectiveness),1,24)[3][1][2],  1)
end

time_step = 10
human_cercariae = [[4],[5],[6],[7],[8]]
env_miracidia = [10, 10, 10, 10, 1]
env_cercariae = 100
contact_rate = 0.01
predisposition  = [0.1, 0.21, 1.4, 3.1, 100]
age_contact_rate = [0.0, 0.06, 0.12, 0.02, 0.06]
vac_status = [0,0,0,0,0]
vaccine_effectiveness = 0.95

@testset "cercariae_uptake" begin
    @test cercariae_uptake(copy( env_miracidia), copy(env_cercariae), copy(time_step), copy(contact_rate),
    community, community_contact_rate, copy(female_worms), copy(male_worms),
        copy(predisposition), copy(age_contact_rate), copy(vac_status), copy(vaccine_effectiveness),1,24)[3][5][2] >0
end

time_step = 10
human_cercariae = [[4],[5],[6],[7],[8]]
env_miracidia = [10, 10, 10, 10, 1]
env_cercariae = 100
contact_rate = 0.01
predisposition  = [0.1, 0.21, 1.4, 3.1, 100]
age_contact_rate = [0.0, 0.06, 0.12, 0.02, 0.06]
vac_status = [0,0,0,0,1]
vaccine_effectiveness = 1

@testset "cercariae_uptake" begin
    @test cercariae_uptake(copy( env_miracidia), copy(env_cercariae), copy(time_step), copy(contact_rate),
    community, community_contact_rate, copy(female_worms), copy(male_worms),
        copy(predisposition), copy(age_contact_rate), copy(vac_status), copy(vaccine_effectiveness),1,1)[3][5][2] == 4
end

female_worms = [[10000,200000],[200000,0]]
male_worms = [[10000,200000],[200000,0]]
worm_stages = 2
average_worm_lifespan = 5.7

@testset "worm_maturity" begin
    @test worm_maturity(female_worms, male_worms, worm_stages, average_worm_lifespan, time_step)[1][1][2] < 200000
end

@testset "worm_maturity" begin
    @test worm_maturity(female_worms, male_worms, worm_stages, average_worm_lifespan, time_step)[1][2][2] > 00
end

@testset "worm_maturity" begin
    @test worm_maturity(female_worms, male_worms, worm_stages, average_worm_lifespan, time_step)[2][1][1] < 10000
end

@testset "worm_maturity" begin
    @test worm_maturity(female_worms, male_worms, worm_stages, average_worm_lifespan, time_step)[2][2][1] < 200000
end

@testset "miracidia_production" begin
    @test miracidia_production([1,2,3],[2,4,1], 10, [1,1,1],[1,1,2], [1,2,3])[1] == 2
end

@testset "miracidia_production" begin
    @test miracidia_production([1,2,3],[2,4,1], 10, [1,1,1],[1,1,1.2], [1,2,3])[end] == 6
end

@testset "miracidia_production" begin
    @test miracidia_production([1,2,3],[2,4,1], 10, [1,1,1],[0,1,0], [1,2,3])[end] == 2
end

@testset "miracidia_production" begin
    @test miracidia_production([1,2,3],[2,4,1], 10, [1,1,1],[0.5,0.5,0.5], [1,2,3])[end] == 6
end

vaccine_info = []
push!(vaccine_info, vaccine_information(0.75, 4, 16, [0,1], 3, 13))
push!(vaccine_info, vaccine_information(0.4, 17, 110, [0,1], 3, 17))

@testset "update_vaccine" begin
    @test update_vaccine(vaccine_info, 1)[1] == 0.4
end

@testset "update_vaccine" begin
    @test update_vaccine(vaccine_info, 1)[2] == 17
end

@testset "update_vaccine" begin
    @test update_vaccine(vaccine_info, 1)[3] ==110
end

@testset "update_vaccine" begin
    @test update_vaccine(vaccine_info, 1)[4] ==17
end

@testset "update_vaccine" begin
    @test update_vaccine(vaccine_info, 1)[5] ==[0,1]
end

@testset "update_vaccine" begin
    @test update_vaccine(vaccine_info, 100)[4] == Inf
end

@testset "death_of_human" begin
    @test death_of_human([2,4], [0,6], [1,0], [0.0002,0.00005], [[2,3,4],[6,3,4]], [15,7],
                                [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                                [0,0], [1,1], 1,
                                [1,1], [1,1])[1] == [4]

end

@testset "death_of_human" begin
    @test death_of_human([2,4], [0,6], [1,0], [0.0002,0.00005], [[2,3,4],[6,3,4]], [15,7],
                                [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                                [0,0], [1,1], 1,
                                [1,1], [1,1])[4] == [0.00005]

end

@testset "death_of_human" begin
    @test death_of_human([120,120], [0,0], [0.4,0.6], [[2,3,4],[6,3,4]], [15,7],
                                [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                                [0,0], [1,1], [0.00001,0.0001], 365,
                                [1,1], [1,1])[1] == []

end
age_death_rate_per_1000 = [6.56, 0.93, 0.3, 0.23, 0.27, 0.38, 0.44, 0.48,0.53, 0.65,
                           0.88, 1.06, 1.44, 2.1, 3.33, 5.29, 8.51, 13.66,
                           21.83, 29.98, 36.98]

contact_rates_by_age = make_age_contact_rate_array(100,"high adult", [], [])

#
# ages, death_ages, gender, predisposition, human_cercariae, eggs, vac_status,
#                         treated, female_worms, male_worms,vaccinated, age_contact_rate,
#                         female_factor, male_factor, contact_rates_by_age,
#                         worm_stages, predis_aggregation, predis_weight,
#                         adherence, death_prob_by_age, ages_for_deaths,
#                         mda_adherence, access, mda_access

death_prob_by_age = [0.0656, 0.0093, 0.003, 0.0023, 0.0027, 0.0038, 0.0044, 0.0048, 0.0053,
                          0.0065, 0.0088, 0.0106, 0.0144, 0.021, 0.0333, 0.0529, 0.0851, 0.1366, 0.2183, 0.2998 , 0.3698, 1]

ages_for_deaths = [1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60,
                        65, 70, 75, 80, 85, 90, 95, 100, 110]

    @testset "birth_of_human" begin
        @test birth_of_human([2,4], [0.44,0], [0,0], [1.1,1], [2,4],[[2,3,4],[6,3,4]], [15,7],
                                [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                                [0,0], [1.1,1], 1, 1, [0.0002,0.00005], 2,1, 1, [1,1],
                            death_prob_by_age, ages_for_deaths,[1,2,3,4], 0.8,[1,1], 0.9)[1] == [2,4,0]
    end
    new_pop = birth_of_human([2,4], [0.44,0], [0,0], [1.1,1],[2,4], [[2,3,4],[6,3,4]], [15,7],
                            [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                            [0,0], [1.1,1], 1, 1, [0.0002,0.00005], 2,1, 1, [1,1],
                        death_prob_by_age, ages_for_deaths,[1,2,3,4,5], 0.8,[1,1], 0.9)
    @testset "birth_of_human" begin
        @test new_pop[1]==[2,4,0]
    end
    @testset "birth_of_human" begin
        @test new_pop[6]==[[2,3,4],[6,3,4],[]]
    end

    @testset "birth_of_human" begin
        @test new_pop[7]==[15,7,0]
    end

@testset "administer_drug" begin
    @test administer_drug([[3,2],[3],[4,10]], [1,2], 1, [1,1,1]) == [[0,0],[0],[4,10]]
end

@testset "administer_drug" begin
    @test administer_drug([[3,2],[3],[4,10]], [1,2], 0, [1,1,1]) == [[3,2],[3],[4,10]]
end

Random.seed!(2525)
big_drug = administer_drug([[3, 2],[1e5],[4,10]], [2], 0.5, [1,1,1])
# I ran it with this seed so know the exact answer, so consider this a regression test.
@testset "administer_drug" begin
    @test big_drug == [[3,2],[50029],[4,10]]
end

@testset "create_mda" begin
     mda_info = create_mda(0, .75, 1, 1, 5, 2, [0,1], [0,1], [0,1], .92)
    @test mda_info[1].coverage == 0
end

out1 = out[]

@testset "create_output" begin
    push!(out1, out(1,2,3,4,5,6,7,8,9,10,11))
    @test out1[1].adult_burden == 3
end

@testset "create_output" begin
    push!(out1, out(1,2,35,4,5,6,7,8,9,10,11))
    @test out1[end].adult_burden == 35
end

@testset "create_mda" begin
     mda_info = create_mda(0, .75, 1, 1, 5, 1, [0,1], [0,1], [0,1], .92)
    @test mda_info[end].time == 5
end

mda_info = create_mda(0, .75, 1, 1, 5, 2, [0,1], [0,1], [0,1], .92)

@testset "update_mda" begin
    @test update_mda(mda_info, 3)[1] == 0
end

@testset "update_mda" begin
    @test update_mda(mda_info, 3)[5] == 3
end

@testset "update_mda" begin
    @test update_mda(mda_info, 30)[5] == Inf
end

  # But also with this big a sample size I think its fairly reasonable to use the 0.0001 and 0.9999 quantiles.
# In R qbinom(0.5, 1e5, c(0.0001, 0.9999))
# So this isn't relying on the seed, and should nearly always be true.
# but if we're repeatedly running it, we still want to use a seed.
@testset "administer_drug" begin
    @test (big_drug[2][1] > 10.0) & (big_drug[2][1] < 99990.0)
end

@testset "mda" begin
    @test isapprox(mda(1, 3, 8, 1,[0,1],
[2,5,4,7,9], [[1],[2],[3],[4],[5]], [[1],[2],[3],[4],[5]],
[[1],[2],[3],[4],[5]], [[1],[2],[3],[4],[5]],0,0,[1,1,1,1,1], [0,0,0,0,0],[0,0,0,0,0])[4],
[[1], [2], [3], [4], [5]])
end

@testset "mda" begin
    @test isapprox(mda(1, 3, 8, 1,[0,1],
[2,5,4,7,9], [[1],[2],[3],[4],[5]], [[1],[2],[3],[4],[5]],
[[1],[2],[3],[4],[5]], [[1],[2],[3],[4],[5]],0,0,[1,1,1,1,1], [1,1,1,1,1],[1,1,1,1,1])[1],
[[1], [0], [0], [0], [5]])
end

@testset "mda" begin
    @test isapprox(mda(1, 3, 8, 1,[0,1],
[2,5,4,7,9], [[1],[2],[3],[4],[5]], [[1],[2],[3],[4],[5]],
[[1],[2],[3],[4],[5]], [[1],[2],[3],[4],[5]],0,0,[1,1,1,1,1], [1,1,1,1,1],[1,1,1,1,1])[2],
[[1], [0], [0], [0], [5]])
end

#=
Run some tests on the population making function
Need to do some setup first though.
=#
 mda_access = 1
N = 1000
max_age = 100
initial_worms = 10
time_step = 10
worm_stages = 2
female_factor = 1
male_factor = 1
initial_miracidia = 1
initial_miracidia_days = trunc(Int,round(41/time_step, digits = 0))
env_cercariae = 0
#const contact_rate = 0.000005
age_death_rate_per_1000 = [6.56, 0.93, 0.3, 0.23, 0.27, 0.38, 0.44, 0.48,0.53, 0.65,
                           0.88, 1.06, 1.44, 2.1, 3.33, 5.29, 8.51, 13.66,
                           21.83, 29.98, 36.98]
predis_aggregation = 0.24
mda_adherence = 0.8
scenario = "high adult"
contact_rates_by_age = make_age_contact_rate_array(max_age, scenario, [], [])

predis_weight = 1

N_communities= 4
community_probs = [1,2,1,3]

@testset "create_pop" begin
@test_throws ErrorException("must provide probabilities for membership of each community") create_population(N, max_age, 100, community_probs, initial_worms, contact_rates_by_age,
worm_stages, female_factor, male_factor,
    initial_miracidia, initial_miracidia_days, predis_aggregation, predis_weight, time_step,
    mda_adherence, mda_access)
end

pop = create_population(N, max_age, N_communities, community_probs, initial_worms, contact_rates_by_age,
worm_stages, female_factor, male_factor,
    initial_miracidia, initial_miracidia_days, predis_aggregation, predis_weight, time_step,
    mda_adherence, mda_access)

pop = create_population(N, max_age, N_communities, community_probs, initial_worms, contact_rates_by_age,
worm_stages, female_factor, male_factor,
    initial_miracidia, initial_miracidia_days, predis_aggregation, predis_weight, time_step,
    mda_adherence, mda_access)

# These should all get give the initial value
@testset "miracidia" begin
    @test all(pop[13] .== initial_miracidia)
end

# Worms should be 0:ininity
# Get first column of worms
mworm1 = [pop[9][i][1] for i=1:length(pop[9])]
fworm1 = [pop[10][i][1] for i=1:length(pop[10])]
@testset "wormsm" begin
    @test all(mworm1 .>= 0)
end
@testset "wormsf" begin
    @test all(fworm1 .>= 0)
end

# All vectors except the last should be length N?
lens = [length(pop[i]) for i=1:(length(pop) - 3)]
push!(lens, length(pop[end]))
push!(lens, length(pop[end-1]))
@testset "N" begin
    @test all(lens .== N)
end

female_worms= [[1,1], [2,1], [2,4]]
male_worms= [[1,1], [2,1], [2,5]]
vaccine_effectiveness = 1
vaccine_coverage = 1
min_age_vaccine = 5
max_age_vaccine = 16
vaccine_gender = [0,1]
ages = [4,10,9]
death_ages = [5,11,66]
human_cercariae = [[1,2],[9,2,4], [ 0,5]]
eggs = [1,2,3]
treated = [1,1,1]
vaccine_duration = 10
vac_status = [0,0,0]
vaccine_round = 1
gender = [1,0,1]
adherence = [1,1,1]
access = [1,1,1]

@testset "vaccinate" begin
@test isapprox(vaccinate(vaccine_coverage, min_age_vaccine, max_age_vaccine, vaccine_effectiveness,
        vaccine_gender, ages, female_worms, male_worms, human_cercariae, eggs,
        treated, vaccine_duration, vac_status, vaccine_round, gender, access)[1],
        [[1,1],[0,0],[0,0]])
end

female_worms= [[1,1], [2,1], [2,4]]
male_worms= [[1,1], [2,1], [2,5]]
access = [1,1,1]

@testset "vaccinate" begin
@test isapprox(vaccinate(vaccine_coverage, min_age_vaccine, max_age_vaccine, vaccine_effectiveness,
            vaccine_gender, ages, female_worms, male_worms, human_cercariae, eggs,
            treated, vaccine_duration, vac_status, vaccine_round, gender,access)[1],
            [[1,1],[0,0],[0,0]])
end

female_worms= [[1,1], [2,1], [2,4]]
male_worms= [[1,1], [2,1], [2,5]]
access = [1,0,0]

@testset "vaccinate" begin
@test isapprox(vaccinate(vaccine_coverage, min_age_vaccine, max_age_vaccine, vaccine_effectiveness,
            vaccine_gender, ages, female_worms, male_worms, human_cercariae, eggs,
            treated, vaccine_duration, vac_status, vaccine_round, gender,access)[1],
            [[1,1],[2,1],[2,4]])
end

time_step = 10
num_time_steps = 1
community = 1,3,1
community_contact_rate = 1,0.5, 1
community_probs = 1,2,1
heavy_burden_threshold = 16
kato_katz_par = 0.87
use_kato_katz = 0
@testset "update_env" begin
    @test update_env(num_time_steps, [1,3,4], [2,5,8],community, community_contact_rate, community_probs, human_cercariae, female_worms, male_worms,
    time_step, 5.7,
    eggs, 0.34, 0.03, 2,
    vac_status, gender, 0.24,predis_weight,
    [0.2, 0.8,1], treated, vaccine_effectiveness,
    0.0005, death_prob_by_age, ages_for_deaths,
    [0,0,0], [0.02,0.04,0.03], env_miracidia,
    env_cercariae, 0.0005, 1, 1,
    1, 1, contact_rates_by_age,
    28*time_step/(1000*365), [], [], [1,1,1], 1,
    [1,1,1], 1,
    1/24,1, 24, heavy_burden_threshold,
    kato_katz_par, use_kato_katz)[1] ==  [1+(num_time_steps*time_step/365),3+(num_time_steps*time_step/365),4+(num_time_steps*time_step/365)]
end

mda_info = create_mda(0, .75, 1, 0, 5, 2, [0,1], [0,1], [0,1], .92)
push!(vaccine_info, vaccine_information(0.75, 4, 16, [0,1], 3, 0))
@testset "update_env" begin
    @test update_env(num_time_steps, [1,3,4], [2,5,8],community, community_contact_rate, community_probs, human_cercariae, female_worms, male_worms,
    time_step, 5.7,
    eggs, 0.34, 0.03, 2,
    vac_status, gender, 0.24,predis_weight,
    [0.2, 0.8,1], treated, vaccine_effectiveness,
    0.0005, death_prob_by_age, ages_for_deaths,
    [0,0,0], [0.02,0.04,0.03], env_miracidia,
    env_cercariae, 0.0005, 1, 1,
    1, 1, contact_rates_by_age,
    28*time_step/(1000*365), mda_info, vaccine_info, [1,1,1], 1,
    [1,1,1], 1,
    1/24,1, 24, heavy_burden_threshold,
    kato_katz_par, use_kato_katz)[1] ==  [1+(num_time_steps*time_step/365),3+(num_time_steps*time_step/365),4+(num_time_steps*time_step/365)]
end

spec_ages = 7639, 7082, 6524, 5674, 4725, 4147, 3928, 3362,
            2636, 1970, 1468, 1166, 943, 718, 455, 244
ages_per_index = 5

@testset "gen_age_dist" begin
    @test generate_age_distribution(spec_ages, ages_per_index)[end] == 1
end

@testset "gen_age_dist" begin
    @test generate_age_distribution(spec_ages, ages_per_index)[1] == spec_ages[1] / (sum(spec_ages)*5)
end

Random.seed!(25251)
@testset "gen_age_dist" begin
    @test specified_age_distribution(5, spec_ages, ages_per_index)[1] == 27
end

 mda_access = 1
N = 1000
max_age = 100
initial_worms = 10
time_step = 10
worm_stages = 2
female_factor = 1
male_factor = 1
initial_miracidia = 1
initial_miracidia_days = trunc(Int,round(41/time_step, digits = 0))
env_cercariae = 0
#const contact_rate = 0.000005
age_death_rate_per_1000 = [6.56, 0.93, 0.3, 0.23, 0.27, 0.38, 0.44, 0.48,0.53, 0.65,
                           0.88, 1.06, 1.44, 2.1, 3.33, 5.29, 8.51, 13.66,
                           21.83, 29.98, 36.98]
predis_aggregation = 0.24
mda_adherence = 0.8
scenario = "high adult"
contact_rates_by_age = make_age_contact_rate_array(max_age, scenario, [], [])

predis_weight = 1

N_communities= 4
community_probs = [1,2,1,3]

pop = create_population_specified_ages(N, N_communities, community_probs, initial_worms, contact_rates_by_age,
        worm_stages, female_factor, male_factor,initial_miracidia,
        initial_miracidia_days, predis_aggregation, predis_weight,
        time_step,
        spec_ages, ages_per_index, death_prob_by_age, ages_for_deaths,
        mda_adherence, mda_access)
        # These should all get give the initial value
@testset "miracidia" begin
    @test all(pop[14] .== initial_miracidia)
end

        # Worms should be 0:ininity
        # Get first column of worms
mworm1 = [pop[9][i][1] for i=1:length(pop[9])]
fworm1 = [pop[10][i][1] for i=1:length(pop[10])]
@testset "wormsm" begin
    @test all(mworm1 .>= 0)
end
@testset "wormsf" begin
    @test all(fworm1 .>= 0)
end

        # All vectors except the last should be length N?
lens = [length(pop[i]) for i=1:(length(pop) - 3)]
push!(lens, length(pop[end]))
push!(lens, length(pop[end-1]))
@testset "N" begin
    @test all(lens .== N)
end

@testset "death_of_human_s_p" begin
    @test death_of_human_save_predis([2,4], [0,6], [1,0], [0.0002,0.00005], [[2,3,4],[6,3,4]], [15,7],
                                [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                                [0,0], [1,1], 1,
                                [1,1], [1,1])[1] == [4]

end

a = death_of_human([2,4], [0,6], [1,0], [0.0002,0.00005], [[2,3,4],[6,3,4]], [15,7],
                                [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                                [0,0], [1,1], 1,
                                [1,1], [1,1])

b = death_of_human_save_predis([2,4], [0,6], [1,0], [0.0002,0.00005], [[2,3,4],[6,3,4]], [15,7],
                           [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                           [0,0], [1,1], 1,
                           [1,1], [1,1])

@testset "death_of_human_s_p" begin
    @test length(a) == length(b)-1
end

@testset "birth_of_human_s_p" begin
    @test birth_of_human_specified_predis([2,4], [0.44,0], [0,0], [1.1,1], [2,4],[[2,3,4],[6,3,4]], [15,7],
                                [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                                [0,0], [1.1,1], 1, 1, [0.0002,0.00005], 2,1, 1, [1,1],
                            death_prob_by_age, ages_for_deaths,[1,2,3,4], 0.8,[1,1], 0.9, 100)[1] == [2,4,0]
end
b1 = birth_of_human_specified_predis([2,4], [0.44,0], [0,0], [1.1,1], [2,4],[[2,3,4],[6,3,4]], [15,7],
                            [0,0], [0,0], [[9,2],[5,3]], [[0,3],[1,12]],
                            [0,0], [1.1,1], 1, 1, [0.0002,0.00005], 2,1, 1, [1,1],
                        death_prob_by_age, ages_for_deaths,[1,2,3,4], 0.8,[1,1], 0.9, 100)

@testset "birth_of_human_s_p" begin
    @test b1[4][end]== 100
end

gamma_k = Gamma(0.87,1/0.87)

@testset "kato_katz" begin
    @test  kato_katz(100, gamma_k) > 0
end

@testset "get_prev" begin
    @test isapprox(get_prevalences(pop[1], pop[7], 0, 16, 0.87, 0).population_burden, [0,0,0])
end

for i in 1:length(pop[7])
    pop[7][i] = 100
end

@testset "get_prev" begin
    @test isapprox(get_prevalences(pop[1], pop[7], 0, 16, 0.87, 0).population_burden, [100,100,100])
end

for i in 1:length(pop[7])
    pop[7][i] = 15
end

@testset "get_prev" begin
    @test isapprox(get_prevalences(pop[1], pop[7], 0, 16, 0.87, 0).population_burden, [100,100,0])
end

for i in 1:length(pop[7])
    pop[7][i] = 2
end

@testset "get_prev" begin
    @test isapprox(get_prevalences(pop[1], pop[7], 0, 16, 0.87, 0).population_burden, [100,0,0])
end

ages , death_ages, gender, predisposition, community,
   human_cercariae, eggs, vac_status,
   treated, female_worms, male_worms, age_contact_rate,
   vaccinated, env_miracidia, adherence, access = create_population_specified_ages(N, N_communities, community_probs, initial_worms, contact_rates_by_age,
                worm_stages, female_factor, male_factor,initial_miracidia,
                initial_miracidia_days, predis_aggregation, predis_weight,
                time_step,
                spec_ages, ages_per_index, death_prob_by_age, ages_for_deaths,
                mda_adherence, mda_access)

filename = "a.jld"

save_population_to_file(filename, ages, gender, predisposition, community, human_cercariae, eggs, vac_status, treated,
        female_worms, male_worms, vaccinated, age_contact_rate, death_ages, env_miracidia, env_cercariae, adherence, access)

@testset "load_data" begin
    aa  = load_population_from_file(filename, N, true)
    @test isapprox(aa[1],ages)
end

@testset "load_data" begin
    aa  = load_population_from_file(filename, 100, false)
    @test length(aa[1]) ==100
end

max_fecundity = 0.34

ages , death_ages, gender, predisposition, community,
   human_cercariae, eggs, vac_status,
   treated, female_worms, male_worms, age_contact_rate,
   vaccinated, env_miracidia, adherence, access = create_population_specified_ages(N, 1, 1, initial_worms, contact_rates_by_age,
                worm_stages, female_factor, male_factor,initial_miracidia,
                initial_miracidia_days, predis_aggregation, predis_weight,
                time_step,
                spec_ages, ages_per_index, death_prob_by_age, ages_for_deaths,
                mda_adherence, mda_access)

                    death_ages = []
                    for i in 1:N
                         push!(death_ages, get_death_age(death_prob_by_age, ages_for_deaths))
                    end
                    mean(death_ages)
                    ages, death_ages = generate_ages_and_deaths(2000, ages, death_ages, death_prob_by_age, ages_for_deaths, time_step)
age_contact_rate = update_contact_rate(ages, age_contact_rate, contact_rates_by_age)

mda_info = []
vaccine_info = []

@testset "update_env_same_pop" begin
ages_equ, death_ages_equ, gender_equ, predisposition_equ, community_equ, human_cercariae_equ, eggs_equ,
vac_status_equ, treated_equ, female_worms_equ, male_worms_equ,
vaccinated_equ, age_contact_rate_equ,
env_miracidia_equ, env_cercariae_equ, adherence_equ,access_equ,
record = update_env_keep_population_same(1, copy(ages), copy(death_ages), copy(community), 1,1,
    copy(human_cercariae), copy(female_worms), copy(male_worms),
    time_step, average_worm_lifespan,
    copy(eggs), max_fecundity, r, worm_stages,
    copy(vac_status), copy(gender), predis_aggregation,predis_weight,
    copy(predisposition), copy(treated), vaccine_effectiveness,
    density_dependent_fecundity, death_prob_by_age, ages_for_deaths,
    copy(vaccinated), copy(age_contact_rate), copy(env_miracidia),
    copy(env_cercariae), contact_rate, env_cercariae_survival_prop, env_miracidia_survival_prop,
    female_factor, male_factor, contact_rates_by_age,
    birth_rate, mda_info, vaccine_info, adherence, mda_adherence, access, mda_access,
    record_frequency, human_cercariae_prop, miracidia_maturity_time, heavy_burden_threshold,
    kato_katz_par, use_kato_katz)

@test env_miracidia_equ[end] > env_miracidia_equ[end-1]
end

@testset "update_env_same_pop_save_pre" begin
ages_equ, death_ages_equ, gender_equ, predisposition_equ, community_equ, human_cercariae_equ, eggs_equ,
vac_status_equ, treated_equ, female_worms_equ, male_worms_equ,
vaccinated_equ, age_contact_rate_equ,
env_miracidia_equ, env_cercariae_equ, adherence_equ,access_equ,
record = update_env_keep_population_same_save_predisposition(1, copy(ages), copy(death_ages), copy(community), 1,1,
    copy(human_cercariae), copy(female_worms), copy(male_worms),
    time_step, average_worm_lifespan,
    copy(eggs), max_fecundity, r, worm_stages,
    copy(vac_status), copy(gender), predis_aggregation,predis_weight,
    copy(predisposition), copy(treated), vaccine_effectiveness,
    density_dependent_fecundity, death_prob_by_age, ages_for_deaths,
    copy(vaccinated), copy(age_contact_rate), copy(env_miracidia),
    copy(env_cercariae), contact_rate, env_cercariae_survival_prop, env_miracidia_survival_prop,
    female_factor, male_factor, contact_rates_by_age,
    birth_rate, mda_info, vaccine_info, adherence, mda_adherence, access, mda_access,
    record_frequency, human_cercariae_prop, miracidia_maturity_time, heavy_burden_threshold,
    kato_katz_par, use_kato_katz)

@test env_miracidia_equ[end] > env_miracidia_equ[end-1]
end

@testset "update_env_to_equ" begin
    ages_equ, gender_equ, predisposition_equ,  human_cercariae_equ, eggs_equ,
    vac_status_equ, treated_equ, female_worms_equ, male_worms_equ,
    vaccinated_equ, env_miracidia_equ, env_cercariae_equ, record_high =
        update_env_to_equilibrium(1, copy(ages), copy(human_cercariae), copy(female_worms), copy(male_worms),
            copy(community),1,
            time_step, average_worm_lifespan,
            copy(eggs), max_fecundity, r, worm_stages,
            copy(vac_status), copy(gender), predis_aggregation,
            copy(predisposition), copy(treated), vaccine_effectiveness,
            density_dependent_fecundity, copy(vaccinated), copy(env_miracidia),
            copy(env_cercariae), contact_rate, env_cercariae_survival_prop, env_miracidia_survival_prop,
            female_factor, male_factor, contact_rates_by_age, record_frequency, copy(age_contact_rate),human_cercariae_prop,
            miracidia_maturity_time, heavy_burden_threshold, kato_katz_par, use_kato_katz)
    @test env_miracidia_equ[end] > env_miracidia_equ[end-1]
end

mda_info = create_mda(0, .75, 1, 1, 5, 2, [0,1], [0,1], [0,1], .92)
vaccine_info = []
push!(vaccine_info, vaccine_information(0.75, 4, 16, [0,1], 3, 0.4))
@testset "update_env_with_mda_no_births_deaths" begin
ages_equ, death_ages_equ, gender_equ, predisposition_equ, community_equ, human_cercariae_equ, eggs_equ,
vac_status_equ, treated_equ, female_worms_equ, male_worms_equ,
vaccinated_equ, age_contact_rate_equ,
env_miracidia_equ, env_cercariae_equ, adherence_equ,access_equ,
record = update_env_with_mda_no_births_deaths(1, copy(ages), copy(death_ages), copy(community), 1, 1,
    copy(human_cercariae), copy(female_worms), copy(male_worms),
    time_step, average_worm_lifespan,
    copy(eggs), max_fecundity, r, worm_stages,
    copy(vac_status), copy(gender), predis_aggregation,predis_weight,
    copy(predisposition), copy(treated), vaccine_effectiveness,
    density_dependent_fecundity, death_prob_by_age, ages_for_deaths,
    copy(vaccinated), copy(age_contact_rate), copy(env_miracidia),
    copy(env_cercariae), contact_rate, env_cercariae_survival_prop, env_miracidia_survival_prop,
    female_factor, male_factor, contact_rates_by_age,
    birth_rate, mda_info, vaccine_info, adherence, mda_adherence, access, mda_access,
    record_frequency, human_cercariae_prop, miracidia_maturity_time, heavy_burden_threshold,
    kato_katz_par, use_kato_katz)

@test env_miracidia_equ[end] > env_miracidia_equ[end-1]
end

@testset "update_env_with_mda_no_births_deaths" begin
ages_equ, death_ages_equ, gender_equ, predisposition_equ, community_equ, human_cercariae_equ, eggs_equ,
vac_status_equ, treated_equ, female_worms_equ, male_worms_equ,
vaccinated_equ, age_contact_rate_equ,
env_miracidia_equ, env_cercariae_equ, adherence_equ,access_equ,
record = update_env_with_mda_no_births_deaths(100, copy(ages), copy(death_ages), copy(community), 1, 1,
    copy(human_cercariae), copy(female_worms), copy(male_worms),
    time_step, average_worm_lifespan,
    copy(eggs), max_fecundity, r, worm_stages,
    copy(vac_status), copy(gender), predis_aggregation,predis_weight,
    copy(predisposition), copy(treated), vaccine_effectiveness,
    density_dependent_fecundity, death_prob_by_age, ages_for_deaths,
    copy(vaccinated), copy(age_contact_rate), copy(env_miracidia),
    copy(env_cercariae), contact_rate, env_cercariae_survival_prop, env_miracidia_survival_prop,
    female_factor, male_factor, contact_rates_by_age,
    birth_rate, mda_info, vaccine_info, adherence, mda_adherence, access, mda_access,
    record_frequency, human_cercariae_prop, miracidia_maturity_time, heavy_burden_threshold,
    kato_katz_par, use_kato_katz)

@test length(ages_equ) == length(ages)
end

times = []
prev = []
sac_prev = []
high_burden = []
high_burden_sac =[]
adult_prev = []
high_adult_burden = []

ages_equ, death_ages_equ, gender_equ, predisposition_equ, community_equ, human_cercariae_equ, eggs_equ,
vac_status_equ, treated_equ, female_worms_equ, male_worms_equ,
vaccinated_equ, age_contact_rate_equ,
env_miracidia_equ, env_cercariae_equ, adherence_equ,access_equ,
record = update_env_with_mda_no_births_deaths(100, copy(ages), copy(death_ages), copy(community), 1, 1,
    copy(human_cercariae), copy(female_worms), copy(male_worms),
    time_step, average_worm_lifespan,
    copy(eggs), max_fecundity, r, worm_stages,
    copy(vac_status), copy(gender), predis_aggregation,predis_weight,
    copy(predisposition), copy(treated), vaccine_effectiveness,
    density_dependent_fecundity, death_prob_by_age, ages_for_deaths,
    copy(vaccinated), copy(age_contact_rate), copy(env_miracidia),
    copy(env_cercariae), contact_rate, env_cercariae_survival_prop, env_miracidia_survival_prop,
    female_factor, male_factor, contact_rates_by_age,
    birth_rate, mda_info, vaccine_info, adherence, mda_adherence, access, mda_access,
    record_frequency, human_cercariae_prop, miracidia_maturity_time, heavy_burden_threshold,
    kato_katz_par, use_kato_katz)

@testset "collect_prevs" begin
    b = collect_prevs(times, prev, sac_prev, high_burden, high_burden_sac, adult_prev, high_adult_burden, record, 1)
    high_ad_burden = (p->p.adult_burden[3]).(record)
    x = randperm(length(high_ad_burden))[1]
    @test b[7][x][1] == high_ad_burden[x]
end

@testset "collect_prevs" begin
    b = collect_prevs(times, prev, sac_prev, high_burden, high_burden_sac, adult_prev, high_adult_burden, record, 2)
    high_ad_burden = (p->p.adult_burden[3]).(record)
    x = randperm(length(high_ad_burden))[1]
    @test ((b[7][x][2] == b[7][x][1]) & (b[7][x][2] == high_ad_burden[x]))
end


text:
using OptimPack
using Base.Test

function rosenbrock_init!{T<:Real}(x0::Array{T,1})
  x0[1:2:end] = -1.2
  x0[2:2:end] =  1.0
  return nothing
end

function rosenbrock_fg!{T<:Real}(x::Array{T,1}, gx::Array{T,1})
  const c1::T = 1
  const c2::T = 2
  const c10::T = 10
  const c200::T = 200
  x1 = x[1:2:end]
  x2 = x[2:2:end]
  t1 = c1 - x1
  t2 = c10*(x2 - x1.*x1)
  g2 = c200*(x2 - x1.*x1)
  gx[1:2:end] = -c2*(x1.*g2 + t1)
  gx[2:2:end] = g2
  return sum(t1.*t1) + sum(t2.*t2)
end

function rosenbrock_test(n::Integer=20, m::Integer=3; single::Bool=false)
  T = (single ? Float32 : Float64)
  x0 = Array(T, n)
  rosenbrock_init!(x0)
  lbfgs(rosenbrock_fg!, x0, m, verb=true)
end

# Run tests in double and single precisions.
for (T, prec) in ((Float64, "double"), (Float32, "single"))

    x0 = Array(T, 20)
    rosenbrock_init!(x0)

    @printf("\nTesting NLCG in %s precision\n", prec)
    x1 = OptimPack.nlcg(rosenbrock_fg!, x0, verb=true)
    @test_approx_eq_eps x1 ones(T,20) 1e-4

    @printf("\nTesting VMLM in %s precision with Oren & Spedicato scaling\n", prec)
    x2 = OptimPack.vmlm(rosenbrock_fg!, x0, verb=true,
                        scaling=OptimPack.SCALING_OREN_SPEDICATO)
    @test_approx_eq_eps x2 ones(T,20) 1e-4

    @printf("\nTesting VMLM in %s precision with Barzilai & Borwein scaling\n", prec)
    x3 = OptimPack.vmlm(rosenbrock_fg!, x0, verb=true,
                    scaling=OptimPack.SCALING_BARZILAI_BORWEIN)
    @test_approx_eq_eps x3 ones(T,20) 1e-4

end

nothing

ext: jl


text:
"""
    set_param!(ref::ComponentReference, name::Symbol, value)

Set a component parameter as `set_param!(reference, name, value)`.
"""
function set_param!(ref::ComponentReference, name::Symbol, value)
    set_param!(ref.model, ref.comp_name, name, value)
end

"""
    set_param!(ref.model, ref.comp_name, name, value)

Set a component parameter as `reference[symbol] = value`.
"""
function Base.setindex!(ref::ComponentReference, value, name::Symbol)
    set_param!(ref.model, ref.comp_name, name, value)
end

"""
    connect_param!(dst::ComponentReference, dst_name::Symbol, src::ComponentReference, src_name::Symbol)

Connect two components as `connect_param!(dst, dst_name, src, src_name)`.
"""
function connect_param!(dst::ComponentReference, dst_name::Symbol, src::ComponentReference, src_name::Symbol)
    connect_param!(dst.model, dst.comp_id, dst_name, src.comp_id, src_name)
end

"""
    connect_param!(dst::ComponentReference, src::ComponentReference, name::Symbol)

Connect two components with the same name as `connect_param!(dst, src, name)`.
"""
function connect_param!(dst::ComponentReference, src::ComponentReference, name::Symbol)
    connect_param!(dst.model, dst.comp_id, name, src.comp_id, name)
end

"""
    Base.getindex(comp_ref::ComponentReference, var_name::Symbol)

Get a variable reference as `comp_ref[var_name]`.
"""
function Base.getindex(comp_ref::ComponentReference, var_name::Symbol)
    VariableReference(comp_ref.model, comp_ref.comp_name, var_name)
end

"""
    Base.setindex!(comp_ref::ComponentReference, var_ref::VariableReference, var_name::Symbol)

Connect two components as `comp_ref[var_name] = var_ref`.
"""
function Base.setindex!(comp_ref::ComponentReference, var_ref::VariableReference, var_name::Symbol)
    comp_ref.model == var_ref.model || error("Can't connect variables defined in different models")

    connect_param!(comp_ref.model, comp_ref.comp_name, var_name, var_ref.comp_name, var_ref.var_name)
end


text:
# syntax: proto2
using Compat
using ProtoBuf
import ProtoBuf.meta

mutable struct RequestHeaderProto <: ProtoType
    methodName::AbstractString
    declaringClassProtocolName::AbstractString
    clientProtocolVersion::UInt64
    RequestHeaderProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RequestHeaderProto
const __req_RequestHeaderProto = Symbol[:methodName,:declaringClassProtocolName,:clientProtocolVersion]
meta(t::Type{RequestHeaderProto}) = meta(t, __req_RequestHeaderProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export RequestHeaderProto


text:
using Distributions
using FreqTables
using Gadfly
using Cairo

#a model of how *sequences* as opposed to
# individual codons evolve
#We implement a random walk model of
#the biomial distribution
#We a force can be added by biasing the probability
#to step to the left. see tag: ##BIAS
# for the relevant line of code.
# When the bias is 0, then the model just reproduces the binomial distribution

srand(1001)
seqLen=6
energs=Array{Float64}(seqLen+1)
energs[1]=0
rang=collect(0:seqLen)
seq=rand(rang)
seq=3
seq2= seqLen - seq
#sequence is (seq,seq2)
numtrial=100000# number of points we sample
relaxTime=2000# time we give the chain to relax
pow=2.0

#outer for
# for ty = collect(1)
rec=[]
for j = 1: numtrial

# we add some random elemnt to relax time to avoid artefacts
numexp=relaxTime +Int64(round(rand()*100))# we add some random elemnt to relax time to avoid artefacts

################
for i = 1:numexp
ppr= Float64((seqLen - seq)^pow) # we assume that seq is k. and on the left of the chain, k=0
ppl =Float64(seq)

pl=ppl/(ppl+ppr)
pr=ppr/(ppl+ppr)

#seq2 is strictly for nothing, I include it anyway
#println(seq, " ", seq2, " ", pl, " " , pr)

if(rand() < pl)
seq = seq -1
seq2 = seq2 + 1
else
seq = seq +1
seq2=seq2 -1
end#if

#println(entropy, "  ", ppr, "   ",ppl, "  seq: ",seq )

end#inner for

#record the last value from last time
push!(rec,seq)
end#outer for

#Calculate the energies:
energs[1]=0;
for i = 2: seqLen+1
l=i-1
energs[i] = -log((seqLen-l +1)^pow/(l)) + energs[i-1]
end
println(energs)
Z=sum(exp.(-energs))
println(Z)
probs=(exp.(-energs))/Z
println(probs)

 outp= [collect(1:numtrial) rec]

#tbl= freqtable(Array{Float64,1}(round.(10000000*(pdf.(dis,rec)))/10000000))
tbl1= freqtable(Array{Float64,1}(rec))
#println(tbl1)
tmp=convert(Array, tbl1)
tmp=tmp/sum(tmp)
eprobs=tmp
println(tmp)
sequ=names(tbl1)[1]#sequences
cub=sum(tmp.*(sequ/seqLen))
#dis=Binomial(seqLen,cub)
dis=Binomial(seqLen,0.14322)
tprobs=pdf.(dis,names(tbl1)[1])

tmp=convert(Array, tbl1)
tmp=tmp/sum(tmp)
#println("Codon usage bias ",cub)

#plot(x=-tprobs,y=-eprobs,Geom.point,Geom.smooth(method=:lm))

tmp=linreg(-eprobs,-tprobs)
println(cub," ",tmp[2])
#end#outer for ty

max_stars_repo_name: jkobject/CodonUsageBias

max_issues_repo_name: jkobject/CodonUsageBias

max_forks_repo_name: jkobject/CodonUsageBias


text:
using CSV
using CodecZstd

# input1 = "../Data/Cortical_SMART/refLoading.csv"
# input2 = "../Data/Cortical_SMART/t_ScaledLogData.csv"
# input3 = "../Data/Cortical_SMART/Data.zst"

# input1="/data/koki/TestData/single_cell/10XGenomics/1M_neurons/refLoading.csv"
# input2="/data/koki/TestData/single_cell/10XGenomics/1M_neurons/t_scale_log_1M_neurons_filtered_gene_bc_matrices_h5.csv"
# input3="/data/koki/TestData/single_cell/10XGenomics/1M_neurons/1M_neurons_filtered_gene_bc_matrices_h5.zst"
# output1="/data2/koki/ICCIPCA/Data/Brain_Chromium/Deconvolution/Eigen_vectors.csv"
# output2="/data2/koki/ICCIPCA/Data/Brain_Chromium/Deconvolution/Eigen_values.csv"

input1 = ARGS[1]
input2 = ARGS[2]
input3 = ARGS[3]
output1 = ARGS[4]
output2 = ARGS[5]

# Return N, M
function nm(input::AbstractString)
    N = zeros(UInt32, 1)
    M = zeros(UInt32, 1)
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, N)
        read!(stream, M)
        close(stream)
    end
    return N[], M[]
end

# Loading
println("Loading...")
# # 8.590286 seconds (6.13 M allocations: 310.771 MiB, 1.01% gc time)
# @time loading = CSV.read(input1, datarow=1)
# 0.533292 seconds (691.03 k allocations: 26.640 MiB, 8.03% gc time)
loading = readcsv(input1)

# Setting
N, M = nm(input3)
N = Int(N)
M = Int(M)
dim = size(loading)[2]
W = zeros(M, dim)

# Split
println("Spliting...")
open(input2 , "r") do f
	counter = 0
    while !eof(f)
        counter += 1
        print("\r", counter)
        xx = readline(f)
        xx = split(xx, ",") # Splited as Character
        # Assume input data is Integer
        x = zeros(Float64, N)
        for i=1:N
            x[i] = parse(Float64, xx[i])
        end
        W[counter, :] = x' * loading
        if counter == M
            break
        end
    end
end

# Eigen value
println("Eigen value...")
σ = Float32[norm(W[:, x]) for x=1:dim]
for n = 1:dim
    W[:, n] ./= σ[n]
end
λ = σ .* σ / M

# Sort by Eigen value
println("Sorting...")
idx = sortperm(λ, rev=true)
W .= W[:, idx]
λ .= λ[idx]

# Output
println("Outputing...")
writedlm(output1, W, ',')
writedlm(output2, λ, ',')

max_stars_repo_name: rikenbit/onlinePCA-experiments

max_issues_repo_name: rikenbit/onlinePCA-experiments

max_forks_repo_name: rikenbit/onlinePCA-experiments


text:
using Documenter

deploydocs(;
    repo="github.com/OutlierDetectionJL/OutlierDetection.jl",
    target = "site",
    push_preview = true,
    deps = Deps.pip("mkdocs==1.1.2", "mkdocs-material==7.1.2", "pymdown-extensions==8.1.1",
                    "pygments==2.8.1", "Jinja2==2.11.3"),
    make = () -> run(`mkdocs build`)
)


text:
using Gtk, Graphics, Cairo, Colors, ColorSchemes

function update_window!(actr, ex)
    refresh!(ex)
    draw_target!(actr, ex)
    draw_focus!(actr, ex)
    bounds = get_min_max(actr)
    for vo in actr.visual_location.iconic_memory
        heat = compute_heat(vo, bounds...)
        draw_object!(ex, vo, heat)
    end
end

function draw_cross!(actr, ex)
    refresh!(ex)
    draw_focus!(actr, ex)
    ex.visible ? sleep(.30/ex.speed) : nothing
    return nothing
end

function draw_object!(ex, vo, heat)
    !vo.visible ? (return) : nothing
    c = ex.canvas
    w = vo.width
    x,y = vo.location
    color = get_color(vo)
    letter = get_text(vo)
    @guarded draw(c) do widget
        ctx = getgc(c)
        circle(ctx, x, y, w/2)
        if vo.attended
            α = get(ColorSchemes.Greys_3, heat)
            set_source_rgba(ctx, α.r, α.g, α.b , .4)
        else
            α = get(ColorSchemes.coolwarm, heat)
            set_source_rgba(ctx, α.r, α.g, α.b , .4)
        end
        fill(ctx)
        select_font_face(ctx, "Arial", Cairo.FONT_SLANT_NORMAL,
             Cairo.FONT_WEIGHT_BOLD);
        set_font_size(ctx, pixels_to_points(vo.width))
        set_source(ctx, Cairo.alphacolor(color,.5))
        extents = text_extents(ctx, letter)
        x′ = x - (extents[3]/2 + extents[1])
        y′ = y - (extents[4]/2 + extents[2])
        move_to(ctx, x′, y′)
        show_text(ctx, letter)
    end
    Gtk.showall(c)
    return nothing
end

function draw_focus!(actr, ex)
    c = ex.canvas
    w = actr.visual_location.visicon[1].width
    x,y = actr.visual.focus
    @guarded draw(c) do widget
        ctx = getgc(c)
        circle(ctx, x, y, w/2)
        set_line_width(ctx, 4);
        set_source_rgba(ctx, 0, 0, 0, 1)
        Cairo.stroke(ctx)
    end
    Gtk.showall(c)
    return nothing
end

function draw_target!(actr, ex)
    target = filter(x->x.target, actr.visual_location.iconic_memory)
    isempty(target) ? (return) : nothing
    c = ex.canvas
    w = target[1].width
    x,y = target[1].location
    @guarded draw(c) do widget
        ctx = getgc(c)
        set_line_width(ctx, 4)
        circle(ctx, x, y, w/2)
        set_source_rgba(ctx, 255, 255, 179, 1)
        Cairo.stroke(ctx)
    end
    Gtk.showall(c)
    return nothing
end

function refresh!(ex)
    c = ex.canvas
    w = ex.array_width
    @guarded draw(c) do widget
        ctx = getgc(c)
        rectangle(ctx, 0, 0, w, w)
        set_source_rgb(ctx, .8, .8, .8)
        fill(ctx)
    end
    Gtk.showall(c)
    return nothing
end

pixels_to_points(pixels) = pixels*3/4

points_to_pixels(points) = points*4/3

get_text(vo) = string(vo.features.shape.value)

get_color(vo) = Colors.parse(Colorant, string(vo.features.color.value))

get_min_max(actr::ACTRV) = get_min_max(actr.visual_location.iconic_memory)

function get_min_max(iconic_memory)
    mn,mx = Inf,-Inf
    for vo in iconic_memory
        vo.activation < mn ? (mn=vo.activation; continue) : nothing
        vo.activation > mx ? (mx=vo.activation; continue) : nothing
    end
    return mn,mx
end

function compute_heat(vo, lb, ub)
    lb == ub ? (return 1.0) : nothing
    return (vo.activation-lb)/(ub-lb)
end

function setup_window(array_width)
	canvas = @GtkCanvas()
    window = GtkWindow(canvas, "ACT-R", array_width, array_width)
    Gtk.visible(window, true)
    @guarded draw(canvas) do widget
        ctx = getgc(canvas)
        rectangle(ctx, 0.0, 0.0, array_width, array_width)
        set_source_rgb(ctx, .8, .8, .8)
        fill(ctx)
    end
	return canvas,window
end


text:
module Utils

include("MapViews.jl")
include("WfnUtils.jl")

using .MapViews
export MapView

using .WfnUtils
export applygate!

end


text:
# Issue #14

using JuliaInterpreter: JuliaInterpreter, pc_expr, evaluate_call!, finish_and_return!, @lookup, enter_call_expr, breakpoints
runframe(frame::Frame, pc=frame.pc[]) = Some{Any}(finish_and_return!(Compiled(), frame))

frame = @make_frame map(x->2x, 1:10)
state = dummy_state(frame)
execute_command(state, Val{:so}(), "so")
@test isnothing(state.frame)
@test state.overall_result == 2 .* [1:10...]

# Issue #12
function complicated_keyword_stuff(args...; kw...)
    args[1] == args[1]
    (args..., kw...)
end
frame = @make_frame complicated_keyword_stuff(1)
state = dummy_state(frame)
execute_command(state, Val{:n}(), "n")
execute_command(state, Val{:so}(), "so")
@test isnothing(state.frame)

@test runframe(JuliaInterpreter.enter_call(complicated_keyword_stuff, 1, 2)) ==
      runframe(@make_frame(complicated_keyword_stuff(1, 2)))
@test runframe(JuliaInterpreter.enter_call(complicated_keyword_stuff, 1, 2; x=7, y=33)) ==
      runframe(@make_frame(complicated_keyword_stuff(1, 2; x=7, y=33)))

# Issue #22
f22() = string(:(a+b))
@test step_through(enter_call_expr(:($f22()))) == "a + b"
f22() = string(QuoteNode(:a))
@test step_through(enter_call_expr(:($f22()))) == ":a"

# Breakpoints
f_outer_break(x) = f_inner_break(x)
f_inner_break(x) = x
JuliaInterpreter.breakpoint(f_inner_break)
frame = @make_frame f_outer_break(2)
state = dummy_state(frame)
execute_command(state, Val{:so}(), "c")
@test state.frame.framecode.scope.name === :f_inner_break
execute_command(state, Val{:so}(), "c")
@test state.frame === nothing
@test state.overall_result == 2

@inline fnothing(x) = 1
frame = @make_frame fnothing(0)
@test chomp(sprint(Debugger.print_next_expr, frame)) == "About to run: return 1"

function f()
    x = 1 + 1
    @info "hello"
end
# https://github.com/JuliaDebug/Debugger.jl/issues/71
frame = Debugger.@make_frame f()
state = dummy_state(frame)
execute_command(state, Val{:n}(), "n")
defline, current_line, body = Debugger.locinfo(state.frame)
@test occursin("handle_message(logger, level, msg", body)

f_unicode() = √

try
    Debugger.set_highlight(Debugger.HIGHLIGHT_SYSTEM_COLORS)
    frame = Debugger.@make_frame f()
    st = chomp(sprint(Debugger.print_status, frame; context = :color => true))
    x_1_plus_1_colored = "\e[39m\e[97mx\e[39m\e[97m \e[39m\e[91m=\e[39m\e[97m \e[39m1\e[97m \e[39m\e[91m+\e[39m\e[97m \e[39m1\e[97m"
    @test occursin(x_1_plus_1_colored, st)

    frame = Debugger.@make_frame f_unicode()
    @info "The eventual warning below is expected:"
    st = chomp(sprint(Debugger.print_status, frame; context = :color => true))
    @test occursin("√", st)
finally
    Debugger.set_highlight(Debugger.HIGHLIGHT_OFF)
end

frame = @make_frame Test.eval(1)
desc = Debugger.locdesc(frame)
@test occursin(Sys.STDLIB, desc)

import InteractiveUtils
@testset "`o` command" begin
    g() = nothing
    LINE = (@__LINE__) - 1
    frame = Debugger.@make_frame g()
    state = dummy_state(frame)
    JuliaInterpreter.breakpoint(InteractiveUtils.edit)
    frame, bp = JuliaInterpreter.@interpret execute_command(state, Val{:o}(), "o")
    JuliaInterpreter.remove()
    locals = JuliaInterpreter.locals(frame)
    @test JuliaInterpreter.Variable(@__FILE__, :file, false) in locals
    @test JuliaInterpreter.Variable(LINE, :line, false) in locals
end

# These are LoadError because the error happens at macro expansion
@test_throws LoadError @macroexpand @enter "foo"
@test_throws LoadError @macroexpand @enter 1
@test_throws LoadError @macroexpand @run [1,2]

# Breakpoints
frame = Debugger.@make_frame sin(1.0)
state = dummy_state(frame)

# add
execute_command(state, Val{:bp}(), "bp add cos")
bp = breakpoints()[1]
@test bp.f === cos
@test bp.sig === nothing
@test bp.condition === nothing
@test bp.line == 0
execute_command(state, Val{:bp}(), "bp rm")
@test length(breakpoints()) == 0
execute_command(state, Val{:bp}(), "bp add cos(x)")
bp = breakpoints()[1]
@test bp.f === cos
@test bp.sig == Tuple{Float64}
JuliaInterpreter.remove()
execute_command(state, Val{:bp}(), "bp add cos(::Float32)")
bp = breakpoints()[1]
@test bp.sig == Tuple{Float32}
JuliaInterpreter.remove()
execute_command(state, Val{:bp}(), """bp add "foo.jl":10 x>3""")
bp = breakpoints()[1]
@test bp.path == "foo.jl"
@test bp.line == 10
@test bp.condition == :(x > 3)
JuliaInterpreter.remove()
execute_command(state, Val{:bp}(), """bp add 10""")
bp = breakpoints()[1]
@test bp.line == 10
@test bp.path == CodeTracking.whereis(@which sin(1.0))[1]

# toggle
JuliaInterpreter.remove()
execute_command(state, Val{:bp}(), """bp add sin""")
execute_command(state, Val{:bp}(), """bp add cos""")
execute_command(state, Val{:bp}(), """bp toggle 1""")
bp = breakpoints()[1]
bp2 = breakpoints()[2]
@test bp.enabled[] == false
@test bp2.enabled[] == true
execute_command(state, Val{:bp}(), """bp toggle""")
@test bp.enabled[] == true
@test bp2.enabled[] == false

# disable / enable
JuliaInterpreter.remove()
execute_command(state, Val{:bp}(), """bp add sin""")
execute_command(state, Val{:bp}(), """bp add cos""")
execute_command(state, Val{:bp}(), """bp disable 1""")
bp = breakpoints()[1]
bp2 = breakpoints()[2]
@test bp.enabled[] == false
execute_command(state, Val{:bp}(), """bp enable 1""")
@test bp.enabled[] == true
execute_command(state, Val{:bp}(), """bp disable""")
@test bp.enabled[] == false
@test bp2.enabled[] == false

JuliaInterpreter.remove()


text:
import Base.isempty

struct Record
    names::Vector{String}
    functions::Dict{String, Function}
    data::Dict{String, AbstractVector}

    Record() = new(String[], Dict(), Dict())
end

Base.isempty(rec::Record) = isempty(rec.names)

function register!(rec, name, func, data_type)
    push!(rec.names, name)
    rec.functions[name] = func
    rec.data[name] = data_type[]
end

function add_record!(rec, g, pop, cfg)
    for name in rec.names
        push!(rec.data[name], rec.functions[name](g, pop, cfg))
    end
end

train_fitness(g, pop, cfg) = @views minimum(getfield.(entry_pool[pop], :train_fitness))
test_fitness(g, pop, cfg) = @views minimum(getfield.(entry_pool[pop], :test_fitness))

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
# 1 / factorial(i)

# how many terms are to be used
const nused_inv_fact = 18

const inv_fact = (
  Double64(1.0,0.0),
  Double64(0.5,0.0),
  Double64(0.16666666666666666,9.25185853854297e-18),
  Double64(0.041666666666666664,2.3129646346357427e-18),
  Double64(0.008333333333333333,1.1564823173178714e-19),
  Double64(0.001388888888888889,-5.300543954373577e-20),
  Double64(0.0001984126984126984,1.7209558293420705e-22),
  Double64(2.48015873015873e-5,2.1511947866775882e-23),
  Double64(2.7557319223985893e-6,-1.858393274046472e-22),
  Double64(2.755731922398589e-7,2.3767714622250297e-23),
  Double64(2.505210838544172e-8,-1.448814070935912e-24),
  Double64(2.08767569878681e-9,-1.20734505911326e-25),
  Double64(1.6059043836821613e-10,1.2585294588752098e-26),
  Double64(1.1470745597729725e-11,2.0655512752830745e-28),
  Double64(7.647163731819816e-13,7.03872877733453e-30),
  Double64(4.779477332387385e-14,4.399205485834081e-31),
  Double64(2.8114572543455206e-15,1.6508842730861433e-31),
  Double64(1.5619206968586225e-16,1.1910679660273754e-32),
  Double64(8.22063524662433e-18,2.2141894119604265e-34),
  Double64(4.110317623312165e-19,1.4412973378659527e-36),
  Double64(1.9572941063391263e-20,-1.3643503830087908e-36),
  Double64(8.896791392450574e-22,-7.911402614872376e-38),
  Double64(3.868170170630684e-23,-8.843177655482344e-40),
  Double64(1.6117375710961184e-24,-3.6846573564509766e-41),
  Double64(6.446950284384474e-26,-1.9330404233703465e-42),
  Double64(2.4795962632247976e-27,-1.2953730964765229e-43),
  Double64(9.183689863795546e-29,1.4303150396787322e-45),
  Double64(3.279889237069838e-30,1.5117542744029879e-46),
  Double64(1.1309962886447716e-31,1.0498015412959506e-47),
  Double64(3.7699876288159054e-33,2.5870347832750324e-49),
  Double64(1.216125041553518e-34,5.586290567888806e-51),
  Double64(3.8003907548547434e-36,1.7457158024652518e-52)
);

const ninv_fact = length(inv_fact)

# [exp(i) for i in 1:64]

const exp_int = (
 Double64(2.718281828459045, 1.4456468917292502e-16),
 Double64(7.38905609893065, -1.7971139497839148e-16),
 Double64(20.085536923187668, -1.8275625525512858e-16),
 Double64(54.598150033144236, 2.8741578015844115e-15),
 Double64(148.4131591025766, 3.4863514900464198e-15),
 Double64(403.4287934927351, 1.2359628024450387e-14),
 Double64(1096.6331584284585, 9.869752640434095e-14),
 Double64(2980.9579870417283, -2.7103295816873633e-14),
 Double64(8103.083927575384, -2.1530877621067177e-13),
 Double64(22026.465794806718, -1.3780134700517372e-12),
 Double64(59874.14171519782, 1.7895764888916994e-12),
 Double64(162754.79141900392, 5.30065881322063e-12),
 Double64(442413.3920089205, 1.2118711752313224e-11),
 Double64(1.2026042841647768e6, -1.5000525764327354e-11),
 Double64(3.2690173724721107e6, -3.075806431120808e-11),
 Double64(8.886110520507872e6, 5.321182483501564e-10),
 Double64(2.41549527535753e7, -7.203995068362157e-10),
 Double64(6.565996913733051e7, 1.4165536846555444e-9),
 Double64(1.7848230096318725e8, 1.333018530234341e-8),
 Double64(4.851651954097903e8, 4.880277289790406e-10),
 Double64(1.3188157344832146e9, 8.043448618843281e-8),
 Double64(3.584912846131592e9, -2.3519384005402157e-7),
 Double64(9.744803446248903e9, -6.74501500127677e-7),
 Double64(2.648912212984347e10, 7.670395527778119e-7),
 Double64(7.200489933738588e10, -6.992440211033874e-6),
 Double64(1.9572960942883878e11, -1.1364989227123904e-5),
 Double64(5.3204824060179865e11, -2.8335783945658822e-5),
 Double64(1.446257064291475e12, 7.602079742299693e-5),
 Double64(3.931334297144042e12, 8.220112058084352e-5),
 Double64(1.0686474581524463e13, -0.0007436345313492586),
 Double64(2.9048849665247426e13, -0.0005501643178883202),
 Double64(7.896296018268069e13, 0.007660978022635108),
 Double64(2.1464357978591606e14, 0.002124297761531261),
 Double64(5.834617425274549e14, 0.006402902734610391),
 Double64(1.5860134523134308e15, -0.02187035537422534),
 Double64(4.311231547115195e15, 0.22711342229285691),
 Double64(1.1719142372802612e16, -0.6912270602088098),
 Double64(3.1855931757113756e16, 0.22032867170129863),
 Double64(8.659340042399374e16, 2.953606932719265),
 Double64(2.3538526683702e17, -14.592100089250966),
 Double64(6.398434935300549e17, 37.22266340351557),
 Double64(1.739274941520501e18, 55.394681303611236),
 Double64(4.727839468229346e18, 257.4744575627443),
 Double64(1.2851600114359308e19, -12.1907003678569),
 Double64(3.4934271057485095e19, 436.0347972334061),
 Double64(9.496119420602448e19, 5929.133649117119),
 Double64(2.5813128861900675e20, -15192.714199784727),
 Double64(7.016735912097631e20, 30185.471599886117),
 Double64(1.9073465724950998e21, -98786.90015904616),
 Double64(5.184705528587072e21, 419031.45332293346),
 Double64(1.4093490824269389e22, -614323.8566876298),
 Double64(3.831008000716577e22, -661524.304512138),
 Double64(1.0413759433029089e23, -7.520901270665062e6),
 Double64(2.830753303274694e23, -4.711377645198593e6),
 Double64(7.694785265142018e23, -3.868399744098706e7),
 Double64(2.091659496012996e24, 5.079641515721467e7),
 Double64(5.685719999335932e24, 2.0801558082063326e8),
 Double64(1.545538935590104e25, 1.2092033491117463e8),
 Double64(4.2012104037905144e25, -1.7624059056928084e9),
 Double64(1.1420073898156842e26, 4.912247462314477e9),
 Double64(3.10429793570192e26, 3.39761293411071e9),
 Double64(8.438356668741454e26, 6.5719328084037315e10),
 Double64(2.29378315946961e27, -7.566162968773138e10),
 Double64(6.235149080811617e27, 1.3899738872492847e11),
);

const n_exp_int = length(exp_int)

# npio32[idx] = (idx-1)*pi / 32

const npio32 = (
Double64(0.0, 0.0),
Double64(0.09817477042468103, 3.827021247335479e-18),
Double64(0.19634954084936207, 7.654042494670958e-18),
Double64(0.2945243112740431, 1.1481063742006436e-17),
Double64(0.39269908169872414, 1.5308084989341915e-17),
Double64(0.4908738521234052, 1.9135106236677395e-17),
Double64(0.5890486225480862, 2.296212748401287e-17),
Double64(0.6872233929727672, 2.678914873134835e-17),
Double64(0.7853981633974483, 3.061616997868383e-17),
Double64(0.8835729338221293, 3.444319122601931e-17),
Double64(0.9817477042468103, 3.827021247335479e-17),
Double64(1.0799224746714915, -6.89250687418254e-17),
Double64(1.1780972450961724, 4.592425496802574e-17),
Double64(1.2762720155208536, -6.127102624715443e-17),
Double64(1.3744467859455345, 5.35782974626967e-17),
Double64(1.4726215563702156, -5.3616983752483455e-17),
Double64(1.5707963267948966, 6.123233995736766e-17),
Double64(1.6689710972195777, -4.5962941257812514e-17),
Double64(1.7671458676442586, 6.888638245203862e-17),
Double64(1.8653206380689398, -3.830889876314156e-17),
Double64(1.9634954084936207, 7.654042494670958e-17),
Double64(2.061670178918302, -3.06548562684706e-17),
Double64(2.159844949342983, -1.378501374836508e-16),
Double64(2.2580197197676637, 1.9904379115123167e-16),
Double64(2.356194490192345, 9.184850993605148e-17),
Double64(2.454369260617026, -1.5346771279128684e-17),
Double64(2.552544031041707, -1.2254205249430885e-16),
Double64(2.650718801466388, 2.1435187614057358e-16),
Double64(2.748893571891069, 1.071565949253934e-16),
Double64(2.84706834231575, -3.868628978676964e-20),
Double64(2.9452431127404313, -1.0723396750496691e-16),
Double64(3.0434178831651124, -2.144292487201471e-16),
Double64(3.141592653589793, 1.2246467991473532e-16),
Double64(3.2397674240144743, 1.5269398699555146e-17),
Double64(3.3379421944391554, -9.192588251562503e-17),
Double64(3.4361169648638366, -1.991211637308052e-16),
Double64(3.5342917352885173, 1.3777276490407724e-16),
Double64(3.6324665057131984, 3.057748368889706e-17),
Double64(3.7306412761378795, -7.661779752628312e-17),
Double64(3.8288160465625607, -1.838130787414633e-16),
Double64(3.9269908169872414, 1.5308084989341916e-16),
Double64(4.025165587411923, -3.9820364117182366e-16),
Double64(4.123340357836604, -6.13097125369412e-17),
Double64(4.221515128261284, 2.755842160979412e-16),
Double64(4.319689898685966, -2.757002749673016e-16),
Double64(4.417864669110647, 6.11936536675809e-17),
Double64(4.516039439535327, 3.9808758230246333e-16),
Double64(4.614214209960009, -1.5319690876277947e-16),
Double64(4.71238898038469, 1.8369701987210297e-16),
Double64(4.810563750809371, -3.675874711931398e-16),
Double64(4.908738521234052, -3.069354255825737e-17),
Double64(5.006913291658733, 3.062003860766251e-16),
Double64(5.105088062083414, -2.450841049886177e-16),
Double64(5.203262832508095, 9.180982364626472e-17),
Double64(5.301437602932776, 4.2870375228114717e-16),
Double64(5.399612373357457, -1.2258073878409563e-16),
Double64(5.497787143782138, 2.143131898507868e-16),
Double64(5.5959619142068195, -3.36971301214456e-16),
Double64(5.6941366846315, -7.737257957353928e-20),
Double64(5.792311455056181, 3.368165560553089e-16),
Double64(5.8904862254808625, -2.1446793500993382e-16),
Double64(5.988660995905543, 1.2242599362494854e-16),
Double64(6.086835766330225, -4.288584974402942e-16),
Double64(6.1850105367549055, -9.19645688054118e-17),
#Double64(6.283185307179586, 2.4492935982947064e-16)
);

const n_npio32 = length(npio32)

# sin_npio32[idx] = sin( (idx-1)*pi / 32 )

const sin_npio32 = (
Double64(0.0, 0.0),
Double64(0.0980171403295606, -1.634582362244256e-18),
Double64(0.19509032201612828, -7.99107906846173e-18),
Double64(0.2902846772544624, -1.892797870777425e-17),
Double64(0.3826834323650898, -1.0050772696461588e-17),
Double64(0.47139673682599764, 6.5166781360690145e-18),
Double64(0.5555702330196022, 4.709410940561677e-17),
Double64(0.6343932841636455, 1.0420901929280035e-17),
Double64(0.7071067811865476, -4.833646656726457e-17),
Double64(0.773010453362737, -3.256590703364977e-17),
Double64(0.8314696123025452, 1.4073856984728037e-18),
Double64(0.881921264348355, -1.9843248405890568e-17),
Double64(0.9238795325112867, 1.7645047084336677e-17),
Double64(0.9569403357322088, 4.05538698618757e-17),
Double64(0.9807852804032304, 1.8546939997825006e-17),
Double64(0.9951847266721969, -4.2486913678304403e-17),
Double64(1.0, -1.1210807766491662e-66),
Double64(0.9951847266721969, -4.248691367830441e-17),
Double64(0.9807852804032304, 1.8546939997825003e-17),
Double64(0.9569403357322088, 4.05538698618757e-17),
Double64(0.9238795325112867, 1.7645047084336674e-17),
Double64(0.881921264348355, -1.9843248405890562e-17),
Double64(0.8314696123025452, 1.4073856984728124e-18),
Double64(0.773010453362737, -3.256590703364977e-17),
Double64(0.7071067811865476, -4.833646656726456e-17),
Double64(0.6343932841636455, 1.0420901929280033e-17),
Double64(0.5555702330196022, 4.7094109405616756e-17),
Double64(0.47139673682599764, 6.516678136069009e-18),
Double64(0.3826834323650898, -1.0050772696461588e-17),
Double64(0.2902846772544624, -1.892797870777425e-17),
Double64(0.19509032201612828, -7.991079068461768e-18),
Double64(0.0980171403295606, -1.6345823622442887e-18),
Double64(-2.9947698097183397e-33, 1.1124542208633653e-49),
Double64(-0.0980171403295606, 1.6345823622442521e-18),
Double64(-0.19509032201612828, 7.991079068461727e-18),
Double64(-0.2902846772544624, 1.892797870777424e-17),
Double64(-0.3826834323650898, 1.0050772696461581e-17),
Double64(-0.47139673682599764, -6.5166781360690145e-18),
Double64(-0.5555702330196022, -4.709410940561676e-17),
Double64(-0.6343932841636455, -1.0420901929280028e-17),
Double64(-0.7071067811865476, 4.833646656726456e-17),
Double64(-0.773010453362737, 3.2565907033649785e-17),
Double64(-0.8314696123025452, -1.407385698472802e-18),
Double64(-0.881921264348355, 1.9843248405890578e-17),
Double64(-0.9238795325112867, -1.7645047084336662e-17),
Double64(-0.9569403357322088, -4.05538698618757e-17),
Double64(-0.9807852804032304, -1.8546939997825006e-17),
Double64(-0.9951847266721969, 4.248691367830441e-17),
Double64(-1.0, 3.068418716632831e-65),
Double64(-0.9951847266721969, 4.248691367830441e-17),
Double64(-0.9807852804032304, -1.8546939997825003e-17),
Double64(-0.9569403357322088, -4.0553869861875694e-17),
Double64(-0.9238795325112867, -1.7645047084336668e-17),
Double64(-0.881921264348355, 1.9843248405890565e-17),
Double64(-0.8314696123025452, -1.407385698472797e-18),
Double64(-0.773010453362737, 3.256590703364977e-17),
Double64(-0.7071067811865476, 4.833646656726457e-17),
Double64(-0.6343932841636455, -1.0420901929280039e-17),
Double64(-0.5555702330196022, -4.709410940561676e-17),
Double64(-0.47139673682599764, -6.516678136069028e-18),
Double64(-0.3826834323650898, 1.0050772696461658e-17),
Double64(-0.2902846772544624, 1.8927978707774248e-17),
Double64(-0.19509032201612828, 7.991079068461796e-18),
Double64(-0.0980171403295606, 1.6345823622442672e-18),
#Double64(5.989539619436679e-33, -2.2249084417267306e-49)
);

const n_sin_npio32 = length(sin_npio32)

# cos_npio32[idx] = cos( (idx-1)*pi / 32 )

const cos_npio32 = (
Double64(1.0, 0.0),
Double64(0.9951847266721969, -4.248691367830441e-17),
Double64(0.9807852804032304, 1.8546939997825006e-17),
Double64(0.9569403357322088, 4.05538698618757e-17),
Double64(0.9238795325112867, 1.7645047084336677e-17),
Double64(0.881921264348355, -1.9843248405890562e-17),
Double64(0.8314696123025452, 1.407385698472803e-18),
Double64(0.773010453362737, -3.256590703364977e-17),
Double64(0.7071067811865476, -4.833646656726457e-17),
Double64(0.6343932841636455, 1.0420901929280033e-17),
Double64(0.5555702330196022, 4.709410940561676e-17),
Double64(0.47139673682599764, 6.5166781360690214e-18),
Double64(0.3826834323650898, -1.0050772696461586e-17),
Double64(0.2902846772544624, -1.8927978707774255e-17),
Double64(0.19509032201612828, -7.99107906846173e-18),
Double64(0.0980171403295606, -1.634582362244275e-18),
Double64(-1.4973849048591698e-33, 5.562271104316826e-50),
Double64(-0.0980171403295606, 1.6345823622442537e-18),
Double64(-0.19509032201612828, 7.991079068461728e-18),
Double64(-0.2902846772544624, 1.892797870777425e-17),
Double64(-0.3826834323650898, 1.0050772696461583e-17),
Double64(-0.47139673682599764, -6.516678136069013e-18),
Double64(-0.5555702330196022, -4.709410940561675e-17),
Double64(-0.6343932841636455, -1.0420901929280036e-17),
Double64(-0.7071067811865476, 4.833646656726457e-17),
Double64(-0.773010453362737, 3.2565907033649766e-17),
Double64(-0.8314696123025452, -1.407385698472808e-18),
Double64(-0.881921264348355, 1.984324840589056e-17),
Double64(-0.9238795325112867, -1.7645047084336677e-17),
Double64(-0.9569403357322088, -4.05538698618757e-17),
Double64(-0.9807852804032304, -1.8546939997825012e-17),
Double64(-0.9951847266721969, 4.2486913678304403e-17),
Double64(-1.0, 4.484323106596665e-66),
Double64(-0.9951847266721969, 4.248691367830441e-17),
Double64(-0.9807852804032304, -1.8546939997825003e-17),
Double64(-0.9569403357322088, -4.0553869861875694e-17),
Double64(-0.9238795325112867, -1.7645047084336674e-17),
Double64(-0.881921264348355, 1.9843248405890562e-17),
Double64(-0.8314696123025452, -1.4073856984728047e-18),
Double64(-0.773010453362737, 3.2565907033649766e-17),
Double64(-0.7071067811865476, 4.8336466567264573e-17),
Double64(-0.6343932841636455, -1.042090192928005e-17),
Double64(-0.5555702330196022, -4.709410940561677e-17),
Double64(-0.47139673682599764, -6.51667813606904e-18),
Double64(-0.3826834323650898, 1.0050772696461555e-17),
Double64(-0.2902846772544624, 1.8927978707774258e-17),
Double64(-0.19509032201612828, 7.991079068461733e-18),
Double64(-0.0980171403295606, 1.6345823622442535e-18),
Double64(-7.8337969295008e-33, 5.173596326540973e-49),
Double64(0.0980171403295606, -1.63458236224422e-18),
Double64(0.19509032201612828, -7.991079068461725e-18),
Double64(0.2902846772544624, -1.8927978707774227e-17),
Double64(0.3826834323650898, -1.0050772696461569e-17),
Double64(0.47139673682599764, 6.516678136069015e-18),
Double64(0.5555702330196022, 4.7094109405616774e-17),
Double64(0.6343932841636455, 1.0420901929280038e-17),
Double64(0.7071067811865476, -4.833646656726457e-17),
Double64(0.773010453362737, -3.256590703364977e-17),
Double64(0.8314696123025452, 1.4073856984728055e-18),
Double64(0.881921264348355, -1.984324840589057e-17),
Double64(0.9238795325112867, 1.7645047084336705e-17),
Double64(0.9569403357322088, 4.05538698618757e-17),
Double64(0.9807852804032304, 1.8546939997825018e-17),
Double64(0.9951847266721969, -4.248691367830441e-17),
#Double64(1.0, -1.7937292426403932e-65)
);

const n_cos_npio32 = length(cos_npio32)

const tan_npio32  = (
    Double64(0.0, 0.0),
    Double64(0.09849140335716425, 5.3100671162822435e-18),
    Double64(0.198912367379658, 8.391794477636538e-19),
    Double64(0.3033466836073424, -1.2766990847826405e-17),
    Double64(0.41421356237309503, 1.4349369327986523e-17),
    Double64(0.5345111359507917, -5.2626469389262167e-17),
    Double64(0.6681786379192989, 4.1042270233610004e-17),
    Double64(0.8206787908286604, -2.2923147594675738e-17),
    Double64(1.0, 0.0),
    Double64(1.2185035255879764, -2.1244788699148428e-17),
    Double64(1.496605762665489, 6.974100888958305e-17),
    Double64(1.8708684117893895, 2.945807107795785e-17),
    Double64(2.414213562373095, 1.2537167179050217e-16),
    Double64(3.2965582089383205, -9.53295779936916e-17),
    Double64(5.027339492125848, 2.95379181037367e-17),
    Double64(10.15317038760886, 5.362306887894472e-16),
    Double64(1.7976931348623157e308, 0.0),
    Double64(-10.15317038760886, -5.362306887894472e-16),
    Double64(-5.027339492125848, -2.95379181037367e-17),
    Double64(-3.2965582089383205, 9.53295779936916e-17),
    Double64(-2.414213562373095, -1.2537167179050217e-16),
    Double64(-1.8708684117893895, -2.945807107795785e-17),
    Double64(-1.496605762665489, -6.974100888958305e-17),
    Double64(-1.2185035255879764, 2.1244788699148428e-17),
    Double64(-1.0, 0.0),
    Double64(-0.8206787908286604, 2.2923147594675738e-17),
    Double64(-0.6681786379192989, -4.1042270233610004e-17),
    Double64(-0.5345111359507917, 5.2626469389262167e-17),
    Double64(-0.41421356237309503, -1.4349369327986523e-17),
    Double64(-0.3033466836073424, 1.2766990847826405e-17),
    Double64(-0.198912367379658, -8.391794477636538e-19),
    Double64(-0.09849140335716425, -5.3100671162822435e-18),
    Double64(0.0, 0.0),
    Double64(0.09849140335716425, 5.3100671162822435e-18),
    Double64(0.198912367379658, 8.391794477636538e-19),
    Double64(0.3033466836073424, -1.2766990847826405e-17),
    Double64(0.41421356237309503, 1.4349369327986523e-17),
    Double64(0.5345111359507917, -5.2626469389262167e-17),
    Double64(0.6681786379192989, 4.1042270233610004e-17),
    Double64(0.8206787908286604, -2.2923147594675738e-17),
    Double64(1.0, 0.0),
    Double64(1.2185035255879764, -2.1244788699148428e-17),
    Double64(1.496605762665489, 6.974100888958305e-17),
    Double64(1.8708684117893895, 2.945807107795785e-17),
    Double64(2.414213562373095, 1.2537167179050217e-16),
    Double64(3.2965582089383205, -9.53295779936916e-17),
    Double64(5.027339492125848, 2.95379181037367e-17),
    Double64(10.15317038760886, 5.362306887894472e-16),
    Double64(-1.7976931348623157e308, 0.0),
    Double64(-10.15317038760886, -5.362306887894472e-16),
    Double64(-5.027339492125848, -2.95379181037367e-17),
    Double64(-3.2965582089383205, 9.53295779936916e-17),
    Double64(-2.414213562373095, -1.2537167179050217e-16),
    Double64(-1.8708684117893895, -2.945807107795785e-17),
    Double64(-1.496605762665489, -6.974100888958305e-17),
    Double64(-1.2185035255879764, 2.1244788699148428e-17),
    Double64(-1.0, 0.0),
    Double64(-0.8206787908286604, 2.2923147594675738e-17),
    Double64(-0.6681786379192989, -4.1042270233610004e-17),
    Double64(-0.5345111359507917, 5.2626469389262167e-17),
    Double64(-0.41421356237309503, -1.4349369327986523e-17),
    Double64(-0.3033466836073424, 1.2766990847826405e-17),
    Double64(-0.198912367379658, -8.391794477636538e-19),
    Double64(-0.09849140335716425, -5.3100671162822435e-18),
    Double64(0.0, 0.0),
#    Double64(0.09849140335716425, 5.3100671162822435e-18)
);

const n_tan_npio32 = length(tan_npio32)

const tan_coeff = const TangentCoeffs = [
	Double64(1.0, 0.0),
	Double64(0.3333333333333333, 1.850371707708594e-17),
	Double64(0.13333333333333333, 1.8503717077085942e-18),
	Double64(0.05396825396825397, -2.5552752154071065e-18),
	Double64(0.021869488536155203, -1.7377829530067485e-19),
	Double64(0.008863235529902197, -7.63300580171831e-19),
	Double64(0.003592128036572481, -1.253823608406629e-19),
	Double64(0.0014558343870513183, -6.214492640136062e-20),
	Double64(0.000590027440945586, 3.478690842383652e-20),
	Double64(0.00023912911424355248, 3.564613898329782e-21),
	Double64(9.691537956929451e-5, -6.2386628755632464e-21),
	Double64(3.927832388331683e-5, 1.3737015743076767e-21),
	Double64(1.5918905069328964e-5, 1.0427554807190543e-21),
	Double64(6.451689215655431e-6, 1.1519922496640058e-22),
	Double64(2.6147711512907546e-6, -9.313685621299801e-23),
	Double64(1.0597268320104654e-6, 2.3670525505213632e-24),
	Double64(4.294911078273806e-7, 1.1643520863702653e-23),
	Double64(1.7406618963571648e-7, -2.4364235361756183e-24),
	Double64(7.054636946400968e-8, 1.81294121906281e-24),
	Double64(2.859136662305254e-8, -6.927061232480201e-25),
	Double64(1.1587644432798853e-8, -6.676141550901535e-25),
	Double64(4.6962953982309016e-9, 5.585087708922479e-26),
	Double64(1.903336833931276e-9, 6.644200217899965e-26),
	Double64(7.713933635359062e-10, 3.330462514587965e-26),
	Double64(3.126339545892087e-10, 1.3431720632001545e-26),
	Double64(1.26705769303054e-10, 3.212977078832784e-27),
	Double64(5.135191408039368e-11, -2.0245650830490024e-27),
	Double64(2.0812146867700473e-11, 8.599281894073044e-28),
	Double64(8.434845419094337e-12, 8.050467205827787e-28),
	Double64(3.4185140868111557e-12, 8.393525824301979e-29),
	Double64(1.385471574294847e-12, -5.907625471272317e-29),
	Double64(5.61510479241468e-13, -1.0849557288985563e-29),
	Double64(2.2757162553728747e-13, 1.4248250730723616e-29),
	Double64(9.223130585139532e-14, -5.829477138214029e-30),
	Double64(3.737994031096739e-14, -1.2136039633979156e-30),
	Double64(1.5149519187148606e-14, -6.234579855436469e-31),
	Double64(6.139868862616814e-15, 1.660718962017884e-31),
	Double64(2.488395122276279e-15, -1.6465031184823015e-31),
	Double64(1.0085085566354096e-15, 6.051385314404599e-32),
	Double64(4.087331226869014e-16, -3.840887770358501e-33),
	Double64(1.6565329513786285e-16, -1.1753361684566294e-32),
	Double64(6.713675175048712e-17, -5.751617343643472e-34),
	Double64(2.7209500613044588e-17, 5.912656012687598e-35),
	Double64(1.102759523372237e-17, 6.889387842184814e-34),
	Double64(4.4693160072374125e-18, -5.963106810322789e-35),
	Double64(1.8113455516997665e-18, 1.6975481764291697e-34),
	Double64(7.3411070113401475e-19, 6.232644202349733e-36),
	Double64(2.975238606536194e-19, 5.973993675861387e-36),
	Double64(1.2058187889305074e-19, -1.1146629551117048e-36),
	Double64(4.8869994781043706e-20, -8.235745703221855e-37),
	Double64(1.980626286324087e-20, 1.0151251301137887e-37),
	Double64(8.027175987339371e-21, 5.155945578317267e-37),
	Double64(3.2532918893703063e-21, -1.8492173686149099e-37),
	Double64(1.3185095398600672e-21, 2.7232745038615834e-38),
	Double64(5.343717888893445e-22, 3.047260525016857e-38),
	Double64(2.165727286213673e-22, 1.4832290335239522e-38),
	Double64(8.777362083427096e-23, 2.1379234735632897e-39),
	Double64(3.557330862200837e-23, 1.700525479402311e-39),
	Double64(1.4417318942624272e-23, -4.05166348193838e-40),
	Double64(5.843119280863155e-24, 1.2167640006197093e-40),
	Double64(2.3681270468016817e-24, -1.719508829736334e-40),
	Double64(9.597657415895207e-25, 8.563087764073255e-41),
	Double64(3.889784038288653e-25, 1.1167780858140198e-42),
	Double64(1.5764700914899153e-25, 4.697215882796148e-42)
];

const n_tan_coeff = length(tan_coeff)

const twokminus1 = [ 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29,
    31, 33, 35, 37, 39, 41, 43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65,
    67, 69, 71, 73, 75, 77, 79, 81, 83, 85, 87, 89, 91, 93, 95, 97, 99, 101,
    103, 105, 107, 109, 111, 113, 115, 117, 119, 121, 123, 125, 127];

# x^-1, x^1, x^3 ...

# bernoulli_even[i] = bernoulli(2*i-2),bernoulli(e) =bernoulli_even[(e+2)>>1]
const bernoulli_even =
[ Double64(1.0, 0.0),                                       # 1:0th
  Double64(0.16666666666666666, 9.25185853854297e-18),      # 2:2nd
  Double64(-0.03333333333333333, -4.625929269271486e-19),   # 3:4th
  Double64(0.023809523809523808, 1.32169407693471e-18),     # 4:6th
  Double64(-0.03333333333333333, -4.625929269271486e-19),
  Double64(0.07575757575757576, -2.10269512239613e-18),
  Double64(-0.2531135531135531, -1.1061562736192037e-17),
  Double64(1.1666666666666667, -7.401486830834377e-17),
  Double64(-7.092156862745098, -3.274069468698501e-16),
  Double64(54.971177944862156, -1.9588897477095493e-16),
  Double64(-529.1242424242424, 6.890111377067638e-16),
  Double64(6192.123188405797, 9.226757844073186e-14),
  Double64(-86580.25311355312, 3.5926706461242705e-12),
  Double64(1.4255171666666667e6, -7.761021455128987e-11),
  Double64(-2.7298231067816094e7, 1.610010519795034e-9),
  Double64(6.015808739006424e8, -2.6635227381825164e-8),
  Double64(-1.5116315767092157e10, 5.011465035232843e-7),
  Double64(4.296146430611667e11, -2.0345052083333332e-5),
  Double64(-1.3711655205088332e13, -0.0007409090879485616),
  Double64(4.883323189735932e14, -0.020833333333333332),
  Double64(-1.9296579341940068e16, -0.14863266814486326),
  Double64(8.416930475736826e17, 55.000553709856035),
  Double64(-4.0338071854059454e19, -949.0768115942029),
  Double64(2.1150748638081993e21, -101647.85460992908),
  Double64(-1.2086626522296526e23, 2.6786766880629174e6),
  Double64(7.500866746076964e24, 2.0081412007575756e8)
];

const n_bernoulli_even = length(bernoulli_even)

#=
bernoulli_evenf
[(1.0, 0.0, 0.0),
   (0.16666666666666666, 9.25185853854297e-18, 5.135813185032629e-34),
   (-0.03333333333333333, -4.625929269271486e-19, -6.419766481290786e-36),
   (0.023809523809523808, 1.32169407693471e-18, 7.336875978618041e-35),
   (-0.03333333333333333, -4.625929269271486e-19, -6.419766481290786e-36),
   (0.07575757575757576, -2.10269512239613e-18, 5.836151346627988e-35),
   (-0.2531135531135531, -1.1061562736192037e-17, 7.607776014751631e-34),
   (1.1666666666666667, -7.401486830834377e-17, -4.108650548026103e-33),
   (-7.092156862745098, -3.274069468698501e-16, 2.0108219152692458e-32),
   (54.971177944862156, -1.9588897477095493e-16, -1.0874022503046528e-32),
   (-529.1242424242424, 6.890111377067638e-16, 5.976218978947059e-34),
   (6192.123188405797, 9.226757844073186e-14, -5.853576572074754e-30),
   (-86580.25311355312, 3.5926706461242705e-12, -1.609669232303282e-28),
   (1.4255171666666667e6, -7.761021455128987e-11, -4.308232357047019e-27),
   (-2.7298231067816094e7, 1.610010519795034e-9, 2.234338436206454e-26),
   (6.015808739006424e8, -2.6635227381825164e-8, -2.310237711361619e-26),
   (-1.5116315767092157e10, 5.011465035232843e-7, -1.9514974765172607e-23),
   (4.296146430611667e11, -2.0345052083333332e-5, -1.1293772630057337e-21),
   (-1.3711655205088332e13, -0.0007409090879485616, -1.4575906467399005e-20),
   (4.883323189735932e14, -0.020833333333333332, -1.1564823173178713e-18),
   (-1.9296579341940068e16, -0.14863266814486326, -1.0970940014219026e-17),
   (8.416930475736826e17, 55.000553709856035, 1.258990450959203e-16),
   (-4.0338071854059454e19, -949.0768115942029, 2.405547580776223e-14),
   (2.1150748638081993e21, -101647.85460992908, 5.366663772163661e-12),
   (-1.2086626522296526e23, 2.6786766880629174e6, 3.3552496116291315e-11),
   (7.500866746076964e24, 2.0081412007575756e8, 1.2643409497810133e-8)]

tan_coeffs = [
  1//1,
  1//3,
  2//15,
  17//315,
  62//2835,
  1382//155925,
  21844//6081075,
  929569//638512875,
  6404582//10854718875,
  443861162//1856156927625,
  18888466084//194896477400625,
  113927491862//2900518163668125,
  58870668456604//3698160658676859375,
  8374643517010684//1298054391195577640625,
  689005380505609448//263505041412702261046875,
  129848163681107301953//122529844256906551386796875,
  1736640792209901647222//4043484860477916195764296875
];

x^1, x^3 ...

 Double64(1.0, 0.0)
 Double64(0.3333333333333333, 1.850371707708594e-17)
 Double64(0.13333333333333333, 1.8503717077085942e-18)
 Double64(0.05396825396825397, -2.5552752154071065e-18)
 Double64(0.021869488536155203, -1.7377829530067485e-19)
 Double64(0.008863235529902197, -7.63300580171831e-19)
 Double64(0.003592128036572481, -1.253823608406629e-19)
 Double64(0.0014558343870513183, -6.214492640136062e-20)
 Double64(0.000590027440945586, 3.478690842383652e-20)
 Double64(0.00023912911424355248, 3.564613898329782e-21)
 Double64(9.691537956929451e-5, -6.2386628755632464e-21)
 Double64(3.927832388331683e-5, 1.3737015743076767e-21)
 Double64(1.5918905069328964e-5, 1.0427554807190543e-21)
 Double64(6.451689215655431e-6, 1.1519922496640058e-22)
 Double64(2.6147711512907546e-6, -9.313685621299801e-23)
 Double64(1.0597268320104654e-6, 2.3670525505213632e-24)
 Double64(4.294911078273806e-7, 1.1643520863702653e-23)

cot_coeffs = [
  1,
  -1//3,
  -1//45,
  -2//945,
  -1//4725,
  -2//93555,
  -1382//638512875,
  -4//18243225,
  -3617//162820783125,
  -87734//38979295480125,
  -349222//1531329465290625,
  -310732//13447856940643125,
  -472728182//201919571963756521875,
  -2631724//11094481976030578125,
  -13571120588//564653660170076273671875,
  -13785346041608//5660878804669082674070015625,
  -7709321041217//31245110285511170603633203125,
  -303257395102//12130454581433748587292890625
]

x^-1, x^1, x^3 ...

  Double64(1.0, 0.0)
 Double64(-0.3333333333333333, -1.850371707708594e-17)
 Double64(-0.022222222222222223, 8.480870326997723e-19)
 Double64(-0.0021164021164021165, 1.4272460344676906e-19)
 Double64(-0.00021164021164021165, 8.851449482249383e-21)
 Double64(-2.1377799155576935e-5, 1.2363213405319781e-21)
 Double64(-2.1644042808063972e-6, 1.4406027968224744e-23)
 Double64(-2.1925947851873778e-7, 4.633975226217266e-25)
 Double64(-2.2214608789979678e-8, -9.8218334027542e-25)
 Double64(-2.2507846516808994e-9, 1.586632784299882e-25)
 Double64(-2.2805151204592183e-10, 5.2026333279551446e-27)
 Double64(-2.3106432599002624e-11, 6.28420623757996e-29)
 Double64(-2.3411706819824882e-12, -1.7365409250001553e-28)
 Double64(-2.3721017400233653e-13, -1.2711970013895883e-29)
 Double64(-2.4034415333307705e-14, -1.5378003549195334e-30)
 Double64(-2.4351954029183367e-15, -1.4005306880921365e-31)
 Double64(-2.4673688045172075e-16, 7.143574177377665e-34)
 Double64(-2.499967277122081e-17, 9.643517068634457e-34)

=#


text:
using ProgressMeter
using .Threads
using Lasso
using NPZ
using Statistics
using PyCall
using JLD
using SparseArrays

py"""
import scipy.sparse as sp

def spload(x):
    return sp.load_npz(x)

def spsave(x, i):
    return sp.save_npz(f"Results/nodal_results_{i}.npz", x)
"""

PyObject(S::SparseMatrixCSC) =
    pyimport("scipy.sparse").csc_matrix((S.nzval, S.rowval .- 1, S.colptr .- 1), shape=size(S))

# G = py"spload"("Data/Gcsc.npz")
# G = SparseMatrixCSC(G.shape[1], G.shape[2], G.indptr .+ 1, G.indices .+ 1, G.data)

G = npzread("Data/nodal_G.npy")

struct InvResult{T,S}
    i::Int
    λmin::T
    dn::T
    m::S
end

function fit_model(G,dv,i)
    dn = std(dv)
    dv ./= dn
    path = fit(LassoPath, G, dv, λminratio = 1e-3, standardize=false, intercept=false, maxncoef=size(G,2))
    λmin = path.λ[argmin(aicc(path))]
    m = coef(path; select=MinAICc())
    return InvResult(i, λmin, dn, m)
end

d = npzread("Data/nodal_wvtdata.npy")
p = Progress(size(d,1) * size(d,2), dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)

results = Array{Any,2}(undef, (size(d,1), size(d,2)))

for j = 1:2
    @threads for i = 1:size(d,2)
        results[j,i] = fit_model(G,d[j,i,:],i)
        next!(p)
    end
end

sres1 = reduce(hcat, [x.m.*x.dn for x in results[1,:]])
py"spsave"(PyObject(sres1),1)

sres2 = reduce(hcat, [x.m.*x.dn for x in results[2,:]])
py"spsave"(PyObject(sres2),2)

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
using Printf

const types = (Complex, Float64, Rational, Int, Bool)

for Tb in types, Te in types
    zb, ze = zero(Tb), zero(Te)
    r = zb ^ ze
    @printf("%10s ^ %-10s = %7s ^ %-7s = %-12s (%s)\n", Tb, Te, zb, ze, r, typeof(r))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using LinearAlgebra
include("shared.jl")

const GRIDN = 24

function checkcovmat(fn, p, v; display_small_piece=false)
  grid1d = range(0.0, 1.0, length=GRIDN)
  pts    = vec(collect.(Iterators.product(grid1d, grid1d)))
  M      = assemble_matrix(fn, pts, [1.0, p, v])
  em     = eigmin(M)
  if display_small_piece # just a debugging thing, really.
    println("Principle 5x5 minor:")
    display(M[1:5, 1:5])
  end
  Mf = cholesky!(M, check=false)
  (issuccess(Mf), em)
end
translate_result(res) = res ? :SUCCESS : :FAILURE

const VRANGE = (0.4, 0.55, 0.755, 0.99, 1.01, 1.55, 2.05, 3.05, 4.05)
const PRANGE = (0.01, 0.1, 1.0, 10.0, 100.0, 1000.0)

for (v, p) in Iterators.product(VRANGE, PRANGE)
  println("\n(v,p) = ($v, $p):")
  print("AMOS:")
  amos_succ = :PLACEHOLDER
  us_succ   = :PLACEHOLDER
  (amos_em, us_em) = (0.0, 0.0)
  try
    (amos_succ, amos_em) = checkcovmat(matern_amos, p, v)
    println("$(translate_result(amos_succ)), $amos_em")
  catch
    amos_succ = :FAILURE
    println(amos_succ)
  end

  print("US:  ")
  try
    (us_succ, us_em) = checkcovmat(matern_us, p, v)
    println("$(translate_result(us_succ)), $us_em")
  catch
    us_succ = :FAILURE
    println(us_succ)
  end
  println("Eigmin difference: $(abs(amos_em-us_em))")
  println("Eigmin rtol:       $(abs(amos_em-us_em)/abs(amos_em))")

  if amos_succ && !us_succ
    println("######################")
    println("!!!!WE FAILED WHERE AMOS SUCCEEDED, INVESTIGATE!!!")
    println("######################")
  end
end


text:
module Scientific_Programming_in_Julia

include("ReverseDiff.jl")
using .ReverseDiff
export track, accum!, σ
export TrackedArray, TrackedMatrix, TrackedVector

include("Ecosystem.jl")

using .Ecosystem
export Grass, Sheep, Wolf, World, Mushroom
export agent_step!, agent_count, world_step!, every_nth

end


text:
#import Base: +
"""
    AbstractConstrains{T}
Abstract type for all possible constructs
"""
abstract type AbstractConstrains{T} end

"""
    Element_Constrains{T}
Abstract Constrains with effect on the element that are allow in the constructs
"""
abstract type Element_Constrains{T} <: AbstractConstrains{T} end

"""
    No_Constrain{T}
If no constrain is given
"""

struct No_Constrain{T} <: Element_Constrains{T}
    con :: T
end

"""
    Construct_Constrains{T}

Constrains that are applicable on a generated construct.
These can only be used in a filter with isn't very effiecent for large design space.
"""
abstract type Construct_Constrains{T} <: Element_Constrains{T} end

"""
    Single_Construct_Constrains{T}

abstract type for one single constrain
"""
abstract type Single_Construct_Constrains{T} <: Construct_Constrains{T} end

"""
    UnOrdered_Constrain{T<:Mod}

an unorder Constrain is a group of modules that can not co-occur in one construct.
These group is given as:  ``Array{T,1} where T <: Mod``
This constrain is unorder because the location of these models in the constructed is not considered

for more information see [`filter_constrain`]@ref

"""
struct UnOrdered_Constrain{T<:Mod} <: Single_Construct_Constrains{T}
    combination::Array{T}
end

"""
    Ordered_Constrain{T<:Mod}

Order Constrains a group of modules that can not co-occur in one construct if these modules are in a specific position.
Order because the location of these models in the constructed is used in the evaluation of the constrain.

These modules are given as ``Array{T} where T <: Mod``
The corresponding positions are ``Array{::Int}``, these are the forbidden indexes for the corresponding module.

for more information see [`filter_constrain`]@ref
"""
struct Ordered_Constrain{T<:Mod} <: Single_Construct_Constrains{T}
    pos::Array{Int}
    combination::Array{T}
end

"""
Compose_Construct_Constrains{T} <: Construct_Constrains{T}

Constrain that are applicable on a generated construct.
These can only be used in a filter with isn't very effiecent for large design space.
"""
struct Compose_Construct_Constrains{T <: Single_Construct_Constrains} <: Construct_Constrains{T}
    construct_con :: Array{T}
end

Base.eltype(K::Compose_Construct_Constrains{T}) where {T} = T

# concentration of multiple constrains can be done using + sign

"""
    +(con1::Single_Construct_Constrains,con2::Single_Construct_Constrains)

Returns a `Compose_Construct_Constrains` containing both constrains `con1` and `con2`.
"""

Base.:+(con1::Single_Construct_Constrains,con2::Single_Construct_Constrains) = Compose_Construct_Constrains([con1,con2])

"""
    Base.:+(con1::Compose_Construct_Constrains{T} where T, con2::N where N <: Single_Construct_Constrains)

Returns a `Compose_Construct_Constrains` containing where `con2` is add to the `con1`
"""

Base.:+(con1::Compose_Construct_Constrains{T} where T, con2::N where N <: Single_Construct_Constrains) = Compose_Construct_Constrains([con1.construct_con;con2])

"""
    Base.:+(con1::Single_Construct_Constrains,con2::Compose_Construct_Constrains)

Returns a `Compose_Construct_Constrains` containing where `con1` is add to the `con2`
"""
Base.:+(con1::Single_Construct_Constrains,con2::Compose_Construct_Constrains) = +(con2,con1)

"""
    Base.:+(con1::Compose_Construct_Constrains{T} where T, con2::Compose_Construct_Constrains{N} where N)

Returns a `Compose_Construct_Constrains` concatenation of `con1` and `con2`
"""
Base.:+(con1::Compose_Construct_Constrains{T} where T, con2::Compose_Construct_Constrains{N} where N) = Compose_Construct_Constrains([con1.construct_con;con2.construct_con;])

##################################
# Code below has currently no application in the package, this will probably not change shortly
# The idea is that for some constraints, we can prevent that unwanted constructs are generated, instead of generational filtering afterwards.
# To allow this a special type was constructed and see how it good be implemented
# No need to check,
#################################

#=
"""
space contrains has effect on how the space is repesanted without real evaluation of the constructs
"""

abstract type Space_Constrains{T} <: Element_Constrains{T} end

"""
Single possition constrain. If allow = 1 it indicaces that only these construcst are allow.
if allow = 0is false thes models aren't allowed on this possition, for non position specific constrain pos = 0
"""

struct Compose_Space_Constrain{Ts <: Space_Constrains{T} where T} <: Space_Constrains{Ts}
    space_con::AbstractArray{Ts}
end

abstract type Single_Space_Constrains{T} <: Space_Constrains{T} end

struct Possition_Constrain{T} <: Single_Space_Constrains{T}
    pos::Int
    mod::T
end

# add two  space Constrains
Base.:+(con1::N  ,con2::N ) where N <: Single_Space_Constrains{T} where T = Compose_Space_Constrain([con1,con2])
# add to compose constrainse
Base.:+(con1::Compose_Space_Constrain,con2::Single_Space_Constrains) = Compose_Space_Constrain([con1.space_con;con2])
Base.:+(con1::Single_Space_Constrains,con2::Compose_Space_Constrain) = Compose_Space_Constrain([con2.space_con;con1])

Base.:+(con1::Compose_Space_Constrain,con2::Compose_Space_Constrain)  = Compose_space_Constrains([con1.space_con;con2.space_con])

Base.:+(con1::Element_Constrains,con2::Element_Constrains...) = +(con1,+(con2...))
"""
Group all elementconstrains
"""

struct Multipe_constrain{T<:Any, Ts<:Space_Constrains, Tc<:Construct_Constrains}  <: Element_Constrains{T}
    space_con:: Ts
    construct_con:: Tc
    function Multipe_constrain(space_con::Space_Constrains{T}, construct_con::Construct_Constrains{V}) where {T, V}
        return new{promote_type(T, V),typeof(space_con),typeof(construct_con)}(space_con,construct_con)
    end
end

Base.:+(con1:: Space_Constrains, con2:: Construct_Constrains) = Multipe_constrain(con1,con2)
Base.:+(con1:: Construct_Constrains, con2:: Space_Constrains) = Multipe_constrain(con2,con1)

Base.:+(con1::Multipe_constrain,con2::Construct_Constrains) = Multipe_constrain(con1.space_con, +(con1.construct_con,con2))
Base.:+(con1::Multipe_constrain,con2::Space_Constrains) = Multipe_constrain(+(con1.space_con,con2),con1.construct_con)

Base.:+(con1::Construct_Constrains,con2::Multipe_constrain) = Multipe_constrain(con2,con1)
Base.:+(con1::Space_Constrains,con2::Multipe_constrain) =  Multipe_constrain(con2,con1)

Base.:+(con1::Multipe_constrain,con2::Multipe_constrain) =  Multipe_constrain(con1.space_con + con2.space_con, con1.construct_con + con2.construct_con)

# fucntion(con1::Multipe_constrain,con2::Space_Constrains)

"""
Group length and constrains to one object
"""
abstract type Full_Constrain{T} <: AbstractConstrains{T} end

"""
Gives an specific lenth to given constrains
"""

struct Len_Constrain{T} <: Full_Constrain{T}
    len::Int
    el_con::Element_Constrains{T}
end

struct Group_Len_Constrain{T} <: Full_Constrain{T}
    group::Tuple{T}
end

group_len_constrain(con1::Len_Constrain,con2::Len_Constrain) = Group_Len_Constrain(Tuple(con1,con2))
group_len_constrain(con1::Len_Constrain,con2::Len_Constrain...) = Group_Len_Constrain(Tuple(con1,Tuple(con2...)))
=#