text:
function AH_eq_coeff(vars::Tuple, ::Outer)
    (
        sigma0, sigma0_x, sigma0_y, sigma0_xx, sigma0_yy, sigma0_xy,
        xi    , xi_x    , xi_y    , xi_xx    , xi_yy, xi_xy,
        B1   , B2   , G   ,  S    , Fx    , Fy    , Sd ,
        B1p  , B2p  , Gp  ,  Sp   , Fxp   , Fyp   , Sdp,
        B1pp , B2pp , Gpp ,  Spp  , Fxpp  , Fypp  ,
        B1_x , B2_x , G_x ,  S_x  , Fx_x  , Fy_x  , Sd_x,
	B1_y , B2_y , G_y ,  S_y  , Fx_y  , Fy_y  , Sd_y,
        B1p_x, B2p_x, Gp_x,  Sp_x , Fxp_x , Fyp_x ,
        B1p_y, B2p_y, Gp_y,  Sp_y , Fxp_y , Fyp_y ,
    ) = vars

    @tilde_outer("B1")
    @tilde_outer("B2")
    @tilde_outer("G")
    @tilde_outer("S")
    @tilde_outer("Fx")
    @tilde_outer("Fy")
    @tilde_outer("Sd")

    @hat_outer("B1")
    @hat_outer("B2")
    @hat_outer("G")
    @hat_outer("S")
    @hat_outer("Fx")
    @hat_outer("Fy")
    @hat_outer("Sd")

    @tilde_outer("B1p")
    @tilde_outer("B2p")
    @tilde_outer("Gp")
    @tilde_outer("Sp")
    @tilde_outer("Fxp")
    @tilde_outer("Fyp")

    @hat_outer("B1p")
    @hat_outer("B2p")
    @hat_outer("Gp")
    @hat_outer("Sp")
    @hat_outer("Fxp")
    @hat_outer("Fyp")

    x0 = exp(B1 + B2)
    x1 = cosh(G)
    x2 = S*x1
    x3 = sinh(G)
    x4 = exp(B2)
    x5 = x3*x4
    x6 = 2*x5
    x7 = -B1
    x8 = exp(B2 + x7)
    x9 = B2h*S
    x10 = Fyp*S
    x11 = Sp*sigma0_y
    x12 = Sp*xi_y
    x13 = 2*B2p
    x14 = S*x13
    x15 = sigma0_y*x14
    x16 = x14*xi_y
    x17 = exp(B1)
    x18 = S*x17
    x19 = Gt*x18
    x20 = Fx*Gp
    x21 = 2*x18*x20
    x22 = 2*sigma0_x
    x23 = Gp*x18
    x24 = x22*x23
    x25 = 2*xi_x
    x26 = x23*x25
    x27 = Fy + xi_y
    x28 = sigma0_y + x27
    x29 = Gp*x28
    x30 = S*(Gh + 2*x29)
    x31 = 2*B1p
    x32 = x13 + x31
    x33 = Fx*S
    x34 = S*x31
    x35 = Fx*Sp
    x36 = B2t*S + Fxp*S - Sp*sigma0_x - Sp*xi_x + St + sigma0_x*x14 + x14*xi_x - x35
    x37 = -Sh
    x38 = Sp*x28
    x39 = x1*x4
    x40 = x28^2
    x41 = x39*x40
    x42 = 3*Spp
    x43 = 3*B2p
    x44 = 3*Sp
    x45 = sigma0_x + xi_x
    x46 = Fx + x45
    x47 = 6*Sp
    x48 = exp(2*B1 + B2)
    x49 = x1*x48
    x50 = x46^2
    x51 = x49*x50
    x52 = Gp*x3
    x53 = Sp*x45
    x54 = St - 2*x35 - 2*x53
    x55 = x17*x3
    x56 = x54*x55
    x57 = B1h + B1p*x28
    x58 = B2p*x28
    x59 = B2h + x58
    x60 = x1*x59
    x61 = Gh + x29
    x62 = x3*x61
    x63 = B2p*x46 + B2t
    x64 = x55*x63
    x65 = Gp*x45 + Gt + x20
    x66 = x1*x17
    x67 = x65*x66
    x68 = x1*x57 - x60 - x62 + x64 + x67
    x69 = S*x68 + x1*(x37 - x38) + x56
    x70 = x4*x69
    x71 = Sh + x38
    x72 = x3*x71
    x73 = St + x35 + x53
    x74 = x66*x73
    x75 = B1p*x46 + B1t
    x76 = x1*x63
    x77 = x3*x65
    x78 = x17*(x1*x75 + x76 + x77)
    x79 = x1*(-Gh - x29) + x3*(-B2h - x58) + x78
    x80 = S*x79 - x72 + x74
    x81 = x0*x80
    x82 = Fxp*x28
    x83 = 2*sigma0_xy + 2*xi_xy
    x84 = x0*x3
    x85 = Fyp*sigma0_y
    x86 = Fyp*x27
    x87 = 6*S
    x88 = Fyp*sigma0_x
    x89 = Fx + xi_x
    x90 = Fyp*x89
    x91 = Fxp*sigma0_x
    x92 = Fxp*x89
    x93 = x0*x1
    x94 = x17*(Sd*x87 + x5*(2*Fyt + x83 + 2*x88 + 2*x90) - x93*(2*Fxt + 2*sigma0_xx + 2*x91 + 2*x92 + 2*xi_xx))
    x95 = -x39*(2*Fyh + 2*sigma0_yy + 2*x85 + 2*x86 + 2*xi_yy) + x84*(2*Fxh + 2*x82 + x83) + x94
    x96 = 2*x28
    x97 = Fyh + sigma0_yy + x85 + x86 + xi_yy
    x98 = sigma0_xy + xi_xy
    x99 = Fxh + x82 + x98
    x100 = 2*x93
    x101 = Fyt + x88 + x90 + x98
    x102 = 2*Gp
    x103 = Fxt + sigma0_xx + x91 + x92 + xi_xx
    x104 = Sph + Spp*x28
    x105 = Gp*x1
    x106 = B2ph + B2pp*x28
    x107 = Gph + Gpp*x28
    x108 = B2pp*x46 + B2pt
    x109 = Fx*Gpp + Gpp*x45 + Gpt
    x110 = 2*Fx
    x111 = x0*(x110 + x22 + x25)
    x112 = Gp*x17
    x113 = x4*(2*Fy + 2*sigma0_y + 2*xi_y)

    axx = -x0*x2

    axy = S*x6

    ayy = -x2*x8

    bx = x4*(x1*(-x17*(B1t*S + sigma0_x*x34 + x32*x33 + x34*xi_x + x36) + x30) + x3*(-Fy*Sp + Fy*x14 + Sh + x10 - x11 - x12 + x15 + x16 - x19 - x21 - x24 - x26 + x9))

    by = x8*(x1*(B1h*S + Fy*(Sp - x14 + x34) + sigma0_y*x34 - x10 + x11 + x12 - x15 - x16 + x19 + x21 + x24 + x26 + x34*xi_y + x37 - x9) + x3*(x17*(x13*x33 + x36) - x30))

    cc = (-B1p*(S*x95 - x111*x80 + x113*x69 + x41*x44 + x44*x51) + Fxp*x46*x47*x49 - 2*Fxp*x81 + 6*Fyp*x38*x39 + 2*Fyp*x70 + Gp*x40*x44*x5 + S*(B1p*x94 + Gp*x100*x99 - Gp*x6*x97 - x13*x39*x97 + x17*(Sd*x47 + Sdp*x87 - x100*x103*(B1p + B2p) - x100*(Fx*Fxpp + Fxpp*x45 + Fxpt) + x101*x102*x39 + x101*x13*x5 - x102*x103*x84 + x6*(Fx*Fypp + Fypp*x45 + Fypt)) + x32*x84*x99 - x39*(2*Fyph + Fypp*x96) + x84*(2*Fxph + Fxpp*x96)) + Sp*x41*x43 + Sp*x51*(6*B1p + x43) + Sp*x95 - x111*(B1p*x74 + Gp*x55*x73 + S*(B1p*x78 - Gp*x60 - Gp*x62 - x1*x107 - x106*x3 + x17*(x1*x108 + x1*(B1pp*x46 + B1pt) + x105*x65 + x109*x3 + x52*x63 + x52*x75)) + Sp*x79 - x104*x3 - x105*x71 + x66*(Fx*Spp + Spp*x45 + Spt)) + x113*(B1p*x56 + Gp*x54*x66 - Gp*x72 + S*(B1p*x64 + B1p*x67 - x1*x106 + x1*(B1ph + B1pp*x28) - x105*x61 - x107*x3 + x108*x55 + x109*x66 + x112*x76 + x112*x77 + x52*x57 - x52*x59) + Sp*x68 - x1*x104 - x55*(Fxp*x44 + Spp*x110 + Spp*x22 + Spp*x25 - Spt)) + x13*x28*x70 - x32*x46*x81 + x41*x42 + x42*x51 + x44*x48*x50*x52)*exp(x7)/2

    return axx, ayy, axy, bx, by, cc
end

function AH_eq_res(vars::Tuple, ::Outer)
    (
        sigma0, sigma0_x, sigma0_y, sigma0_xx, sigma0_yy, sigma0_xy,
        xi    , xi_x    , xi_y    , xi_xx    , xi_yy, xi_xy,
        B1   , B2   , G   ,  S    , Fx    , Fy    , Sd ,
        B1p  , B2p  , Gp  ,  Sp   , Fxp   , Fyp   , Sdp,
        B1pp , B2pp , Gpp ,  Spp  , Fxpp  , Fypp  ,
        B1_x , B2_x , G_x ,  S_x  , Fx_x  , Fy_x  , Sd_x,
        B1_y , B2_y , G_y ,  S_y  , Fx_y  , Fy_y  , Sd_y,
    ) = vars

    @tilde_outer("B1")
    @tilde_outer("B2")
    @tilde_outer("G")
    @tilde_outer("S")
    @tilde_outer("Fx")
    @tilde_outer("Fy")
    @tilde_outer("Sd")

    @hat_outer("B1")
    @hat_outer("B2")
    @hat_outer("G")
    @hat_outer("S")
    @hat_outer("Fx")
    @hat_outer("Fy")
    @hat_outer("Sd")

    x0 = Fy + xi_y
    x1 = sigma0_y + x0
    x2 = cosh(G)
    x3 = exp(B2)
    x4 = x2*x3
    x5 = 3*Sp
    x6 = Fx + xi_x
    x7 = sigma0_x + x6
    x8 = 2*sigma0_y
    x9 = Sp*x1
    x10 = 2*Fx
    x11 = sigma0_x + xi_x
    x12 = Sp*x11
    x13 = exp(B1)
    x14 = sinh(G)
    x15 = x13*x14
    x16 = B2p*x1
    x17 = Gp*x1
    x18 = B2p*x7 + B2t
    x19 = Fx*Gp + Gp*x11 + Gt
    x20 = x13*x2
    x21 = exp(B1 + B2)
    x22 = 2*sigma0_xy + 2*xi_xy
    x23 = 2*sigma0_x
    x24 = 2*x6

    (S*(x13*(6*S*Sd + x14*x3*(Fyp*x23 + Fyp*x24 + 2*Fyt + x22) - x2*x21*(Fxp*x23 + Fxp*x24 + 2*Fxt + 2*sigma0_xx + 2*xi_xx)) + x14*x21*(2*Fxh + 2*Fxp*x1 + x22) - x4*(2*Fyh + 2*Fyp*x0 + Fyp*x8 + 2*sigma0_yy + 2*xi_yy)) + x1^2*x4*x5 + x2*x5*x7^2*exp(2*B1 + B2) - x21*(x10 + x23 + 2*xi_x)*(S*(x13*(x14*x19 + x18*x2 + x2*(B1p*x7 + B1t)) + x14*(-B2h - x16) + x2*(-Gh - x17)) - x14*(Sh + x9) + x20*(Fx*Sp + St + x12)) + x3*(2*Fy + x8 + 2*xi_y)*(S*(-x14*(Gh + x17) + x15*x18 + x19*x20 + x2*(B1h + B1p*x1) - x2*(B2h + x16)) + x15*(-Sp*x10 + St - 2*x12) + x2*(-Sh - x9)))*exp(-B1)/2

end


text:
using RareEventABC
using Plots; pgfplots()
default(guidefont=font(25), tickfont=font(25), legendfont=font(25), size=(500,500));
using StatPlots
using ProgressMeter
using DataFrames

##ABAKALIKI DATA
size_obs = 30;
popsize = 120;
robs = cumsum([0,13,7,2,3,0,0,1,4,5,3,2,0,2,0,5,3,1,4,0,1,1,1,2,0,1,5,0,5,5]);
duration_obs = maximum(robs);
robs = [robs; [Inf for i in 1:90]];

##MODELS
using RareEventABC
using Distributions

##Markoivan SIR binned to 5 day periods
sir_binned = RareEventABC.SEIRDetails(robs, popsize, θ->RareEventABC.Exponential(1.), θ->RareEventABC.DiracDelta(0.0), θ->Distributions.Exponential(θ[2]), x->floor(x/5.)*5.);

##Gamma infectious period
sir_gamma_I = RareEventABC.SEIRDetails(robs, popsize, θ->RareEventABC.Exponential(1.), θ->RareEventABC.DiracDelta(0.0), θ->Distributions.Gamma(θ[2], θ[3]), x->x);

##Weibull pressure thresholds
sir_weibull_P = RareEventABC.SEIRDetails(robs, popsize, θ->RareEventABC.Weibull(θ[2], 1.), θ->RareEventABC.DiracDelta(0.0), θ->Distributions.Exponential(θ[3]), x->x);

prior_2 = MvIndependent([Exponential(10.), Exponential(10.)]); ##Corresponds to Kypraios et al 2016 prior (n.b. Paper specified Exp(0.1) under rate parameterisation. Julia uses scale parameterisation.)
prior_3 = MvIndependent([Exponential(10.), Exponential(10.), Exponential(10.)]); ##Corresponds to Kypraios et al 2016 prior (n.b. Paper specified Exp(0.1) under rate parameterisation. Julia uses scale parameterisation.)

eps = 15.0;

priors = Distribution[prior_2, prior_3, prior_3];
models = RareEventABC.SEIRDetails[sir_binned, sir_gamma_I, sir_weibull_P];

##Do a pilot run, display some diagnostics, and return suggestions for main run: schedule, start point and random walk increment distribution
function pilot_run(model::RareEventABC.SEIRDetails, prior::Distribution, increment::Distribution, pars_init::Array{Float64, 1}, nparticles::Int; niterations=200)
    npars = length(pars_init)
    ##Run pilot
    schedule = doSMC(model, pars_init, eps, nparticles=nparticles, d=size_obs, quiet=true).schedule
    schedule = unique(schedule)
    pars, lls = doPMMH(model, schedule, niterations, nparticles, prior, pars_init, increment)
    ##Plot output
    plot(lls, line=:steppost) |> display
    plot(pars, line=:steppost, layout=npars) |> display
    ##scatter(pars[:,1], pars[:,2]) |> display
    ##Diagnostics
    pilot_mean = vec(mean(pars, 1))
    print("Mean state $(round(pilot_mean, 2))\n")
    post_cov = cov(pars[100:niterations,:]) * (2.562^2 / 2)
    print("Recommended RW standard deviations $(round(sqrt(diag(post_cov)), 2))\n")
    print("Estimating log likelihood standard deviation...\n")
    llests = @showprogress Float64[doSMC(model, pilot_mean, schedule, nparticles=nparticles, d=size_obs, quiet=true).llest for i in 1:50];
    print("Estimated log likelihood standard deviation $(round(std(llests), 2)) or $(round(std(llests[llests .> -Inf]), 2)) (removing $(sum(isinf(llests))) zero estimates)\n")
    return schedule, pilot_mean, MvNormal(zeros(npars), post_cov)
end

##Pilot runs. Not automated across models as tuning parameters are tweaked manually.
schedules = Array(Array{Float64, 1}, 3);
starts = Array(Array{Float64, 1}, 3);
increments = Array(Distribution, 3);

srand(1);
schedules[1], starts[1], increments[1] = pilot_run(sir_binned, prior_2,
                                                     MvNormal(zeros(2), [0.05, 4.0]), ##RW increment
                                                     [0.1, 10.0], ##Initial parameters
                                                     400); ##Number of particles

srand(1);
schedules[2], starts[2], increments[2] = pilot_run(sir_gamma_I, prior_3,
                                                     MvNormal(zeros(3), [0.02, 4.0, 2.0]), ##RW increment
                                                     [0.1, 7.0, 2.5], ##Initial parameters
                                                     300); ##Number of particles

srand(1);
schedules[3], starts[3], increments[3] = pilot_run(sir_weibull_P, prior_3,
                                                     MvNormal(zeros(3), [0.05, 0.2, 2.0]), ##RW increment
                                                     [0.1, 1.0, 14.0], ##Initial parameters
                                                     200); ##Number of particles

main_nparticles = Int[400,300,200];

function main_run(model::RareEventABC.SEIRDetails, prior::Distribution, schedule::Array{Float64, 1}, increment::Distribution, pars_init::Array{Float64, 1}, nparticles::Int; niterations=2000)
    npars = length(pars_init)
    ##PMMH run
    pars, lls = doPMMH(model, schedule, niterations, nparticles, prior, pars_init, increment);
    ##Plot output
    plot(lls, line=:steppost) |> display
    plot(pars, line=:steppost, layout=npars) |> display
    ##Diagnostics
    print("Parameter means $(round(mean(pars, 1),2)) and standard deviations $(round(std(pars, 1),2))\n")
    return pars, lls
end

main_pars = Array(Array{Float64, 2}, 3);
main_lls = Array(Float64, (2000, 3));

srand(1);
for i in 1:3
    main_pars[i], main_lls[:,i] = main_run(models[i], priors[i], schedules[i], increments[i], vec(starts[i]), main_nparticles[i])
end

##PARAMETER ESTIMATES
for i in 1:3
    print("Means: $(round(mean(main_pars[i], 1), 3)) \n Sdev: $(round(std(main_pars[i], 1), 3)) \n")
end

##R0 ESTIMATES
R0_1 = main_pars[1][:,1] .* main_pars[1][:,2];
mean(R0_1), std(R0_1)

R0_2 = main_pars[2][:,1] .* main_pars[2][:,2] .* main_pars[2][:,3];
mean(R0_2), std(R0_2)

##ESTIMATE OF OTHER QUANTITIES
temp = main_pars[1][:,1]; ##Mean and standard deviation of pressure thresholds
mean(temp), std(temp)
temp = main_pars[1][:,2]; ##Mean and standard deviation of infectious period
mean(temp), std(temp)
temp = main_pars[2][:,1]; ##Mean and standard deviation of pressure thresholds
mean(temp), std(temp)
temp = main_pars[2][:,2] .* main_pars[2][:,3]; ##Mean of infectious period
mean(temp), std(temp)
temp = sqrt(main_pars[2][:,2]) .* main_pars[2][:,3]; ##Standard deviation of infectious period
mean(temp), std(temp)
temp = [mean(Weibull(main_pars[3][i,2],main_pars[3][i,1])) for i in 1:2000]; ##Mean of pressure thresholds
mean(temp), std(temp)
temp = [std(Weibull(main_pars[3][i,2],main_pars[3][i,1])) for i in 1:2000]; ##Standard deviation of pressure thresholds
mean(temp), std(temp)
temp = main_pars[3][:,3]; ##Mean and standard deviation of infectious period
mean(temp), std(temp)

##PLOT LOG LIKELIHOODS
df_ll = vcat(
             DataFrame(iteration=1:2000, loglikelihood=main_lls[:,2], model="Gamma infectious period"),
             DataFrame(iteration=1:2000, loglikelihood=main_lls[:,3], model="Weibull pressure thresholds"),
             DataFrame(iteration=1:2000, loglikelihood=main_lls[:,1], model="5 day bins")
);
p1 = plot(df_ll, :iteration, :loglikelihood, group=:model,
     xaxis=("Iteration"), yaxis=("Log likelihood"), legend=:bottomright);

savefig(p1, "Abakaliki_lls.pdf")

##PLOT EXAMPLE EPIDEMICS
function plot_example_sims(model::RareEventABC.SEIRDetails, schedule::Array{Float64, 1}, nparticles::Int, pars::Array{Float64, 2}; title="")
    ##Plot some example simulations
    rData = DataFrame(Day = robs[1:30], R = 1:30);
    simData = DataFrame(Day=Float64[], R=Int[], Rep=Int[]);
    i = 1
    while i<=10
        p = doSMC(model, vec(pars[100*i,:]), schedule, nparticles=nparticles, quiet=true);
        hits = RareEventABC.getDistances(model, p) .<= p.ε
        if sum(hits) > 0
            j = sample((1:p.nparticles)[hits])
            x = p.epidemics[j].compartmentTimes[:,3]
            x = x[!isinf(x)]; x = sort(x); x .-= x[1]
            simData = vcat(simData, DataFrame(Day = x, R=1:length(x), Rep=string(i)))
            i += 1
        end
    end
    default(guidefont=font(25), tickfont=font(25), legendfont=font(25), size=(500,500));
    p1=plot(rData, :Day, :R, line=:steppost, linewidth=3, linecolour=:black, title=title, legend=false);
    p1=plot!(simData, :Day, :R, group=:Rep, line=:steppost);
    display(p1)
    return p1
end

p2 = plot_example_sims(models[1], schedules[1], main_nparticles[1], main_pars[1], title="5 day bins");
p3 = plot_example_sims(models[2], schedules[2], main_nparticles[2], main_pars[2], title="Gamma infectious period");
p4 = plot_example_sims(models[3], schedules[3], main_nparticles[3], main_pars[3], title="Weibull pressure thresholds");
p5 = plot(p2, p3, p4, layout=(1,3), size=(1000, 500));
savefig(p5, "example_epidemics.pdf")

ext: jl


text:
struct EmptyBridge{T} <: MOIB.Constraint.AbstractBridge
end

function MOIB.Constraint.bridge_constraint(
    ::Type{EmptyBridge{T}}, model::MOI.ModelLike, f::MOI.AbstractVectorFunction,
    s::SOS.EmptyCone) where {T}
    @assert MOI.output_dimension(f) == MOI.dimension(s)
    return EmptyBridge{T}()
end

function MOI.supports_constraint(
    ::Type{<:EmptyBridge}, ::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:SOS.EmptyCone})
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:EmptyBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:EmptyBridge})
    return Tuple{DataType, DataType}[]
end
function MOIB.Constraint.concrete_bridge_type(
    ::Type{<:EmptyBridge{T}}, ::Type{<:MOI.AbstractVectorFunction},
    ::Type{SOS.EmptyCone}) where T
    return EmptyBridge{T}
end

# Indices
function MOI.delete(model::MOI.ModelLike, bridge::EmptyBridge) end

function MOI.get(::MOI.ModelLike,
                 ::Union{MOI.ConstraintDual, MOI.ConstraintPrimal},
                 bridge::EmptyBridge{T}) where T
    return T[]
end


text:
##
OSBORNEA

CERI651CLS
DEVGLA1
GAUSS2LS
NELSONLS
OSBORNEA

##
S = [1 2 3 4 5; 6 7 8 9 10; 11 12 13 14 15; 16 17 18 19 20]
##
m = 8
for i = 0:m
    print(S[:,mod(i,m)+1])
end
##
S


text:
# ---
# title: 167. Two Sum II - Input array is sorted
# id: problem167
# author: zhwang
# date: 2022-01-19
# difficulty: Easy
# categories: Array, Two Pointers, Binary Search
# link: <https://leetcode.com/problems/two-sum-ii-input-array-is-sorted/description/>
# hidden: true
# ---
#
# Given an array of integers that is already **_sorted in ascending order_** ,
# find two numbers such that they add up to a specific target number.
#
# The function twoSum should return indices of the two numbers such that they
# add up to the target, where index1 must be less than index2.
#
# **Note:**
#
#   * Your returned answers (both index1 and index2) are not zero-based.
#   * You may assume that each input would have _exactly_ one solution and you may not use the _same_ element twice.
#
#
#
# **Example 1:**
#
#
#
#     Input: numbers = [2,7,11,15], target = 9
#     Output: [1,2]
#     Explanation: The sum of 2 and 7 is 9. Therefore index1 = 1, index2 = 2.
#
#
# **Example 2:**
#
#
#
#     Input: numbers = [2,3,4], target = 6
#     Output: [1,3]
#
#
# **Example 3:**
#
#
#
#     Input: numbers = [-1,0], target = -1
#     Output: [1,2]
#
#
#
#
# **Constraints:**
#
#   * `2 <= nums.length <= 3 * 104`
#   * `-1000 <= nums[i] <= 1000`
#   * `nums` is sorted in **increasing order**.
#   * `-1000 <= target <= 1000`
#
#
## @lc code=start
using LeetCode

## using double pointers
function two_sum_167(numbers::Vector{Int}, target::Int)::Vector{Int}
    left, right = 1, length(numbers)
    while left < right
        num = numbers[left] + numbers[right]
        if num > target
            right -= 1
        elseif num < target
            left += 1
        else
            return [left, right]
        end
    end
    [-1, -1]
end

## @lc code=end


text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random
path1=realpath(dirname(@__FILE__)*"/../..")
include(string(path1,"/src/ANNBN.jl"))

# Polynomial
function ff(x)
    -x[1]^3+5x[1]^2+5x[1]*x[2]
end

i_train=1000
vars=2
xx_train=0.8*rand(i_train,vars).+0.1
yy_train=zeros(i_train); for i=1:i_train yy_train[i]=ff(xx_train[i,:]) end
# mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=0.8;yy_train.+=0.1
i_test=1000
xx_test=0.8*rand(i_test,vars).+0.1
yy_test=zeros(i_test); for i=1:i_test yy_test[i]=ff(xx_test[i,:]) end
# yy_test.-=mi;yy_test./=ma;yy_test.*=0.8;yy_test.+=0.1

neurons=Int64(floor(i_train/(vars+1)))
neurons=100
inds_all,n_per_part=ann_by_parts.___clustering(neurons,xx_train,200)
cc1=1.0
# Use # 2 r^4
a_all,a_layer1,layer1=ann_by_parts.train_layer_1_rbf(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train,cc1)
predl1=ann_by_parts.predict_new_rbf(a_all,a_layer1,xx_train,i_train,neurons,n_per_part,inds_all,xx_train,vars,cc1)
predl1_test=ann_by_parts.predict_new_rbf(a_all,a_layer1,xx_test,i_test,neurons,n_per_part,inds_all,xx_train,vars,cc1)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))

##### APPROXIMATION OF DERIVATIVES
# first variable
yy_train_deriv=-3xx_train[:,1].^2+10xx_train[:,1]+5xx_train[:,2]
yy_train_deriv_deriv=-6xx_train[:,1]+10ones(i_train)
yy_test_deriv=-3xx_test[:,1].^2+10xx_test[:,1]+5xx_test[:,2]
yy_test_deriv_deriv=-6xx_test[:,1]+10ones(i_test)

predl1_deriv=ann_by_parts.predict_new_rbf_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,1)
mean(abs.(predl1_deriv-yy_train_deriv))
predl1_deriv_deriv=ann_by_parts.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,1)
mean(abs.(predl1_deriv_deriv-yy_train_deriv_deriv))
predl1_test_deriv=ann_by_parts.predict_new_rbf_deriv(a_all,a_layer1,xx_test,i_test,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,1)
mean(abs.(predl1_test_deriv-yy_test_deriv))
predl1_test_deriv_deriv=ann_by_parts.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_test,i_test,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,1)
mean(abs.(predl1_test_deriv_deriv-yy_test_deriv_deriv))

# second variable
yy_train_deriv=5xx_train[:,1]
yy_train_deriv_deriv=zeros(i_train)
yy_test_deriv=5xx_test[:,1]
yy_test_deriv_deriv=zeros(i_test)

predl1_deriv=ann_by_parts.predict_new_rbf_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,2)
mean(abs.(predl1_deriv-yy_train_deriv))
predl1_deriv_deriv=ann_by_parts.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,2)
mean(abs.(predl1_deriv_deriv-yy_train_deriv_deriv))

predl1_test_deriv=ann_by_parts.predict_new_rbf_deriv(a_all,a_layer1,xx_test,i_test,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,2)
mean(abs.(predl1_test_deriv-yy_test_deriv))
predl1_test_deriv_deriv=ann_by_parts.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_test,i_test,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,2)
mean(abs.(predl1_test_deriv_deriv-yy_test_deriv_deriv))


text:
abstract type GenericMomentBasedEstimator <: MathProgBase.AbstractNLPEvaluator end
abstract type Constraint end
abstract type Weighting end

abstract type MomentFunction end

abstract type AnaMomFun <: MomentFunction end

struct FADMomFun{F1, F2, K} <: MomentFunction
    g::F1            ## Moment Function
    s::F2            ## Smoothed moment function
    kern::K
end

struct AnaGradMomFun{F1, F2, F3, F4, F5, K} <: AnaMomFun
    g::F1            ## Moment Function
    s::F2            ## Smoothed moment function
    Dsn::F3
    Dws::F4
    Dsl::F5
    kern::K
end

struct AnaFullMomFun{F1, F2, F3, F4, F5, F6, K} <: AnaMomFun
    g::F1            ## Moment fun
    s::F2            ## Smoothed moment fun
    Dsn::F3
    Dws::F4
    Dsl::F5
    Hwsl::F6
    kern::K
end

mutable struct MomentBasedEstimatorResults
    status::Symbol
    objval::Float64
    coef::Array{Float64, 1}
    H::Array{Float64,2}      ## Hessian of the objective function
end

mutable struct GMMEstimator{M<:MomentFunction, V<:IterationManager, S<:Constraint, T<:Weighting} <: GenericMomentBasedEstimator
    mf::M
    c::S
    x0::Array{Float64, 1}
    lb::Array{Float64, 1}
    ub::Array{Float64, 1}
    glb::Array{Float64, 1}
    gub::Array{Float64, 1}
    mgr::V
    ist::IterationState
    W::Array{Array{Float64,2},1}
    wtg::T
    gele::Int64
    hele::Int64
    nobs::Int64
    npar::Int64
    nmom::Int64
end

mutable struct MDEstimator{M<:MomentFunction, V<:Divergence, S<:Constraint, T<:Weighting} <: GenericMomentBasedEstimator
    mf::M
    c::S
    x0::Array{Float64, 1}
    lb::Array{Float64, 1}
    ub::Array{Float64, 1}
    glb::Array{Float64, 1}
    gub::Array{Float64, 1}
    wlb::Array{Float64, 1}
    wub::Array{Float64, 1}
    div::V
    wtg::T
    gele::Int64
    hele::Int64
    nobs::Int64
    npar::Int64
    nmom::Int64
end

mutable struct MomentBasedEstimatorOptions
    ##options
    ##maybe optimization options?
end

struct Unweighted <: Weighting end

struct Weighted <: Weighting
    wtg::StatsBase.Weights{Float64}
end

struct Unconstrained <: Constraint end

mutable struct Constrained <: Constraint
    h::Function
    hlb::Array{Float64, 1}
    hub::Array{Float64, 1}
    nc::Int64  ## Number of constraints: row of h(θ)
end

mutable struct MomentBasedEstimator{T<:GenericMomentBasedEstimator, S<:MathProgBase.AbstractMathProgSolver, M<:MathProgBase.AbstractMathProgModel}
    e::T
    r::MomentBasedEstimatorResults
    s::S
    m::M
    status::Vector{Symbol}
end

## Basic MomentBasedEstimator constructor
function MomentBasedEstimator(e::GenericMomentBasedEstimator)
    MomentBasedEstimator(e, MomentBasedEstimatorResults(
                                                        :Unsolved, 0.0,
                                                        Array{Float64}(undef, npar(e)),
                                                        Array{Float64}(undef, npar(e), npar(e))),
                         DEFAULT_SOLVER(e),
                         MathProgBase.NonlinearModel(DEFAULT_SOLVER(e)),
                         [:Uninitialized])
end

setW0(mgr::TwoStepGMM, m::Int64) = [Array{Float64}(undef, m, m) for i=1:2]
setW0(mgr::OneStepGMM, m::Int64) = [Array{Float64}(undef, m, m) for i=1:1]
setW0(mgr::IterativeGMM, m::Int64) = [Array{Float64}(undef, m, m) for i=1:mgr.maxiter+1]

function make_fad_mom_fun(g::Function,
                          kernel::SmoothingKernel = IdentitySmoother())
    FADMomFun(g, θ -> smooth(g(θ), kernel), kernel)
end

function make_ana_mom_fun(::Type{GMMEstimator}, g::Function, ∇g::Function)
    AnaGradMomFun(g, g, ∇g, identity, identity, IdentitySmoother())
end

function make_ana_mom_fun(::Type{MDEstimator}, g::Function, ∇g::Tuple{Function, Function, Function})
    AnaGradMomFun(g, g, ∇g..., IdentitySmoother())
end

function make_ana_mom_fun(::Type{MDEstimator}, g::Function, ∇g::Tuple{Function, Function, Function, Function})
    AnaFullMomFun(g, g, ∇g..., IdentitySmoother())
end

mutable struct DEFAULT_SOLVER{T <: GenericMomentBasedEstimator}
    s::MathProgBase.AbstractMathProgSolver
end

function DEFAULT_SOLVER(::GMMEstimator)
  IpoptSolver(hessian_approximation = "limited-memory", print_level=2, sb = "yes")
end

DEFAULT_SOLVER(::MDEstimator)  = IpoptSolver(print_level=2, sb = "yes")


text:
@testset "convexhull.jl" begin

end


text:
##
## rand for Flint-Finite fields
##
## fq_nmod has no base ring
function rand(R::FqNmodFiniteField)
  #gen is not a generator for the group!
  r = R(0)
  for i=0:degree(R)
    r *= gen(R)
    r += rand(1:characteristic(R))
  end

  return r
end

function rand(R::FqFiniteField)
  #gen is not a generator for the group!
  r = R(0)
  for i=0:degree(R)
    r *= gen(R)
    r += rand(1:characteristic(R))
  end

  return r
end

################################################################################
#
#  Iterators for finite fields
#
################################################################################

# FqNmodFiniteField

Base.start(F::FqNmodFiniteField) = zeros(UInt, degree(F))

function Base.next(F::FqNmodFiniteField, i::Array{UInt, 1})
  d = F()
  ccall((:fq_nmod_init2, :libflint), Void,
        (Ptr{fq_nmod}, Ptr{FqNmodFiniteField}), &d, &F)
  for j in 1:length(i)
         ccall((:nmod_poly_set_coeff_ui, :libflint), Void,
              (Ptr{fq_nmod}, Int, UInt), &d, j - 1, i[j])
  end

  i[1] = i[1] + 1
  for j in 1:(length(i) - 1)
    if i[j] == F.n
      i[j] = 0
      i[j + 1] = i[j + 1] + 1
    end
  end
  return d, i
end

Base.length(F::FqNmodFiniteField) = BigInt(F.n)^degree(F)

function Base.done(F::FqNmodFiniteField, i::Array{UInt, 1})
  for j in 1:length(i) - 1
    if i[j] != 0
      return false
    end
  end
  if i[end] != F.n
    return false
  end
  return true
end

Base.eltype(::FqNmodFiniteField) = fq_nmod

# FqFiniteField

Base.start(F::FqFiniteField) = zeros(FlintZZ, degree(F))

function Base.next(F::FqFiniteField, i::Array{fmpz, 1})
  d = F()
  ccall((:fq_init2, :libflint), Void,
        (Ptr{fq}, Ptr{FqFiniteField}), &d, &F)
  for j in 1:length(i)
         ccall((:fmpz_mod_poly_set_coeff_fmpz, :libflint), Void,
               (Ptr{fq}, Int, Ptr{fmpz}), &d, j - 1, &(i[j]))
  end

  i[1] = i[1] + 1
  for j in 1:(length(i) - 1)
    if i[j] == characteristic(F)
      i[j] = 0
      i[j + 1] = i[j + 1] + 1
    end
  end
  return d, i
end

Base.length(F::FqFiniteField) = BigInt(characteristic(F))^degree(F)

function Base.done(F::FqFiniteField, i::Array{fmpz, 1})
  for j in 1:length(i) - 1
    if i[j] != 0
      return false
    end
  end
  if i[end] != characteristic(F)
    return false
  end
  return true
end

Base.eltype(::FqFiniteField) = fq


text:
# %% Modules

include("DASTDP.jl")
include("NeuronModel.jl")

using .DASTDP
using .NeuronModel
using Statistics

# %% Constants

const S = [rand(1:800,50) for _ in 1:100] # Stimulis Sk
const T = 3600

# %% functions

function stimuli_fire(s_index::Int64,I::Array{Float64})
    ind = S[s_index]
    I[ind] = I[ind] .+ 200  #super threshold current
    return I
end

function means(s::Array{Float64,2})
    mean_cs = mean(s[S[1],:])
    mean_us = mean([mean(s[S[index],:]) for index in 2:50])
    return mean_cs,mean_us
end

function reward(rew::Array{Int64},s_index::Int64,time::Int64)
    if s_index == 1
        append!(rew,time+rand(1:1000))
    end
    return rew
end

# %% Network structure

mutable struct NeuralNet
    v::Array{Float64}
    u::Array{Float64}
    s::Array{Float64,2}
    sd::Array{Float64,2}
    STDP::Array{Float64,2}
    firings::Array{Int64,2}
    DA::Float64
    rew::Array{Int64}
    s_del::Array{Int64}
    I::Array{Float64}
    shist::Array{Float64,2}

    function NeuralNet()
        v = -65.0*ones(N)
        u = 0.2*v
        s = vcat(0.5 .* ones(Ne,M),-0.5 .* ones(Ni,M))
        sd = 0.0 .* zeros(N,M)
        STDP = 0.0 .* zeros(N,1001+D)
        firings = [-D 0]
        DA = 0.0
        rew = []
        s_del = [0,0]
        I = Float64[]
        shist = zeros(1000*T, 2)
        new(v,u,s,sd,STDP,firings,DA,rew,s_del,I,shist)
    end
end

# %% main loop

net = NeuralNet()

for sec in 0:T-1
    for msec in 1:1000
        net.I = 13*(rand(N).-0.5)
        time = 1000*sec+msec
        fired = findall(x->x>=thresh,net.v)
        net.v,net.u = izhikevicmodel_fire(net.v,net.u,fired)
        net.STDP = STDP_fire(net.STDP,fired,msec)
        net.sd = LTP(net.STDP,net.sd,fired,msec)
        net.firings = vcat(net.firings,hcat(msec.*ones(length(fired)),fired))
        net.I,net.sd = LTD(net.STDP,net.sd,net.s,net.firings,net.I,msec)
        if net.s_del[1] == net.s_del[2]
            s_index = rand(1:100)
            net.I = stimuli_fire(s_index,net.I)
            net.s_del = [0,rand(100:300)]
            net.rew = reward(net.rew,s_index,time)
        end
        net.v,net.u = izhikevicmodel_step(net.v,net.u,net.I)
        net.STDP,net.DA = DA_STDP_step(net.STDP,net.DA,msec)
        net.s,net.sd = synweight_step(net.sd,net.s,net.DA,msec)
        net.DA = DA_inc(net.rew,net.DA,time)
        net.shist[time,:] .= means(net.s)
        net.s_del[1] += 1
    end
    net.STDP,net.firings = time_reset(net.STDP,net.firings)
    if sec%100==0
        print("\rsec = $sec")
    end
end

# %% Plot means
using Plots
gr()
x1 = 0.001.*collect(1:length(net.shist[:,1]))
y1 = net.shist[:,1]
x2 = x1
y2 = net.shist[:,2]
fig = plot()
plot!(x1,y1,color="blue",label="from S1", legend = true)
plot!(x2,y2,color="green",label="mean", legend = true)
xlabel!("Time (sec)")
ylabel!("synaptic weigth (mV)")

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
int_rules_9_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 9.1*Integrand*simplification*rules =#
    #= Antiderivative((~u')*((~v)+(~w))^(~p'),~x) := Antiderivative(u*w^p,x) <-- FreeQ(p,x) &&  EqQ(v,0) =#
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(p'), ~x) => Antiderivative(~u * (~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~n, ~p], ~x) && EqQ(~a, 0)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p'), ~x) => Antiderivative(~u * (~a) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~n, ~p], ~x) && EqQ(~b, 0)
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(j')) ^ ~(p'), ~x) => Antiderivative(~u * (~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~j, 2 * ~n) && EqQ(~a, 0))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(j')) ^ ~(p'), ~x) => Antiderivative(~u * (~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~j, 2 * ~n) && EqQ(~b, 0))
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(j')) ^ ~(p'), ~x) => Antiderivative(~u * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && (EqQ(~j, 2 * ~n) && EqQ(~c, 0))
    @apply_utils Antiderivative(~(u') * (~(a') * ~v + ~(b') * ~v + ~(w')) ^ ~(p'), ~x) => Antiderivative(~u * ((~a + ~b) * ~v + ~w) ^ ~p, ~x) <-- FreeQ([~a, ~b], ~x) && Not(FreeQ(~v, ~x))
    @apply_utils Antiderivative(~(u') * (~Px) ^ ~p, ~x) => Antiderivative(~u * (~Px) ^ Simplify(~p), ~x) <-- PolyQ(~Px, ~x) && (Not(RationalQ(~p)) && (FreeQ(~p, ~x) && RationalQ(Simplify(~p))))
    @apply_utils Antiderivative(~a, ~x) => ~a * ~x <-- FreeQ(~a, ~x)
    @apply_utils Antiderivative(~a * (~b + ~(c') * ~x), ~x) => (~a * (~b + ~c * ~x) ^ 2) / (2 * ~c) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(-(~u), ~x) => Identity(-1) * Antiderivative(~u, ~x)
    @apply_utils Antiderivative(Complex(0, ~a) * ~u, ~x) => Complex(Identity(0), ~a) * Antiderivative(~u, ~x) <-- FreeQ(~a, ~x) && EqQ((~a) ^ 2, 1)
    @apply_utils Antiderivative(~a * ~u, ~x) => ~a * Antiderivative(~u, ~x) <-- FreeQ(~a, ~x) && Not(MatchQ(~u, ~b * ~v <-- FreeQ(b, ~x)))
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := ShowStep("", "Antiderivative(a*u + b*v + ...,x)", "a*Integrate(u,x) + b*Integrate(v,x) + ...", Hold( IntSum(u, x))) <-- SimplifyFlag && SumQ(u), Antiderivative((~u), ~x) := IntSum(u, x) <-- SumQ(u))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~u, ~x) => Antiderivative(ExpandIntegrand((~c * ~x) ^ ~m * ~u, ~x), ~x) <-- FreeQ([~c, ~m], ~x) && (SumQ(~u) && (Not(LinearQ(~u, ~x)) && Not(MatchQ(~u, ~a + ~(b') * ~v <-- FreeQ([a, b], ~x) && InverseFunctionQ(v)))))
    @apply_utils Antiderivative(~(u') * (~(a') * (~x) ^ ~n) ^ ~m, ~x) => (((~a) ^ IntPart(~m) * (~a * (~x) ^ ~n) ^ FracPart(~m)) / (~x) ^ (~n * FracPart(~m))) * Antiderivative(~u * (~x) ^ (~m * ~n), ~x) <-- FreeQ([~a, ~m, ~n], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative(~(u') * (~v) ^ ~(m') * (~b * ~v) ^ ~n, ~x) => (1 / (~b) ^ ~m) * Antiderivative(~u * (~b * ~v) ^ (~m + ~n), ~x) <-- FreeQ([~b, ~n], ~x) && IntegerQ(~m)
    @apply_utils Antiderivative(~(u') * (~(a') * ~v) ^ ~m * (~(b') * ~v) ^ ~n, ~x) => (((~a) ^ (~m + 1 / 2) * (~b) ^ (~n - 1 / 2) * sqrt(~b * ~v)) / sqrt(~a * ~v)) * Antiderivative(~u * (~v) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~m], ~x) && (Not(IntegerQ(~m)) && (IGtQ(~n + 1 / 2, 0) && IntegerQ(~m + ~n)))
    #= Antiderivative((~u')*((~a')*(~v))^(~m)*((~b')*(~v))^(~n),~x) :=  b^(n-1/2)*sqrt(b*v)/(a^(n-1/2)*sqrt(a*v))*Antiderivative(u*(a*v)^(m+n),x) <--  FreeQ([a,b,m],x) && Not(IntegerQ(m)) && IGtQ(n+1/2,0) &&  Not(IntegerQ(m+n)) =#
    @apply_utils Antiderivative(~(u') * (~(a') * ~v) ^ ~m * (~(b') * ~v) ^ ~n, ~x) => (((~a) ^ (~m - 1 / 2) * (~b) ^ (~n + 1 / 2) * sqrt(~a * ~v)) / sqrt(~b * ~v)) * Antiderivative(~u * (~v) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~m], ~x) && (Not(IntegerQ(~m)) && (ILtQ(~n - 1 / 2, 0) && IntegerQ(~m + ~n)))
    #= Antiderivative((~u')*((~a')*(~v))^(~m)*((~b')*(~v))^(~n),~x) :=  b^(n+1/2)*sqrt(a*v)/(a^(n+1/2)*sqrt(b*v))*Antiderivative(u*(a*v)^(m+n),x) <--  FreeQ([a,b,m],x) && Not(IntegerQ(m)) && ILtQ(n-1/2,0) &&  Not(IntegerQ(m+n)) =#
    @apply_utils Antiderivative(~(u') * (~(a') * ~v) ^ ~m * (~(b') * ~v) ^ ~n, ~x) => (((~a) ^ (~m + ~n) * (~b * ~v) ^ ~n) / (~a * ~v) ^ ~n) * Antiderivative(~u * (~v) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && (Not(IntegerQ(~n)) && IntegerQ(~m + ~n)))
    @apply_utils Antiderivative(~(u') * (~(a') * ~v) ^ ~m * (~(b') * ~v) ^ ~n, ~x) => (((~b) ^ IntPart(~n) * (~b * ~v) ^ FracPart(~n)) / ((~a) ^ IntPart(~n) * (~a * ~v) ^ FracPart(~n))) * Antiderivative(~u * (~a * ~v) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (Not(IntegerQ(~m)) && (Not(IntegerQ(~n)) && Not(IntegerQ(~m + ~n))))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * ~v) ^ ~(m') * (~c + ~(d') * ~v) ^ ~(n'), ~x) => (~b / ~d) ^ ~m * Antiderivative(~u * (~c + ~d * ~v) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (IntegerQ(~m) && (Not(IntegerQ(~n)) || SimplerQ(~c + ~d * ~x, ~a + ~b * ~x))))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * ~v) ^ ~m * (~c + ~(d') * ~v) ^ ~n, ~x) => (~b / ~d) ^ ~m * Antiderivative(~u * (~c + ~d * ~v) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (GtQ(~b / ~d, 0) && Not(IntegerQ(~m) || IntegerQ(~n))))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * ~v) ^ ~m * (~c + ~(d') * ~v) ^ ~n, ~x) => ((~a + ~b * ~v) ^ ~m / (~c + ~d * ~v) ^ ~m) * Antiderivative(~u * (~c + ~d * ~v) ^ (~m + ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && Not(IntegerQ(~m) || (IntegerQ(~n) || GtQ(~b / ~d, 0))))
    #= Antiderivative((~u')*((~a')*(~v))^(~m)*((~b')*(~v)+(~c')*(~v)^2),~x) :=  1/a*Antiderivative(u*(a*v)^(m+1)*(b+c*v),x) <-- FreeQ([a,b,c],x) && LeQ(m,-1) =#
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * ~v) ^ ~m * (~(A') + ~(B') * ~v + ~(C') * (~v) ^ 2), ~x) => (1 / (~b) ^ 2) * Antiderivative(~u * (~a + ~b * ~v) ^ (~m + 1) * Simp((~b * ~B - ~a * ~C) + ~b * ~C * ~v, ~x), ~x) <-- FreeQ([~a, ~b, ~A, ~B, ~C], ~x) && (EqQ((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C, 0) && LeQ(~m, -1))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(m') * (~c + ~(d') * (~x) ^ ~(q')) ^ ~(p'), ~x) => (~d / ~a) ^ ~p * Antiderivative((~u * (~a + ~b * (~x) ^ ~n) ^ (~m + ~p)) / (~x) ^ (~n * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~q, -(~n)) && (IntegerQ(~p) && (EqQ(~a * ~c - ~b * ~d, 0) && Not(IntegerQ(~m) && NegQ(~n)))))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(m') * (~c + ~(d') * (~x) ^ ~j) ^ ~(p'), ~x) => (-((~b) ^ 2) / ~d) ^ ~m * Antiderivative(~u * (~a - ~b * (~x) ^ ~n) ^ -(~m), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~j, 2 * ~n) && (EqQ(~p, -(~m)) && (EqQ((~b) ^ 2 * ~c + (~a) ^ 2 * ~d, 0) && (GtQ(~a, 0) && LtQ(~d, 0)))))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p'), ~x) => Antiderivative(~u * Cancel((~b / 2 + ~c * ~x) ^ (2 * ~p) / (~c) ^ ~p), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / (~c) ^ ~p) * Antiderivative(~u * (~b / 2 + ~c * (~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(~p)))
end

max_stars_repo_path: src/rules/9 Miscellaneous/9.1 Integrand simplification rules.jl


text:
"""
$(SIGNATURES)

Makes named example neural networks available as a [`NeuralNet`](@ref) object.

# Arguments
* `name::String`: Name of example neural network. Options:
    * `'MNIST.n1'`:
        * Architecture: Two fully connected layers with 40 and 20 units.
        * Training: Trained regularly with no attempt to increase robustness.
    * `'MNIST.WK17a_linf0.1_authors'`.
        * Architecture: Two convolutional layers (stride length 2) with 16 and
          32 filters respectively (size 4 × 4 in both layers), followed by a
          fully-connected layer with 100 units.
        * Training: Network trained to be robust to attacks with \$l_\\infty\$ norm
          at most 0.1 via method in [Provable defenses against adversarial examples
          via the convex outer adversarial polytope](https://arxiv.org/abs/1711.00851).
          Is MNIST network for which results are reported in that paper.
    * `'MNIST.RSL18a_linf0.1_authors'`.
        * Architecture: One fully connected layer with 500 units.
        * Training: Network trained to be robust to attacks with \$l_\\infty\$ norm
          at most 0.1 via method in [Certified Defenses against Adversarial Examples
          ](https://arxiv.org/abs/1801.09344).
          Is MNIST network for which results are reported in that paper.
"""
function get_example_network_params(name::String)::NeuralNet
    if name == "MNIST.n1"
        param_dict = prep_data_file(joinpath("weights", "mnist"), "n1.mat") |> matread
        fc1 = get_matrix_params(param_dict, "fc1", (784, 40))
        fc2 = get_matrix_params(param_dict, "fc2", (40, 20))
        logits = get_matrix_params(param_dict, "logits", (20, 10))

        nn = Sequential([Flatten(4), fc1, ReLU(interval_arithmetic), fc2, ReLU(), logits], name)
        return nn
    elseif name == "MNIST.WK17a_linf0.1_authors"
        param_dict =
            prep_data_file(
                joinpath("weights", "mnist", "WK17a", "linf0.1"),
                "master_seed_1_epochs_100.mat",
            ) |> matread
        conv1 = get_conv_params(param_dict, "conv1", (4, 4, 1, 16), expected_stride = 2)
        conv2 = get_conv_params(param_dict, "conv2", (4, 4, 16, 32), expected_stride = 2)
        fc1 = get_matrix_params(param_dict, "fc1", (1568, 100))
        logits = get_matrix_params(param_dict, "logits", (100, 10))

        nn = Sequential(
            [
                conv1,
                ReLU(interval_arithmetic),
                conv2,
                ReLU(),
                Flatten([1, 3, 2, 4]),
                fc1,
                ReLU(),
                logits,
            ],
            name,
        )
        return nn
    elseif name == "MNIST.RSL18a_linf0.1_authors"
        param_dict =
            prep_data_file(joinpath("weights", "mnist", "RSL18a", "linf0.1"), "two-layer.mat") |>
            matread
        fc1 = get_matrix_params(param_dict, "fc1", (784, 500))
        logits = get_matrix_params(param_dict, "logits", (500, 10))

        nn = Sequential([Flatten([1, 3, 2, 4]), fc1, ReLU(interval_arithmetic), logits], name)
        return nn
    else
        throw(DomainError("No example network named $name."))
    end
end

# TODO (vtjeng): Add mnist networks Ragunathan/Steinhardt/Liang.
# TODO (vtjeng): Make network naming case insensitive.


text:
export  excelInitPopPlot

function excelInitPopPlot( mpSim::MPsim, filename::AbstractString,
    sheetname::AbstractString="Initial population plots" )

    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )

    if !ispath( filename )
        @warn string( "File '", filename,
            "' not found. Can't generate plots/reports." )
        return
    end  # if !ispath( filename )

    makeCompositionPlot = nothing
    plotsToMake = Vector{Tuple}()

    XLSX.openxlsx( filename ) do xf
        if !XLSX.hassheet( xf, sheetname )
            @warn string( "Sheet '", sheetname,
                "' not found in file. Can't generate plots/reports." )
            return
        end  # if !XLSX.hassheet( xf, sheetname )

        xs = xf[sheetname]
        makeCompositionPlot = lowercase.( xs[CR( 5, 2, 7, 2 )][:] ) .== "yes"

        nToProcess = xs["B12"]

        for ii in (1:nToProcess) .+ 14
            # The age resolution and type of the plot.
            timeRes = xs[ii, 1]
            ageType = lowercase( string( xs[ii, 2] ) )

            if !( timeRes isa Real ) || ( timeRes <= 0 ) ||
                !haskey( ageTypes, ageType )
                continue
            end  # if !( timeRes isa Real ) || ...

            # Show plot in browser.
            showPlot = lowercase( string( xs[ii, 3] ) ) == "yes"

            # Save plot.
            savePlot = lowercase( string( xs[ii, 4] ) ) == "yes"

            # Generate report.
            makeReport = lowercase( string( xs[ii, 5] ) ) == "yes"

            push!( plotsToMake,
                (timeRes, ageTypes[ageType], showPlot, savePlot, makeReport) )
        end  # for ii in 1:nToProcess
    end  # XLSX.openxlsx( filename ) do xf

    plotdir = filename[1:(end-5)]

    if !ispath( plotdir )
        mkpath( plotdir )
    end  # if !ispath( plotdir )

    dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
    isFirst = true
    reportname = string( plotdir, "/initial population report ", dateStr )

    if any( makeCompositionPlot )
        tStart = now()
        report = initPopReport( mpSim )
        reportGenerationTime = ( now() - tStart ).value / 1000.0

        if !isempty( report )
            makePlot = any( makeCompositionPlot[1:2] )
            makeReport = makeCompositionPlot[3]

            if makePlot
                plotname = string( plotdir, "/initial population (", dateStr,
                    ")" )
                generateInitialPopPlot( report, makeCompositionPlot[1],
                    makeCompositionPlot[2], plotname )
            end  # makePlot

            if makeReport
                writeInitialPopData( report, mpSim, 12.0, reportGenerationTime,
                    reportname, isFirst )
                isFirst = false
            end  # if makeReport
        end  # if !isempty( report )
    end  # if any( makeCompositionPlot )

    for plotkey in plotsToMake
        if !any( [plotkey[3], plotkey[4], plotkey[5]] )
            continue
        end  # if !any( [plotkey[3], plotkey[4], plotkey[5]] )

        tStart = now()
        report = initPopAgeReport( mpSim, plotkey[1], plotkey[2] )
        reportGenerationTime = ( now() - tStart ).value / 1000.0

        if isempty( report )
            continue
        end  # if isempty( report )

        dateStr = string( " (", Dates.format( now(), "yyyymmdd HHMMSS" ), ")" )
        makePlot = plotkey[3] || plotkey[4]
        makeReport = plotkey[5]

        if makePlot
            plotname = string( plotdir, "/initial population ",
                ageTypesS[plotkey[2]], dateStr )
            generateInitialPopAgePlot( report, 12.0, plotkey[3], plotkey[4],
                plotname )
        end  # if makePlot

        if makeReport
            writeInitialPopAgeData( report, mpSim, plotkey[2], 12.0,
                reportGenerationTime, reportname, isFirst )
            isFirst = false
        end  # if makeReport
    end  # for plotkey in plotsToMake

end  # excelInitPopPlot( mpSim, filename, sheetname )


text:
#!/bin/env julia

import Pkg
Pkg.activate(;temp=true)
Pkg.add("PackageCompiler")
using PackageCompiler

Pkg.activate(".")

create_sysimage(["CairoMakie", "Dagger"]; sysimage_path="sysimage.so", precompile_execution_file="precompile.jl")

max_stars_repo_name: JuliaLabs/DaggerDiffusion

max_issues_repo_name: JuliaLabs/DaggerDiffusion

max_forks_repo_name: JuliaLabs/DaggerDiffusion


text:
"""
x = coord_desc(A, B; ...)

Solves non-negative least-squares problem by coordinate descent (CD).

Optional arguments:
    ε: tolerance for stopping (small number times sqrt[m*n] by default)

References:
    Reference has been lost...
"""

function coord_desc(A, B; ε = sqrt(size(A,2)*size(B,2))*1e-15, max_iter = 1000)
    rows,vars = size(A)
    ATA,ATB = transpose(A) * A, transpose(A) * B
    a = zeros(vars)
    μ = -ATB
    Haf = similar(ATB)
    @inbounds for iter in 1:max_iter
        Haf .= ATA * a .- ATB
        all(>=(-ε), Haf) && break
        for v in 1:vars
            initial = a[v]
            a[v,1]  = max(a[v] - μ[v] / ATA[v,v], 0.0)
            ∇       = a[v] - initial
            μ      .+= ∇ .* @views ATA[:,v]
        end
    end
    return a
end

# variables = 100
# X    = randn(200, variables)
# beta = rand(variables, 1)
# Y    = X * beta
# r = NNLS_CD(X, Y; ϵ=1e-9, max_iters=300)
# sum(abs,beta .- r)


text:
# See https://en.wikipedia.org/wiki/Prime_number

using Primes

expressions = [
    :(primes(0, 100))
    :(primes(1000000, 1000100))
    :(isprime(6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151))
    :(nextprime(6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151))
    :(totient(271))
    :(totient(2718281))
    :(totient(27182818284))
    :(totient(271828182845904))
    :(totient(2718281828459045235360287471352662497757247))
]

for expr in expressions
    println("[] :julia-answer \"\"\"", escape_string(string(expr)), " = ", eval(expr), "\"\"\".")
end

max_stars_repo_name: pmoura/eye

max_issues_repo_name: pmoura/eye

max_forks_repo_name: pmoura/eye


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 71c5500c-00cc-11eb-3162-a10318958fc2
using Flux

# ╔═╡ 703d8bca-00cd-11eb-1b57-af64a482bd8a
using LinearAlgebra

# ╔═╡ 78ac553e-00cd-11eb-26b7-f5bdbe229dd7
using Gadfly

# ╔═╡ 9ae94d2a-00cd-11eb-37b9-1dc46a16a49d
using Statistics

# ╔═╡ a23e157c-00cd-11eb-268c-d992c466fdd8
using Distributions

# ╔═╡ bd0aafe6-00cd-11eb-21ca-fd6c71a68f2f
using Random

# ╔═╡ 34ebec92-00d2-11eb-3f06-c30d260fddda
using DataFrames

# ╔═╡ 0eb7608c-00cc-11eb-0abe-bd2bca0d2b3e
md"# Sigmod regression

use to clssification

The loss function : cross entropy function
"

# ╔═╡ afdc663e-00cd-11eb-04b7-899986c22aee
Random.seed!(123)

# ╔═╡ c37e9a4a-00cd-11eb-32ce-b514071f2623
GroupOne = rand(MvNormal([10.0, 10.0], 10.0 * Matrix{Float64}(I, 2, 2)), 100)

# ╔═╡ cc6f4438-00cd-11eb-01c8-578b14b406b6
GroupTwo = rand(MvNormal([0.0, 0.0], 10 * Matrix{Float64}(I, 2, 2)), 100)

# ╔═╡ df38d270-00cd-11eb-2899-419fbf880361
X = Flux.normalise(hcat(GroupOne, GroupTwo), dims=2)

# ╔═╡ 2a87f08e-00d2-11eb-2889-1f2385032787
dt = DataFrame()

# ╔═╡ 44233b5c-00d2-11eb-3bc6-47248934072c
dt.x = X[1, :]

# ╔═╡ 50de7302-00d2-11eb-183a-752f712cbb38
dt.y = X[2, :]

# ╔═╡ 7177531a-00d2-11eb-1f6c-53d12ebffb5c
plot(dt, x=:x, y=:y, color=:label)

# ╔═╡ ffd9738e-00cd-11eb-0681-4d71a42001b1
begin
	y = []
	for i in 1:200
		if 1<=i<=100
			push!(y, "A")
		else
			push!(y, "B")
		end
	end
end

# ╔═╡ 57182a38-00d2-11eb-1aa3-f39349334c9c
dt.label = y

# ╔═╡ 258db6e4-00ce-11eb-08c7-0ff693c39c95
labels = Flux.onehotbatch(y, [:"A", :"B"])

# ╔═╡ 656fb6e0-00ce-11eb-014b-6d34bd579e77
train_data = Flux.Data.DataLoader((X, labels), batchsize=10, shuffle=true)

# ╔═╡ 808a8b6c-00ce-11eb-27f7-fb9c543b53ff
accuracy(x, y, model) = mean(Flux.onecold(model(x)) .== Flux.onecold(y))

# ╔═╡ a24cb676-00ce-11eb-36cd-0f7f0e955ec1
m = Chain(Dense(2,2))

# ╔═╡ a7bb0414-00ce-11eb-35b6-e596f46224d3
loss(x,y) = Flux.logitcrossentropy(m(x),y)

# ╔═╡ ac248b60-00ce-11eb-3589-e5522bb8f0ea
optimiser = Descent(0.3)

# ╔═╡ b0f07cf8-00ce-11eb-37a0-59880c3055b4
Flux.train!(loss, Flux.params(m), train_data, optimiser)

# ╔═╡ b43b19c2-00ce-11eb-032f-11a2e2c82b91
accuracy_score = accuracy(X, labels, m)

# ╔═╡ Cell order:
# ╟─0eb7608c-00cc-11eb-0abe-bd2bca0d2b3e
# ╠═71c5500c-00cc-11eb-3162-a10318958fc2
# ╠═703d8bca-00cd-11eb-1b57-af64a482bd8a
# ╠═78ac553e-00cd-11eb-26b7-f5bdbe229dd7
# ╠═9ae94d2a-00cd-11eb-37b9-1dc46a16a49d
# ╠═a23e157c-00cd-11eb-268c-d992c466fdd8
# ╠═bd0aafe6-00cd-11eb-21ca-fd6c71a68f2f
# ╠═34ebec92-00d2-11eb-3f06-c30d260fddda
# ╠═afdc663e-00cd-11eb-04b7-899986c22aee
# ╠═c37e9a4a-00cd-11eb-32ce-b514071f2623
# ╠═cc6f4438-00cd-11eb-01c8-578b14b406b6
# ╠═df38d270-00cd-11eb-2899-419fbf880361
# ╠═2a87f08e-00d2-11eb-2889-1f2385032787
# ╠═44233b5c-00d2-11eb-3bc6-47248934072c
# ╠═50de7302-00d2-11eb-183a-752f712cbb38
# ╠═57182a38-00d2-11eb-1aa3-f39349334c9c
# ╠═7177531a-00d2-11eb-1f6c-53d12ebffb5c
# ╠═ffd9738e-00cd-11eb-0681-4d71a42001b1
# ╠═258db6e4-00ce-11eb-08c7-0ff693c39c95
# ╠═656fb6e0-00ce-11eb-014b-6d34bd579e77
# ╠═808a8b6c-00ce-11eb-27f7-fb9c543b53ff
# ╠═a24cb676-00ce-11eb-36cd-0f7f0e955ec1
# ╠═a7bb0414-00ce-11eb-35b6-e596f46224d3
# ╠═ac248b60-00ce-11eb-3589-e5522bb8f0ea
# ╠═b0f07cf8-00ce-11eb-37a0-59880c3055b4
# ╠═b43b19c2-00ce-11eb-032f-11a2e2c82b91

max_stars_repo_name: ZhouZhuofei/juliaDL

max_issues_repo_name: ZhouZhuofei/juliaDL

max_forks_repo_name: ZhouZhuofei/juliaDL


text:
mutable struct AlexNet <: Photon.Layer
    features::Sequential
    output::Dense

    function AlexNet(; classes = 1000)
        features = Sequential(
            Conv2D(64, (11, 11), relu, strides = 4, padding = 2),
            MaxPool2D(3, strides = 2),
            Conv2D(192, (5, 5), relu, padding = 2),
            MaxPool2D(3, strides = 2),
            Conv2D(384, (3, 3), relu, padding = 1),
            Conv2D(256, (3, 3), relu, padding = 1),
            Conv2D(256, (3, 3), relu, padding = 1),
            MaxPool2D(3, strides = 2),
            Dense(4096, relu),
            Dropout(0.5),
            Dense(4096, relu),
            Dropout(0.5),
        )

        output = Dense(classes)
        new(features, output)
    end
end

function (m::AlexNet)(X)
    X = m.features(X)
    m.output(X)
end

@info "Loaded AlexNet model"


text:
"""
Compute preimages of monotonic sequences
"""

using IntervalArithmetic
using .Contractors

"""
Strict-precedes on intervals, extended to non-intervals
"""
lt(a::Interval, b::Interval) = strictprecedes(a, b)
lt(a, b::Interval) = a < b.lo
lt(a::Interval, b) = a.hi < b
lt(a, b) = a < b

gt(a, b) = lt(b, a)

"""
Computes the pair (skip, last), where `skip` is the number of elements
at the beginning of `seq.v` that do not intersect with `X`, and `last` is the
last element after `skip` to intersect with `X`.

If `X` falls entirely between `v[i]` and `v[i+1]`, then (i, i) is returned.
"""
function skipandlast(seq, X)
    cmp = ifelse(seq.increasing, lt, gt)
    return (searchsortedfirst(seq.v, X, lt=cmp) - 1, searchsortedlast(seq.v, X, lt=cmp))
end

"""
Type used to store monotonic sequences of preimages. `skip` is a certain number of initial elements that are skipped
with respect to an original reference array: for instance, preimages of a certain vector `y = [y(1), y(2), y(3), y(4), y(5)]`
may only contain f^{-1}(y(3)) and f^{-1}(y(4)), so we set skip=2 and construct a v of length 2.
`increasing` tells if the monotonic sequence `v` is increasing (true) or decreasing (false)
"""
struct PointSequence{T<:AbstractVector}
    v::T
    skip::Int
    increasing::Bool
end
PointSequence(v, skip=0, increasing=unique_increasing(v[begin], v[end])) = PointSequence{typeof(v)}(v, skip, increasing)

"""
Type used to represent a "branch" of a dynamic. The branch is represented by a monotonic map `f` with domain `X=(a,b)` with a≤b (where typically a,b are intervals).
`Y=(f(a),f(b))` and `increasing` may be provided (for instance if we know that `Y=(0,1)`), otherwise they are computed automatically.
"""
struct Branch{T,S}
    f::T
    X::Tuple{S, S}
    Y::Tuple{S, S}
    increasing::Bool
end
Branch(f, X, Y=(f(Interval(X[1])), f(Interval(X[2]))), increasing=unique_increasing(Y[1], Y[2])) = Branch{typeof(f), typeof(X[1])}(f, X, Y, increasing)

"""
Construct preimages of a monotonic array y under a monotonic function f in a domain X.

In general, there may be a certain number of points in y that have no preimage at the beginning and the end of the sequence, because
they fall out of the interval R = [f(X.lo), f(X.hi)]. In the worst case, no point has a preimage, because y[i] < R < y[i+1] for some
i (or vice versa with orientations).

We return a PointSequence where `v` is a vector of preimages, and `skip` is the number of elements of y that have no preimage
at the *beginning* of the interval. So for instance if f and y are increasing but y[i] < R < y[i+1], we return skip=i and v = [].
The number of intervals that have no image at the *end* of the interval is then `length(y) - length(v) - skip`

Rationale: this is more efficient than returning a vector with a lot of empty intervals at the beginning/end, and callers
need to know how many of the empty values come from the *beginning* of the array rather than its end.

Fills the array by using a bisection strategy to save computations: if y ∈ [a,b], then f⁻¹(y) ∈ [f⁻¹(a),f⁻¹(b)] (paying attention to orientation).
So we can fill v by filling in first entries `v[k+1]` with higher dyadic valuation of k.

Currently this works only for 1-based 1-dimensional arrays y.

TODO: Idea: replace this with a function that returns an *increasing* vector v that is guaranteed to start/end with v = (a,...,b), and a sequence of indices `locs` like k:m or m:-1:k
such that v[i], v[i+1] is inside `locs[i]` for all `i`. This constructs "slices" more explicitly, and composes much better. Unclear: how to keep 'full-branched' info in this framework?
"""
function preimages(seq, branch, ϵ = 0.0)
    v_increasing = !(seq.increasing ⊻ branch.increasing)
    (skip, last) = skipandlast(seq, hull(branch.Y...))

    n = last - skip

    v = fill((-∞..∞)::typeof(Interval(branch.X[1])), n)
    if n == 0
        return PointSequence(v, seq.skip+skip, v_increasing)
    end
    v[1] = preimage(seq.v[skip+1], branch.f, hull(branch.X...), ϵ)
    if n == 1
        return PointSequence(v, seq.skip+skip, v_increasing)
    end
    v[end] = preimage(seq.v[skip+n], branch.f, hull(branch.X...), ϵ)
    stride = prevpow(2, n-1)
    while stride >= 1
        # fill in v[i] using v[i-stride] and v[i+stride]
        for i = 1+stride:2*stride:n-1
            X = hull(v[i-stride], v[min(i+stride, n)]) #TODO: this hull() could be replaced with the proper [a.lo..b.hi], since we know orientations
            v[i] = preimage(seq.v[skip+i], branch.f, X, ϵ)
        end
        stride = stride ÷ 2
    end
    return PointSequence(v, seq.skip+skip, v_increasing)
end


text:
function Homogenization(U, Chis, N, E, nu)
    function homogenizeP()
        function MacroP(dir::Int64, U::Matrix{Float64}, elID::Int64, coords::Matrix{Float64}, E::Float64, nu::Float64)
            gps = gaussPoints1D(2)
            if dir == 1
                P = zeros(Float64, 4, 1)
                for gp in eachrow(gps)
                    B, detJ, PK1_1D, A_2D = constitutiveEqs(U, coords, 1/2, gp[2], E, nu)
                    P += gp[1]*PK1_1D*(1/2)
                end
                P11, P21 = P[1], P[4]
                return P11, P21
            elseif dir == 2
                P = zeros(Float64, 4, 1)
                for gp in eachrow(gps)
                    B, detJ, PK1_1D, A_2D = constitutiveEqs(U, coords, gp[2], 1/2, E, nu)
                    P += gp[1]*PK1_1D*(1/2)
                end
                P22, P12 = P[2], P[3]
                return P22, P12
            end
        end

        P11, P21 = 0, 0
        for i in N:N:N*N
            Ce = vec(DOFs[i, :, :]')
            P11_i, P21_i = MacroP(1, U[Ce, :], i, elCoords[i, :, :], E, nu)
            P11 += P11_i
            P21 += P21_i
        end

        P22, P12 = 0, 0
        for i in N*(N-1)+1:1:N*N
            Ce = vec(DOFs[i, :, :]')
            P22_i, P12_i = MacroP(2, U[Ce, :], i, elCoords[i, :, :], E, nu)
            P22 += P22_i
            P12 += P12_i
        end
        return L/(L*L)*hcat(P11, P22, P12, P21)
    end

    function homogenizeA()
        function MacroA(dir::Int64, U::Matrix{Float64}, elID::Int64, coords::Matrix{Float64}, E::Float64, nu::Float64)
            gps = gaussPoints1D(2)
            if dir == 1
                Ael = zeros(Float64, 4, 4)
                Ce = vec(DOFs[elID, :, :]')
                for i in 1:4
                    for j in 1:4
                        for gp in eachrow(gps)
                            B, detJ, PK1_1D, A_2D = constitutiveEqs(U, coords, 1/2, gp[2], E, nu)
                            Ql = B*Chis[i][Ce,:]
                            Qr = B*Chis[j][Ce,:]
                            Ael[i,j] += gp[1]*(Ql'*A_2D*Qr)[1]*(dx/2)
                        end
                    end
                end
                return Ael
            elseif dir == 2
                Ael = zeros(Float64, 4, 4)
                Ce = vec(DOFs[elID, :, :]')
                for i in 1:4
                    for j in 1:4
                        for gp in eachrow(gps)
                            B, detJ, PK1_1D, A_2D = constitutiveEqs(U, coords, gp[2], 1/2, E, nu)
                            Ql = B*Chis[i][Ce,:]
                            Qr = B*Chis[j][Ce,:]
                            Ael[i,j] += gp[1]*(Ql'*A_2D*Qr)[1]*(dx/2)
                        end
                    end
                end
                return Ael
            end
        end

        A14 = zeros(Float64, 4, 4)
        for i in N:N:N*N
            Ce = vec(DOFs[i, :, :]')
            A_i = MacroA(1, U[Ce, :], i, elCoords[i, :, :], E, nu)
            A14 += A_i
        end

        A23 = zeros(Float64, 4, 4)
        for i in N*(N-1)+1:1:N*N
            Ce = vec(DOFs[i, :, :]')
            A_i = MacroA(2, U[Ce, :], i, elCoords[i, :, :], E, nu)
            A23 += A_i
        end

        return L/(L*L)*vcat(A14[1,:]', A23[2,:]', A23[3,:]', A14[4,:]')
    end
    return homogenizeP(), homogenizeA()
end

max_stars_repo_name: likemaestro/Nonlinear-Homogenization

max_issues_repo_name: likemaestro/Nonlinear-Homogenization

max_forks_repo_name: likemaestro/Nonlinear-Homogenization


text:
import Convex
export variablesin

hash(x::Convex.AbstractExpr) = x.id_hash

# TODO memoize variables in expressions

# id2var maps variable ids to variables
function variablesin(x::Variable, id2var = Dict{UInt64, Variable}())
  myid = hash(x)
  id2var[myid] = x
  return Set{UInt64}([myid])
end
variablesin(x::Constant, id2var = Dict{UInt64, Variable}()) = Set{UInt64}()
function variablesin(x::Convex.AbstractExpr, id2var = Dict{UInt64, Variable}())
  myvars = Set{UInt64}()
  for child in x.children
    union!(myvars, variablesin(child, id2var))
  end
  myvars
end
function variablesin(x::Constraint, id2var = Dict{UInt64, Variable}())
  return union(variablesin(x.rhs, id2var), variablesin(x.lhs, id2var))
end
function variablesin(x, id2var = Dict{UInt64, Variable}())
  warn("why are you asking me this? i'm just a $(typeof(x))"); Set{UInt64}()
end

# Minkowski or outer product of sets
function outerproduct(s::Set{T}, t::Set{T}) where T
  results = Set{@compat Tuple{T,T}}()
  for si in s
    for ti in t
      union!(results, [(si, ti)])
    end
  end
  return results
end


text:
using Grid

include("derivative_numer.jl")

## Interpolation ##
# Low-level components
dims = (4, 6)
z = zeros(dims)
ic = InterpGridCoefs(z, InterpQuadratic)
x = [2,3.5]
set_position(ic, BCnearest, true, x)
@assert ic.coord1d == Vector{Int}[[1,2,3],[3,4,5]]
@assert ic.coef1d == Vector{Float64}[[1/8,3/4,1/8],[1/2,1/2,0]]
@assert ic.wrap == false
z[ic.offset+ic.offset_base+1] = ic.coef
@assert z == [0 0 1/16 1/16 0 0; 0 0 3/8 3/8 0 0; 0 0 1/16 1/16 0 0; 0 0 0 0 0 0]
# Test that the generic coef algorithm gives the same result as the
# dimension-specialized versions
for n_dims = 1:4
    dims = fill(3, n_dims)
    s = similar(dims)
    s[1] = 1
    for idim = 1:n_dims-1
        s[idim+1] = dims[idim]*s[idim]
    end
    ic = InterpGridCoefs(Float64, InterpQuadratic, dims, s)
    x = rand(n_dims) + 1.5
    x[x .== 2.5] = 2.49999
    set_position(ic, BCnearest, true, x)
    @assert ic.wrap == false
    @assert abs(ic.coef[iceil(length(ic.coef)/2)] - prod(3/4-(x-2).^2)) < eps()
    index = ic.offset+ic.offset_base+1
    indices_generic = similar(index)
    coefs_generic = similar(ic.coef)
    Grid.interp_index_coef_generic(indices_generic, coefs_generic, ic.coord1d, ic.coef1d, s)
    @assert index == indices_generic
    @assert ic.coef == coefs_generic
end

# On-grid values
A = randn(4,10)
const EPS = sqrt(eps())
for bc in (BCnil, BCnan, BCna, BCreflect, BCperiodic, BCnearest, BCfill)
    for it in (InterpNearest, InterpLinear, InterpQuadratic)
        ig = bc == BCfill ? InterpGrid(A, 0.0, it) : InterpGrid(A, bc, it)
        for i = 1:size(A,1)
            for j = 1:size(A,2)
                @assert abs(ig[i,j] - A[i,j]) < EPS
            end
        end
    end
end

# Quadratic interpolation
c = 2.3
a = 8.1
o = 1.6
qfunc = x -> a*(x-c).^2 + o
xg = Float64[1:5]
y = qfunc(xg)
yc = copy(y)
interp_invert!(yc, BCnan, InterpQuadratic)
ic = InterpGridCoefs(yc, InterpQuadratic)

x = 1.8
set_position(ic, BCnil, true, [x])
yi = interp(ic, yc)
@assert abs(yi - qfunc(x)) < 100*eps()
set_gradient_coordinate(ic, 1)
g = interp(ic, yc)
@assert abs(g - 2*a*(x-c)) < 100*eps()

# High-level interface
yi = InterpGrid(y, BCnil, InterpQuadratic)
@assert abs(yi[x] - qfunc(x)) < 100*eps()
v,g = valgrad(yi, x)
@assert abs(v - qfunc(x)) < 100*eps()
@assert abs(g - 2*a*(x-c)) < 100*eps()

# Test derivatives with other boundary conditions
Eps = sqrt(eps())
func = x -> yi[x]
gnum = derivative_numer(func, x)
@assert abs(g-gnum) < Eps*(abs(g)+abs(gnum))
for BC in (BCnan, BCna, BCnearest, BCperiodic, BCreflect, 0)
    yi = InterpGrid(y, BC, InterpQuadratic)
    func = x -> yi[x]
    v,g = valgrad(yi, x)
    gnum = derivative_numer(func, x)
    @assert abs(g-gnum) < Eps*(abs(g)+abs(gnum))
end

# Test derivatives of 2d interpolation
y = rand(7,8)
yi = InterpGrid(y, BCreflect, InterpQuadratic)
x = [2.2, 3.1]
v,g = valgrad(yi, x...)
func = (x1,x2) -> yi[x1,x2]
gnum = derivative_numer(func, tuple(x...), 1)
@assert abs(g[1]-gnum) < Eps*(abs(g[1])+abs(gnum))
gnum = derivative_numer(func, tuple(x...), 2)
@assert abs(g[2]-gnum) < Eps*(abs(g[2])+abs(gnum))

#### Interpolation on irregularly-spaced grids ####
x = [100.0,110.0,150.0]
y = rand(3)
iu = InterpIrregular(x, y, -200, InterpNearest)
@assert iu[99] == -200
@assert iu[101] == y[1]
@assert iu[106] == y[2]
@assert iu[149] == y[3]
@assert iu[150.1] == -200
iu = InterpIrregular(x, y, BCna, InterpLinear)
@assert isnan(iu[99])
@assert abs(iu[101] - (0.9*y[1] + 0.1*y[2])) < Eps
@assert abs(iu[106] - (0.4*y[1] + 0.6*y[2])) < Eps
@assert abs(iu[149] - (y[2]/40 + (39/40)*y[3])) < Eps
@assert isnan(iu[150.1])

#### Restrict/prolong ####
# 1d
a = zeros(5)
a[3] = 1
p = prolong(a, [9])
@assert all(p .== [0,0,0,0.5,1,0.5,0,0,0])
r = restrict([0,0,0,0,1.0,0,0,0,0], 1)
@assert all(r .== [0,0,1,0,0])
r = restrict([0,1.0,0,0,0,0,0,0,0], 1)
@assert all(r .== [0.5,0.5,0,0,0])
r = restrict(p, 1)
@assert all(r .== [0,0.25,1.5,0.25,0])
p = prolong(a, [8])
@assert all(p .== [0,0,0.25,0.75,0.75,0.25,0,0])
a[3] = 0
a[1] = 1
p = prolong(a, [9])
@assert all(p .== [1,0.5,0,0,0,0,0,0,0])
p = prolong(a, [8])
@assert all(p .== [0.75,0.25,0,0,0,0,0,0])
a = [0,0,1.0,0]
p = prolong(a, [7])
@assert all(p .== [0,0,0,0.5,1,0.5,0])
p = prolong(a, [6])
@assert all(p .== [0,0,0.25,0.75,0.75,0.25])
r = restrict([0,1.0,0,0,0,0,0,0], 1)
@assert all(r .== [0.25,0.75,0,0,0])
r = restrict([0,0,1.0,0,0,0,0,0], 1)
@assert all(r .== [0,0.75,0.25,0,0])
# Check with non-BLAS types
a = Vector{Float64}[[0.0],[0.0],[1.0],[0.0],[0.0]]
p = prolong(a, [9])
@assert all(p .== Vector{Float64}[[0.0],[0.0],[0.0],[0.5],[1.0],[0.5],[0.0],[0.0],[0.0]])
p = prolong(a, [8])
@assert all(p .== Vector{Float64}[[0.0],[0.0],[0.25],[0.75],[0.75],[0.25],[0.0],[0.0]])
a = Vector{Float64}[[0],[0],[1.0],[0]]
p = prolong(a, [7])
@assert all(p .== Vector{Float64}[[0],[0],[0],[0.5],[1],[0.5],[0]])
p = prolong(a, [6])
@assert all(p .== Vector{Float64}[[0],[0],[0.25],[0.75],[0.75],[0.25]])
# 2d
a = zeros(5,5)
a[3,3] = 1
p = prolong(a, 1, 9)
@assert all(p .== [zeros(9,2) [0,0,0,0.5,1,0.5,0,0,0] zeros(9,2)])
p = prolong(a, 2, 9)
@assert all(p .== [zeros(9,2) [0,0,0,0.5,1,0.5,0,0,0] zeros(9,2)]')
p = prolong(a, 1, 8)
@assert all(p .== [zeros(8,2) [0,0,0.25,0.75,0.75,0.25,0,0] zeros(8,2)])
p = prolong(a, 2, 8)
@assert all(p .== [zeros(8,2) [0,0,0.25,0.75,0.75,0.25,0,0] zeros(8,2)]')
r = restrict(a, [true,false])
@assert all(r .== [zeros(3,2) [0,1.0,0] zeros(3,2)])
r = restrict(a, [false,true])
@assert all(r .== [zeros(3,2) [0,1.0,0] zeros(3,2)]')
a = zeros(6,6)
a[3,3] = 1
r = restrict(a, [true,false])
@assert all(r .== [zeros(4,2) [0,0.75,0.25,0] zeros(4,3)])
r = restrict(a, [false,true])
@assert all(r .== [zeros(4,2) [0,0.75,0.25,0] zeros(4,3)]')
# Non-BLAS
a = fill([0.0], 5, 5)
a[3,3] = [1.0]
p = prolong(a, 1, 9)
ptest = fill([0.0], 9, 5)
ptest[6,3] = ptest[4,3] = [0.5]
ptest[5,3] = [1.0]
@assert all(p .== ptest)
p = prolong(a, 2, 9)
@assert all(p .== ptest')
r = restrict(a, [true,false])
rtest = fill([0.0], 3, 5)
rtest[2,3] = [1.0]
@assert all(r .== rtest)
r = restrict(a, [false,true])
@assert all(r .== rtest')
a = fill([0.0], 5, 5)
a[2,3] = [1.0]
rtest = fill([0.0], 3, 5)
rtest[1,3] = rtest[2,3] = [0.5]
r = restrict(a, [true,false])
@assert all(r .== rtest)
a = fill([0.0], 4, 4)
a[2,3] = [1.0]
r = restrict(a, [true,false])
rtest = fill([0.0], 3, 4)
rtest[1,3] = [0.25]
rtest[2,3] = [0.75]
@assert all(r .== rtest)
r = restrict(a, [false,true])
rtest = fill([0.0], 4, 3)
rtest[2,2] = [0.75]
rtest[2,3] = [0.25]
@assert all(r .== rtest)


text:
@testset "73.set-matrix-zeroes.jl" begin
    @test set_zeroes([[1, 1, 1], [1, 0, 1], [1, 1, 1]]) == [[1, 0, 1], [0, 0, 0], [1, 0, 1]]
    @test set_zeroes([[0, 1, 2, 0], [3, 4, 5, 2], [1, 3, 1, 5]]) ==
          [[0, 0, 0, 0], [0, 4, 5, 0], [0, 3, 1, 0]]
end


text:
module DistDict

using ComputeFramework
import ComputeFramework: split_range, PartitionScheme, alignfirst, DomainBranch, project, isempty, intersect, domain, cat_data

immutable DictDomain <: Domain
    hashrange::UnitRange
end

Base.isempty(d::DictDomain) = isempty(d.hashrange)
Base.intersect(d::DictDomain, d2::DictDomain) =
    DictDomain(intersect(d.hashrange, d2.hashrange))
project(d::DictDomain, d2::DictDomain) = d
Base.getindex(bigd::DictDomain, d::DictDomain) = d
alignfirst(a::DictDomain) = a

function Base.getindex(dict::Dict, d::DictDomain)
    res = Dict()
    for k in keys(dict)
        if hash(k) in d.hashrange
            res[k] = dict[k]
        end
    end
    res
end

immutable DictPartition <: PartitionScheme
    bins::Int
end

function ComputeFramework.partition(p::DictPartition, d::DictDomain)
    rs = collect(1:UInt(typemax(UInt)/p.bins+1):typemax(UInt))
    ranges = map(UnitRange, rs[1:end-1], rs[2:end])
    ComputeFramework.DomainBranch(d, map(DictDomain, ranges))
end

domain(x::Dict) = DictDomain(zero(UInt):(typemax(UInt)-1))

cat_data(x::DictPartition, dmn::DomainBranch, parts::AbstractArray) = merge(parts...)

Base.map(f, d::Dict) = Dict(Pair[k=>f(v) for (k,v) in d])

end # module


text:
export ¬, ∧, ∨, →, ⟷
export AbstractExpression, LogicalSymbol, LogicalExpression, istree, isnode, name, metadata, variables, operation, operations, arguments, parents, left, right, isassociative, iscommutative
export isunary, isbinary, argument_count
export @symbols, @unique_symbols

"""
    AbstractExpression

The abstract type which all expression trees are built from.
"""
abstract type AbstractExpression end

"""
    LogicalSymbol(name::Symbol, metadata::Any=nothing)

Represent a logical symbolically with a provided name. By default the attached metadata is set to `nothing`.
"""
struct LogicalSymbol <: AbstractExpression
    name::Symbol
    metadata::Any
end
LogicalSymbol(name::Symbol) = LogicalSymbol(name, nothing)
"""
    name(sym::LogicalSymbol)

The name of the symbol provided at instantiation. Equivalent to `symbol(sym)`.
"""
name(sym::LogicalSymbol) = sym.name
"""
    metadata(sym::LogicalSymbol)

The metadata of the symbol provided at instantiation, if any. Returns `nothing` if none was provided.
"""
metadata(sym::LogicalSymbol) = sym.metadata

Base.show(io::IO, sym::LogicalSymbol) = print(io, string(sym.name))
Base.hash(sym::LogicalSymbol, h::UInt) = hash(sym.name, hash(sym.metadata, h))
Base.:(==)(sym1::LogicalSymbol, sym2::LogicalSymbol) = sym1.name == sym2.name && isequal(metadata(sym1), metadata(sym2))
Base.isless(sym1::LogicalSymbol, sym2::LogicalSymbol) = Base.isless(sym1.name, sym2.name)

# fake copy methods since symbols are immutable
Base.copy(sym::LogicalSymbol) = sym
Base.deepcopy(sym::LogicalSymbol) = LogicalSymbol(name(sym), deepcopy(metadata(sym)))

# convenience macros
"""
Define any number of [`LogicalSymbols`](@ref Deductive.LogicalSymbol) with the names provided.

# Examples
```julia-repl
julia> @symbols a  # defines symbol `a` in the current scope

julia> @symbols b c d  # defines symbols `a`, `b`, and `c`

julia> @symbols α β γ  # defines symbols `α`, `β`, and `γ`
```
"""
macro symbols(syms...)
    definitions = [:(
        $(esc(sym)) = LogicalSymbol($(esc(Symbol))($(string(sym))))
    ) for sym ∈ syms]
    quote
        $(definitions...)
        nothing
    end
end

_unique_index = 1
macro unique_symbols(syms...)
    global _unique_index

    unique_names = [Symbol("u" * subscript_number(_unique_index + i)) for i ∈ 0:(length(syms) - 1)]
    definitions = [:(
        $(esc(sym)) = LogicalSymbol($(esc(Symbol))($(string(unique_name))))
    ) for (sym, unique_name) ∈ zip(syms, unique_names)]

    _unique_index += length(syms)

    quote
        $(definitions...)
        nothing
    end
end

"""
    LogicalOperation(bool_fn::Function, name::Symbol, argument_count::Int, associative::Bool, commutative::Bool)

Defines a logical operator which can be used to form expressions with. Usually you don't need to define your own
operators, just use the ones built-in to the package if possible.
"""
struct LogicalOperation
    bool_fn::Function
    name::Symbol
    argument_count::Int
    associative::Bool
    commutative::Bool
end
LogicalOperation(bool_fn::Function, name::Symbol, argument_count::Int) = LogicalOperation(bool_fn, name, argument_count, false, false)
name(op::LogicalOperation) = op.name
"""
    argument_count(op::LogicalOperation)

The number of arguments which the given operation expects to receive.
"""
argument_count(op::LogicalOperation) = op.argument_count
"""
    isunary(op::LogicalOperation)

True if the given operation only receives one argument, false otherwise.
"""
isunary(op::LogicalOperation) = argument_count(op) == 1
"""
    isbinary(op::LogicalOperation)

True if the given operation receives two arguments, false otherwise.
"""
isbinary(op::LogicalOperation) = argument_count(op) == 2
"""
    isassociative(op::LogicalOperation)

True if the given operation is associative, false otherwise.
"""
isassociative(op::LogicalOperation) = op.associative
"""
    iscommutative(op::LogicalOperation)

True if the given operation is commutative, false otherwise.
"""
iscommutative(op::LogicalOperation) = op.commutative
Base.show(io::IO, op::LogicalOperation) = print(io, string(op.name))
Base.hash(op::LogicalOperation, h::UInt) = hash(op.name, hash(argument_count(op), hash(op.associative, hash(op.commutative, h))))
function Base.:(==)(op1::LogicalOperation, op2::LogicalOperation)
    op1.name == op2.name && argument_count(op1) == argument_count(op2) && isassociative(op1) == isassociative(op2) && iscommutative(op1) == iscommutative(op2)
end

"""
    (op::LogicalOperation)(args::AbstractExpression...)

This function is the foundation of expression building, as it allows instances of `LogicalOperations` to be treated as real
Julia functions. When called this function returns a `LogicalExpression` representing the application of the given operator
on the provided arguments.
"""
function (op::LogicalOperation)(args::AbstractExpression...)
    if length(args) != argument_count(op)
        throw(ErrorException("Invalid argument count $(length(args)). Expected $(argument_count(op)) arguments."))
    end
    return LogicalExpression(AbstractExpression[args...], op)
end
(op::LogicalOperation)(args::Bool...) = op.bool_fn(args...)
(op::LogicalOperation)(args::BitVector) = op.bool_fn(args...)

recursivevariables(args::Vector{AbstractExpression}) = reduce(∪, variables.(args))
recursiveoperations(args::Vector{AbstractExpression}, rootop::LogicalOperation) = reduce(∪, operations.(args)) ∪ Set([rootop])
"""
    LogicalExpression(arguments::Vector{AbstractExpression}, operation::LogicalOperation)

Constructs an expression with given arguments and logical operation. Please refrain from using this syntax, instead using
the "operators as functions" syntax, where a `LogicalOperation` instance can be called to produce a `LogicalExpression`.
"""
mutable struct LogicalExpression <: AbstractExpression
    arguments::Vector{AbstractExpression}
    operation::LogicalOperation

    parents::Set{LogicalExpression}

    cached_variables::Set{LogicalSymbol}
    cached_operations::Set{LogicalOperation}
    cached_variables_valid::Bool
    cached_operations_valid::Bool

    function LogicalExpression(arguments::Vector{AbstractExpression}, operation::LogicalOperation)
        expr = new(arguments, operation, Set{LogicalExpression}(), recursivevariables(arguments), recursiveoperations(arguments, operation), true, true)
        add_to_parents!(expr)
        expr
    end
end
"""
    operation(expr::LogicalExpression)

The operation which is performed on the arguments of the given `LogicalExpression`.
"""
operation(expr::LogicalExpression) = getfield(expr, :operation)
"""
    parents(expr::LogicalExpression)

The parent expressions of a given subexpression. If an expression is never used within another, it will have no parents.
In most expressions there will only be one parent, but it is possible for an expression assigned to a variable to have
multiple parents by using it as a "named subexpression".
"""
parents(expr::LogicalExpression) = getfield(expr, :parents)

"""
    add_to_parents!(expr::LogicalExpression)

Adds an expression to the `parents` set of all its arguments. For internal use only.
"""
add_to_parents!(expr::LogicalExpression) = add_to_parents!(expr, arguments(expr))
function add_to_parents!(expr::LogicalExpression, relevant_args::Vector{AbstractExpression})
    for arg ∈ relevant_args
        if arg isa LogicalExpression
            push!(parents(arg), expr)
        end
    end
end

"""
    remove_from_parents!(expr::LogicalExpression)

Removes an expression from the `parents` set of all its arguments. For internal use only.
"""
remove_from_parents!(expr::LogicalExpression) = remove_from_parents!(expr, arguments(expr))
function remove_from_parents!(expr::LogicalExpression, relevant_args::Vector{AbstractExpression})
    for arg ∈ relevant_args
        if arg isa LogicalExpression
            delete!(parents(arg), expr)
        end
    end
end
metadata(::LogicalExpression) = nothing
"""
    left(expr::LogicalExpression)

If the expression is binary, this method returns the left-hand operand.
"""
left(expr::LogicalExpression) = isbinary(operation(expr)) ? arguments(expr)[1] : throw(ErrorException("Operation $(operation(expr)) is not binary"))
"""
    right(expr::LogicalExpression)

If the expression is binary, this method returns the right-hand operand.
"""
right(expr::LogicalExpression) = isbinary(operation(expr)) ? arguments(expr)[2] : throw(ErrorException("Operation $(operation(expr)) is not binary"))
"""
    isassociative(expr::LogicalExpression)

Checks if the entire expression is associative based on the associative property of its constituent operations.
"""
isassociative(expr::LogicalExpression) = length(operations(expr)) == 1 && isassociative(operation(expr))
"""
    iscommutative(expr::LogicalExpression)

Checks if the entire expression is commutative based on the commutative property of its constituent operations.
"""
iscommutative(expr::LogicalExpression) = length(operations(expr)) == 1 && iscommutative(operation(expr))
Base.hash(expr::LogicalExpression, h::UInt) = hash(arguments(expr), hash(operation(expr), h))
Base.:(==)(expr1::LogicalExpression, expr2::LogicalExpression) = operation(expr1) == operation(expr2) && all(arguments(expr1) .== arguments(expr2))

# copy methods
Base.copy(expr::LogicalExpression) = LogicalExpression(Vector{AbstractExpression}(arguments(expr)), operation(expr))
Base.deepcopy(expr::LogicalExpression) = LogicalExpression(Vector{AbstractExpression}(deepcopy.(arguments(expr))), operation(expr))

# mutability methods
function Base.setproperty!(expr::LogicalExpression, name::Symbol, x)
    if name == :arguments
        remove_from_parents!(expr)
        setfield!(expr, name, x)
        add_to_parents!(expr)
        expr.cached_variables_valid = false
    elseif name == :operation
        @assert argument_count(operation(expr)) == argument_count(x) "cannot assign operation with $(argument_count(x)) arguments to an expression with $(argument_count(operation(expr)))"
        setfield!(expr, name, x)
        expr.cached_operations_valid = false
    elseif name == :cached_variables_valid || name == :cached_operations_valid
        setfield!(expr, name, x)
        if !x  # cache invalidation, so propagate the invalidation up the expression tree
            for parent ∈ parents(expr)
                setproperty!(parent, name, x)
            end
        end
    elseif name == :cached_variables || name == :cached_operations
        throw(ErrorException("type LogicalExpression has protected field `$(name)`"))
    else
        throw(ErrorException("type LogicalExpression has no field `$(name)`"))
    end
end
function Base.getproperty(expr::LogicalExpression, name::Symbol)
    if name == :arguments
        return FakeVector(expr, name, getfield(expr, name))
    end

    throw(ErrorException("use method `$(name)(my_expr)` instead"))
end
# see utils.jl#FakeVector
function setvectorindex!(expr::LogicalExpression, name::Symbol, x, index::Int)
    if name == :arguments
        remove_from_parents!(expr, arguments(expr)[index:index])
        getfield(expr, :arguments)[index] = x
        add_to_parents!(expr, arguments(expr)[index:index])
        expr.cached_variables_valid = false
    else
        @warn "unexpected `setvectorindex!` call for field `$(name)`"
    end
end
# this method is probably not optimized but it doesn't matter since expressions usually have two or fewer arguments, and this method
# will hardly ever be used regardless
function setvectorindex!(expr::LogicalExpression, name::Symbol, x, index::Union{UnitRange{Int}, StepRange{Int, Int}})
    li = 1  # linear index, to account for StepRange
    for i ∈ index
        setvectorindex!(expr, name, x[li], i)
        li += 1
    end
end

function Base.show(io::IO, expr::LogicalExpression)
    showparens(expr) = (expr isa LogicalExpression) && !isunary(operation(expr))

    if isunary(operation(expr))
        arg = first(arguments(expr))
        show(io, operation(expr))
        if showparens(arg)
            print(io, "(")
        end
        print(io, arg)
        if showparens(arg)
            print(io, ")")
        end
    elseif isbinary(operation(expr))
        args = arguments(expr)

        if showparens(args[1])
            print(io, "(")
        end
        show(io, arguments(expr)[1])
        if showparens(args[1])
            print(io, ")")
        end

        print(io, " ")
        show(io, operation(expr))
        print(io, " ")

        if showparens(args[2])
            print(io, "(")
        end
        show(io, arguments(expr)[2])
        if showparens(args[2])
            print(io, ")")
        end
    end
end

# Standard AbstractExpression interface methods
"""
    istree(expr::T) where {T <: AbstractExpression}

Returns true when `expr` is a `LogicalExpression` and false otherwise.
"""
istree(::LogicalSymbol) = false
istree(::LogicalExpression) = true

"""
    isnode(expr::T) where {T <: AbstractExpression}

Returns true when `expr` is a `LogicalSymbol` and false otherwise.
"""
isnode(::LogicalSymbol) = true
isnode(::LogicalExpression) = false

"""
    variables(expr::T) where {T <: AbstractExpression}

Returns variables present in an entire expression tree. When the first argument is a `LogicalSymbol`, this is singleton
set containing only the provided `LogicalSymbol`. When the argument is a `LogicalExpression`, this is a set containing the union
of all the variables present in each argument.
"""
variables(sym::LogicalSymbol) = Set{LogicalSymbol}(LogicalSymbol[sym])
function variables(expr::LogicalExpression)
    if !getfield(expr, :cached_variables_valid)
        setfield!(expr, :cached_variables, recursivevariables(arguments(expr)))
        expr.cached_variables_valid = true
    end
    getfield(expr, :cached_variables)
end

"""
    operations(expr::T) where {T <: AbstractExpression}

Returns the operations present in an entire expression tree. When given a `LogicalSymbol`, this is an empty set, whereas
when given a `LogicalExpression`, this is a set containing the expression's own `LogicalOperation` and the operations set
of each of its arguments.
"""
operations(::LogicalSymbol) = Set{LogicalOperation}([])
function operations(expr::LogicalExpression)
    if !getfield(expr, :cached_operations_valid)
        setfield!(expr, :cached_operations, recursiveoperations(arguments(expr), operation(expr)))
        expr.cached_operations_valid = true
    end
    getfield(expr, :cached_operations)
end

"""
    arguments(expr::T) where {T <: AbstractExpression}

Returns the arguments which an `AbstractExpression` contains. `LogicalSymbols` contain no arguments and `LogicalExpressions`
can contain any number of arguments.
"""
arguments(::LogicalSymbol) = AbstractExpression[]
arguments(expr::LogicalExpression) = getfield(expr, :arguments)

# ASSOCIATION (perhaps move this elsewhere eventually?)
"""
    associative_ordering(expr::LogicalExpression)

Descends the expression tree with a left-side-first depth first search. Each symbol encountered is added to a list in
the order it appears in this search and is returned by this function.
"""
associative_ordering(expr::LogicalExpression) = reduce(vcat, associative_ordering.(arguments(expr)))
associative_ordering(sym::LogicalSymbol) = [sym]

"""
    isequal_associative(expr1::LogicalExpression, expr2::LogicalExpression)

Checks whether an expression is equal with another ignoring the associative ordering of each expression.
"""
function isequal_associative(expr1::LogicalExpression, expr2::LogicalExpression)
    length(operations(expr1)) == 1 && isassociative(first(operations(expr1))) && operations(expr1) == operations(expr2) && associative_ordering(expr1) == associative_ordering(expr2)
end
isequal_associative(sym1::LogicalSymbol, sym2::LogicalSymbol) = isequal(sym1, sym2)
isequal_associative(::AbstractExpression, ::AbstractExpression) = false

_associative_tree_count_cache = Dict()
"""
    associative_tree_count(nodes::Int)

Function which calculates how many ways to arrange parenthesis in an expression there are with a given node count.
"""
function associative_tree_count(nodes::Int)
    if nodes == 0
        return 1
    end
    if nodes == 1
        return 1
    end

    if haskey(_associative_tree_count_cache, nodes)
        return _associative_tree_count_cache[nodes]
    end

    agg = 0
    for i ∈ 0:(nodes - 1)
        agg += associative_tree_count(nodes - i - 1) * associative_tree_count(i)
    end

    _associative_tree_count_cache[nodes] = agg

    agg
end

# OPERATORS
"""
    ¬(x)

Logical negation operator, typed with `\\neg`. Boolean equivalent is the `!` operator.
"""
const ¬ = LogicalOperation(x -> !x, :¬, 1, false, false)

# binary operators
"""
    x ∧ y

Logical conjunction operator, typed with `\\wedge`. Boolean equivalent is the `&&` operator.
"""
const ∧ = LogicalOperation((x, y) -> x && y, :∧, 2, true, true)
"""
    x ∨ y

Logical disjunction operator, typed with `\\vee`. Boolean equivalent is the `||` operator.
"""
const ∨ = LogicalOperation((x, y) -> x || y, :∨, 2, true, true)
"""
    x → y

Logical implication operator, typed with `\\rightarrow`.
"""
const → = LogicalOperation((x, y) -> (¬x ∨ y), :→, 2, false, false)
"""
    x ⟷ y

Logical equivalence operator, typed with `\\longleftrightarrow`.
"""
const ⟷ = LogicalOperation((x, y) -> (x ∧ y) ∨ (¬x ∧ ¬y), :⟷, 2, true, true)

"""
    random_expression(subchance=0.5, symbol_alphabet::Vector{LogicalSymbol}=LogicalSymbol.([:a, :b, :c]), operator_alphabet::Vector{LogicalOperation}=[∧, ∨, →, ⟷, ¬])

Creates a random expression from a symbol and operator alphabet. The `subchance` parameter describes how likely it is that
an argument will be a subexpression as opposed to an expression.
"""
function random_expression(subchance=0.5, symbol_alphabet::Vector{LogicalSymbol}=LogicalSymbol.([:a, :b, :c]), operator_alphabet::Vector{LogicalOperation}=[∧, ∨, →, ⟷, ¬])
    if rand() > subchance
        return rand(symbol_alphabet)
    else
        rand_op = rand(operator_alphabet)
        return LogicalExpression(AbstractExpression[random_expression(subchance, symbol_alphabet, operator_alphabet) for _ ∈ 1:argument_count(rand_op)], rand_op)
    end
end