text:
using GibbsSeaWater
using Compat
using Base.Test

testTolerances = @compat Dict(:abs_pressure_from_p => 2.300031483173370e-004,
                              :adiabatic_lapse_rate_from_ct => 5.699743845487985e-019,
                              :adiabatic_lapse_rate_from_t => 5.699738675609156e-019,
                              :adiabatic_lapse_rate_t_exact => 5.699738675609156e-019,
                              :alpha => 8.264713918662917e-015,
                              :alpha_ct_exact => 8.257102480598889e-015,
                              :alpha_ct_exact_rab => 8.257102480598889e-015,
                              :alpha_on_beta => 1.052907760978883e-11,
                              :alpha_rab => 8.264713918662917e-015,
                              :alpha_wrt_ct_t_exact => 8.257094010269417e-015,
                              :alpha_wrt_pt_t_exact => 8.599068316130290e-015,
                              :alpha_wrt_t_exact => 8.594856868316542e-015,
                              :beta => 1.846179459308317e-015,
                              :beta_const_ct_t_exact => 1.847372081698051e-015,
                              :beta_const_pt_t_exact => 1.805738718274608e-015,
                              :beta_const_t_exact => 1.804871356536619e-015,
                              :beta_ct_exact => 1.847372081698051e-015,
                              :beta_ct_exact_rab => 1.847372081698051e-015,
                              :beta_rab => 1.846179459308317e-015,
                              :brinesa_ct => 1.300080043620255e-010,
                              :brinesa_t => 1.300293206440983e-010,
                              :cabbeling => 1.634722766223964e-16,
                              :c_from_sp => 6.163816124171717e-010,
                              :chem_potential_salt_t_exact => 4.805315256817266e-007,
                              :chem_potential_t_exact => 1.317864928296331e-009,
                              :chem_potential_water_t_exact => 4.811790859093890e-007,
                              :cp_t_exact => 2.813976607285440e-009,
                              :ct_freezing => 2.257127817983928e-011,
                              :ct_from_entropy => 6.261107188265669e-010,
                              :ct_from_pt => 6.261107188265669e-010,
                              :ct_from_rho => 6.298552790440226e-010,
                              :ct_from_rho_exact => 6.280096442878858e-010,
                              :ct_from_t => 6.261107188265669e-010,
                              :ct_maxdensity => 6.688338771709823e-011,
                              :ct_maxdensity_exact => 5.839062566792563e-011,
                              :ct_pt => 2.964295475749168e-013,
                              :ct_pt_pt => 5.551115123125783e-014,
                              :ct_sa => 1.006231122729906e-012,
                              :ct_sa_pt => 1.457167719820518e-014,
                              :ct_sa_sa => 1.431146867680866e-014,
                              :deltasa_from_sp => 6.963318810448982e-013,
                              :deltasa_atlas => 6.945514042372425e-013,
                              :depth_from_z => 2.287223921371151e-008,
                              :distance => 4.470348358154297e-008,
                              :drho_dct => 8.315403920988729e-12,
                              :drho_dp => 1.782157321023048e-18,
                              :drho_dsa => 1.943112337698949e-12,
                              :dynamic_enthalpy => 2.288754734930485e-007,
                              :dynamic_enthalpy_ct_exact => 2.288797986693680e-007,
                              :dynamic_enthalpy_t_exact => 2.288943505845964e-007,
                              :enthalpy => 2.499356924090534e-006,
                              :enthalpy_ct_exact => 2.499349648132920e-006,
                              :enthalpy_diff => 1.154347728515859e-010,
                              :enthalpy_diff_ct_exact => 7.275957614183426e-011,
                              :enthalpy_t_exact => 2.499349648132920e-006,
                              :entropy_from_ct => 9.028163105995191e-009,
                              :entropy_from_pt => 9.028163105995191e-009,
                              :entropy_from_t => 9.028163105995191e-009,
                              :entropy_t_exact => 9.028163105995191e-009,
                              :eta_ct => 3.137579085432662e-011,
                              :eta_ct_ct => 2.381358998881922e-013,
                              :eta_sa => 4.653527563291959e-012,
                              :eta_sa_ct => 4.981922535098049e-014,
                              :eta_sa_sa => 6.995931611797346e-013,
                              :f => 1.000000000000000e-015,
                              :fdelta => 2.702939055302528e-014,
                              :g_ct => 2.854610769986721e-009,
                              :geo_strf_cunningham => 9.926097845891491e-007,
                              :geo_strf_dyn_height => 9.969444363377988e-007,
                              :geo_strf_dyn_height_pc => 4.210555459849275e-008,
                              :geo_strf_dyn_height_pc_p_mid => 1.000000000000000e-015,
                              :geo_strf_isopycnal => 4.297791695861974e-007,
                              :geo_strf_isopycnal_pc => 8.540336438045415e-009,
                              :geo_strf_isopycnal_pc_p_mid => 1.000000000000000e-015,
                              :geo_strf_montgomery => 9.915690188933013e-007,
                              :geo_strf_velocity => 8.024344404222727e-009,
                              :geo_strf_velocity_mid_lat => 1.000000000000000e-015,
                              :geo_strf_velocity_mid_long => 1.000000000000000e-015,
                              :grav => 5.329070518200751e-014,
                              :h_ct => 3.160494088660926e-010,
                              :h_ct_ct => 1.135647131889073e-011,
                              :helmholtz_energy_t_exact => 2.440137905068696e-007,
                              :h_p => 2.818925648462312e-016,
                              :h_sa => 2.371365326325758e-010,
                              :h_sa_ct => 2.188554892867956e-012,
                              :h_sa_sa => 7.264189250122399e-013,
                              :internal_energy => 2.499342372175306e-006,
                              :internal_energy_ct_exact => 2.499335096217692e-006,
                              :internal_energy_t_exact => 2.499335096217692e-006,
                              :ionic_strength_from_sa => 2.768674178810215e-012,
                              :ipvfn2 => 3.474816878679121e-009,
                              :isochoric_heat_cap_t_exact => 1.614353095646948e-009,
                              :isopycnal_slope_ratio => 3.781384094736495e-010,
                              :kappa => 1.717743939542931e-21,
                              :kappa_const_t_exact => 1.697064424229105e-021,
                              :kappa_t_exact => 1.712677458291044e-021,
                              :latentheat_evap_ct => 1.455657184123993e-006,
                              :latentheat_evap_t => 1.443084329366684e-006,
                              :latentheat_melting => 6.286427378654480e-008,
                              :molality_from_sa => 4.446665258228677e-012,
                              :n2 => 1.578186366313350e-014,
                              :ntpptct => 5.024869409453459e-013,
                              :osmotic_coefficient_t_exact => 3.583799923490005e-013,
                              :osmotic_pressure_t_exact => 1.023465756588848e-009,
                              :p_from_abs_pressure => 2.300066626048647e-008,
                              :p_from_z => 2.301931090187281e-008,
                              :p_mid_g_ct => 2.300021151313558e-008,
                              :p_mid_ipvfn2 => 2.300021151313558e-008,
                              :p_mid_n2 => 2.300021151313558e-008,
                              :p_mid_tursr => 2.300021151313558e-008,
                              :pot_enthalpy_from_pt => 2.499356924090534e-006,
                              :pot_rho_t_exact => 2.929709808086045e-010,
                              :pt0_from_t => 6.054037271496782e-010,
                              :pt_ct => 2.733369086627135e-013,
                              :pt_ct_ct => 4.440892098500626e-014,
                              :pt_from_ct => 6.054037271496782e-010,
                              :pt_from_entropy => 6.054072798633570e-010,
                              :pt_from_t => 6.054037271496782e-010,
                              :pt_sa => 9.670458878119348e-013,
                              :pt_sa_ct => 1.199127602768968e-014,
                              :pt_sa_sa => 2.081668171172169e-014,
                              :r_from_sp => 1.436317731418058e-011,
                              :rho => 2.945625965367071e-010,
                              :rho_ct_exact => 2.944489096989855e-010,
                              :rho_ct_exact_rab => 2.944489096989855e-010,
                              :rho_rab => 2.944489096989855e-010,
                              :rho_t_exact => 2.944489096989855e-010,
                              :rsubrho => 1.709803143512545e-008,
                              :sa_from_rho => 1.308677610722953e-010,
                              :sa_from_rho_ct_exact => 1.306119656874216e-010,
                              :sa_from_rho_t_exact => 1.304769625676272e-010,
                              :sa_from_sp => 1.300080043620255e-010,
                              :sa_from_sp_baltic => 1.300080043620255e-010,
                              :sa_from_sstar => 1.300222152167407e-010,
                              :sa_sa_sstar_from_sp => 1.300008989346679e-010,
                              :sigma0 => 2.933120413217694e-010,
                              :sigma0_ct_exact => 2.929709808086045e-010,
                              :sigma0_pt0_exact => 2.929709808086045e-010,
                              :sigma1 => 2.999058779096231e-010,
                              :sigma1_ct_exact => 2.994511305587366e-010,
                              :sigma2 => 3.060449671465904e-010,
                              :sigma2_ct_exact => 3.055902197957039e-010,
                              :sigma3 => 3.119566827081144e-010,
                              :sigma3_ct_exact => 3.117293090326712e-010,
                              :sigma4 => 3.180957719450817e-010,
                              :sigma4_ct_exact => 3.176410245941952e-010,
                              :sound_speed => 2.596152626210824e-009,
                              :sound_speed_ct_exact => 2.590240910649300e-009,
                              :sound_speed_t_exact => 2.590240910649300e-009,
                              :specvol_anom => 2.810252031082428e-016,
                              :specvol_anom_ct_exact => 2.805915222392486e-016,
                              :specvol_anom_t_exact => 2.805915222392486e-016,
                              :specvol => 2.821094052807283e-016,
                              :specvol_ct_exact => 2.818925648462312e-016,
                              :specvol_t_exact => 2.818925648462312e-016,
                              :sp_from_c => 1.297193463756230e-010,
                              :sp_from_r => 2.681299626772216e-012,
                              :sp_from_sa => 1.297113527698457e-010,
                              :sp_from_sa_baltic => 1.297113527698457e-010,
                              :sp_from_sr => 1.297122409482654e-010,
                              :sp_from_sstar => 1.297122409482654e-010,
                              :sp_salinometer => 1.297131291266851e-010,
                              :sr_from_sp => 1.303233077010191e-010,
                              :sstar_from_sa => 1.300008989346679e-010,
                              :sstar_from_sp => 1.300008989346679e-010,
                              :sstar_sa_sstar_from_sp => 1.300008989346679e-010,
                              :steric_height => 1.017674460257467e-007,
                              :t90_from_t48 => 5.997407015456702e-010,
                              :t90_from_t68 => 5.998579410970706e-010,
                              :t_freezing => 2.157829470661454e-011,
                              :t_from_ct => 6.000142604989378e-010,
                              :t_from_rho_exact => 6.032974120273593e-010,
                              :thermobaric => 4.890907320111513e-23,
                              :t_maxdensity_exact => 6.274447628129565e-011,
                              :tu => 2.190718007000214e-008,
                              :z_from_depth => 2.287223921371151e-008,
                              :z_from_p => 2.287223921371151e-008)

# Test data
sp =  35.5e0
sa = 35.7e0
sstar = 35.5e0
sr = 35.5e0
t = 15e0
ct = 20e0
pt = 15e0
p = 300e0
p_bs = 50e0
p_ref = 100e0
c = 43.6e0
lon = 260e0
long_bs = 20e0
lat = 16e0
lat_bs = 60e0
saturation_fraction = 0.5e0

function run_gsw_test(funcsym, expected, args...)
    tol = testTolerances[funcsym]
    res = eval(funcsym)(args...)
    @test abs(res-expected) < tol
end

# Tests copied from gsw_check_functions.c

# Practical Salinity, PSS-78
run_gsw_test(:sp_from_c, 35.500961780774482e0, c,t,p);
run_gsw_test(:c_from_sp, 43.598945605280484e0, sp,t,p);

# Absolute Salinity, Preformed Salinity and Conservative Temperature
run_gsw_test(:sa_from_sp, 35.671358392019094e0, sp,p,lon,lat);
run_gsw_test(:sstar_from_sp, 35.866946753006239e0, sa,p,lon,lat);
run_gsw_test(:ct_from_t,  14.930280459895560e0, sa, t, p);

# Other conversions between temperatures, salinities, entropy, pressure and height
run_gsw_test(:deltasa_from_sp, 3.96067773336028495e-3, sp,p,lon,lat);
run_gsw_test(:sr_from_sp, 35.667397714285734e0, sp);
run_gsw_test(:sp_from_sr, 35.333387933015295e0, sr);
run_gsw_test(:sp_from_sa, 35.528504019167094e0, sa,p,lon,lat);
run_gsw_test(:sstar_from_sa, 35.694648791860907e0, sa,p,lon,lat);
run_gsw_test(:sp_from_sstar, 35.334761242083573e0, sstar,p,lon,lat);
run_gsw_test(:sa_from_sstar, 35.505322027120805e0, sstar,p,lon,lat);
run_gsw_test(:pt_from_ct, 20.023899375975017e0, sa,ct);
run_gsw_test(:t_from_ct, 20.079820359223014e0, sa,ct,p);
run_gsw_test(:ct_from_pt, 14.976021403957613e0, sa,pt);
run_gsw_test(:pt0_from_t, 14.954241363902305e0, sa,t,p);
run_gsw_test(:pt_from_t, 14.969381237883740e0, sa,t,p,p_ref);
run_gsw_test(:z_from_p, -2.980161553316402e2, p,lat);
run_gsw_test(:entropy_from_t, 212.30166821093002e0, sa,t,p);
run_gsw_test(:adiabatic_lapse_rate_from_ct, 1.877941744191155e-8, sa,ct,p);

# Density and enthalpy, based on the 48-term expression for density
run_gsw_test(:rho, 1026.4562376198473e0, sa,ct,p);
run_gsw_test(:alpha, 2.62460550806784356e-4, sa,ct,p);
run_gsw_test(:beta, 7.29314455934463365e-4, sa,ct,p);
run_gsw_test(:alpha_on_beta, 0.359872958325632e0, sa,ct,p);

let
    # WIP: rho_first_derivatives expects a pointer, so need to pass an array for now
    drho_dsa = Float64[0.0]
    drho_dct = Float64[0.0]
    drho_dp = Float64[0.0]

    rho_first_derivatives(sa, ct, p, drho_dsa, drho_dct, drho_dp);
    drho_dsa_error = abs(drho_dsa[1] - 0.748609372480258e0);
    drho_dct_error = abs(drho_dct[1] + 0.269404269504765e0);
    drho_dp_error = abs(drho_dp[1] - 4.287533235942749e-7);
    @test drho_dsa_error < testTolerances[:drho_dsa]
    @test drho_dct_error < testTolerances[:drho_dct]
    @test drho_dp_error < testTolerances[:drho_dp]
end

run_gsw_test(:specvol, 9.74225654586897711e-4, sa,ct,p);
run_gsw_test(:specvol_anom, 2.90948181201264571e-6, sa,ct,p);
run_gsw_test(:sigma0, 25.165674636323047e0, sa,ct);
run_gsw_test(:sigma1, 29.434338510752923e0, sa,ct);
run_gsw_test(:sigma2, 33.609842926904093e0, sa,ct);
run_gsw_test(:sigma3, 37.695147569371784e0, sa,ct);
run_gsw_test(:sigma4, 41.693064726656303e0, sa,ct);
run_gsw_test(:sound_speed, 1527.2011773569989e0, sa,ct,p);
run_gsw_test(:kappa, 4.177024873349404e-010, sa,ct,p);
run_gsw_test(:cabbeling, 9.463053321129075e-6, sa,ct,p);
run_gsw_test(:thermobaric, 1.739078662082863e-12, sa,ct,p);
run_gsw_test(:internal_energy, 79740.482561720783e0, sa,ct,p);
run_gsw_test(:enthalpy, 82761.872939932495e0, sa,ct,p);
run_gsw_test(:dynamic_enthalpy, 2924.5137975399025e0, sa,ct,p);
ρ = rho(sa,ct,p);
run_gsw_test(:sa_from_rho, sa, ρ,ct,p);

# Water column properties, based on the 48-term expression for density
let
    nz = 3
    sa_profile = [35.5e0, 35.7e0, 35.6e0]
    ct_profile = [12.5e0, 15e0, 10e0]
    p_profile = [00e0, 50e0, 100e0]
    lat_profile = [10e0, 10e0, 10e0]

    n2 = zeros(Float64, 2)
    n2_error = zeros(Float64, 2)
    p_mid_n2 = zeros(Float64, 2)
    p_mid_n2_error = zeros(Float64, 2)

    nsquared(sa_profile, ct_profile, p_profile, lat_profile, nz, n2, p_mid_n2);
    n2_error[1] = abs(n2[1] + 0.070960392693051e-3)
    n2_error[2] = abs(n2[2] - 0.175435821615983e-3)
    p_mid_n2_error[1] = abs(p_mid_n2[1] - 25e0);
    p_mid_n2_error[2] = abs(p_mid_n2[2] - 75e0);
    @test n2_error[1] < testTolerances[:n2]
    @test n2_error[1] < testTolerances[:n2]
    @test p_mid_n2_error[1] < testTolerances[:p_mid_n2]
    @test p_mid_n2_error[2] < testTolerances[:p_mid_n2]
end

let
    sa_profile = [35.5e0, 35.7e0, 35.6e0]
    ct_profile = [12.5e0, 15e0, 10e0]
    p_profile = [00e0, 50e0, 100e0]
    nz = 3

    tu = zeros(Float64, 2)
    tu_error = zeros(Float64, 2)
    rsubrho = zeros(Float64, 2)
    rsubrho_error = zeros(Float64, 2)
    p_mid_tursr = zeros(Float64, 2)
    p_mid_tursr_error = zeros(Float64, 2)

    turner_rsubrho(sa_profile, ct_profile, p_profile, nz, tu, rsubrho, p_mid_tursr);
    tu_error[1] = abs(tu[1] + 1.187243981606485e2);
    tu_error[2] = abs(tu[2] - 0.494158257088517e2);
    rsubrho_error[1] = abs(rsubrho[1] - 3.425146897090065e0);
    rsubrho_error[2] = abs(rsubrho[2] - 12.949399443139164e0);
    p_mid_tursr_error[1] = abs(p_mid_tursr[1] - 25e0);
    p_mid_tursr_error[2] = abs(p_mid_tursr[2] - 75e0);
    @test tu_error[1] < testTolerances[:tu]
    @test tu_error[2] < testTolerances[:tu]
    @test rsubrho_error[1] < testTolerances[:rsubrho]
    @test rsubrho_error[2] < testTolerances[:rsubrho]
    @test p_mid_tursr_error[1] < testTolerances[:p_mid_tursr]
    @test p_mid_tursr_error[2] < testTolerances[:p_mid_tursr]
end

let
    sa_profile = [35.5e0, 35.7e0, 35.6e0]
    ct_profile = [12.5e0, 15e0, 10e0]
    p_profile = [00e0, 50e0, 100e0]
    nz = 3

    ipvfn2 = zeros(Float64, 2)
    ipvfn2_error = zeros(Float64, 2)
    p_mid_ipvfn2 = zeros(Float64, 2)
    p_mid_ipvfn2_error = zeros(Float64, 2)

    ipv_vs_fnsquared_ratio(sa_profile, ct_profile, p_profile, nz, ipvfn2, p_mid_ipvfn2)
    ipvfn2_error[1] = abs(ipvfn2[1] - 0.996783975249010e0);
    ipvfn2_error[2] = abs(ipvfn2[2] - 0.992112251478320e0);
    p_mid_ipvfn2_error[1] = abs(p_mid_ipvfn2[1] - 25e0);
    p_mid_ipvfn2_error[2] = abs(p_mid_ipvfn2[2] - 75e0);
    @test ipvfn2_error[1] < testTolerances[:ipvfn2]
    @test ipvfn2_error[2] < testTolerances[:ipvfn2]
    @test p_mid_ipvfn2_error[1] < testTolerances[:p_mid_ipvfn2]
    @test p_mid_ipvfn2_error[2] < testTolerances[:p_mid_ipvfn2]
end

# Freezing temperatures
run_gsw_test(:ct_freezing, -2.1801450326174852e0, sa,p,saturation_fraction);
run_gsw_test(:t_freezing, -2.1765521998023516e0, sa,p,saturation_fraction);

# Isobaric melting enthalpy and isobaric evaporation enthalpy
run_gsw_test(:latentheat_melting, 329330.54839618353e0, sa,p);
run_gsw_test(:latentheat_evap_ct, 2450871.0228523901e0, sa,ct);
run_gsw_test(:latentheat_evap_t, 2462848.2895522709e0, sa,t);

# Basic thermodynamic properties in terms of in-situ t, based on the exact Gibbs function
run_gsw_test(:rho_t_exact, 1027.7128170207150e0, sa,t,p);
run_gsw_test(:pot_rho_t_exact, 1026.8362655887486e0, sa,t,p,p_ref);
run_gsw_test(:alpha_wrt_t_exact, 2.19066952410728916e-4, sa,t,p);
run_gsw_test(:beta_const_t_exact, 7.44744841648729426e-4, sa,t,p);
run_gsw_test(:specvol_t_exact, 9.73034473676164815e-4, sa,t,p);
run_gsw_test(:sound_speed_t_exact, 1512.2053940303056e0, sa,t,p);
run_gsw_test(:kappa_t_exact, 4.25506953386609075e-010, sa,t,p);
run_gsw_test(:enthalpy_t_exact, 62520.680485510929e0, sa,t,p);
run_gsw_test(:cp_t_exact, 3982.7832563441461e0, sa,t,p);

# Library functions of the GSW toolbox
run_gsw_test(:deltasa_atlas, 3.87660373016291727e-3, p,lon,lat);
run_gsw_test(:fdelta, 1.49916256924158942e-004, p,lon,lat);
run_gsw_test(:sa_from_sp_baltic, 35.666154857142850e0, sp,long_bs,lat_bs);
run_gsw_test(:sp_from_sa_baltic, 35.533769845749660e0, sa,long_bs,lat_bs);


text:
for F in (:typemax, :typemin, :zero, :one)
  @eval begin
    @pure function $F(::Type{M}) where {M<:MarkableInteger}
        Unmarked($F(itype(M)))
    end
    @pure function $F(x::M) where {M<:MarkableInteger}
        Unmarked($F(itype(M)))
    end
  end
end

for F in (:iszero, :isone, :iseven, :isodd)
  @eval begin
    function $F(x::M) where {M<:MarkableInteger}
        xx = itype(x)
        return $F(xx)
    end
  end
end

for F in (:leading_zeros, :leading_ones, :trailing_zeros, :trailing_ones,
          :count_zeros, :count_ones, :abs, :(-))
  @eval begin
    function $F(x::M) where {M<:MarkableInteger}
        xx = itype(x) >> 1
        zz = $F(xx)
        z = Unmarked(zz)
        z |= any_lsbits(x)
        return z
    end
  end
end

for F in (:(|), :(&), :(⊻))
  @eval begin
    function $F(x::M, y::M) where {M<:MarkableInteger}
        xx = ityped(x) >> 1
        yy = ityped(y) >> 1
        zz = $F(xx, yy)
        z = zz << 1
        !(zz === z >> 1) && throw(DomainError("$x, $y"))
        z |= any_lsbits(x, y)
        return reinterpret(M, z)
    end
    $F(x::M1, y::M2) where {M1<:MarkableInteger, M2<:MarkableInteger} =
         $F(promote(x,y)...,)
    $F(x::M, y::I) where {M<:MarkableInteger, I<:Union{Signed,Unsigned}} =
         $F(promote(x,y)...,)
    $F(x::I, y::M) where {M<:MarkableInteger, I<:Union{Signed,Unsigned}} =
         $F(promote(x,y)...,)
  end
end

for F in(:flipsign, :copysign,
         :(+), :(-), :(*), :(^), :div, :cld, :fld, :rem, :mod,
         :checked_add, :checked_cld, :checked_div, :checked_fld,
         :checked_mod, :checked_mul, :checked_rem, :checked_sub)
  @eval begin
    function $F(x::M, y::M) where {M<:MarkableInteger}
        xx = ityped(x) >> 1
        yy = ityped(y) >> 1
        zz = $F(xx, yy)
        z = zz << 1
        !(zz === z >> 1) && throw(DomainError("$x, $y"))
        z |= any_lsbits(x, y)
        return reinterpret(M, z)
    end
    $F(x::M1, y::M2) where {M1<:MarkableInteger, M2<:MarkableInteger} =
         $F(promote(x,y)...,)
    $F(x::M, y::I) where {M<:MarkableInteger, I<:Union{Signed,Unsigned}} =
         $F(promote(x,y)...,)
    $F(x::I, y::M) where {M<:MarkableInteger, I<:Union{Signed,Unsigned}} =
         $F(promote(x,y)...,)
  end
end

string(x::T) where {T<:MarkableSigned} = string(ityped(x) >> 1)
string(x::T) where {T<:MarkableUnsigned} = string(ityped(x) >> 1)

show(io::IO, x::T) where {T<:MarkableInteger} = print(io, string(x))
show(x::T) where {T<:MarkableInteger} = show(StdOut, x)

ext: jl


text:
# AUTO GENERATED FILE - DO NOT EDIT

export group

"""
    group(;kwargs...)
    group(children::Any;kwargs...)
    group(children_maker::Function;kwargs...)

A Group component.
Compose elements and components in flex container. For more information, see: https://mantine.dev/core/group/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `align` (a value equal to: "stretch", "center", "flex-end", "flex-start"; optional): Defines align-items css property
- `className` (String; optional): Often used with CSS to style elements with common properties
- `direction` (a value equal to: "row", "column"; optional): Defines flex-direction property, row for horizontal, column for vertical
- `grow` (Bool; optional): Defines flex-grow property for each element, true -> 1, false -> 0
- `noWrap` (Bool; optional): Defined flex-wrap property
- `position` (a value equal to: "right", "center", "left", "apart"; optional): Defines justify-content property
- `spacing` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Space between elements
- `style` (Dict; optional): Inline style override
"""
function group(; kwargs...)
        available_props = Symbol[:children, :id, :align, :className, :direction, :grow, :noWrap, :position, :spacing, :style]
        wild_props = Symbol[]
        return Component("group", "Group", "dash_mantine_components", available_props, wild_props; kwargs...)
end

group(children::Any; kwargs...) = group(;kwargs..., children = children)
group(children_maker::Function; kwargs...) = group(children_maker(); kwargs...)


text:
@testset "abc exponential" begin
    gs_p = randn()^2

    data = Dict(
        :y => rand(Exponential(gs_p), 1000)
    )

    mod = Model(
        y=Stochastic(1, (l) -> Exponential(l), false),
        l=Stochastic(() -> LogNormal(), true)
    )

    inits = [
        Dict{Symbol,Any}(
          :y => data[:y],
          :l => rand(LogNormal()),
            )
        for i in 1:3
    ]

    f = x -> [mean(x), std(x)]
    scheme = [ABC(:l, 0.01, f, 0.01, nsim = 3, maxdraw=15)]

    setsamplers!(mod, scheme)

    sim = mcmc(mod, data, inits, 100000, burnin=25000, thin=10, chains=2)

    while maximum(gelmandiag(sim).value[1,2,:]) > 1.1
        @show maximum(gelmandiag(sim).value[1,2,:])
        sim = mcmc(sim, 10000)
    end

    m = mean(sim.value[:,1,:])
    sd = std(sim.value[:,1,:])

    @test m - sd <= gs_p <= m + sd
end

@testset "abc Normal" begin
    gs_m = randn()
    gs_s = randn()^2
    data = Dict(
        :y => rand(Normal(gs_m, gs_s), 1000)
    )

    mod = Model(
        y=Stochastic(1, (m, s) -> Normal(m, s), false),
        m=Stochastic(() -> Normal(), true),
        s=Stochastic(() -> Exponential(), true)
    )

    inits = [
        Dict{Symbol,Any}(
          :y => data[:y],
          :m => gs_m,
          :s => rand(),
        )
        for i in 1:3
    ]

    f = x -> [mean(x), std(x)]
    scheme = [ABC(:m, 0.01, f, 0.001, nsim = 3, maxdraw=15),
              ABC(:s, 0.01, f, 0.001, nsim = 3, maxdraw=15)]

    setsamplers!(mod, scheme)

    sim = mcmc(mod, data, inits, 100000, burnin=25000, thin=10, chains=2)

    while maximum(gelmandiag(sim).value[1:2,2,:]) > 1.1
        @show maximum(gelmandiag(sim).value[1:2,2,:])
        sim = mcmc(sim, 10000)
    end

    m = mean(sim.value[:,1,:])
    sd = std(sim.value[:,1,:])

    @test m - sd <= gs_s <= m + sd

    m = mean(sim.value[:,2,:])
    sd = std(sim.value[:,2,:])

    @test m - sd <= gs_m <= m + sd
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
#general validation strategy:
#   check all used subsequences individually, adding them to ValidationState if they pass testing
#   check all sequence transitions
#      iterate chronologically through sequences
#      check whether the ValidationState has an entry for the next sequence.
#      If so, proceed directly to next sequence
#      If not, keep the last min_samps-1 samples and read the first min_samps-1 samples of the next sequence.  Check that sequence
#           note that if there are not min_samps-1 samples in the NEXT sequence then we can't add it to the dictionary.  in that case:
#                keep reading more sequences until we get enough samples, do the check
#                assuming check passes, proceed to the next sequence without adding anything to the dictionary
#      Assuming that the combined sequence passes validation, add the transition to ValidationState
#

mutable struct ValidationState
    sequences::Set{String}
    transitions::Dict{String, Set{String}}
end

valid_sequences(vs::ValidationState) = vs.sequences
valid_transitions(vs::ValidationState) = vs.transitions

set_validated!(vs::ValidationState, seq_name::String) = push!(valid_sequences(vs), seq_name)
function set_validated!(vs::ValidationState, trans_from::String, trans_to::String)
    if !haskey(valid_transitions(vs), trans_from)
        valid_transitions(vs)[trans_from] = Set{String}()
    end
    push!(valid_transitions(vs)[trans_from], trans_to)
end

is_validated(vs::ValidationState, seq_name::String) = in(seq_name, valid_sequences(vs))
is_validated(vs::ValidationState, trans_from::String, trans_to::String) = haskey(valid_transitions(vs), trans_from) && in(trans_to, valid_transitions(vs)[trans_from])

ValidationState() = ValidationState(Set{String}(), Dict{String, Set{String}}())

function window_validate!(vs::ValidationState, val_func::Function, window_sz::Int, sig::ImagineSignal)
    seq_nms = sequence_names(sig)
    lookup = sequence_lookup(sig)
    seq_lens = Dict()
    seqs_used = unique(seq_nms)
    for nm in seqs_used #precompute compressed vector lengths
        seq_lens[nm] = full_length(lookup[nm])
    end
    #check each sequence used
    for nm in seqs_used
        if !is_validated(vs, nm) && seq_lens[nm] >= window_sz
            samps = get_samples(sig, nm; sampmap=:raw).data
            isgood = val_func(samps, window_sz)
            if !isgood
                print("This sequence is invalid: $nm\n")
                error() #The caller should catch this to give a more specific error message
            else
                set_validated!(vs, nm)
            end
        end
    end
    #check all transitions between sequences
    for i = 1:(length(seq_nms)-1)
        this_nm = seq_nms[i]
        next_nm = seq_nms[i+1]
        if !is_validated(vs, this_nm, next_nm)
            this_len = seq_lens[this_nm]
            @assert this_len >= 1
            next_len = seq_lens[next_nm]
            #NOTE: below is inefficient because we read the whole sequence.  If it becomes an issue, we should define getindex for RLEVectors
            this_samps = view(get_samples(sig, this_nm; sampmap=:raw).data, max(1, this_len - window_sz + 1):this_len)
            next_samps = -1;
            record_validation = next_len >= (window_sz-1) #only record validation if we have enough samples in the next sequence to do a check
            out_of_seqs = false #useful for if we don't have enough samples to check the last sequence transition
            if !record_validation #keep reading more samples until we have enough.  This is pretty inefficient, but it should happen very rarely.
                nexti = i+1
                enough = false
                next_samps = get_samples(sig, next_nm; sampmap=:raw).data
                while !enough
                    needed = (window_sz-1) - length(next_samps)
                    next_next = get_samples(sig, seq_nms[nexti]; sampmap=:raw).data
                    if length(next_next) >= needed
                        enough = true
                        next_samps = append!(next_samps, next_next[1:needed])
                    elseif nexti > length(seq_nms) #we've run out of sequences
                        out_of_seqs = true
                        break
                    else
                        next_samps = append!(next_samps, next_next)
                        nexti+=1
                    end
                end
            else
                next_samps = view(get_samples(sig, next_nm; sampmap=:raw).data, 1:window_sz)
            end
            samps_to_check = out_of_seqs ? get_samples(sig, length(sig) - 2*window_sz + 2, length(sig)).data : vcat(this_samps, next_samps)
            isgood = val_func(samps_to_check, window_sz)
            if !isgood
                print("Transition with error: from $this_nm to $next_nm\n")
                error() #The caller should catch this to give a more specific error message
            elseif record_validation
                set_validated!(vs, this_nm, next_nm)
            end
        end
    end
    return vs
end

window_validate(val_func::Function, window_sz::Int, sig::ImagineSignal) = window_validate!(ValidationState(), val_func, window_sz, sig)

max_forks_repo_name: HolyLab/ImagineInterface


text:
"""
    is_observable_fixed(ego::VehicleState, car::VehicleState, env::OccludedEnv)
    is_observable_fixed(x::Float64, y::Float64, ego::VehicleState, env::OccludedEnv)
Occlusion checker considering only fixed obstacles in the environment:
- is_observable_fixed(ego::VehicleState, car::VehicleState, env::OccludedEnv) check if car is observable from ego
- is_observable_fixed(x::Float64, y::Float64, ego::VehicleState, env::OccludedEnv) check if x, y is observable from ego
"""
function is_observable_fixed(ego::VehicleState, car::VehicleState, env::OccludedEnv)
    m = length(env.obstacles)
    front = ego.posG + polar(VehicleDef().length/2, ego.posG.θ)
    angle = atan(car.posG.y - front.y, car.posG.x - front.x)
    ray = Projectile(VecSE2(front.x, front.y, angle), 1.0)
    for i = 1:m
        if is_colliding(ray, env.obstacles[i], car.posG)
            return false
        end
    end
    return true
end

# check if a point in space is observable
function is_observable_fixed(x::Float64, y::Float64, ego::VehicleState, env::OccludedEnv)
    m = length(env.obstacles)
    front = ego.posG + polar(VehicleDef().length/2, ego.posG.θ)
    angle = atan(y - front.y, x - front.x)
    ray = Projectile(VecSE2(front.x, front.y, angle), 1.0)
    for i = 1:m
        if is_colliding(ray, env.obstacles[i], VecSE2(x,y,0.))
            return false
        end
    end
    return true
end

"""
    is_observable_dyna(ego::Vehicle, car::Vehicle, scene::Scene)
Occlusion checker verifying is `car` is observable from `ego`
It considers only other vehicles (not pedestrians, not fixed obstacles)
"""
function is_observable_dyna(ego::Vehicle, car::Vehicle, scene::Scene)
    front = ego.state.posG + polar(ego.def.length/2, ego.state.posG.θ)
    angle = atan(car.state.posG.y - front.y, car.state.posG.x - front.x)
    ray = Projectile(VecSE2(front.x, front.y, angle), 1.0)
    for veh in scene
        if veh.id == car.id || veh.id == ego.id || veh.def == PEDESTRIAN_DEF
            continue
        end
        box = get_oriented_bounding_box(veh)
        if is_colliding(ray, box, car.state.posG)
            return false
        end
    end
    return true
end

"""
    is_observable(ego::Vehicle, car::Vehicle, scene::Scene, env::OccludedEnv)
Check if car is observable from ego, by considering both fixed and dynamic obstacles
"""
function is_observable(ego::Vehicle, car::Vehicle, scene::Scene, env::OccludedEnv)
    return is_observable_dyna(ego, car, scene) && is_observable_fixed(ego.state, car.state, env)
end

function AutomotiveDrivingModels.is_colliding(P::Projectile, poly::ConvexPolygon, A::VecSE2)
    # collides if at least one of the segments collides with the ray
    point_time = norm(VecE2(A - P.pos))
    for i in 1 : length(poly)
        seg = get_edge(poly, i)
        obs_time = get_intersection_time(P, seg)
        if !isinf(obs_time) && obs_time < point_time
            return true
        end
    end
    false
end

# check if x,y is in veh
function AutomotiveDrivingModels.is_colliding(x::Float64, y::Float64, veh::VehicleState)
    v = VecSE2(x,y, 0.)
    poly = get_oriented_bounding_box(Vehicle(veh, VehicleDef(), -1))
    return contains(poly, v)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using Sockets
conn = Sockets.connect(8080)
write(conn,"Do you Hear me?\n")

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using CalibrationAnalysis
using Plots
using StableRNGs

using Test

@testset "CalibrationAnalysis.jl" begin
    # We only check that a simple example runs without errors
    # Otherwise we rely on tests in CalibrationErrors, CalibrationTests, and
    # ReliabilityDiagrams
    rng = StableRNG(8317)
    n = 100
    confidence = rand(rng, n)
    outcomes = rand(rng, Bool, n)

    # ReliabilityDiagrams.jl
    reliabilityplot(confidence, outcomes) # Plots
    reliability(confidence, outcomes) # Makie

    # CalibrationErrors.jl
    ece = ECE(UniformBinning(10))
    @test ece(confidence, outcomes) > 0

    kernel = ExponentialKernel() ⊗ WhiteKernel()
    skce = SKCE(kernel)
    skceb = SKCE(kernel; unbiased=false)
    @test skce(confidence, outcomes) < skceb(confidence, outcomes)
    @test skceb(confidence, outcomes) > 0

    # CalibrationTests.jl
    @test pvalue(StableRNG(1447), ConsistencyTest(ece, confidence, outcomes)) < 0.05
    @test pvalue(StableRNG(9483), AsymptoticSKCETest(kernel, confidence, outcomes)) < 0.05
end


text:
abstract Nucleotide2bit

## 2-bit DNA Seq Type ##

type Nucleotide2bitSeq <: Nucleotide2bit
  b1::BitVector
  b2::BitVector
  Nucleotide2bitSeq(b1::BitVector,b2::BitVector) = new(b1,b2)
end

Nucleotide2bitSeq(len::Int) = Nucleotide2bitSeq(BitArray(len),BitArray(len))

## Nucleotide2bitBase Type ##

immutable Nucleotide2bitBase <: Nucleotide2bit
  b1::Bool
  b2::Bool
  Nucleotide2bitBase(b1::Bool,b2::Bool) = new(b1,b2)
end

## Selection ##

getindex(s::Nucleotide2bitSeq,ind::Int)  = Nucleotide2bitBase(getindex(s.b1,ind),getindex(s.b2,ind))

getindex(s::Nucleotide2bitSeq,ind...)    = Nucleotide2bitSeq(getindex(s.b1,ind...),getindex(s.b2,ind...))

## Assignation ##

function setindex!(s::Nucleotide2bitSeq,x::Nucleotide2bitBase,ind...)
  setindex!(s.b1,x.b1,ind...)
  setindex!(s.b2,x.b2,ind...)
end

function setindex!(s::Nucleotide2bitSeq,x,ind...)
  inseq = convert(Nucleotide2bitBase,x)
  setindex!(s.b1,inseq.b1,ind...)
  setindex!(s.b2,inseq.b2,ind...)
end

## Copy ##

copy(seq::Nucleotide2bitSeq) = Nucleotide2bitSeq(copy(seq.b1),copy(seq.b2))

## length ##

length(s::Nucleotide2bitSeq) = length(s.b1)

## Convertions ##

##	A	C	T	G
## b1	0	0	1	1
## b2	0	1	0	1

const _convert_to_Char = ['A' 'C';'T' 'G']

function convert{T<:Union{Integer, Char}}(::Type{Array{T,1}},s::Nucleotide2bitSeq)
  len = length(s)
  out = Array(T,len)
  for i in 1:len
    out[i] = _convert_to_Char[s.b1[i]+1,s.b2[i]+1]
  end
  out
end

function convert{T<:Union{Integer, Char}}(::Type{T},s::Nucleotide2bitBase)
    convert(T,_convert_to_Char[s.b1+1,s.b2+1])
end

const _convert_to_base2_index = zeros(UInt8,256)
_convert_to_base2_index[Int('A')] = 1
_convert_to_base2_index[Int('a')] = 1
_convert_to_base2_index[Int('C')] = 2
_convert_to_base2_index[Int('c')] = 2
_convert_to_base2_index[Int('T')] = 3
_convert_to_base2_index[Int('t')] = 3
_convert_to_base2_index[Int('U')] = 3
_convert_to_base2_index[Int('u')] = 3
_convert_to_base2_index[Int('G')] = 4
_convert_to_base2_index[Int('g')] = 4

const _convert_to_base2 = Array(Nucleotide2bitBase,4)

_convert_to_base2[1] = Nucleotide2bitBase(false,false)
_convert_to_base2[2] = Nucleotide2bitBase(false,true)
_convert_to_base2[3] = Nucleotide2bitBase(true,false)
_convert_to_base2[4] = Nucleotide2bitBase(true,true)

function convert{T<:Union{Integer, Char}}(::Type{Nucleotide2bitSeq},s::Vector{T})
  len = length(s)
  seq = Nucleotide2bitSeq(len)
  for i in 1:len
    seq[i] = _convert_to_base2[ _convert_to_base2_index[ s[i] ] ]
  end
  seq
end

function convert{T<:Union{Integer, Char}}(::Type{Nucleotide2bitBase},s::T)
  _convert_to_base2[ _convert_to_base2_index[ s ] ]
end

nucleotide2bit{T<:Union{Integer, Char}}(x::Vector{T}) = convert(Nucleotide2bitSeq,x)
nucleotide2bit{T<:Union{Integer, Char}}(x::T) = convert(Nucleotide2bitBase,x)

macro nt2_str(s);  :(nucleotide2bit(@b_str($s))); end

convert(::Type{ASCIIString}, seq::Nucleotide2bitSeq) = ASCIIString(convert(Vector{UInt8},seq))
convert(::Type{Nucleotide2bitSeq}, str::ASCIIString) = nucleotide2bit(str.data)

bytestring(seq::Nucleotide2bitSeq) = bytestring(convert(Vector{UInt8},seq))

nucleotide2bit(x::ASCIIString) = convert(Nucleotide2bitSeq,x)

nucleotide(seq::Nucleotide2bitSeq)  = convert(Vector{Nucleotide},seq)

## == ##

promote_rule{T<:Union{Integer, Char},B<:Nucleotide2bit}(::Type{B}, ::Type{T} ) = T
promote_rule{T<:Union{Integer, Char},B<:Nucleotide2bit}(::Type{T}, ::Type{B} ) = T

promote_rule{T<:Union{Integer, Char}}(::Type{Nucleotide2bitSeq}, ::Type{Vector{T}} ) = Vector{T}
promote_rule{T<:Union{Integer, Char}}(::Type{Vector{T}}, ::Type{Nucleotide2bitSeq} ) = Vector{T}

=={T<:Nucleotide2bit}(S1::T, S2::T) = ==(S1.b1,S2.b1) && ==(S1.b2,S2.b2)

=={T<:Union{Integer, Char},B<:Nucleotide2bit}(S1::B, S2::T) = ==(promote(S1,S2)...)
=={T<:Union{Integer, Char},B<:Nucleotide2bit}(S1::T, S2::B) = ==(promote(S1,S2)...)

=={T<:Union{Integer, Char}}(S1::Nucleotide2bitSeq, S2::Vector{T}) = ==(promote(S1,S2)...)
=={T<:Union{Integer, Char}}(S1::Vector{T}, S2::Nucleotide2bitSeq) = ==(promote(S1,S2)...)

# Julia 0.2 compatibility
isequal{T<:Nucleotide2bit}(S1::T, S2::T) = isequal(S1.b1,S2.b1) && isequal(S1.b2,S2.b2)

isequal{T<:Union{Integer, Char},B<:Nucleotide2bit}(S1::B, S2::T) = isequal(promote(S1,S2)...)
isequal{T<:Union{Integer, Char},B<:Nucleotide2bit}(S1::T, S2::B) = isequal(promote(S1,S2)...)

isequal{T<:Union{Integer, Char}}(S1::Nucleotide2bitSeq, S2::Vector{T}) = isequal(promote(S1,S2)...)
isequal{T<:Union{Integer, Char}}(S1::Vector{T}, S2::Nucleotide2bitSeq) = isequal(promote(S1,S2)...)

## Loop ##

isempty(s::Nucleotide2bitSeq) = length(s) == 0

start(s::Nucleotide2bitSeq) = 1
next(s::Nucleotide2bitSeq,i) = (s[i],i+1)
done(s::Nucleotide2bitSeq,i) = (i > length(s))

## Write and Show ##

function show(io::IO,seq::Nucleotide2bitSeq)
  len = length(seq)
  print(len)
  println(" bp Nucleotide2bitSeq:")
  if len > 0
    screen = Base.tty_size()[1] > 6 ? Base.tty_size()[1] - 5 : Base.tty_size()[1]
    if length(seq) <= screen
      print(" ")
      print(Char(seq[1]))
      if(len>1)
	for i in 2:len
	  print("\n")
	  print(" ")
	  print(Char(seq[i]))
	end
      end
    else
      partlen = Int((screen)/2) - 1
      for i in 1:partlen
	print(" ")
	println(Char(seq[i]))
      end
      print(" \u22ee")
      for i in (len-partlen):len
	print("\n")
	print(" ")
	print(Char(seq[i]))
      end
    end
  end
end

function write(io::IO,seq::Nucleotide2bitSeq)
  for bp in seq
    write(io,Char(bp))
  end
end

show(io::IO,x::Nucleotide2bitBase) = (write(io,UInt8(x)); nothing)

## Complement ##

complement!{T<:Nucleotide2bit}(seq::T) = (seq.b1 = ~seq.b1 ; seq)
complement{T<:Nucleotide2bit}(seq::T) = complement!(copy(seq))

reversecomplement!{T<:Nucleotide2bit}(seq::T) = (reverse!(seq.b1); seq.b1 = ~seq.b1; reverse!(seq.b2); seq )
reversecomplement{T<:Nucleotide2bit}(seq::T) = reversecomplement!(copy(seq))

## Faster using bit level parallelism:

##	A	C	T	G
## b1	0	0	1	1
## b2	0	1	0	1

## A: ~ (b1 | b2)
isadenine{T<:Nucleotide2bit}(seq::T)	= ~ ( seq.b1 | seq.b2 )
## G: b1 & b2
isguanine{T<:Nucleotide2bit}(seq::T)	= seq.b1 & seq.b2
## C: ~b1 & b2
iscytosine{T<:Nucleotide2bit}(seq::T)	= ~seq.b1 & seq.b2
## T: b1 & ~b2
isthymine{T<:Nucleotide2bit}(seq::T)	= seq.b1 & ~seq.b2
## AT: ~b2
isweak{T<:Nucleotide2bit}(seq::T)		= ~seq.b2
## CG: b2
isstrong{T<:Nucleotide2bit}(seq::T)	= seq.b2
## CT: b1 $ b2
ispyrimidine{T<:Nucleotide2bit}(seq::T)	= seq.b1 $ seq.b2
## AG: ~( b1 $ b2 )
ispurine{T<:Nucleotide2bit}(seq::T)	= ~( seq.b1 $ seq.b2 )

function percentGC(seq::Nucleotide2bitSeq)
  len = length(seq)
  c2 = seq.b2.chunks
  chunks_len = length(c2)
  sum_ones = 0
  for i in 1:chunks_len
    sum_ones += count_ones(c2[i])
  end
  sum_ones/len
end

# Hash as Nucleotide
hash(seq::Nucleotide2bitSeq) = hash(nucleotide(seq))
hash(base::Nucleotide2bitBase) = hash(nucleotide(base))

# Hash Nucleotide2bitSeq as Vector{Nucleotide}
hash(seq::Nucleotide2bitSeq) = hash(nucleotide(seq))

==(S1::Nucleotide2bitSeq, S2::Vector{Nucleotide}) = ==(nucleotide(S1),S2)
==(S1::Vector{Nucleotide}, S2::Nucleotide2bitSeq) = ==(S1,nucleotide(S2))

# Julia 0.2 compatibility
isequal(S1::Nucleotide2bitSeq, S2::Vector{Nucleotide}) = isequal(nucleotide(S1),S2)
isequal(S1::Vector{Nucleotide}, S2::Nucleotide2bitSeq) = isequal(S1,nucleotide(S2))


text:
# %% Saltzman iceage models
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

function sm91(state, param, t)
    u, ν, p, r, s, q = param
    x, y, z = state
    dx = -x -y -u*z # - ν*astronomical(t/10) # ensure time is in kyr for interpolation
    dy = -p*z +r*y -s*y^2 - y^3
    dz = -q*(x+z)
    return SVector(dx, dy, dz)
end

# Astronomical forcing:
using DelimitedFiles, Statistics
X = readdlm(projectdir("exercise_data", "16.csv"))
tsol = X[:, 1] # time in kyr
Isol = X[:, 5]
# Transform simulation to last 500 kyr and reverse time axis
# L = length(Isol)
# Isol = Isol[1:L÷2]
# tsol = tsol[1:L÷2]
Isol = standardize(Isol)
tsol .+= abs(tsol[1])
using Interpolations
astronomical = LinearInterpolation(tsol, Isol)
quasiperiodic(t) = cospi(2t/22.1) - cospi(2t/44.1)
# Default parameter values
p = 1.0
q = 2.5
r = 1.3
s = 0.6
u = 0.2
ν = 0.0
t₀ = 10.0 # timescale in kyr
params = [u, ν, p, r, s, q]
u₀ = [-0.2, 1.0, 0.0]

ds = ContinuousDynamicalSystem(sm91, u₀, params)

fig, axs = subplots(1,3)

# Unforced solution
T = 50.0
Δt = 0.01
tr = trajectory(ds, T; Δt)
t = (0:Δt:T) .* t₀
for i in 1:3;
    x = tr[:, i]
    axs[1].plot(t, x, ls = LINESTYLES[i], color = "C$(i-1)")
    imax = argmax(x)
    label = "\$"*("x", "y", "z")[i]*"\$"
    axs[1].text(t[imax], x[imax], label; color = "C$(i-1)")
end
axs[1].set_xlabel("\$t\$ (kyr, arbtr.)")
axs[1].set_ylim(-1.5, 1.5)
fig.tight_layout(pad=0.3)

# # Forced solution
# set_parameter!(ds, 2, 0.5)
# tr = trajectory(ds, 120; Δt)
# axs[3].plot(tr[:, 2])

# transform tsol into same time units as above
# axs[2].plot(tsol, Isol; c = "C3")
# axs[2].set_xlabel("\$t\$ (kyr)")
# axs[2].text(200, 410, "\$F(t)\$", color = "C3")

# Van der poll version
function vdp(u, p, t)
    x, y = u
    α, β, γ, τ, A = p
    dx = -y - β + γ*astronomical_forcing(A, t)
    # dx = -y - β + γ*sinpi(2t/41)
    dy = -α*(y^3/3 - y - x)
    return SVector(dx/τ, dy/τ)
end

include("12.7_astronomical_forcing.jl")

α = 11.11
β = 0.25
γ = 0.75
τ = 40.09
p = (α, β, γ, τ, A)

u₀ = [0.1, -0.2]
T = 600.0
Δt = 0.1
t0 = 500.0
t = (0:Δt:T) .- t0

F = astronomical_forcing.(Ref(A), t)
axs[2].plot(t, F; color = "C3")
axs[2].text(-400, -5, "\$F(t)\$", color = "C3")
# %%
axs[3].clear()

ds = ContinuousDynamicalSystem(vdp, u₀, p; t0 = -t0)
tr = trajectory(ds, T; Δt)
axs[3].plot(t, standardize(tr[:, 1]))

# Plot vostok
V = readdlm(projectdir("exercise_data", "17.csv"))
Vt = V[:, 1] ./ 1000 # in ky
VT = V[:, 2]
Vt = -reverse(Vt)
VT = reverse(VT)
# Interpolate Vostok data for better plot clarity
using Dierckx
# interpolation object:
spl = Spline1D(Vt, VT; k=3, bc="extrapolate")
Vti = Vt[1]:(10Δt):Vt[end]
VTi = spl(Vti)

axs[3].plot(Vti, standardize(VTi);
    color = "C2", zorder = 99, alpha = 0.8, ls = "-.", lw = 2,
    # marker = "o", mew = 1, mfc = "none", alpha = 0.2, ms = 5,
)

axs[3].text(38.6, 1.61, "\$x\$"; color = "C0")
axs[3].text(-120, 2.6, "Vostok"; color = "C2")

axs[2].set_xlabel("\$t\$ (kyr)")
axs[3].set_xlabel("\$t\$ (kyr)")

add_identifiers!(fig)
wsave(plotsdir("12", "iceages"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
#Tests des packages
include("../ALAS/include.jl")

using MPCCmod
using OutputRelaxationmod
#package pour plot
using PyPlot

using NLPModels
using JuMP
using NLPModels
using Ipopt
using NLopt

path_to_data=""
include(string(path_to_data,"Portfolio.jl"))

#Algo choices:
solvers=["Ipopt"]
allrelaxations=["NL","SS","KS","ButterflyTheta1","ButterflyTheta1_C1","ButterflyTheta1_C2","Thetal0","Theta2l0","Thetal0scaled","Theta2l0scaled"];
thetal0=["Thetal0","Theta2l0","Thetal0scaled","Theta2l0scaled"]
relaxations=["ButterflyALAS"]
#problem infeasible with NL

pbs200=["size200/orl200-05-a","size200/orl200-005-a","size200/orl200-05-b","size200/orl200-005-b","size200/orl200-05-c","size200/orl200-005-c","size200/orl200-05-d","size200/orl200-005-d","size200/orl200-05-e","size200/orl200-005-e","size200/orl200-05-f","size200/orl200-005-f","size200/orl200-05-g","size200/orl200-005-g","size200/orl200-05-h","size200/orl200-005-h","size200/orl200-05-i","size200/orl200-005-i","size200/orl200-05-j","size200/orl200-005-j","size200/pard200_a","size200/pard200_b","size200/pard200_c","size200/pard200_d","size200/pard200_e","size200/pard200_f","size200/pard200_g","size200/pard200_h","size200/pard200_i","size200/pard200_j"]
pbs300=["size300/orl300_05_a","size300/orl300_005_a","size300/orl300_05_b","size300/orl300_005_b","size300/orl300_05_c","size300/orl300_005_c","size300/orl300_05_d","size300/orl300_005_d","size300/orl300_05_e","size300/orl300_005_e","size300/orl300_05_f","size300/orl300_005_f","size300/orl300_05_g","size300/orl300_005_g","size300/orl300_05_h","size300/orl300_005_h","size300/orl300_05_i","size300/orl300_005_i","size300/orl300_05_j","size300/orl300_005_j","size300/pard300_a","size300/pard300_b","size300/pard300_c","size300/pard300_d","size300/pard300_e","size300/pard300_f","size300/pard300_g","size300/pard300_h","size300/pard300_i","size300/pard300_j"]
pbs400=["size400/orl400_05_a","size400/orl400_005_a","size400/orl400_05_b","size400/orl400_005_b","size400/orl400_05_c","size400/orl400_005_c","size400/orl400_05_d","size400/orl400_005_d","size400/orl400_05_e","size400/orl400_005_e","size400/orl400_05_f","size400/orl400_005_f","size400/orl400_05_g","size400/orl400_005_g","size400/orl400_05_h","size400/orl400_005_h","size400/orl400_05_i","size400/orl400_005_i","size400/orl400_05_j","size400/orl400_005_j","size400/pard400_a","size400/pard400_b","size400/pard400_c","size400/pard400_d","size400/pard400_e","size400/pard400_f","size400/pard400_g","size400/pard400_h","size400/pard400_i","size400/pard400_j"]
pbs=vcat(pbs200,pbs300,pbs400)

names=["VaR","CVaR","RVaR","RCVar"];
beta=[0.9 0.95 0.99];
cb=[1.2816 1.6449 2.3263;1.7550 2.0627 2.6652;1.3333 2.0647 4.9247;3.0000 4.3589 9.9499];

#sorties
npb=length(pbs);npb=1;
nrelax=length(relaxations);nrelax=1;
nsolv=length(solvers);nsolv=1;
nmod=size(cb,1);nmod=1;
nbeta=length(beta);nbeta=1;

obj=zeros(npb,nrelax,nsolv,nmod);
tm=zeros(npb,nrelax,nsolv,nmod);
feas=zeros(npb,nrelax,nsolv,nmod);

prec=1e-6
r0=0.1
s0=0.1
#t0=sqrt(0.1) #t=o(r) normalement
t0=0.01
sig_r=0.1
sig_s=0.1
sig_t=0.01

for mo=1:nmod
 for b=1:nbeta
 cbmo=cb[mo,b]
  for p=1:npb
	pb=pbs[p]
	@printf("%f %s \n",p,pb)

	#nous les vecteurs l et u
	bds=readdlm(string("data/mv/",pb,".bds"));
	u=bds[:,2];

	# nous donne la matrice Q
	Q=readdlm(string("data/mv/",pb,".mat"))
	n=size(Q,1);

	txt=readdlm(string(path_to_data,"data/mv/",pb,".txt"));
	mu=txt[:,1];

	kappa=20;
#n=5;Q=Q[1:n,1:n];mu=mu[1:n];u=u[1:n];kappa=2;
	#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	#%CHOIX DU POINT INITIAL
	x0=prec*ones(n);y0=ones(n);z0=vcat(x0,y0)
	#%CHOIX DE LA GESTION DES PARAMETRES (r,s,t,tb)
	nlp,G,H=Portfolio(cbmo,mu,Q,u,kappa,z0,prec)

	#CREATION DU MPCCmod:
	omcc=MPCCmod.MPCC(nlp,G,H)

	#APPEL DU SOLVER:
@time xb,fb,orb,nb_eval = MPCCsolve.solve(omcc,r0,sig_r,s0,sig_s,t0,sig_t)

	try
 	oa=orb.inner_output_alas[1]
 	@printf("%s	",orb.solve_message)
	end
	@show nb_eval

  end
 end
end


text:
module Buoyancy

export
    BuoyancyTracer, SeawaterBuoyancy, buoyancy_perturbation,
    LinearEquationOfState, RoquetIdealizedNonlinearEquationOfState,
    ∂x_b, ∂y_b, ∂z_b, buoyancy_perturbation, buoyancy_frequency_squared

# Physical constants
# https://en.wikipedia.org/wiki/Gravitational_acceleration#Gravity_model_for_Earth (30 Oct 2019)
const g_Earth = 9.80665

"""
    AbstractBuoyancy{EOS}

Abstract supertype for buoyancy models.
"""
abstract type AbstractBuoyancy{EOS} end

"""
    AbstractEquationOfState

Abstract supertype for equations of state.
"""
abstract type AbstractEquationOfState end

"""
    AbstractNonlinearEquationOfState

Abstract supertype for nonlinar equations of state.
"""
abstract type AbstractNonlinearEquationOfState <: AbstractEquationOfState end

include("buoyancy_utils.jl")
include("no_buoyancy.jl")
include("buoyancy_tracer.jl")
include("seawater_buoyancy.jl")
include("linear_equation_of_state.jl")
include("nonlinear_equation_of_state.jl")
include("roquet_idealized_nonlinear_eos.jl")
include("show_buoyancy.jl")

end


text:
# EXCLUDE FROM TESTING
# this file doesn't have an entry point, see `verify.jl` instead

# Fast parallel reduction for Kepler hardware
# - uses shuffle and shared memory to reduce efficiently
# - support for large arrays
#
# Based on devblogs.nvidia.com/parallelforall/faster-parallel-reductions-kepler/

using CUDAdrv, CUDAnative, CuArrays

"""Reduce a value across a warp"""
@inline function reduce_warp(op::F, val::T)::T where {F <: Function,T}
    offset = CUDAnative.warpsize() ÷ 2
    # TODO: this can be unrolled if warpsize is known...
    while offset > 0
        val = op(val, shfl_down(val, offset))
        offset ÷= 2
    end
    return val
end

"""Reduce a value across a block, using shared memory for communication"""
@inline function reduce_block(op::F, val::T)::T where {F <: Function,T}
    # shared mem for 32 partial sums
    shared = @cuStaticSharedMem(T, 32)

    wid, lane = fldmod1(threadIdx().x, CUDAnative.warpsize())

    # each warp performs partial reduction
    val = reduce_warp(op, val)

    # write reduced value to shared memory
    if lane == 1
        @inbounds shared[wid] = val
    end

    # wait for all partial reductions
    sync_threads()

    # read from shared memory only if that warp existed
    @inbounds val = (threadIdx().x <= fld(blockDim().x, CUDAnative.warpsize())) ? shared[lane] : zero(T)

    # final reduce within first warp
    if wid == 1
        val = reduce_warp(op, val)
    end

    return val
end

"""Reduce an array across a complete grid"""
function reduce_grid(op::F, input::CuDeviceVector{T}, output::CuDeviceVector{T},
                     len::Integer) where {F <: Function,T}
    # TODO: neutral element depends on the operator (see Base's 2 and 3 argument `reduce`)
    val = zero(T)

    # reduce multiple elements per thread (grid-stride loop)
    # TODO: step range (see JuliaGPU/CUDAnative.jl#12)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    step = blockDim().x * gridDim().x
    while i <= len
        @inbounds val = op(val, input[i])
        i += step
    end

    val = reduce_block(op, val)

    if threadIdx().x == 1
@inbounds output[blockIdx().x] = val
end

    return
end

"""
Reduce a large array.

Kepler-specific implementation, ie. you need sm_30 or higher to run this code.
"""
function gpu_reduce(op::Function, input::CuArray{T}, output::CuArray{T}) where {T}
    len = length(input)

    function get_config(kernel)
        fun = kernel.fun
        config = launch_configuration(fun)

        # pad to cover as many elements as possible (`cld`), but keep in mind
        # how we launch this kernel for the final reduction step `min(..., threads)`
        blocks = min(cld(len, config.threads), config.threads)

        # the output array must have a size equal to or larger than the number of blocks
        # in the grid because each block writes to a unique location within the array.
        if length(output) < blocks
            throw(ArgumentError("output array too small, should be at least $blocks elements"))
        end

        return (threads = config.threads, blocks = blocks)
    end

    # NOTE: manual expansion of @cuda because we need the result of the dynamic launch config
    # @cuda config=get_config reduce_grid(op, input, output, len)
    args = (op, input, output, len)
    GC.@preserve args begin
        kernel_args = cudaconvert.(args)
        kernel_tt = Tuple{Core.Typeof.(kernel_args)...}
        kernel = cufunction(reduce_grid, kernel_tt)
        kernel_config = get_config(kernel)
        kernel(kernel_args...; kernel_config...)
    end

    @cuda threads = kernel_config.threads reduce_grid(op, output, output, kernel_config.blocks)
end

# FURTHER IMPROVEMENTS:
# - use atomic memory operations
# - dynamic block/grid size based on device capabilities
# - vectorized memory access
#   devblogs.nvidia.com/parallelforall/cuda-pro-tip-increase-performance-with-vectorized-memory-access/


text:
using Jecco
using Plots
using Printf
pyplot()

ts = OpenPMDTimeSeries("./data")

for it in ts.iterations
    phi, grid=get_field(ts, it=it, field="phi c=1");
    ti = ts.current_t
    heatmap(phi[1,:,:],
            aspect_ratio=:equal,
            clims=(-0.8, 0.8),
            # title="t =  $ti",
            title="t =  $(@sprintf("%.2f", ti))",
            )

    padit = lpad(it, 4, "0")
    savefig("figs/phi_$padit.png")
end


text:
function get_sub_matrices(cur_matrix::FourBlockMatrix)
  A, B, C, D = map(
    cur_field -> getfield(cur_matrix, cur_field),
    [:A, :B, :C, :D]
  )

  [A, B, C, D]
end


text:
struct NmmSystem
    pot::Function
    mass::Float64
    beta::Float64

    function NmmSystem(pot::Function, mass::Float64, temperature::Float64)
        beta = 1 / (kB * temperature)
        new(pot, mass, beta)
    end
end


text:
using ExtendedKronigPennyMatrix
using LinearAlgebra
using PyPlot

function plot_KP(v0, rho)
   clf()
   cm=get_cmap("tab10")

   pot=KronigPennyPotential(v0, rho)
   pf = get_potential(pot)
   a=1
   xs=-2a:a/100:3a
   plot(xs, pf.(xs), "k")

   nmax=60
   Ka=0
   model=KronigPennyModel(pot, nmax, Ka)

   for Ka in (-18:18)/18*π
      update!(model, Ka=Ka)
      ev = eigvals(model.hnm)
      for i in 1:5
         plot(Ka/ π, ev[i], ".", color=cm(i-1))
      end
   end
   xlim(-1,1)
   ylim(-2,32)
   xlabel(L"$Ka / \pi$")
   ylabel(L"Energy / $E_0$")
   title( L"$\rho =$"*string(rho))
end

function main()
   #
   plot_KP(      0, 0.5)
   ylim(-0.5,25)
   savefig("Pavelich_Fig4a.png")
   #
   plot_KP(     10, 0.5)
   ylim(-0.5,30)
   savefig("Pavelich_Fig4b.png")
   savefig("KP_10_05.png")
   #
   plot_KP(     10, 0.8)
   ylim(-0.5,30)
   savefig("KP_10_08.png")
   #
   plot_KP(20.5607, 0.5)
   ylim(-0.5,40)
   savefig("Pavelich_Fig6.png")
   #
   plot_KP(10.8775, 0.8)
   ylim(-0.5,30)
   savefig("Pavelich_Fig7.png")
end

main()


text:
module Deps # Bukdu

import HTTP
import MbedTLS
import URIParser

const Request = HTTP.Messages.Request

end # module Bukdu.Deps

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
using LinearAlgebra
using LinearAlgebra: I

using ForwardDiff
using StaticArrays

import IntervalArithmetic: where_bisect

export Bisection, Newton, Krawczyk, GradientContractor, TrivialContractor

struct Newton end
struct Krawczyk end
struct Bisection end

function _contract(f, derivative, region::T, mid_point, ::Newton) where {T}
    m = convert(T, mid(region, mid_point))
    return m - (f(m) / derivative(region))
end

function _contract(f, jacobian, region::R, mid_point, ::Newton) where {T, R <: AbstractVector{T}}
    m = T.(mid.(region, mid_point))
    J = jacobian(region)
    fm = f(m)
    try
        return convert(R, m .- (J \ fm))
    catch
        return region .± Inf
    end
end

function _contract(f, derivative, region::T, mid_point, ::Krawczyk) where T
    m = convert(T, mid(region, mid_point))
    Y = 1 / derivative(m)
    return m - Y * f(m) + (1 - Y * derivative(region)) * (region - m)
end

function _contract(f, jacobian, region::T, mid_point, ::Krawczyk) where T <: AbstractVector
    m = mid.(region, mid_point)
    mm = convert(T, m)
    J = jacobian(region)
    fmm = f(mm)
    try
        Y = mid.(inv(jacobian(mm)))
        return m - Y * fmm + (I - Y * J) * (region - m)
    catch
        return region .± Inf
    end
end

"""
    GradientContractor(f, method, jacobian)
    GradientContractor(f, method, region)

Contractor for `f` based on `method` (`Newton()` or `Krawczyk()`), either using
the explicitely supplied function `jacobian(region)` or a jacobian method for
`region` using automatic differentiation of `f` from `ForwardDiff.jl`.

    (contractor::GradientContractor)(region, mid_point = IntervalArithmetic.where_bisect)

Calling instances of `GradientContractor` returns a contraction of `region`.
"""
struct GradientContractor{T}
    f
    method::T
    jacobian
end

function GradientContractor(f, method, ::R) where R <: AbstractVector
    function grad(region)
        jac = similar(region, Size(length(region), length(region)))
        ForwardDiff.jacobian!(jac, f, region)
        return jac
    end
    return GradientContractor(f, method, grad)
end

function GradientContractor(f, method, ::R) where R <: Number
    derivative(region) = ForwardDiff.derivative(f, region)
    return GradientContractor(f, method, derivative)
end

function (contractor::GradientContractor{T})(region, mid_point = where_bisect) where T
    contraction = _contract(
        contractor.f,
        contractor.jacobian,
        region,
        mid_point,
        contractor.method
    )
    if _isempty(contraction)
        throw(DomainError(region))
    end
    return contraction
end

"""
    TrivialContractor()

Trivial contractor used for benchmarking which solely relies on bisections.

    (contractor::TrivialContractor)(region) = region

Calling instances of `TrivialContractor` returns `region` itself.
"""
struct TrivialContractor end

(contractor::TrivialContractor)(region) = region


text:
using FastGroupBy,RCall

const M=100_000_000; const K=100;
srand(1);
svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M);
@time radixsort!(svec1); #18, 18.6, 17.24
issorted(svec1);

R"""
memory.limit(2e9)
N=$M; K=$K
set.seed(1)
library(data.table)
id3 = sample(sprintf("i%07d",1:(N/K)), N, TRUE)
pt = proc.time()
system.time(sort(id3, method="radix"))
data.table::timetaken(pt) # 18.9, 20.6 22.1
"""

# Roughly on par with R for length 8 id strings

# using longer string to compare speed
srand(1);
svec1 = rand(["id"*dec(k,10) for k in 1:M÷K], M);
@time radixsort!(svec1); # 26.89
issorted(svec1);

R"""
memory.limit(2e9)
N=$M; K=$K
set.seed(1)
library(data.table)
id3 = sample(sprintf("i%010d",1:(N/K)), N, TRUE)
pt = proc.time()
system.time(sort(id3, method="radix"))
data.table::timetaken(pt) # 23.1
"""

# Roughly on par with R for length 10 id strings; but is about 3 seconds slower

# sorting strings of variable length upto length 8
const M=10_000_000; const K=100;
srand(1);
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M);
@time sort!(svec1); # 50 seconds
issorted(svec1)

using RCall
R"""
memory.limit(2e9)
N=$M; K=$K
set.seed(1)
library(data.table)
bs = sapply(1:(N/K), function(i) rawToChar(sample(as.raw(32:126), sample(8), replace=T)))
id3 = sample(bs, N, replace = T)
pt = proc.time()
system.time(sort(id3, method="radix"))
data.table::timetaken(pt) #24.2
"""

# using longer string to compare speed
srand(1);
svec1 = rand([string(rand(Char.(32:126), rand(1:16))...) for k in 1:M÷K], M);
# using FastGroupBy.radixsort! to sort strings of length 16
@time radixsort!(svec1); # 41 seconds on 100m
issorted(svec1)

R"""
memory.limit(2e9)
N=$M; K=$K
set.seed(1)
library(data.table)
bs = sapply(1:(N/K), function(i) rawToChar(sample(as.raw(32:126), sample(8), replace=T)))
id3 = sample(bs, N, replace = T)
pt = proc.time()
system.time(sort(id3, method="radix"))
data.table::timetaken(pt) #24.2
"""

# Slower than R for length of 16; 1.7x times slower


text:
"""
    BlockingChequerboardSwitch <: BlockingType

A decorator indicating a switch (in a chequerboard pattern) that changes the
locations of delimiters of blocks.
"""
struct BlockingChequerboardSwitch <: BlockingType
    block_half_len::Int64
end


text:
unset = Unset()

struct FuncArg
    name
    type
    default
end

struct FuncHeader
    name  :: Any
    args  :: Any
    kwargs :: Any
    ret   :: Any
    fresh :: Any
end

FuncHeader() = FuncHeader(unset, unset, unset, unset, unset)

is_func_header(a::FuncHeader) = a.args != unset

function func_arg(@nospecialize(ex))::FuncArg
    @match ex begin
        :(:: $ty) => @with func_arg(gensym("_")).type = ty
        :($var :: $ty) => @with func_arg(var).type = ty
        Expr(:kw, var, default) => @with func_arg(var).default = default
        Expr(:(=), var, default) => @with func_arg(var).default = default
        var => FuncArg(var, unset, unset)
    end
end

function func_header(@nospecialize(ex))::FuncHeader
    @match ex begin
        :($hd::$ret) => @with func_header(hd).ret = ret

        :($f($(args...); $(kwargs...))) =>
            begin inter = @with func_header(f).args = map(func_arg, args)
                  @with inter.kwargs = map(func_arg, kwargs)
            end
        :($f($(args...))) => @with func_header(f).args = map(func_arg, args)

        :($f where {$(args...)}) => @with func_header(f).fresh = args

        Expr(:tuple, Expr(:parameters, kwargs...), args...) =>
            begin inter = @with FuncHeader().args = map(func_arg, args)
                  @with inter.kwargs = map(func_arg, kwargs)
            end
        Expr(:tuple, args...) => @with FuncHeader().args = map(func_arg, args)

        f => @with FuncHeader().name = f
    end
end

function of_args(::Unset)
    Argument[]
end

function of_args(args::AbstractArray{FuncArg})
    ret = Argument[]
    for (i, each) in enumerate(args)
        name = each.name === unset ? gensym("_$i") : each.name
        type = each.type === unset ? nothing : each.type
        arg = Argument(name, type, each.default)
        push!(ret, arg)
    end
    ret
end

extract_tvar(var :: Union{Symbol, Expr})::Symbol =
    @match var begin
        :($a <: $_) => a
        :($a >: $_) => a
        :($_ >: $a >: $_) => a
        :($_ <: $a <: $_) => a
        a::Symbol         => a
    end


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

import Test
Test.@testset "Package StruckVMEDevices" begin

include("MemRegisters/runtests.jl")
include("VMEGateways/runtests.jl")
include("SIS3316Digitizers/runtests.jl")

end # testset


text:
using Test, NLSolvers, LinearAlgebra
@testset "truncated SPR solver" begin
m = NLSolvers.TCG()

H = [0.3 0.0; 0.0 0.9]
g = [0.2,0.4]

m(g, H, 0.7, rand(2), 1)

# Since -H is negative definite the solution is guaranteed
# to be at the boundary unless g = 0
for Δ in range(0, 100; step=0.1)
    @test norm(m(g, -H, Δ, rand(2), 1).p, 2) ≈ Δ
end

# Small gradient entries
g .= [1e-12, 1e-9]
for Δ in range(0, 100; step=0.1)
    @test norm(m(g, -H, Δ, rand(2), 1).p, 2) ≈ Δ
end

# Mixed gradient entries
g .= [1e12, 1e-9]
for Δ in range(0, 100; step=0.1)
    @test norm(m(g, -H, Δ, rand(2), 1).p, 2) ≈ Δ
end

# Zero case
g = [0.0, 0.0]
for Δ in range(0, 100; step=0.1)
    @test norm(m(g, -H, Δ, rand(2), 1).p, 2) == 0
end
end


text:
include("misc.jl")

function leastSquares(X,y)

	# Find regression weights minimizing squared error
	w = (X'*X)\(X'*y)

	# Make linear prediction function
	predict(Xhat) = Xhat*w

	# Return model
	return GenericModel(predict)
end

max_stars_repo_name: nicolasigor/codeAppareanceLatex

max_issues_repo_name: nicolasigor/codeAppareanceLatex

max_forks_repo_name: nicolasigor/codeAppareanceLatex


text:
ONE_HOT_TO_DNA = Dict([1,0,0,0] => "A", [0,1,0,0] => "C", [0,0,1,0] => "G", [0,0,0,1] => "T");

convert_one_hot(x) = ONE_HOT_TO_DNA[x];

function one_hot_to_dna(seq)
    dna = ""
    for b in 1:4:length(seq)
         dna = string(dna, convert_one_hot(seq[b:b+3]));
    end

    dna
end

"""
encode_8x3x8(N::Int, alpha::Float64, ofile::String)

This will run a very basic 8x3x8 encoder to verify that feed forward and back-
propagation setps work. It will -un the encoder for N iterations with a learning
rate of ALPHA and output updated weights to the file specified at OFILE.
"""
function encode_8x3x8(N, alpha, ofile; write_output = false, verbose = true)

    input = one_hot_8x3x8(8)

    # Add a bias term to the input
    input = hcat(ones(8, 1), input)

    labels = one_hot_8x3x8(8)

    # draw random weights for first layer to hidden layaer
    w1 = rand(9, 3)
    # draw random weights for hidden layer to output layer
    w2 = rand(4, 8)

    for j in 1:N
        for i in 1:size(input, 1)

            inp = input[i,:]
            lab = labels[i,:]

            hidden, pred = feedforward(inp, w1, w2)

            u1, u2 = backpropogate(pred, lab, hidden, w1, w2)
            u1 = u1[1:3]

            # Now update weights
            w1 -= alpha * (u1.' * inp.').'
            w2 -= alpha * (u2.' * hidden).'

        end
    end

    # Report final predictions
    output = zeros(8,8)
    raw = zeros(8, 8)
    h_output = zeros(8, 3)

    sserr = 0.0
    for i in 1:size(input, 1)

        inp = input[i,:]
        lab = labels[i,:]
        hidden, pred = feedforward(inp, w1, w2)

        raw[i,:] = pred
        output[i,:] = threshold(pred)
        h_output[i,:] = hidden[1:3]

        sserr += sum((pred - lab.').^2)
    end

    if verbose
        println(string("Sum of Squared Error: ", sserr))
    end

    if write_output
        write_results(output, h_output, ofile)
    end

    raw, output, h_output, sserr

end

function nn_predict_on_data(tdata, labels, test_data, hl_size, N, alpha, ofile, write_output, C)
    # Add bias term to training data
    tdata = hcat(ones(size(tdata, 1), 1), tdata);

    # draw random weights for input layer to hidden layer
    w1 = rand(size(tdata, 2), hl_size);
    # draw random weights for hidden layer to output layer
    w2 = rand(hl_size+1, size(labels, 2));

    # Cross-Validate C times
    total_accuracy = 0.0
    for c in 1:C

        println(string("Cross Validation Round: ", c))

        # Select 90% as training and 10% as testing
        train_ii = rand(1:size(tdata, 1), convert(Int, round(0.9 * size(tdata, 1))));
        test_ii = setdiff(1:size(tdata, 1), train_ii);

        training_input = tdata[train_ii,:];
        training_labels = labels[train_ii, :];

        test_input = tdata[test_ii,:];
        test_labels = labels[test_ii,:];

        for j in 1:N

            if j % 100 == 0
                println(string("Epoch: ", j));
            end

            for i in 1:size(training_input, 1)

                inp = training_input[i,:];
                lab = training_labels[i,:];

                hidden, pred = feedforward(inp, w1, w2);

                u1, u2 = backpropogate(pred, lab, hidden, w1, w2);
                u1 = u1[1:hl_size];

                # Now update weights
                w1 -= alpha * (u1.' * inp.').';
                w2 -= alpha * (u2.' * hidden).';
            end

        end

        # Now test on the held out data to see how we did
        predictions = zeros(size(test_labels));
        prob_predictions = zeros(size(test_labels))
        dna_seqs = Array{String, 1}(size(test_labels, 1))
        for i in 1:size(test_input, 1)

            inp = test_input[i,:];
            lab = test_labels[i,:];
            hidden, pred = feedforward(inp, w1, w2);

            prob_predictions[i,:] = pred;
            predictions[i,:] = threshold(pred);

            # Remove bias term from the input for converting back to dna sequences
            dna_seqs[i] = one_hot_to_dna(inp[2:end]);

        end

        accuracy = calc_accuracy(test_labels, predictions);
        println(string("Accuracy: ", accuracy))

        total_accuracy += accuracy;
    end

    println(string("~~~~~~~Final Cross-Validation Accuracy~~~~~~\n", total_accuracy / C));

    if write_output

        # Add bias term for test_data to pass through neural network
        test_data = hcat(ones(size(test_data, 1), 1), test_data);

        prob_predictions = zeros(size(test_data, 1))
        dna_seqs = Array{String, 1}(size(test_data, 1))

        for i in 1:size(test_data, 1)

            inp = test_data[i,:];
            hidden, pred = feedforward(inp, w1, w2);

            prob_predictions[i,:] = pred;

            # Remove bias term from the input for converting back to dna sequences
            dna_seqs[i] = one_hot_to_dna(inp[2:end]);

        end

        write_seq_predictions(dna_seqs, prob_predictions, ofile);

    end

    total_accuracy / C

end

function train_nn(tdata, labels, hl_size, N, alpha, ofile, write_output::Bool, C::Int)

    # Add bias term to training data
    tdata = hcat(ones(size(tdata, 1), 1), tdata);

    # Cross-Validate C times
    total_accuracy = 0.0
    for c in 1:C

        println(string("Cross Validation Round: ", c))

        # Select 90% as training and 10% as testing
        train_ii = rand(1:size(tdata, 1), convert(Int, round(0.9 * size(tdata, 1))));
        test_ii = setdiff(1:size(tdata, 1), train_ii);

        training_input = tdata[train_ii,:];
        training_labels = labels[train_ii, :];

        test_input = tdata[test_ii,:];
        test_labels = labels[test_ii,:];

        # draw random weights for input layer to hidden layer
        w1 = rand(size(training_input, 2), hl_size);
        # draw random weights for hidden layer to output layer
        w2 = rand(hl_size+1, size(training_labels, 2));

        for j in 1:N

            if j % 100 == 0
                println(string("Epoch: ", j));
            end

            for i in 1:size(training_input, 1)

                inp = training_input[i,:];
                lab = training_labels[i,:];

                hidden, pred = feedforward(inp, w1, w2);

                u1, u2 = backpropogate(pred, lab, hidden, w1, w2);
                u1 = u1[1:hl_size];

                # Now update weights
                w1 -= alpha * (u1.' * inp.').';
                w2 -= alpha * (u2.' * hidden).';
            end

        end

        # Now test on the held out data to see how we did
        predictions = zeros(size(test_labels));
        prob_predictions = zeros(size(test_labels))
        dna_seqs = Array{String, 1}(size(test_labels, 1))
        for i in 1:size(test_input, 1)

            inp = test_input[i,:];
            lab = test_labels[i,:];
            hidden, pred = feedforward(inp, w1, w2);

            prob_predictions[i,:] = pred;
            predictions[i,:] = threshold(pred);

            # Remove bias term from the input for converting back to dna sequences
            dna_seqs[i] = one_hot_to_dna(inp[2:end]);

        end

        # if we want, we can just write out the last set of dna_seqs and probability
        # that they're binding sites
        if c == C && write_output
            write_seq_predictions(dna_seqs, prob_predictions, ofile);
        end

        accuracy = calc_accuracy(test_labels, predictions);
        println(string("Accuracy: ", accuracy))

        total_accuracy += accuracy;
    end

    println(string("~~~~~~~Final Accuracy~~~~~~\n", total_accuracy / C));

    total_accuracy / C

end

function nn_3layer(input, labels, hl_size, N, alpha, ofile; verbose = false)

    # Add bias term to input layer
    input = hcat(ones(size(input, 1), 1), input);

    # draw random weights for first layer to hidden layaer
    w1 = rand(size(input, 2), hl_size);
    # draw random weights for hidden layer to output layer
    w2 = rand(hl_size+1, size(labels, 2));

    for j in 1:N

        if j % 100 == 0
            println(string("Epoch: ", j));
        end

        for i in 1:size(input, 1)

            inp = input[i,:];
            lab = labels[i,:];

            hidden, pred = feedforward(inp, w1, w2);

            u1, u2 = backpropogate(pred, lab, hidden, w1, w2);
            u1 = u1[1:hl_size];

            # Now update weights
            w1 -= alpha * (u1.' * inp.').';
            w2 -= alpha * (u2.' * hidden).';
        end

    end

    output = zeros(size(labels));
    raw = zeros(size(labels));
    h_output = zeros(size(labels, 1), hl_size);

    sserr = 0.0
    for i in 1:size(input, 1)

        inp = input[i,:];
        lab = labels[i,:];
        hidden, pred = feedforward(inp, w1, w2);

        raw[i,:] = pred;
        output[i,:] = threshold(pred);
        h_output[i,:] = hidden[1:hl_size];

        sserr += sum((pred - lab.').^2);
    end

    accuracy = calc_accuracy(labels, output);

    println(string("Sum of Squared Error: ", sserr));
    println(string("Final Prediction Accuracy: ", accuracy));
    raw, output, h_output, sserr, accuracy

end

function calc_accuracy(labels, output)

    accur = 0.0
    for i in 1:size(labels, 1)

        if labels[i,:] == output[i,:]
            accur += 1
        end
    end

    accur / size(labels, 1)
end

function one_hot_8x3x8(n)

    eye(n)

end

# Define sigmoid activation function
sigmoid(x) = 1 / (1 + exp(-x))

"""
feed_forward_encoder(input::Array{Float64, 2}, w1::Array{Float64, 2}, w2::Array{Float64, 2})

Feed an INPUT array through a simple 8x3x8 encoder with input weights W1 and
hidden layer weights W2. Report back the hidden predictions as well as the
output predictions.
"""
function feedforward(x, w1, w2)

    # Feed to first layer
    h = x.' * (w1)
    hp = sigmoid.(h)
    hp = hcat(hp, 1)

    pred = sigmoid.(hp * w2)
    hp, pred

end

function backpropogate(predictions, label, hl, w1, w2)

    u2 = predictions .* (1 - predictions) .* (predictions - label.')

    u1 = hl .* (1 - hl) .* (w2 * u2.').'

    u1, u2

end

function threshold(arr)

    for i in 1:length(arr)
        if arr[i] > 0.5
            arr[i] = 1
        else
            arr[i] = 0
        end
    end

    arr

end

function write_results(output, preds, ofile)

    open(ofile, "w") do f

        for i in 1:size(output, 1)
            write(f, string(output[i,:], "\t", preds[i,:], "\n"))
        end
    end
end

function write_seq_predictions(output, preds, ofile)

    open(ofile, "w") do f

        for i in 1:size(output, 1)
            write(f, string(output[i,:][1], "\t", preds[i,:][1], "\n"))
        end
    end
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    demorect()

Demonstrate of rectangluar curve.

# Examples:
```julia
julia> demorect()
```
"""
function demorect()
crv = nrbtform(nrbrect(2.0,1.0),vecrotz(deg2rad(35)));
nrbplot(crv,4,
        title = "Construction and rotation of a rectangle.",
        framestyle=:box,
        legend = false);
end # demorect


text:
"""
    perfect_cube(number)

Check if a number is a perfect cube or not.

# Example
```jula
perfect_cube(27) # returns true
perfect_cube(4)  # returns false
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function perfect_cube(number)
    val = number^(1/3)
    return (val * val * val) == number
end

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia


text:
const ROOT_PATH       = pwd()
const DOC_ROOT_PATH   = ROOT_PATH     * "/public"
const CONFIG_PATH     = ROOT_PATH     * "/config"
const ENV_PATH        = CONFIG_PATH   * "/env"
const APP_PATH        = ROOT_PATH     * "/app"
const RESOURCES_PATH  = APP_PATH      * "/resources"
const TEST_PATH       = ROOT_PATH     * "/test"
const TEST_PATH_UNIT  = TEST_PATH     * "/unit"
const LIB_PATH        = ROOT_PATH     * "/lib"
const HELPERS_PATH    = APP_PATH      * "/helpers"
const LOG_PATH        = ROOT_PATH     * "/log"
const LAYOUTS_PATH    = APP_PATH      * "/layouts"
const TASKS_PATH      = ROOT_PATH     * "/tasks"
const BUILD_PATH      = ROOT_PATH     * "/build"

const GENIE_AUTHORIZATOR_FILE_NAME      = "authorization.yml"
const GENIE_DB_CONFIG_FILE_NAME         = "database.yml"

const GENIE_CONTROLLER_FILE_POSTFIX     = "Controller.jl"
const GENIE_CHANNEL_FILE_POSTFIX        = "Channel.jl"

const ROUTES_FILE_NAME  = joinpath(CONFIG_PATH, "routes.jl")

const PARAMS_REQUEST_KEY    = :REQUEST
const PARAMS_RESPONSE_KEY   = :RESPONSE
const PARAMS_SESSION_KEY    = :SESSION
const PARAMS_FLASH_KEY      = :FLASH
const PARAMS_ACL_KEY        = :ACL
const PARAMS_WS_CLIENT      = :WS_CLIENT

const TEST_FILE_IDENTIFIER = "_test.jl"

# Used to store log info during app bootstrap, when the logger itself is not available.
# The queue is automatically emptied by the logger upon load.
const GENIE_LOG_QUEUE = Vector{Tuple{String,Symbol}}()


text:
################################################################
"""
    timesum(snoop)

Calculates and prints the total time measured by a snoop macro.

It is used inside @snoopi_bench. Julia can cache inference results so to measure the effect of adding _precompile_() sentences generated by snoopi to your package, use the [`@snoopi_bench`](@ref). This benchmark measures inference time taken during loading and running of a package.

# Examples
```julia
using SnoopCompile
data = @snoopi begin
    include(joinpath(dirname(dirname(pathof(MatLang))),"test","runtests.jl"))
end;
println(timesum(data));
```

## Manual Benchmark (withtout using [`@snoopi_bench`](@ref))
- dev your package

- comment the precompile part of your package (`include()` and `_precompile_()`)
- run the following benchmark
- restart Julia

- uncomment the precompile part of your package (`include()` and `_precompile_()`)
- run the following benchmark
- restart Julia

### Benchmark
```julia
using SnoopCompile

println("Package load time:")
loadSnoop = @snoopi using MatLang

timesum(loadSnoop)

println("Running Examples/Tests:")
runSnoop = @snoopi begin
    using MatLang
    include(joinpath(dirname(dirname(pathof(MatLang))),"test","runtests.jl"))
end

timesum(runSnoop)
```
"""
function timesum(snoop::Vector{Tuple{Float64, Core.MethodInstance}})
    if isempty(snoop)
        return 0.0
    else
        return sum(first, snoop)
    end
end

################################################################
"""
    @snoopi_bench(packageName::String, snoopScript::Expr)
    @snoopi_bench(packageName::String)

Performs an infertime benchmark by activating and deactivating the _precompile_()
# Examples
Benchmarking the load infer time
```julia
println("loading infer benchmark")

@snoopi_bench "MatLang" using MatLang
```

Benchmarking the example infer time
```julia
println("examples infer benchmark")

@snoopi_bench "MatLang" begin
    using MatLang
    examplePath = joinpath(dirname(dirname(pathof(MatLang))), "examples")
    # include(joinpath(examplePath,"Language_Fundamentals", "usage_Entering_Commands.jl"))
    include(joinpath(examplePath,"Language_Fundamentals", "usage_Matrices_and_Arrays.jl"))
    include(joinpath(examplePath,"Language_Fundamentals", "Data_Types", "usage_Numeric_Types.jl"))
end
```
"""
macro snoopi_bench(packageName::String, snoopScript::Expr)

    ################################################################
    packagePath = joinpath(pwd(),"src","$packageName.jl")
    precompilePath, precompileFolder = precompile_pather(packageName)

    juliaCode = """
    using SnoopCompile; data = @snoopi begin
        $(string(snoopScript));
    end;
    println(timesum(data));
    """
    juliaCmd = `julia --project=@. -e "$juliaCode"`
    quote
        packageSym = Symbol($packageName)
        ################################################################
        using SnoopCompile
        println("""*******************
        Benchmark Started
        *******************
        """)
        ################################################################
        println("""Precompile Deactivated Benchmark
        ------------------------
        """)
        precompile_deactivator($packagePath, $precompilePath);
        ### Log the compiles
        run($juliaCmd)
        ################################################################
        println("""Precompile Activated Benchmark
        ------------------------
        """)
        precompile_activator($packagePath, $precompilePath)
        ### Log the compiles
        run($juliaCmd)
        println("""*******************
        Benchmark Finished
        *******************
        """)
    end

end

"""
    @snoopi_bench packageName::String

Benchmarking the infer time of the tests:
```julia
@snoopi_bench "MatLang"
```
"""
macro snoopi_bench(packageName::String)
    package = Symbol(packageName)
    snoopScript = :(
        using $(package);
        runtestpath = joinpath(dirname(dirname(pathof($(package)))), "test", "runtests.jl");
        include(runtestpath);
    )
    return quote
        @snoopi_bench $packageName $(snoopScript)
    end
end


text:
---------------------------------------------
# the Rosenbrock-Lerner polynomial
N=60
b=15
@polyvar x[1:N]
f_R=10*sum((x[3:N-1].+x[2:N-2].-x[1:N-3].^2).^2)+sum((ntuple(i->1,N-3).-x[1:N-3].-x[4:N]).^2)
f_Q=sum([x[i]*i/j*x[j] for i=1:b,j=1:b])
f=f_R+f_Q
---------------------------------------------

---------------------------------------------
# the Broyden tridiagonal polynomial
l=5
p=20
n=l*p
@polyvar x[1:n]
f=((3-2*x[1])*x[1]-2*x[2]+1)^2
for i=2:n-1
    global f+=((3-2*x[i])*x[i]-x[i-1]-2*x[i+1]+1)^2
end
f+=((3-2*x[n])*x[n]-x[n-1]+1)^2
pop=[f]
for i=1:l
    push!(pop,1-sum(x[(i-1)*p+1:i*p].^2))
end
---------------------------------------------

---------------------------------------------
# the Chained singular polynomial
l=2
p=15
n=l*p
@polyvar x[1:n]
f=(x[1]+10*x[2])^2+5*(x[3]-x[4])^2+(x[2]-2*x[3])^4+10*(x[1]-10*x[4])^4
for i=3:2:n-3
    global f+=(x[i]+10*x[i+1])^2+5*(x[i+2]-x[i+3])^2+(x[i+1]-2*x[i+2])^4+10*(x[i]-10*x[i+3])^4
end
pop=[f]
for i=1:l
    push!(pop,1-sum(x[(i-1)*p+1:i*p].^2))
end
---------------------------------------------

---------------------------------------------
# the Broyden banded polynomial
n=400
@polyvar x[1:n]
jset=Array{Any}(undef, n)
f=x[1]
for i=1:n
    jset[i]=max(1,i-5):min(n,i+1)
    jset[i]=setdiff(jset[i],i)
    f0=(1+x[jset[i][1]])*x[jset[i][1]]
    ljset=length(jset[i])
    for j=2:ljset
        f0+=(1+x[jset[i][j]])*x[jset[i][j]]
    end
    global f+=(x[i]*(2+5*x[i]^2)+1-f0)^2
end
f-=x[1]
---------------------------------------------

---------------------------------------------
# the chained Wood polynomial
l=50
p=20
n=l*p
@polyvar x[1:n]
f=x[1]+1
for i=1:2:n-3
    global f+=100*(x[i+1]-x[i]^2)^2+(1-x[i])^2+90*(x[i+3]-x[i+2]^2)^2+(1-x[i+2])^2+10*(x[i+1]+x[i+3]-2)^2+0.1*(x[i+1]-x[i+3])^2
end
f-=x[1]
pop=[f]
for i=1:l
    push!(pop,1-sum(x[(i-1)*p+1:i*p].^2))
end
---------------------------------------------

---------------------------------------------
# the generalized Rosenbrock polynomial
l=50
p=20
n=l*p
@polyvar x[1:n]
f=x[1]+1
for i=2:n
    global f+=100*(x[i]-x[i-1]^2)^2+(1-x[i])^2
end
f-=x[1]
pop=[f]
for i=1:l
    push!(pop,1-sum(x[(i-1)*p+1:i*p].^2))
end
---------------------------------------------


text:
#!/usr/bin/env julia
using JSON
using AWSSDK.SQS
using AWSCore
using ProgressMeter
using ArgParse

const AWS_CREDENTIAL = AWSCore.aws_config()

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--jsonfile", "-j"
            help = "the id list was contained in the consensus json file"
            arg_type = String
        "--sqsqueue", "-q"
            help = "AWS SQS queue name"
            arg_type = String
            default = "skeleton"
        "--skip", "-s"
            help="skip some huge neurons which will explode memory. example: 34,45,78"
            arg_type = String
            default = ""
        "--skip-dir", "-d"
            help="skip existing file"
            arg_type=String
            default = ""
    end
    return parse_args(s)
end

function main()

    args = parse_commandline()
    @show args

    skipNeuronIdSet = map(Meta.parse, split(args["skip"], ",")) |> Set{Int}
    println("skipping some neuron ids: ", skipNeuronIdSet)

    d = JSON.parsefile(args["jsonfile"] |> expanduser)
    neuronIdSet = Set{Int}()
    for v in d |> values
        for neuronIdStr in keys(v)
            neuronId = Meta.parse(neuronIdStr)
            if !(neuronId in skipNeuronIdSet)
                push!(neuronIdSet, neuronId)
            end
        end
    end

    messageList = Vector{String}()
    skipDir = expanduser(args["skip-dir"])
    for id in neuronIdSet
        if !isempty(args["skip-dir"]) &&
            (isfile(joinpath(skipDir, "$(id).swc")) ||
            isfile(joinpath(skipDir, "$(id).swc.bin")) ||
            isfile(joinpath(skipDir, "$(id).bin")))
            # this neuron was already skeletonized
            continue
        end
        push!(messageList, string(id))
    end

    # put to AWS SQS queue
    queueUrl = SQS.get_queue_url( QueueName=args["sqsqueue"] )["QueueUrl"]
    println("SQS queue url: ", queueUrl)

#    for message in messageList
#        @show message
#        SQS.send_message( QueueUrl=queueUrl, MessageBody=message)
#    end

    @showprogress 1 "producing tasks..." for i in 1:10:length(messageList)
        messages = messageList[i:min(i+9, length(messageList))]
        messageBatch = map((x,y)->["Id"=>string(i+x-1), "MessageBody"=>y],
                                                1:length(messages), messages)
        SQS.send_message_batch(AWS_CREDENTIAL; QueueUrl=queueUrl,
                               SendMessageBatchRequestEntry=messageBatch)
    end

    nothing
end

main()