text:
using JLD2
using FileIO

using Glob
using LinearAlgebra

include("simulations-bcel-DM-diss.jl")

function give_alpha_set(n,αmax)
    return rand(n).*αmax
end

function takememory_alpha_set(β,v,α)

    prob=probmemory_alpha_set(β,v,α)
    wahrscheinlikeit = rand()
    result=false
    if wahrscheinlikeit < prob
        result=true
    end
    result
end

function probmemory_alpha_set(β,v,α)
    prob = 1- exp(-β *mean(energymemory.(v,α)))
    prob
end

function onestep_dritf_real_alpha_set(β,δ,oldv,αs,n)
    pos=0
    if (takememory_alpha_set(β,oldv,αs))
        probs=energymemory.(oldv,αs)
        d = Categorical(probs./sum(probs))
        pos=rand(d)
        α=αs[pos]
        free=freeenergy(β,oldv,α)
        diss=dissipation(β,oldv,α)
        unet=free -diss
        oldv=newv_drift(δ,n)
        n+=1

    else
        unet = 0
        free=0
        diss=0
        oldv=newv_drift(δ,1)

        n=2
    end

    return (unet,oldv,free,diss,n,pos)
end

function livesycle_drift_real_alpha_set(β,δ,αs,meanencounters,reps)
    meanunet=0 #initializing for all the rounds for fist cost
    meanfree=0
    meandiss=0
    naiveuses = 0
    for round in 1:reps
        encounters = meanencounters *1
        localnetu= 0
        netu = 0
        localfree= 0
        free = 0
        localdiss= 0
        diss = 0

        v=newv_drift(δ,1) #start with memory produced at t= 0
        naiveuses += 1/reps/encounters
        n=2
        for encounter in 2:encounters  #start at 1 because initialization is allways at 0 in the first round without any prior memory
             naiveuse =0    # Used naive for the first virus
            (netu,v,free,diss,n,dummi)=onestep_dritf_real_alpha_set(β,δ,v,αs,n)
            localnetu+=netu
            localfree+=free
            localdiss+=diss
        end

        meanunet+= localnetu/reps/encounters
        meanfree+= localfree/reps/encounters
        meandiss+= localdiss/reps/encounters
    end

    return (meanunet*meanencounters,meanfree*meanencounters,meandiss*meanencounters)
end

function livesycle_drift_real_alpha_set_deltauniform(β,δmax,αs,meanencounters,reps)
    meanunet=0 #initializing for all the rounds for fist cost
    meanfree=0
    meandiss=0
    for round in 1:reps
        δ=rand()*δmax
    localnetu,localfree,localdiss=livesycle_drift_real_alpha_set(β,δ,αs,meanencounters,1)

        meanunet+= localnetu/reps
        meanfree+= localfree/reps
        meandiss+= localdiss/reps
    end

    return (meanunet,meanfree,meandiss)
end

function alpharange(α,amin,amax) # puts alpha in to the range between alpha_min and alpha_max
    return max(amin,min(α,amax))
end
function optimization_set_alphas(β,δmax,αs,meanencounters,reps,steps,gradients,λ1,λ2)  ## does "steps" optimization steps starting from αs; for  αs==zeros(l) gives rand alpha of langth(l)
    amax=4
    amin=0.001
    l=length(αs)
    if (αs[1]==0)
        αs=rand(l)*amax
    end
    αs=alpharange.(αs,amin,amax)
    astart=αs.*1
    Estart,dunmi1,dumi2=livesycle_drift_real_alpha_set_deltauniform(β,δmax,αs,meanencounters,reps)
    for ii in 1:steps
        E0,dunmi1,dumi2=livesycle_drift_real_alpha_set_deltauniform(β,δmax,αs,meanencounters,reps)
        gradient=zeros(l)
        for jj in 1:gradients
            alocal=alpharange.(αs.+(λ1*randn(l)),amin,amax)
            E1,dunmi1,dumi2=livesycle_drift_real_alpha_set_deltauniform(β,δmax,alocal,meanencounters,reps)
            if (norm((alocal.-αs),2) >0)
                gradient=gradient.+(E1-E0)*(alocal.-αs)/norm((alocal.-αs),2)/gradients
            end
        end
        αs=alpharange.(αs.+λ2*gradient,amin,amax)
    end
    Efinal,dunmi1,dumi2=livesycle_drift_real_alpha_set_deltauniform(β,δmax,αs,meanencounters,reps)
    return αs,astart,Estart,Efinal
end

function simulation_alphaopt_rep(β,δmax,alphas,meanencounters,reps,steps,gradients,λ,λ2,rep)

  ##αs==zeros(l) gives rand alpha of langth(l)
parameters=(β,δmax,alphas,meanencounters,reps,steps,gradients,λ,λ2,rep)

    filename = "Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)/reps_$(reps)/beta_$(β)-Delmax_$(δmax)-asteps_$(alphas)-grad_$(gradients)_steps_$(steps)_l1_$(λ)_l2_$(λ2)_rep_$(rep).jld2"
    println("I want to do  $(filename)")
    if isfile(filename)
        println("Die daten $(filename) exestieren schon ich hoehre auf")
    else

        if !isdir("Data/DM_bcell_drift_alphaopt")
            mkdir("Data/DM_bcell_drift_alphaopt")
        end

        if !isdir("Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)")
            mkdir("Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)")
        end
        if !isdir("Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)/reps_$(reps)")
            mkdir("Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)/reps_$(reps)")
        end

        afin,astat,Estart,efin=optimization_set_alphas(β,δmax,zeros(alphas),meanencounters,reps,steps,gradients,λ,λ2)

        save(filename,"Astart",  astat,"Afin",afin,"Estart",Estart,"Efin",efin,"parameters",parameters)
        println("I'm saving the file $(filename)")
    end
end

function simulation_alphaopt_rep_nextsteps(β,δmax,alphas,meanencounters,reps,steps,gradients,λ,λ2,rep,startstep)

  ##αs==zeros(l) gives rand alpha of langth(l)
parameters=(β,δmax,alphas,meanencounters,reps,steps,gradients,λ,λ2,rep)

    filename = "Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)/reps_$(reps)/beta_$(β)-Delmax_$(δmax)-asteps_$(alphas)-grad_$(gradients)_steps_$(startstep)_l1_$(λ)_l2_$(λ2)_rep_$(rep).jld2"
    if isfile(filename)

        astart_original= load(filename,"Astart")
        afin_original= load(filename,"Afin")

        println("Loaded data and will start simulation")

        filename = "Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)/reps_$(reps)/beta_$(β)-Delmax_$(δmax)-asteps_$(alphas)-grad_$(gradients)_steps_$(steps)_l1_$(λ)_l2_$(λ2)_rep_$(rep).jld2"
        println("I want to do  $(filename)")
        if isfile(filename)
            println("Die daten $(filename) exestieren schon ich hoehre auf")
        else

            if !isdir("Data/DM_bcell_drift_alphaopt")
                mkdir("Data/DM_bcell_drift_alphaopt")
            end

            if !isdir("Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)")
                mkdir("Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)")
            end
            if !isdir("Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)/reps_$(reps)")
                mkdir("Data/DM_bcell_drift_alphaopt/meanencounters_$(meanencounters)/reps_$(reps)")
            end
            rounds=round(Int64,steps-startstep)
            #println("ich starte mit $((β,δmax,afin_original,meanencounters,reps,rounds,gradients,λ,λ2))")
            afin,astat,Estart,efin=optimization_set_alphas(β,δmax,afin_original,meanencounters,reps,rounds,gradients,λ,λ2)
            #println("ich bin fertig mit mit $((β,δmax,afin_original,meanencounters,reps,rounds,gradients,λ,λ2))")
        #filename = "/Users/oskar/Dropbox/PHD/Simulations/b-cell/Simulation/src/julia/data/test.jld"
            save(filename,"Astart",  astart_original,"Afin",afin,"Estart",Estart,"Efin",efin,"parameters",parameters)
            println("I'm saving the file $(filename)")
        end

    else
        println("Data $(filename) allready exists")
    end

end

# Single alpha

max_stars_repo_name: StatPhysBio/ImmuneMemoryDM

max_issues_repo_name: StatPhysBio/ImmuneMemoryDM

max_forks_repo_name: StatPhysBio/ImmuneMemoryDM


text:
using SpikingRNNs ; global const S=SpikingRNNs
using Test
using LinearAlgebra,SparseArrays,Statistics,Distributions
using Random ; Random.seed!(0)

# utility test functions

# connections
function onesparsemat(w::Real)
  mat=Matrix{Float64}(undef,1,1) ; mat[1,1]=w
  return sparse(mat)
end

# testing utility function
function sparse_thingies(w::SparseMatrixCSC)
  return nonzeros(w),rowvals(w),SparseArrays.getcolptr(w)
end

##

@testset "Connectivity matrices" begin
  m,n = (800,1000)
  μtest = 3.0
  ptest = 0.1
  σtest = 0.5
  wtest = S.sparse_wmat_lognorm(m,n,ptest,μtest,σtest)
  wtestvals = nonzeros(wtest)
  @test isapprox(mean(wtestvals),μtest;atol=0.1)
  @test isapprox(std(wtestvals),σtest;atol=0.1)
  @test tr(Diagonal(abs.(wtest))) == 0.0 # no autpses
  m,n = (80,100)
  μtest = -3.0
  ptest = 0.333
  wtest = S.sparse_wmat_lognorm(m,n,ptest,μtest,σtest;rowsum=123.0)
  @test all(isapprox.(sum(wtest;dims=2),-123.0;atol=0.001))

  wtest = S.make_sparse_weights(m,n,ptest,-0.123)
  wvals = nonzeros(wtest)
  @test all( isapprox.(wvals,-0.123) )

end

@testset "2D rate model" begin

  neuron_e = S.NTReLU(1.,1.)
  neuron_i = S.NTReLU(1.,1.)
  pse  = S.PSRate(neuron_e,1)
  psi  = S.PSRate(neuron_i,1)

  (w_ee,w_ie,w_ei,w_ii) = let w = 20. , k = 1.1
    onesparsemat.((w,w,-k*w,-k*w))
  end

  conn_ee = S.BaseConnection(w_ee)
  conn_ei = S.BaseConnection(w_ei)
  conn_ie = S.BaseConnection(w_ie)
  conn_ii = S.BaseConnection(w_ii)

  pop_e = S.Population(pse,(conn_ee,pse),(conn_ei,psi))
  pop_i = S.Population(psi,(conn_ie,pse),(conn_ii,psi))

  # initial conditions
  pse.state_now[1] = S.ioinv(10.0,pse)
  psi.state_now[1] = S.ioinv(5.0,pse)

  dt = 1E-4
  T = 5.0
  times = 0:dt:T
  ntimes = length(times)
  mynetwork = S.RecurrentNetwork(dt,(pop_e,pop_i))

  e_out = Vector{Float64}(undef,ntimes)
  i_out = Vector{Float64}(undef,ntimes)

  for (k,t) in enumerate(times)
    e_out[k] = S.iofunction(pse.state_now[1],pse)
    i_out[k] = S.iofunction(psi.state_now[1],psi)
    # rate model with constant input  does not really depend on absolute time (first argument)
    S.dynamics_step!(0.0,mynetwork)
  end

  @test all(e_out .>= 0.0)
  @test all(isfinite.(e_out))
  @test all(i_out .>= 0.0)
  @test all(isfinite.(i_out))
  # I expect aplification, let's say that 1/3 of elements are above the starting rate
  # but the last is below
  @test i_out[end]<i_out[1]
  @test e_out[end]<e_out[1]
  @test count(i_out .> i_out[1])/length(i_out) > 0.333
  @test count(e_out .> e_out[1])/length(e_out) > 0.333

end

@testset "2D rate model with input" begin

  neuron_e_and_i = S.NTReLU(1.,1.)
  pse  = S.PSRate(neuron_e_and_i,1)
  psi  = S.PSRate(neuron_e_and_i,1)

  (w_ee,w_ie,w_ei,w_ii) = (2.,2.5,-3.,-0.5)
  wmat = [w_ee  w_ei
          w_ie w_ii]
  make_connection(w)=S.BaseConnection(onesparsemat(w))
  conn_ee = make_connection(w_ee)
  conn_ei = make_connection(w_ei)
  conn_ie = make_connection(w_ie)
  conn_ii = make_connection(w_ii)

  # inputs
  h_e = 50.33
  h_i = 2.8
  h_vec = [h_e,h_i]

  fpoint = - inv(Matrix(wmat)-I)*h_vec

  ## input connection!
  in_state_e = S.PSSimpleInput(S.InputSimpleOffset(h_e))
  in_state_i = S.PSSimpleInput(S.InputSimpleOffset(h_i))
  pop_e = S.Population(pse,(conn_ee,pse),(conn_ei,psi),
      (S.FakeConnection(),in_state_e))
  pop_i = S.Population(psi,(conn_ie,pse),(conn_ii,psi),
      (S.FakeConnection(),in_state_i))
  ##
  dt = 1E-2
  T = 60.0
  times = 0:dt:T
  ntimes = length(times)
  mynetwork = S.RecurrentNetwork(dt,(pop_e,pop_i))

  ei_out = Matrix{Float64}(undef,2,ntimes)
  # initial conditions
  pse.state_now .= S.ioinv(10.0,pse)
  psi.state_now .= S.ioinv(10.0,psi)

  for (k,t) in enumerate(times)
    ei_out[1,k] = S.iofunction(pse.state_now[1],neuron_e_and_i)
    ei_out[2,k] = S.iofunction(psi.state_now[1],neuron_e_and_i)
    # rate model with constant input  does not really depend on absolute time (first argument)
    S.dynamics_step!(mynetwork)
  end
  @test all(isapprox.(ei_out[:,end],fpoint;atol=1E-1))

end

@testset "Spike inputs in single LIF neuron" begin
  dt = 1E-3
  Ttot = 10.0
  # One LIF neuron
  myτ = 0.2
  vth = 10.
  v_r = -5.0
  τrefr= 0.3 # refractoriness
  τpcd = 0.2 # post synaptic current decay
  myinput = 0.0 # constant input to E neuron
  ps_e = S.PSLIF(myτ,vth,v_r,τrefr,τpcd,1)

  # one static input
  in_state_e = S.PSSimpleInput(S.InputSimpleOffset(myinput))
  # connection will be FakeConnection()

  # let's produce a couple of trains
  train1 = let rat = 1.0
    sort(rand(Uniform(0.05,Ttot),round(Integer,rat*Ttot) ))
  end
  train2 = let rat = 0.5
    sort(rand(Uniform(0.05,Ttot),round(Integer,rat*Ttot) ))
  end
  # input population
  ps_train_in=S.PSFixedSpiketrain([train1,train2],myτ)

  # and connection object
  conn_e_in = let w_intrain2e = sparse([eps() Inf ; ])
    S.ConnSpikeTransfer(w_intrain2e)
  end

  # connected populations
  # two populations: the input population (unconnected)
  # and the E neuron connected to input
  pop_in = S.UnconnectedPopulation(ps_train_in)
  pop_e = S.Population(ps_e,(conn_e_in,ps_train_in),
    (S.FakeConnection(),in_state_e))

  # that's it, let's make the network
  myntw = S.RecurrentNetwork(dt,pop_in,pop_e)

  # record spiketimes and internal potential
  krec = 1
  rec_state_e = S.RecStateNow(ps_e,krec,dt,Ttot)
  rec_spikes_e = S.RecSpikes(ps_e,100.0,Ttot)
  rec_spikes_in = S.RecSpikes(ps_train_in,100.0,Ttot)

  ## Run

  times = (0:myntw.dt:Ttot)
  nt = length(times)
  # clean up
  S.reset!.([rec_state_e,rec_spikes_e,rec_spikes_in])
  S.reset!.([ps_e,ps_train_in])
  # initial conditions
  ps_e.state_now[1] = 0.0

  for (k,t) in enumerate(times)
    rec_state_e(t,k,myntw)
    rec_spikes_e(t,k,myntw)
    rec_spikes_in(t,k,myntw)
    S.dynamics_step!(t,myntw)
  end
  # this is useful for visualization only
  S.add_fake_spikes!(1.5vth,rec_state_e,rec_spikes_e)
  ##
  train1_sim,train2_sim = let (spkt,spkneu) = S.get_spiketimes_spikeneurons(rec_spikes_in)
    spkt[spkneu .== 1],spkt[spkneu .== 2]
  end
  train_e_sim = let (spkt,spkneu) = S.get_spiketimes_spikeneurons(rec_spikes_e)
    spkt
  end
  @test all(isapprox.(train1,train1_sim;atol=1.1dt))
  @test all(isapprox.(train2,train2_sim;atol=1.1dt))
  @test all(isapprox.(train_e_sim,sort(vcat(train1_sim,train2_sim));atol=1.1dt))
  # but the error between train_e_sim and train1, train2 is up to 2dt!
  ## part 2, expected period of LIF neuron
  dt = 5E-4
  Ttot = 3.0
  myτ = 0.1
  vth = 12.
  v_r = -6.123
  τrefr = 0.0
  τpcd = 1E10
  myinput = 14.0
  ps_e = S.PSLIF(myτ,vth,v_r,τrefr,τpcd,1)
  # create static input
  in_state_e = S.PSSimpleInput(S.InputSimpleOffset(myinput))
  # only one population: E with input
  pop_e = S.Population(ps_e,(S.FakeConnection(),in_state_e))
  # that's it, let's make the network
  myntw = S.RecurrentNetwork(dt,pop_e)

  times = (0:myntw.dt:Ttot)
  nt = length(times)
  S.expected_period_norefr(ps_e.neurontype,myinput)
  # spike recorder
  rec_spikes = let exp_freq = inv(S.expected_period_norefr(ps_e.neurontype,myinput))
    S.RecSpikes(ps_e,1.5*exp_freq,Ttot)
  end
  # reset and run
  S.reset!(rec_spikes)
  S.reset!(ps_e)
  # initial conditions
  ps_e.state_now[1] = v_r

  # run!
  for (k,t) in enumerate(times)
    rec_spikes(t,k,myntw)
    S.dynamics_step!(t,myntw)
  end
  ##
  spkt,_ = S.get_spiketimes_spikeneurons(rec_spikes)
  # period of first spike
  @test isapprox(S.expected_period_norefr(ps_e.neurontype,myinput),spkt[1] ;
    atol = 0.02)
end

@testset "Exact spike time generation" begin
  therate = 123.4
  sgentest = S.SGPoisson(therate)

  t_final = let t_current = 0.0
    for i in 1:1_000
      t_current=S._get_spiketime_update(t_current,sgentest,2)
    end
    t_current
  end
  ratenum = 1_000/t_final
  @test(isapprox(ratenum,therate;atol=5.0))
end

@testset "Inputs faster than Euler Δt" begin
  # first test , Poissoin fixed input, faster than Euler dt
  # I expect that doubling the rate will double the firing
  # careful:
  # this would not work with an exponential spike generation kernel

  dt = 0.1E-3
  myτe = 20E-3 # seconds
  τrefr= 0.1E-3 # refractoriness
  vth_e = -20.0   # mV
  Cap = 300.0 #capacitance mF
  v_rest_e = -60.0
  v_rev_e = 0.0
  v_leak_e = v_rest_e
  v_reset_e = v_rest_e

  # synaptic kernel
  taueplus = 6E-3 #e synapse decay time
  taueminus = 1E-3 #e synapse rise time

  Ne = 200
  # if spike genrator kernel is absent, then inputs scale linearly
  nt_e = let sker = S.SKExpDiff(taueplus,taueminus)
    sgen = S.SpikeGenNone()
    S.NTLIFConductance(sker,sgen,myτe,Cap,
      vth_e,v_reset_e,v_rest_e,τrefr,v_rev_e)
  end
  ps_e = S.PSLIFConductance(nt_e,Ne)

  nt_in = let in_rate = 10E3
    sker = S.SKExpDiff(taueplus,taueminus)
    sgen = S.SGPoisson(in_rate)
    S.NTInputConductance(sgen,sker,v_rev_e)
  end
  in_weight = 7.0
  ps_in = S.PSInputPoissonConductanceExact(nt_in,in_weight,Ne)

  pop_e = S.Population(ps_e,(S.FakeConnection(),ps_in))
  ntw = S.RecurrentNetwork(dt,pop_e)

  Ttot = 2.0
  # record spiketimes and internal potential
  n_e_rec = Ne
  t_wup = 0.0
  rec_spikes_e = S.RecSpikes(ps_e,500.0,Ttot;idx_save=collect(1:n_e_rec),Tstart=t_wup)

  times = (0:ntw.dt:Ttot)
  nt = length(times)
  S.reset!.([ps_e,ps_in,rec_spikes_e])
  # initial conditions
  ps_e.state_now .= v_reset_e
  for (k,t) in enumerate(times)
    rec_spikes_e(t,k,ntw)
    S.dynamics_step!(t,ntw)
  end
# _ = let plt=plot(),ts=rec_state_e.times,
#   neu = 1 , vs = rec_state_e.state_now[neu,:]
#   plot!(plt,ts,vs;leg=false,linewidth=2)
# end
  rates_e = let rdic=S.get_mean_rates(rec_spikes_e)
    ret = fill(0.0,n_e_rec)
    for (k,v) in pairs(rdic)
      ret[k] = v
    end
    ret
  end
  rates_input1 = mean(rates_e)

  # now double the input rate, see if output doubles too

  nt_in = let in_rate = 2*10E3
    sker = S.SKExpDiff(taueplus,taueminus)
    sgen = S.SGPoisson(in_rate)
    S.NTInputConductance(sgen,sker,v_rev_e)
  end
  ps_in = S.PSInputPoissonConductanceExact(nt_in,in_weight,Ne)
  pop_e = S.Population(ps_e,(S.FakeConnection(),ps_in))
  ntw = S.RecurrentNetwork(dt,pop_e)

  rec_spikes_e = S.RecSpikes(ps_e,500.0,Ttot;idx_save=collect(1:n_e_rec),Tstart=t_wup)

  S.reset!.([ps_e,ps_in,rec_spikes_e])
  # initial conditions
  ps_e.state_now .= v_reset_e

  for (k,t) in enumerate(times)
    rec_spikes_e(t,k,ntw)
    S.dynamics_step!(t,ntw)
  end

  rates_e = let rdic=S.get_mean_rates(rec_spikes_e)
    ret = fill(0.0,n_e_rec)
    for (k,v) in pairs(rdic)
      ret[k] = v
    end
    ret
  end
  rates_input2 = mean(rates_e)
  @test isapprox(2*rates_input1,rates_input2;rtol=0.2)
end

@testset "Homeostatic plasticity, easy implementation" begin

  # incoming connections : rows
  # Npost neurons have Npre incoming connections
  Npost = 13
  Npre = 17
  sum_max = Npre/2 - 0.1
  wtest_in_start = sparse(rand(Npost,Npre))
  sum_start = sum(wtest_in_start;dims=2)
  plast_method = S.HeterosynapticAdditive(-Inf,10.0)
  plast_target = S.HeterosynapticIncoming(sum_max)
  wtest_in = copy(wtest_in_start)
  S._apply_easy_het_plasticity!(wtest_in,plast_method,plast_target)
  sum_end = sum(wtest_in;dims=2)
  @test all(sum_end .<= sum_max+1E-6)
  idx_good = findall(<=(sum_max),sum_start)
  @test all(sum_end[idx_good] .== sum_start[idx_good])

  # outgoing connections, columns

  Npost = 11
  Npre = 19
  sum_max = Npost/2 - 0.1
  wtest_in_start = sparse(rand(Npost,Npre))
  sum_start = sum(wtest_in_start;dims=2)
  plast_method = S.HeterosynapticAdditive(-Inf,10.0)
  plast_target = S.HeterosynapticOutgoing(sum_max)
  wtest_in = copy(wtest_in_start)
  S._apply_easy_het_plasticity!(wtest_in,plast_method,plast_target)
  sum_end = sum(wtest_in;dims=1)
  @test all(sum_end .<= sum_max+1E-6)
  idx_good = findall(<=(sum_max),sum_start)
  @test all(sum_end[idx_good] .== sum_start[idx_good])

  # Internals for other method
  N = 400
  sp = 0.4
  testmat = sprand(N,N,sp)
  meansum = N*sp*0.5
  maxsum = meansum - 10.0

  hup = S.HetUpperLimit(maxsum,-Inf,Inf,0.0)
  hprecise = S.HetStrictSum(maxsum,-Inf,Inf,1E-4)
  hetincoming = S.HetIncoming()
  hetoutgoing = S.HetOutgoing()
  hetboth = S.HetBoth()

  idxtry = sample((1:N),200;replace=false)
  # outgoing, along columns
  tochange,sumsvals,nels =  S._hetplast_check(idxtry,testmat,
      S.HetOutgoing(),hup)
  testmatfix = copy(testmat)
  S._apply_hetplast_spiketriggered!(tochange,idxtry,sumsvals,nels,
    testmatfix, S.HetOutgoing(),S.HetAdditive(),
    hup)
  @test all(sum(testmatfix;dims=1)[idxtry] .< (maxsum + 1E-5))

  # incoming, along rows
  tochange,sumsvals,nels =  S._hetplast_check(idxtry,testmat,
      S.HetIncoming(),hprecise)

  testmatfix = copy(testmat)
  S._apply_hetplast_spiketriggered!(tochange,idxtry,sumsvals,nels,
    testmatfix, S.HetIncoming(),S.HetAdditive(),
    hprecise)
  @test all(isapprox.(sum(testmatfix;dims=2)[idxtry],maxsum;atol=1E-3))

  # both, exact sum, for a single neuron
  _lowval = 1.234
  hprecise_low = S.HetStrictSum(_lowval,-Inf,Inf,1E-4)
  # idxtry = [div(N,3),]
  idxtry = [div(N,3)]

  tochange,sumsvals,nels  =  S._hetplast_check(idxtry,testmat,
      S.HetBoth(),hprecise_low)
  testmatfix = copy(testmat)
  S._apply_hetplast_spiketriggered!(tochange,idxtry,sumsvals,nels,
    testmatfix, S.HetBoth(),S.HetAdditive(),
    hprecise_low)

  @test isapprox(sum(testmatfix[:,idxtry]),_lowval;atol=1E-4)
  @test isapprox(sum(testmatfix[idxtry,:]),_lowval;atol=1E-4)
  # TO-DO  a test with both, but where only row/col passes the limit

  # more subroutines
  N = 400
  mattest = sparse(fill(1.0,N,N))
  maxsum = N/2

  constr = S.HetStrictSum(maxsum,-Inf,Inf,1E-4)

  allocrows = fill(NaN,N)
  alloccols = fill(NaN,N)

  S._het_plasticity_fix_rows!(allocrows,fill(0,N),
      mattest,constr,S.HetAdditive(),S.HetBoth())
  S._het_plasticity_fix_cols!(alloccols,fill(0,N),
      mattest,constr,S.HetAdditive(),S.HetBoth())
  # apply the fix
  @test all( isapprox.(allocrows,-maxsum/N;atol=1E-4))
  @test all( isapprox.(alloccols,-maxsum/N;atol=1E-4))
  mattestfix = copy(mattest)
  S._het_plasticity_apply_fix!(
        allocrows,alloccols,
        mattestfix,
        constr,S.HetAdditive(),S.HetBoth())

  @test all(isapprox.(sum(mattestfix;dims=1),maxsum;atol=1E-3))
  @test all(isapprox.(sum(mattestfix;dims=2),maxsum;atol=1E-3))

end

@testset "Structural plasticity" begin
  Ntot = 1000
  Ttot = 5.0
  dt = 0.1E-3
  mydensity = 0.3
  rates = rand(Uniform(20.,50.),Ntot)

  # test 1 : only death, connections expected to drop of half
  # if death freq is log(2)/Tot  (because half-life is log(2)/freq )

  # dummy neurons
  nt_all = let sker = S.SKExp(Inf)
    spkgen = S.SGPoissonMulti(rates)
    S.NTInputConductance(spkgen,sker,-Inf)
  end
  ps_all = S.PSInputConductance(nt_all,Ntot)

  ## initial weights

  weights_start = S.sparse_constant_wmat(Ntot,Ntot,mydensity,1.0;no_autapses=true)

  ## plasticity rule
  plast_struct = let  νdeath = log(2)/(Ttot),
    Δt = 30E-3,
    ρ = 0.0
    syngen = S.SynapticGenerationConstant(1.11)
    srt_type = S.StructuralPlasticityPlain(;
      connection_density=ρ,
      death_rate = νdeath)
    S.PlasticityStructural(srt_type,syngen,Δt)
  end
  ##
  conn_all=S.ConnectionPlasticityTest(copy(weights_start),plast_struct)
  pop_all=S.Population(ps_all,(conn_all,ps_all))
  myntw = S.RecurrentNetwork(dt,pop_all)
  ## Run
  times = (0:myntw.dt:Ttot)
  nt = length(times)
  # reset weights
  copy!(conn_all.weights,weights_start)
  wstart = copy(weights_start)
  density_start = nnz(wstart)/Ntot^2
  @test isapprox(density_start,mydensity;atol=0.01)
  # clean up
  S.reset!(ps_all)
  S.reset!(conn_all)

  for (k,t) in enumerate(times)
    S.dynamics_step!(t,myntw)
  end
  wend = conn_all.weights
  density_end =  nnz(wend)/Ntot^2
  @test isapprox(density_end*2,mydensity;atol=0.03)

  # test 2 : 1/10 of connections survive, but density is preserved
  # so I use νdeath = log(10)/(Ttot)

  nt_all = let sker = S.SKExp(Inf)
    spkgen = S.SGPoissonMulti(rates)
    S.NTInputConductance(spkgen,sker,-Inf)
  end
  ps_all = S.PSInputConductance(nt_all,Ntot)

  ## initial weights
  weights_start = S.sparse_constant_wmat(Ntot,Ntot,mydensity,1.0;no_autapses=true)

  ## plasticity rule
  plast_struct = let  νdeath = log(10)/(Ttot),
    Δt = 30E-3,
    ρ = mydensity,
    syngen = S.SynapticGenerationConstant(2.0)
    srt_type = S.StructuralPlasticityPlain(;
      connection_density=ρ,
      death_rate = νdeath)
    S.PlasticityStructural(srt_type,syngen,Δt)
  end

  conn_all=S.ConnectionPlasticityTest(copy(weights_start),plast_struct)
  pop_all=S.Population(ps_all,(conn_all,ps_all))
  myntw = S.RecurrentNetwork(dt,pop_all)

  times = (0:myntw.dt:Ttot)
  nt = length(times)
  # reset weights
  copy!(conn_all.weights,weights_start)
  wstart = copy(weights_start)
  density_start = nnz(wstart)/Ntot^2
  # clean up
  S.reset!(ps_all)
  S.reset!(conn_all)

  for (k,t) in enumerate(times)
    S.dynamics_step!(t,myntw)
  end
  wend = conn_all.weights
  density_end =  nnz(wend)/Ntot^2
  @test isapprox(density_end,density_start;atol=0.01)
  n_start = count(==(1.0),nonzeros(wstart))
  @test n_start == nnz(wstart)
  n_kept = count(==(1.0),nonzeros(wend))
  density_survivors = n_kept/Ntot^2
  @test isapprox(density_survivors*10,density_start;atol=0.02)
end

ext: jl


text:
"""
	DataCell{T,S,U}(x::T, y::S, f::U, layer::Int, tinfo::String, oinfo::String)

Constructs the basic object for holding data, the `DataCell`. The constructor is parametric so it can be
easily extended to new data types.

# Main `DataCell` constructors
  * `DataCell(x::AbstractArray, tinfo::String="")` creates an `unlabeled` object.
  * `DataCell(x::AbstractArray, y::AbstractVector, tinfo::String="")` creates a `labeled` object. Useful for classification.
  * `DataCell(x::AbstractArray, y::AbstractMatrix, tinfo::String="")` creates a `multi-labeled` object. Useful for classification/regression

For the creation of data cells, the `datacell` function is quite well suited.
"""
#############
# DataCells #
#############
struct DataCell{T,S,U} <: AbstractCell{T,S,U}
	x::T                                                                                # Data
    	y::S                                                                                # Labels/probabilities etc.
    	f::U                                                                                # Usually a function
	layer::Int										# Layer
    	tinfo::String										# Title or name
    	oinfo::String                                                                        	# Other information
end

DataCell(x::T where T<:AbstractArray, tinfo::String = "") = DataCell(x, nothing, nothing, 0, tinfo, oinfoglobal)
DataCell(x::T where T<:AbstractArray, y::S where S<:Void, tinfo::String = "") = DataCell(x, nothing, nothing, 0, tinfo, oinfoglobal)
DataCell(x::T where T<:AbstractArray, y::S where S<:AbstractArray{Void}, tinfo::String = "") = DataCell(x, nothing, nothing, 0, tinfo, oinfoglobal)
DataCell(x::T where T<:AbstractArray, y::S where S<:AbstractArray, tinfo::String = "") = DataCell(x, checktargets(x,y), nothing, 0, tinfo, oinfoglobal)

####################
# DataCell aliases #
####################
const DataElement{T<:Union{Number, AbstractString, Symbol}} = Union{T,Nullable{T}}
const FixedElement{T<:Union{Number, AbstractString, Symbol}} = T
const NullableElement{T<:Union{Number, AbstractString, Symbol}} = Nullable{T}
const CellDataVec{T<:AbstractVector} = DataCell{T, Void, Void}
const CellDataVecVec{T<:AbstractVector, S<:AbstractVector} = DataCell{T, S, Void}
const CellDataVecMat{T<:AbstractVector, S<:AbstractMatrix} = DataCell{T, S, Void}
const CellDataMat{T<:AbstractMatrix} = DataCell{T, Void, Void}
const CellDataMatVec{T<:AbstractMatrix, S<:AbstractVector} =  DataCell{T, S, Void}
const CellDataMatMat{T<:AbstractMatrix, S<:AbstractMatrix} =  DataCell{T, S, Void}
const CellDataVectors{T<:AbstractVector, S<:Union{AbstractArray, Void}} = DataCell{T, S, Void}
const CellDataMatrices{T<:AbstractMatrix, S<: Union{AbstractArray, Void}} = DataCell{T, S, Void}
const CellDataU{T<:AbstractArray} = DataCell{T, Void, Void} 					# Data cell, unlabeled
const CellDataL{T<:AbstractArray, S<:AbstractVector} = DataCell{T, S, Void}			# Data cell labeled
const CellDataLL{T<:AbstractArray, S<:AbstractMatrix} = DataCell{T, S, Void}			# Data cell multiple labels
const CellData{T<:AbstractArray, S} = DataCell{T, S, Void}
#const TrainableData{U <:Real, T<:AbstractMatrix{U}, S<:AbstractArray} = DataCell{T, S, Void}

##############################################################################################################################
# Operators [Data cells]												     #
##############################################################################################################################

# Vertical datacell concatenation (e.g. data concatenation) and auxiliary methods
dcat(x::T where T<:PTuple{<:Void}) = nothing
dcat(::Void) = nothing
dcat(::Void,::Void) = nothing
dcat(x,::Void) = dcat(x)
dcat(::Void,x) = dcat(x)
dcat(x::AbstractArray) = getobs(x)
dcat(x::Array) = x
dcat(x::AbstractVector, y::AbstractVector) = vcat(mat(x, ObsDim.Constant{2}()), mat(y,ObsDim.Constant{2}()))
dcat(x::AbstractMatrix, y::AbstractVector) = vcat(dcat(x), mat(y, ObsDim.Constant{2}()))
dcat(x::AbstractVector, y::AbstractMatrix) = vcat(mat(x, ObsDim.Constant{2}()), dcat(y))
dcat(x::AbstractMatrix, y::AbstractMatrix) = vcat(dcat(x), dcat(y))
dcat(x,y,z...) = dcat(dcat(x,y),z...)
dcat(x) = dcat(x...)
vcat(c::T...) where T<:CellData = begin
	@assert all(labx == laby for labx in gety!.(c), laby in gety!.(c)) "[vcat] 'y' fields have to be identical for all DataCells."
	datacell(dcat(getx!.(c)), gety!(c[1]))
end

# Horizontal datacell concatenation (e.g. observation concatenation) and auxiliary methods
ocat(x::T where T<:PTuple{<:Void}) = nothing
ocat(::Void) = nothing
ocat(::Void,::Void) = nothing
ocat(x,::Void) = ocat(x)
ocat(::Void,x) = ocat(x)
ocat(x::AbstractArray) = getobs(x)
ocat(x::Array) = x
ocat(x::AbstractVector, y::AbstractVector) = vcat(ocat(x), ocat(y))
ocat(x::AbstractArray, y::AbstractArray) = hcat(ocat(x), ocat(y))
ocat(x,y,z...) = ocat(ocat(x,y),z...)
ocat(x) = ocat(x...)
hcat(c::T...) where T<:CellData = datacell(ocat(getx!.(c)), ocat(gety!.(c)))

# Piping operators
|>(x::T where T<:AbstractArray, c::S where S<:CellData) = begin
    	info("[operators] Creating new data cell with Array contents..." )
    	datacell(x)
end

|>(x::T where T<:Tuple, c::S where S<:CellData) = begin 					# Generally used in the case: (data, labels ) |> datacell([])
	if (nobs(c) > 0)
		info("[operators] Creating new data cell appending Tuple contents..." )
		datacell(dcat(x[1],getx!(c)), dcat(x[2], gety!(c)))
	else
		info("[operators] Creating new data cell with Tuple contents..." )
		datacell(x...)
	end
end

|>(x1::T where T<:CellDataU, x2::S where S<:CellDataU) = begin                 			# Generally used to glue together two DataCells
	@assert(nobs(x1) == nobs(x2))
	datacell(dcat(getx!(x1),getx!(x2)))
end

|>(x1::T where T<:CellDataU, x2::S where S<:CellData) = begin                 			# Generally used to glue together two DataCells
	@assert(nobs(x1) == nobs(x2))
	datacell(dcat(getx!(x1),getx!(x2)), gety!(x2))
end

|>(x1::T where T<:CellData, x2::S where S<:CellDataU) = begin                 			# Generally used to glue together two DataCells
	@assert(nobs(x1) == nobs(x2))
	datacell(dcat(getx!(x1),getx!(x2)), gety!(x1))
end

|>(x1::T where T<:CellData, x2::S where S<:CellData) = begin                 			# Generally used to glue together two DataCells
	@assert(nobs(x1) == nobs(x2))
	if isequal(gety!(x1), gety!(x2))
		datacell(dcat(getx!(x1),getx!(x2)), gety!(x1))
	else
		datacell(dcat(getx!(x1),getx!(x2)), dcat(gety!(x1),gety!(x2)))
	end
end

###########################
# Indexing for data cells #
###########################

# Generic setindex! for data cells
setindex!(c::T where T<:CellData, data, inds...) = setindex!(getx!(c), data, inds...)

# Generic indexing operator, no argument
getindex(c::CellData) = j4pr.datacell(eltype(getx!(c))[])
getindex(c::CellData, ::Colon) = c
getindex(c::CellData, ::Colon, ::Colon) = c

# Single integer indexing (select sample)
getindex(c::Union{CellDataVec, CellDataVecVec}, i::Int64) = DataCell([getx!(c)[i]], [gety!(c)[i]] )
getindex(c::CellDataVecMat, i::Int64) = DataCell([getx!(c)[i]], gety!(c)[:,i:i] )
getindex(c::Union{CellDataMat, CellDataMatVec}, i::Int64) = DataCell(getx!(c)[:,i:i], [gety!(c)[i]] )
getindex(c::CellDataMatMat, i::Int64) = DataCell(getx!(c)[:,i:i], gety!(c)[:,i:i] )

# Integer + Colon indexing (select variable) - only for matrix data
getindex(c::CellDataMatrices , i::Int64, ::Colon) = DataCell(getx!(c)[i:i,:], gety!(c) )

# Colon + Integer indexing (select sample) - only for matrix data
getindex(c::Union{CellDataMat, CellDataMatVec}, ::Colon, i::Int64) = DataCell(getx!(c)[:,i:i], [gety!(c)[i]])
getindex(c::CellDataMatMat, ::Colon, i::Int64) = DataCell(getx!(c)[:,i:i], gety!(c)[:,i:i])

# Range/Vector indexing (select samples) - only for vector data
getindex(c::Union{CellDataVec, CellDataVecVec}, i::Union{UnitRange{Int64}, Vector{Int64}}) = DataCell(getx!(c)[i], gety!(c)[i])
getindex(c::CellDataVecMat, i::Union{UnitRange{Int64}, Vector{Int64}}) = DataCell(getx!(c)[i], gety!(c)[:,i])

# Range/Vector + Colon indexing (select variables) - only for matrix data
getindex(c::Union{CellDataMat, CellDataMatVec}, i::Union{UnitRange{Int64}, Vector{Int64}}, ::Colon) = DataCell(getx!(c)[i,:], gety!(c)[:])
getindex(c::CellDataMatMat, i::Union{UnitRange{Int64}, Vector{Int64}}, ::Colon) = DataCell(getx!(c)[i,:], gety!(c)[:,:])

# Colon + Range/Vector indexing (select samples) - only for matrix data
getindex(c::Union{CellDataMat, CellDataMatVec}, ::Colon, i::Union{UnitRange{Int64}, Vector{Int64}}) = DataCell(getx!(c)[:,i], gety!(c)[i])
getindex(c::CellDataMatMat, ::Colon, i::Union{UnitRange{Int64}, Vector{Int64}}) = DataCell(getx!(c)[:,i], gety!(c)[:,i])

# Range/Vector + Range/Vector indexing (select samples and variables) - only for matrix data
getindex(c::Union{CellDataMat, CellDataMatVec}, i::Union{UnitRange{Int64}, Vector{Int64}}, j::Union{UnitRange{Int64}, Vector{Int64}}) = DataCell(getx!(c)[i,j], gety!(c)[j])
getindex(c::CellDataMatMat, i::Union{UnitRange{Int64}, Vector{Int64}}, j::Union{UnitRange{Int64}, Vector{Int64}}) = DataCell(getx!(c)[i,j], gety!(c)[:,j])

# Getindex for void arguments (needed to get indexes for CellDataU objects)
getindex(c::T where T<:Void, args...) = nothing

######################################
# Iteration interface for data cells #
######################################
start(c::CellData) = 1
next(c::CellData, state) = (datasubset(c, state), state+1)
done(c::CellData, state) = state > nobs(c)
endof(c::CellData) = nobs(c)
eltype(c::CellData) = typeof(c)
length(c::CellData) = nobs(c)

"""
	datacell([x,y,name=""])

Calls the DataCell constructor for creating `data cells`. It is the easyest way of creating such
objects from `AbstractArrays`. Consult `methods(datacell)` for a full list of available methods.
If no arguments are provided, a `FunctionCell` that expects an input compatible with the `DataCell`
constructor is returned.
"""
# Tier 1 - low level
datacell(x::T where T<:DataElement; name = "") = DataCell([x], name)
datacell(x::T where T<:DataElement, y::S where S<:DataElement; name = "") = DataCell([x], [y], name)
datacell(x::T where T<:DataElement, y::S where S<:AbstractArray; name = "") = DataCell([x], y, name)
datacell(x::T where T<:AbstractArray, y::S where S<:DataElement; name = "") = DataCell(x, [y], name)
datacell(x::T where T<:AbstractArray; name = "") = DataCell(x, name)
datacell(x::T where T<:AbstractArray, ::Void; name = "") =	DataCell(x, name)
datacell(x::T where T<:AbstractArray, y::S where S<:AbstractArray; name = "") = DataCell(x, y, name)
datacell(t :: Tuple{T} where T<:AbstractArray; name = "") = DataCell(t[1], name)
datacell(t :: Tuple{T,S} where {T<:AbstractArray, S<:AbstractArray}; name = "") = DataCell(t[1], t[2], name)
datacell(t :: Tuple{T,S} where {T<:AbstractArray,S}; name = "") = DataCell(t[1], name)

# Tier 2 - j4pr
datacell(x::T where T<:CellData; name = x.tinfo) = DataCell(getx!(x), gety!(x), name)
datacell(x::T where T<:CellData, y::S where S<:CellData; name = "") = datacell(getx!(x), getx!(y), name)

# DataCell creator
datacell() = FunctionCell((args...)->j4pr.DataCell(args...), (), "DataCell creator")


text:
using DotEnv
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

println("Testing DotEnv")

const dir = dirname(@__FILE__)

@testset "basic" begin
    #basic input
    str = "BASIC=basic"
    file = joinpath(dir, ".env.override")
    file2 = joinpath(dir, ".env")

    #iobuffer, string, file
    @test DotEnv.parse(str) == Dict("BASIC"=>"basic")
    @test DotEnv.parse(read(file)) == Dict("CUSTOMVAL123"=>"yes","USER"=>"replaced value")
    @test DotEnv.config(file) == Dict("CUSTOMVAL123"=>"yes","USER"=>"replaced value")

    #should trigger a warning too, but I cant test that
    @test DotEnv.config("inexistentfile.env") == nothing

    #length of returned values
    @test length(DotEnv.config(file2)) === 10

    #shouldn't replace ENV vars
    previous_value = ENV["USER"]
    cfg = DotEnv.config(file)

    @test ENV["USER"] != cfg["USER"]
    @test ENV["USER"] == previous_value

    #appropiately laoded into ENV if CUSTOM_VAL is non existent
    @test ENV["CUSTOMVAL123"] == "yes"

    #test alias
    @test DotEnv.load(file) == DotEnv.config(file)
end

@testset "parsing" begin

    #comment
    @test DotEnv.parse("#HIMOM") == Dict()

    #spaces without quotes
    @test begin
        p = DotEnv.parse("TEST=hi  the  re")
        count(c -> c == ' ', collect(p["TEST"])) == 4
    end

    #single quotes
    @test DotEnv.parse("TEST=''")["TEST"] == ""
    @test DotEnv.parse("TEST='something'")["TEST"] == "something"

    #double quotes
    @test DotEnv.parse("TEST=\"\"")["TEST"] == ""
    @test DotEnv.parse("TEST=\"something\"")["TEST"] == "something"

    #inner quotes are mantained
    @test DotEnv.parse("TEST=\"\"json\"\"")["TEST"] == "\"json\""
    @test DotEnv.parse("TEST=\"'json'\"")["TEST"] == "'json'"
    @test DotEnv.parse("TEST=\"\"\"")["TEST"] == "\""
    @test DotEnv.parse("TEST=\"'\"")["TEST"] == "'"

    #line breaks
    @test DotEnv.parse("TEST=\"\\n\"")["TEST"] == "" #It's null because of final trim
    @test DotEnv.parse("TEST=\"\\n\\n\\nsomething\"")["TEST"] == "something"
    @test DotEnv.parse("TEST=\"something\\nsomething\"")["TEST"] == "something\nsomething"
    @test DotEnv.parse("TEST=\"something\\n\\nsomething\"")["TEST"] == "something\n\nsomething"
    @test DotEnv.parse("TEST='\\n'")["TEST"] == "\\n"
    @test DotEnv.parse("TEST=\\n")["TEST"] == "\\n"

    #empty vars
    @test DotEnv.parse("TEST=")["TEST"] == ""

    #trim spaces with and without quotes
    @test DotEnv.parse("TEST='  something  '")["TEST"] == "something"
    @test DotEnv.parse("TEST=\"  something  \"")["TEST"] == "something"
    @test DotEnv.parse("TEST=  something  ")["TEST"] == "something"
    @test DotEnv.parse("TEST=    ")["TEST"] == ""
end


text:
module Testing

using ..NodalDG
using Test

VX = [0.0, 1.0, 2.0, 3.0]
EToV = [1 2; 2 3; 3 4]

g = NodalDG.Grid1D(NodalDG.Mesh1D(VX, EToV), NodalDG.gaussLobatto1D(2))

println(g.vmapM, g.vmapP)

# Test based on p58/9 in book. Note errata on website
EToV = [1 2; 2 3; 3 5; 5 4]

end


text:
##
path = @__DIR__()
using DelimitedFiles

##
mutable struct IntCode
    debug::Bool
    state::Union{Int,Nothing}
    pc::Int
    tape::Vector{Int}
    IntCode(tape::Vector{Int}; debug = false) = new(debug, 0, 0, copy(tape))
end
function debug_info(ic::IntCode, i::Int, res)
    println("IntCode ($(ic.state)): pc=$(ic.pc); i=$i; res=$res")
end
function debug_failure(ic::IntCode, i::Int, res)
    println("IntCode ($(ic.state)) FAILED: pc=$(ic.pc); i=$i; res=$res")
end
function debug_failure(ic::IntCode, msg::AbstractString)
    println("IntCode ($(ic.state)) FAILED: pc=$(ic.pc) -- $msg")
end

function add(ic::IntCode)
    address1 = ic.tape[ic.pc+2]+1
    address2 = ic.tape[ic.pc+3]+1
    address3 = ic.tape[ic.pc+4]+1
    x = ic.tape[address1]
    y = ic.tape[address2]
    res = x + y
    ic.tape[address3] = res
    if ic.debug
        println("IntCode ($(ic.state)) opcode[$(ic.pc+1)]: (1, $address1, $address2, $address3)")
    end
    res
end
function mul(ic::IntCode)
    address1 = ic.tape[ic.pc+2]+1
    address2 = ic.tape[ic.pc+3]+1
    address3 = ic.tape[ic.pc+4]+1
    x = ic.tape[address1]
    y = ic.tape[address2]
    res = x * y
    ic.tape[address3] = res
    if ic.debug
        println("IntCode ($(ic.state)) opcode[$(ic.pc+1)]: (2, $address1, $address2, $address3)")
    end
    res
end
function stop(ic::IntCode)
    if ic.debug
        println("IntCode ($(ic.state)) opcode[$(ic.pc+1)]: (99, -, -, -)")
    end
    ic.state = -1
end
function fail(ic::IntCode)
    opcode = ic.tape[ic.pc+1]
    if ic.debug
        println("IntCode ($(ic.state)) opcode[$(ic.pc+1)]: ($opcode, -, -, -)")
    end
    ic.state = nothing
end

function (ic::IntCode)()
    i = ic.tape[ic.pc+1]
    res = ic(i)
    i, res
end
function (ic::IntCode)(i::Int)
    i == 1 && return add(ic)
    i == 2 && return mul(ic)
    i == 99 && return stop(ic)
    fail(ic)
end

function increment(ic::IntCode)
    if ic.pc >= length(ic.tape)
        ic.state = nothing
        if ic.debug; debug_failure(ic, "pc past code") end
    end
    ic.pc += 4
end

isrunning(ic::IntCode) = !isnothing(ic.state) && ic.state != -1
issuccess(ic::IntCode) = !isnothing(ic.state)

function evaluate_intcode(ic::IntCode)::Union{Nothing,Some{Vector{Int}}}
    while isrunning(ic)
        i, res = ic()
        if ic.debug; debug_info(ic, i, res) end
        if !issuccess(ic); debug_failure(ic, i, res) end
        increment(ic)
    end
    !issuccess(ic) && return nothing
    Some(ic.tape)
end

function set_input(ic::IntCode, noun::Int, verb::Int)
    ic.tape[2] = noun
    ic.tape[3] = verb
end

function get_output(ic::IntCode)
    ic.tape |> first
end

##
test_cases = (
    ([1,0,0,0,99],                    [2,0,0,0,99]),
    ([2,3,0,3,99],                    [2,3,0,6,99]),
    ([2,4,4,5,99,0],                  [2,4,4,5,99,9801]),
    ([1,1,1,4,99,5,6,0,99],           [30,1,1,4,2,5,6,0,99]),
    ([1,9,10,3,2,3,11,0,99,30,40,50], [3500,9,10,70,2,3,11,0,99,30,40,50]),
)
function test_intcode(tests)
    println("Run test for IntCode")
    res = all(
        map(tests) do (i,o)
            ic = IntCode(i; debug=true)
            println("\nTesting $i => $o...")
            evaluate_intcode(ic)
            res = ic.tape == o
            if res
                println("PASSED")
            else
                println("FAILED: $(ic.tape)")
            end
            res
        end
    )
    if res
        println("\nAll tests PASSED!")
    else
        println("\nTest FAILED")
    end
end

##
test_intcode(test_cases)

##
init_data = readdlm(joinpath(path, "input.txt"), ',', Int) |> vec;

##
noun = 12
verb = 2
output = 19690720
println("noun=$noun; verb=$verb; res=$(100noun+verb)")
data = copy(init_data)
data[2] = noun
data[3] = verb
ic = IntCode(data)
evaluate_intcode(ic)
println("IntCode result: $(first(ic.tape))")

##
function find_noun_verb(init_data::Vector{Int}, output::Int)
    nouns = 0:99
    verbs = 0:99
    found = nothing
    for noun in nouns, verb in verbs
        ic = IntCode(init_data)
        set_input(ic, noun, verb)
        evaluate_intcode(ic)
        if get_output(ic) == output
            found = Some(100noun + verb)
            break
        end
    end
    return found
end

##
res = find_noun_verb(init_data, output)
if isnothing(res)
    println("Found nothing")
else
    println("Found: $res")
end


text:
using Test
using LinearAlgebra
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using HDGElasticity

function allapprox(v1,v2)
    return all(v1 .≈ v2)
end

function allapprox(v1,v2,atol)
    @assert length(v1) == length(v2)
    return all([isapprox(v1[i],v2[i],atol=atol) for i = 1:length(v1)])
end

vbasis = TensorProductBasis(2,1)
sbasis = TensorProductBasis(1,1)
squad = tensor_product_quadrature(1,2)

LH = zeros(1,1)
NED = [ones(1,1)]
linemap = HDGElasticity.LineMap([-1.,-1.],[1.,-1.])
vfunc(x) = [(1.0-x[1])*(1.0-x[2])/4.0]
sfunc(x) = [(1.0-x[1])/2.0]
HDGElasticity.LHop!(LH,vfunc,sfunc,squad,linemap,NED,1.0,1)
@test allapprox(LH,[2/3])

fill!(LH,0.0)
HDGElasticity.LHop!(LH,vfunc,sfunc,squad,linemap,NED,2.0,1)
@test allapprox(LH,[4/3])

imap = InterpolatingPolynomial(2,sbasis)
update!(imap,[0.,-1.,-1.,0.])
vf(x) = [x[1] + x[2]]
sf(x) = [x]
scale = ones(length(squad))
LH = zeros(1,1)
normals = 1.0/sqrt(2)*ones(2,length(squad))
HDGElasticity.LHop!(LH,vf,sf,squad,imap,normals,NED,scale,1)
@test allapprox(LH,[0.0],1e-15)

facemaps = HDGElasticity.reference_cell_facemaps(2)
LH = zeros(4,2)
testjac = 3.0
HDGElasticity.LHop!(LH,vbasis,sbasis,squad,facemaps[1],NED,testjac,1)
LHtest =  testjac*[2/3  1/3
                   0.0  0.0
                   1/3  2/3
                   0.0  0.0]
@test allapprox(LH,LHtest)

LH = zeros(4,2)
HDGElasticity.LHop!(LH,vbasis,sbasis,squad,facemaps[2],NED,1.0,1)
LHtest = [0.0   0.0
          0.0   0.0
          2/3   1/3
          1/3   2/3]
@test allapprox(LH,LHtest)

LH = zeros(4,2)
HDGElasticity.LHop!(LH,vbasis,sbasis,squad,facemaps[3],NED,1.0,1)
LHtest =  [0.0  0.0
           2/3  1/3
           0.0  0.0
           1/3  2/3]
@test allapprox(LH,LHtest)

LH = zeros(4,2)
HDGElasticity.LHop!(LH,vbasis,sbasis,squad,facemaps[4],NED,1.0,1)
LHtest =  [2/3  1/3
           1/3  2/3
           0.0  0.0
           0.0  0.0]
@test allapprox(LH,LHtest)

normals = repeat([1.,1.]/sqrt(2.),inner=(1,length(squad)))
Dhalf = diagm(ones(3))
cellmap = HDGElasticity.CellMap([-1.,-1.],[1.,1.])
imap = InterpolatingPolynomial(2,sbasis)
update!(imap,[0.,-1.,-1.,0.])
LH = HDGElasticity.LHop_on_interface(vbasis,sbasis,squad,imap,normals,Dhalf,cellmap)
@test size(LH) == (12,4)

Dhalf = diagm(ones(3))
cellmap = HDGElasticity.CellMap([0.,0.],[2.,1.])
facequads = [squad,squad,squad,squad]
normals = HDGElasticity.reference_normals()
LH = HDGElasticity.LHop(vbasis,sbasis,facequads,facemaps,normals,Dhalf,cellmap)
@test all([size(L) == (12,4) for L in LH])

UH = zeros(1,1)
vfunc(x) = [(1.0-x[1])*(1.0-x[2])/4.0]
sfunc(x) = [(1.0-x[1])/2.0]
linemap = HDGElasticity.LineMap([-1.,-1.],[1.,-1.])
HDGElasticity.UHop!(UH,vfunc,sfunc,squad,linemap,1.0,1)
@test allapprox(UH,[2/3])

UH = zeros(1,1)
linemap = HDGElasticity.LineMap([0.,-1.],[1.,-1.])
HDGElasticity.UHop!(UH,vfunc,sfunc,squad,linemap,1.0,1)
@test allapprox(UH,[1/6])

UH = zeros(4,2)
testjac = 2.0
linemap = HDGElasticity.LineMap([1.0,-1.0],[1.0,1.0])
HDGElasticity.UHop!(UH,vbasis,sbasis,squad,linemap,testjac,1)
UHtest = testjac*[0.0  0.0
                   0.0  0.0
                   2/3  1/3
                   1/3  2/3]
@test allapprox(UH,UHtest)

vf(x) = [x[1] + x[2]]
sf(x) = [x]
scale = ones(length(squad))
imap = InterpolatingPolynomial(2,1,1)
update!(imap,[0.,-1.,-1.,0.])
UH = reshape([0.0],1,1)
HDGElasticity.UHop!(UH,vf,sf,squad,imap,scale,1)
@test allapprox(UH,[0.0],1e-15)

UH = HDGElasticity.UHop(vbasis,sbasis,facequads,facemaps,1.0,cellmap)
@test all([size(i) == (8,4) for i in UH])

cellmap = HDGElasticity.CellMap([-1.,-1.],[1.,1.])
update!(imap,[0.,-1.,-1.,0.])
normals = reshape(1/sqrt(2)*ones(2*length(squad)),2,:)
UH = HDGElasticity.UHop_on_interface(vbasis,sbasis,squad,imap,normals,1.,cellmap)
@test size(UH) == (8,4)

normals = HDGElasticity.reference_normals()
lhc = HDGElasticity.local_hybrid_operator(vbasis,sbasis,facequads,
    facemaps,normals,Dhalf,1.,cellmap)
@test length(lhc) == 4
@test all([size(lhc[i]) == (20,4) for i = 1:4])

update!(imap,[0.,-1.,-1.,0.])
normals = 1/sqrt(2)*ones(2,length(squad))
lhc = HDGElasticity.local_hybrid_operator_on_interface(vbasis,sbasis,squad,
    imap,normals,Dhalf,1.0,cellmap)
@test size(lhc) == (20,4)

function plane_distance_function(coords,n,x0)
    return [n'*(coords[:,idx]-x0) for idx in 1:size(coords)[2]]
end

mesh = UniformMesh([0.,0.],[2.,1.],[1,1])
poly = InterpolatingPolynomial(1,2,1)
coords = HDGElasticity.nodal_coordinates(mesh,poly.basis)
coeffs = reshape(plane_distance_function(coords,[1.,1.]/sqrt(2.),[0.5,0.]),4,1)
dgmesh = HDGElasticity.DGMesh(mesh,coeffs,poly)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,1,4,coeffs,poly)
Dhalf = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(1.,2.))
lhc = HDGElasticity.local_hybrid_operator(ufs.vbasis,ufs.sbasis,
    ufs.fquads[2,1],dgmesh.facemaps,ufs.fnormals,Dhalf,1.0,cellmap)
@test length(lhc) == 2
@test all([size(lhc[i]) == (20,4) for i = 1:2])


text:
function ComputeFCFV(mesh, sex, sey, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, tau)
    ae = zeros(mesh.nel)
    be = zeros(mesh.nel,2)
    ze = zeros(mesh.nel,2,2)

    # Assemble FCFV elements
    @tturbo for iel=1:mesh.nel

        be[iel,1] += mesh.vole[iel]*sex[iel]
        be[iel,2] += mesh.vole[iel]*sey[iel]

        for ifac=1:mesh.nf_el

            nodei = mesh.e2f[iel,ifac]
            bc    = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau,dAi,mesh.vole[iel],mesh.type,mesh.ke[iel])                              # Stabilisation parameter for the face

            # Assemble
            ze[iel,1,1] += (bc==1) * dAi*ni_x*VxDir[nodei] # Dirichlet
            ze[iel,1,2] += (bc==1) * dAi*ni_x*VyDir[nodei] # Dirichlet
            ze[iel,2,1] += (bc==1) * dAi*ni_y*VxDir[nodei] # Dirichlet
            ze[iel,2,2] += (bc==1) * dAi*ni_y*VyDir[nodei] # Dirichlet
            be[iel,1]   += (bc==1) * dAi*taui*VxDir[nodei] # Dirichlet
            be[iel,2]   += (bc==1) * dAi*taui*VyDir[nodei] # Dirichlet
            ae[iel]     +=           dAi*taui

        end
    end
    return ae, be, ze
end

#--------------------------------------------------------------------#

function ComputeElementValues(mesh, Vxh, Vyh, Pe, ae, be, ze, VxDir, VyDir, tau)
    Vxe         = zeros(mesh.nel);
    Vye         = zeros(mesh.nel);
    Txxe        = zeros(mesh.nel);
    Tyye        = zeros(mesh.nel);
    Txye        = zeros(mesh.nel);

    @tturbo for iel=1:mesh.nel

        Vxe[iel]  =  be[iel,1]/ae[iel]
        Vye[iel]  =  be[iel,2]/ae[iel]
        Txxe[iel] =  mesh.ke[iel]/mesh.vole[iel]*ze[iel,1,1]
        Tyye[iel] =  mesh.ke[iel]/mesh.vole[iel]*ze[iel,2,2]
        Txye[iel] =  mesh.ke[iel]/mesh.vole[iel]*0.5*(ze[iel,1,2]+ze[iel,2,1])

        for ifac=1:mesh.nf_el

            # Face
            nodei = mesh.e2f[iel,ifac]
            bc    = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau,dAi,mesh.vole[iel],mesh.type,mesh.ke[iel])      # Stabilisation parameter for the face

            # Assemble
            Vxe[iel]  += (bc!=1) *  dAi*taui*Vxh[nodei]/ae[iel]
            Vye[iel]  += (bc!=1) *  dAi*taui*Vyh[nodei]/ae[iel]
            Txxe[iel] += (bc!=1) *  mesh.ke[iel]/mesh.vole[iel]*dAi*ni_x*Vxh[nodei]
            Tyye[iel] += (bc!=1) *  mesh.ke[iel]/mesh.vole[iel]*dAi*ni_y*Vyh[nodei]
            Txye[iel] += (bc!=1) *  mesh.ke[iel]*0.5*( 1.0/mesh.vole[iel]*dAi*( ni_x*Vyh[nodei] + ni_y*Vxh[nodei] ) )
         end
        Txxe[iel] *= 2.0
        Tyye[iel] *= 2.0
        Txye[iel] *= 2.0
    end
    return Vxe, Vye, Txxe, Tyye, Txye
end

#--------------------------------------------------------------------#

function ElementAssemblyLoop(mesh, ae, be, ze, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, gbar, tau)
    # Assemble element matrices and rhs
    f   = zeros(mesh.nf)
    Kuu = zeros(mesh.nel, 2*mesh.nf_el, 2*mesh.nf_el)
    fu  = zeros(mesh.nel, 2*mesh.nf_el)
    Kup = zeros(mesh.nel, 2*mesh.nf_el);
    fp  = zeros(mesh.nel)

    @tturbo for iel=1:mesh.nel

        for ifac=1:mesh.nf_el

            nodei = mesh.e2f[iel,ifac]
            bci   = mesh.bc[nodei]
            dAi   = mesh.dA[iel,ifac]
            ni_x  = mesh.n_x[iel,ifac]
            ni_y  = mesh.n_y[iel,ifac]
            taui  = StabParam(tau,dAi,mesh.vole[iel], mesh.type,mesh.ke[iel])

            for jfac=1:mesh.nf_el

                nodej = mesh.e2f[iel,jfac]
                bcj   = mesh.bc[nodej]
                dAj   = mesh.dA[iel,jfac]
                nj_x  = mesh.n_x[iel,jfac]
                nj_y  = mesh.n_y[iel,jfac]
                tauj  = StabParam(tau,dAj,mesh.vole[iel], mesh.type,mesh.ke[iel])

                # Delta
                del = 0.0 + (ifac==jfac)*1.0

                # Element matrix
                nitnj = ni_x*nj_x + ni_y*nj_y;
                Ke_ij =-dAi * (1.0/ae[iel] * dAj * taui*tauj - mesh.ke[iel]/mesh.vole[iel]*dAj*nitnj - taui*del)
                yes   = (bci!=1) & (bcj!=1)
                Kuu[iel,ifac,           jfac           ] = yes * Ke_ij
                Kuu[iel,ifac+mesh.nf_el,jfac+mesh.nf_el] = yes * Ke_ij
                # Kuu[iel,ifac,           jfac           ] = (bci!=1) * (bcj!=1) * Ke_ij
                # Kuu[iel,ifac+mesh.nf_el,jfac+mesh.nf_el] = (bci!=1) * (bcj!=1) * Ke_ij
            end
            # RHS
            Xi      = 0.0 + (bci==2)*1.0
            Ji      = 0.0 + (bci==3)*1.0
            tix     = ni_x*SxxNeu[nodei] + ni_y*SxyNeu[nodei]
            tiy     = ni_x*SyxNeu[nodei] + ni_y*SyyNeu[nodei]
            nitze_x = ni_x*ze[iel,1,1] + ni_y*ze[iel,2,1]
            nitze_y = ni_x*ze[iel,1,2] + ni_y*ze[iel,2,2]
            feix    = (bci!=1) * -dAi * (mesh.ke[iel]/mesh.vole[iel]*nitze_x - tix*Xi - gbar[iel,ifac,1]*Ji - 1.0/ae[iel]*be[iel,1]*taui)
            feiy    = (bci!=1) * -dAi * (mesh.ke[iel]/mesh.vole[iel]*nitze_y - tiy*Xi - gbar[iel,ifac,2]*Ji - 1.0/ae[iel]*be[iel,2]*taui)
            # up block
            Kup[iel,ifac]                            -= (bci!=1) * dAi*ni_x;
            Kup[iel,ifac+mesh.nf_el]                 -= (bci!=1) * dAi*ni_y;
            # Dirichlet nodes - uu block
            Kuu[iel,ifac,ifac]                       += (bci==1) * 1e6
            Kuu[iel,ifac+mesh.nf_el,ifac+mesh.nf_el] += (bci==1) * 1e6
            fu[iel,ifac]                             += (bci!=1) * feix + (bci==1) * VxDir[nodei] * 1e6
            fu[iel,ifac+mesh.nf_el]                  += (bci!=1) * feiy + (bci==1) * VyDir[nodei] * 1e6
            # Dirichlet nodes - pressure RHS
            fp[iel]                                  += (bci==1) * -dAi*(VxDir[nodei]*ni_x + VyDir[nodei]*ni_y)
        end
    end
    return Kuu, fu, Kup, fp
end

#--------------------------------------------------------------------#

function CreateTripletsSparse(mesh, Kuu_v, fu_v, Kup_v)
    # Create triplets and assemble sparse matrix fo Kuu
    e2fu = mesh.e2f
    e2fv = mesh.e2f .+ mesh.nf
    e2f  = hcat(e2fu, e2fv)
    idof = 1:mesh.nf_el*2
    ii   = repeat(idof, 1, length(idof))'
    ij   = repeat(idof, 1, length(idof))
    Ki   = e2f[:,ii]
    Kif  = e2f[:,ii[1,:]]
    Kj   = e2f[:,ij]
    @time Kuu  = sparse(Ki[:], Kj[:], Kuu_v[:], mesh.nf*2, mesh.nf*2)
    # file = matopen(string(@__DIR__,"/results/matrix_uu.mat"), "w" )
    # write(file, "Ki",       Ki[:] )
    # write(file, "Kj",    Kj[:] )
    # write(file, "Kuu",  Kuu_v[:] )
    # write(file, "nrow",  mesh.nf*2 )
    # write(file, "ncol",  mesh.nf*2 )
    # close(file)
    @time fu   = sparse(Kif[:], ones(size(Kif[:])), fu_v[:], mesh.nf*2, 1)
    u   = Array(fu)
    droptol!(Kuu, 1e-6)
    # Create triplets and assemble sparse matrix fo Kup
    idof = 1:mesh.nf_el*2
    ii   = repeat(idof, 1, mesh.nel)'
    ij   = repeat(1:mesh.nel, 1, length(idof))
    Ki   = e2f
    Kj   = ij
    @time Kup  = sparse(Ki[:], Kj[:], Kup_v[:], mesh.nf*2, mesh.nel  )
    # file = matopen(string(@__DIR__,"/results/matrix_up.mat"), "w" )
    # write(file, "Ki",       Ki[:] )
    # write(file, "Kj",    Kj[:] )
    # write(file, "Kup",  Kup_v[:] )
    # write(file, "nrow",  mesh.nf*2 )
    # write(file, "ncol",  mesh.nel )
    # close(file)
    return Kuu, fu, Kup
end

#--------------------------------------------------------------------#


text:
using Dates

"""
Runs a Monte Carlo simulation with the FUND model over it's distirbutional parameters, and calculates
the Social Cost of Carbon for the specified `years` and discount `rates`.
`trials`: the number of trials to run.
`ntimesteps`: how many timesteps to run.
`output_dir`: an output directory; if none provided, will create and use "output/SCC yyyy-mm-dd HH-MM-SS MCtrials".
`save_trials`: whether or not to generate and save the MC trial values up front to a file.
"""
function run_fund_scc_mcs(trials = 10000; years = [2020], rates = [0.03], ntimesteps = MimiFUND.default_nsteps + 1, output_dir = nothing, save_trials = false)

    # Set up output directories
    output_dir = output_dir == nothing ? joinpath(@__DIR__, "../../output/", "SCC $(Dates.format(now(), "yyyy-mm-dd HH-MM-SS")) MC$trials") : output_dir
    mkpath("$output_dir/results")

    # Write header in SCC output file
    scc_file = joinpath(output_dir, "scc.csv")
    open(scc_file, "w") do f
        write(f, "trial, rate, year, scc\n")
    end

    # Scenario set up
    scenario_args = Any[
        :rate           => rates
        :emissionyear   => years
    ]

    mcs = getmcs()

    # Generate trials
    if save_trials
        filename = joinpath(@__DIR__, "$output_dir/trials.csv")
        generate_trials!(mcs, trials; filename=filename)
    else
        generate_trials!(mcs, trials)
    end

    # Get FUND marginal model
    mm = create_marginal_FUND_model()
    set_models!(mcs, mm)

    # Define scenario function
    function _scenario_func(mcs::Simulation, tup::Tuple)

        # Unpack scenario tuple argument
        (rate, emissionyear) = tup

        # Get models
        base, marginal = mcs.models

        # Perturb emissons in the marginal model
        perturb_marginal_emissions!(marginal, emissionyear)

    end

    # Define post trial function
    function scc_calculation(mcs::Simulation, trialnum::Int, ntimesteps::Int, tup::Tuple)

        # Unpack scenario tuple argument
        (rate, emissionyear) = tup

        # Get marginal damages
        base, marginal = mcs.models
        marginaldamages = (marginal[:impactaggregation, :loss] - base[:impactaggregation, :loss]) / 10000000.0

        # Calculate discount factor
        T = ntimesteps == typemax(Int) ? length(Mimi.dimension(base, :time)) : ntimesteps
        discount_factor = zeros(T)
        idx = getindexfromyear(emissionyear)
        discount_factor[idx:T] = [1 / ((1 + rate) ^ t) for t in 0:T-idx]

        # Sum discounted global damages to scc
        scc = sum(sum(marginaldamages, dims = 2)[2:T] .* discount_factor[2:end])

        # Write output
        open(scc_file, "a") do f
            write(f, "$trialnum, $rate, $emissionyear, $scc\n")
        end

        return nothing
    end

    # Run monte carlo trials
    run_sim(mcs;
        models_to_run = 2,
        ntimesteps = ntimesteps,
        output_dir = "$output_dir/results",
        scenario_args = scenario_args,
        scenario_func = _scenario_func,
        post_trial_func = scc_calculation
        )

    return nothing
end


text:
using TensorNetworkCodes
using Test

@testset "surface_code" begin
    surface = surface_code(3)
    @test verify(surface)
    @test find_distance_logicals(surface)[1] == 4
end

@testset "almost_rotated_surface_code" begin
    almost_rot = almost_rotated_surface_code(3,five_qubit_code(),[2,2])
    @test verify(almost_rot)
    @test find_distance_logicals(almost_rot)[1] == 3
end

@testset "rotated_surface_code" begin
    rot_surface = rotated_surface_code(3)
    @test verify(rot_surface)
    @test find_distance_logicals(rot_surface)[1] == 3
end


text:
using Downloads
using ZipFile
url = "https://m-selig.ae.illinois.edu/ads/archives/coord_seligFmt.zip"
archive = joinpath(@__DIR__, "coord_seligFmt.zip")
reader = ZipFile.Reader(archive)
output = open(joinpath(@__DIR__, "../src/database.jl"), "w")
println(output, "#################### AUTO GENERATED! DO NOT MODIFY ####################")
println(output, "const database = AirfoilData[")
index = []
asciionly(s) =  join(map(m->getfield(m, :match), eachmatch(r"[0-9a-z ]", s)))

function formatpoints(points)
    x0, y0 = parse.(Float64, split(points[1]))
    iox = IOBuffer()
    ioy = IOBuffer()
    print(iox, '[')
    print(ioy, '[')
    if x0>2 && y0>2 # Lednicer's format
        np1 = round(Int, x0)
        np2 = round(Int, y0)
        @assert np1+np2 == length(points)-1
        for i in np1:-1:2
            x, y = split(points[i+1])
            print(iox, x, ',')
            print(ioy, y, ',')
        end
        for i in 1:np2
            x, y = split(points[i+1+np1])
            print(iox, x, ',')
            print(ioy, y, ',')
        end
    else # original format
        for l in points
            x, y = split(l)
            print(iox, x, ',')
            print(ioy, y, ',')
        end
    end
    print(iox, "],")
    write(iox, take!(ioy))
    print(iox, ']')
    String(take!(iox))
end

for f in reader.files
    dir, fname = splitdir(f.name)
    name, ext = splitext(fname)
    name = uppercase(name)
    lines = readlines(f)
    desc = strip(lines[1])
    points = filter(l->!isempty(l), strip.(@view lines[2:end]))
    @info "    $name: $desc"
    println(output, "    AirfoilData(\"", name, "\",\"", desc, "\",", formatpoints(points), "),")
    push!(index, lowercase(name) * " " * asciionly(lowercase(desc)))
end

println(output, "]")
println(output)
println(output, "const index = String[")
for l in index
    println(output, "    \"", l, "\",")
end
println(output, "]")
println(output)
println(output, "const db_size = $(length(index))")
close(output)


text:
export scatter!, scatter

## Scatter API
#   - Scatter:
#     - scatter(op, src, idx)
#     - scatter!(op, dst, src, idx)
#   - Scatter destination backpropagation
#     - ∇scatter_dst!
#   - Scatter source backpropagation
#     - ∇scatter_src
#     - ∇scatter_src!
#

function _check_dims(Ndst, Nsrc, N, Nidx)
    @assert Ndst - N == Nsrc - Nidx "Incompatible input shapes of (dst, src, idx) = ($Ndst, $Nsrc, $Nidx)."
    dims = Ndst - N
    if dims < 0
        throw(ArgumentError("dims must be non-negative but got dims=$dims."))
    end
    return dims
end

typelength(::Type{<:Number}) = 1
typelength(::Type{<:NTuple{M}}) where M = M

"""
    scatter!(op, dst, src, idx)

Scatter operation, which scatters data in `src` and assigns to `dst` according to `idx`.
With the data going to the same place, specified aggregate operation is applied on to reduce
data. For each index `k` in `idx`, accumulate values in `dst` according to

    dst[:, ..., idx[k]...] = (op).(dst[:, ..., idx[k]...], src[:, ..., k...])

# Arguments
- `op`: operations to be applied on `dst` and `src`, e.g. `+`, `-`, `*`, `/`, `max`, `min`
and `mean`.
- `dst`: the destination for `src` to aggregate to. This argument will be mutated.
- `src`: the source data for aggregating.
- `idx`: the mapping for aggregation from source (index) to destination (value).
The index of `idx` is corresponding to the index of `src` and the dimensions of `idx` must
aligned with the last few dimensions of `src`. The value of `idx` is corresponding to the
index of `dst` and the value of `idx` must indicate the last few dimensions of `dst`.
Once the dimensions match, arrays are aligned automatically. The value of `idx` can be
`Int` or `Tuple` type.
"""
function scatter!(op,
                  dst::AbstractArray{Tdst,Ndst},
                  src::AbstractArray{Tsrc,Nsrc},
                  idx::AbstractArray{Tidx,Nidx}) where {Tdst,Tsrc,Tidx<:IntOrIntTuple,Ndst,Nsrc,Nidx}
    M = typelength(Tidx)
    dims = _check_dims(Ndst, Nsrc, M, Nidx)
    scatter!(op, dst, src, idx, Val(dims))
end

function scatter!(op, dst::AbstractArray{Tdst}, src::AbstractArray{Tsrc}, idx::AbstractArray{<:IntOrIntTuple},
                  dims::Val{N}) where {Tdst,Tsrc,N}
    colons = Base.ntuple(_->Colon(), dims)
    for k in CartesianIndices(idx)
        dst_v = view(dst, colons..., idx[k]...)
        src_v = view(src, colons..., k)
        dst_v .= (op).(dst_v, src_v)
    end
    dst
end

function scatter!(op::typeof(mean),
                  dst::AbstractArray{Tdst,Ndst},
                  src::AbstractArray{Tsrc,Nsrc},
                  idx::AbstractArray{<:IntOrIntTuple,Nidx}) where {Tdst,Tsrc,Ndst,Nsrc,Nidx}
    Ns = scatter!(+, zero(dst), one.(src), idx)
    dst_ = scatter!(+, zero(dst), src, idx)
    dst .+= safe_div.(dst_, Ns)
    return dst
end

"""
    scatter(op, src, idx)

Scatter operation, which applies specified operation on `src` according to `idx`
and gives an new array `dst`.
For each index `k` in `idx`, accumulate values in `dst` according to

    dst[:, ..., idx[k]...] = (op).(src[:, ..., k...])

# Arguments
- `op`: operations to be applied on `dst` and `src`, e.g. `+`, `-`, `*`, `/`, `max` and `min`.
- `src`: the source data for aggregating.
- `idx`: the mapping for aggregation from source (index) to destination (value).
The index of `idx` is corresponding to the index of `src` and the value of `idx` is
corresponding to the index of `dst`. The value of `idx` can be `Int` or `Tuple` type.
"""
function scatter end

for op in [+, -]
    @eval function scatter(op::typeof($op),
                           src::AbstractArray{T,Nsrc},
                           idx::AbstractArray{<:IntOrIntTuple,Nidx}) where {T,Nsrc,Nidx}
        dims = Nsrc - Nidx
        dstsize = (size(src)[1:dims]..., maximum_dims(idx)...)
        dst = similar(src, T, dstsize)
        fill!(dst, Base.reduce_empty(+, T))
        scatter!(op, dst, src, idx)
    end
end

for op in [*, /]
    @eval function scatter(op::typeof($op),
                           src::AbstractArray{T,Nsrc},
                           idx::AbstractArray{<:IntOrIntTuple,Nidx}) where {T,Nsrc,Nidx}
        dims = Nsrc - Nidx
        dstsize = (size(src)[1:dims]..., maximum_dims(idx)...)
        dst = similar(src, T, dstsize)
        fill!(dst, Base.reduce_empty(*, T))
        scatter!(op, dst, src, idx)
    end
end

function scatter(op::typeof(max),
                 src::AbstractArray{T,Nsrc},
                 idx::AbstractArray{<:IntOrIntTuple,Nidx}) where {T,Nsrc,Nidx}
    dims = Nsrc - Nidx
    dstsize = (size(src)[1:dims]..., maximum_dims(idx)...)
    dst = similar(src, T, dstsize)
    fill!(dst, typemin(T))
    scatter!(op, dst, src, idx)
end

function scatter(op::typeof(min),
                 src::AbstractArray{T,Nsrc},
                 idx::AbstractArray{<:IntOrIntTuple,Nidx}) where {T,Nsrc,Nidx}
    dims = Nsrc - Nidx
    dstsize = (size(src)[1:dims]..., maximum_dims(idx)...)
    dst = similar(src, T, dstsize)
    fill!(dst, typemax(T))
    scatter!(op, dst, src, idx)
end

function scatter(op::typeof(mean),
                 src::AbstractArray{T,Nsrc},
                 idx::AbstractArray{<:IntOrIntTuple,Nidx}) where {T,Nsrc,Nidx}
    FT = float(T)
    dims = Nsrc - Nidx
    dstsize = (size(src)[1:dims]..., maximum_dims(idx)...)
    dst = similar(src, T, dstsize)
    fill!(dst, Base.reduce_empty(+, FT))
    scatter!(op, dst, src, idx)
end


text:
mutable struct RHS_IIF_Scalar{F,uType,tType,aType,P} <: Function
  f::F
  tmp::uType
  t::tType
  dt::tType
  a::aType
  p::P
end

function (f::RHS_IIF_Scalar)(resid,u)
  resid[1] = first(u) - f.tmp - (f.a*f.dt)*first(f.f.f2(first(u),f.p,f.t+f.dt,))
end

function initialize!(integrator,cache::Union{GenericIIF1ConstantCache,GenericIIF2ConstantCache})
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  A = integrator.f.f1
  cache.uhold[1] = integrator.f.f2(integrator.uprev,integrator.p,integrator.t)
  integrator.fsalfirst = integrator.f.f1(integrator.uprev,integrator.p,integrator.t) .+ cache.uhold[1]

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsalfirst
end

function perform_step!(integrator,cache::Union{GenericIIF1ConstantCache,GenericIIF2ConstantCache},repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack uhold,rhs,nl_rhs = cache
  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg

  # If adaptive, this should be computed after and cached
  A = integrator.f.f1
  if typeof(cache) <: GenericIIF1ConstantCache
    rhs.tmp = exp(A*dt)*(uprev)
  elseif typeof(cache) <: GenericIIF2ConstantCache
    @muladd rhs.tmp = exp(A*dt)*(uprev + 0.5dt*uhold[1]) # This uhold only works for non-adaptive
  end

  if integrator.success_iter > 0 && !integrator.reeval_fsal
    uhold[1] = current_extrapolant(t+dt,integrator)
  end # else uhold is previous value.

  rhs.t = t
  rhs.dt = dt
  nlres = alg.nlsolve(nl_rhs,uhold)
  uhold[1] = integrator.f.f2(nlres[1],integrator.p,t+dt)
  u = nlres[1]
  integrator.fsallast = A*u + uhold[1]
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

mutable struct RHS_IIF{F,uType,tType,aType,DiffCacheType,P} <: Function
  f::F
  tmp::uType
  t::tType
  dt::tType
  a::aType
  dual_cache::DiffCacheType
  p::P
end
function (f::RHS_IIF)(resid,u)
  _du = get_du(f.dual_cache, eltype(u))
  du = reinterpret(eltype(u),_du)
  f.f.f2(du,u,f.p,f.t+f.dt)
  @. resid = u - f.tmp - (f.a*f.dt)*du
end

function initialize!(integrator,cache::Union{GenericIIF1Cache,GenericIIF2Cache})
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  A = integrator.f.f1
  integrator.f.f2(cache.rtmp1,integrator.uprev,integrator.p,integrator.t)
  mul!(cache.k,A,integrator.uprev)
  @. integrator.fsalfirst = cache.k + cache.rtmp1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator,cache::Union{GenericIIF1Cache,GenericIIF2Cache},repeat_step=false)
  @unpack rtmp1,tmp,k = cache
  @unpack rhs,nl_rhs = cache
  @unpack t,dt,uprev,u,f,p = integrator
  alg = typeof(integrator.alg) <: CompositeAlgorithm ? integrator.alg.algs[integrator.cache.current] : integrator.alg

  @. k = uprev
  if typeof(cache) <: GenericIIF2Cache
    @muladd @. k = k + 0.5dt*rtmp1
  end

  mul!(tmp,cache.expA,k)

  if integrator.success_iter > 0 && !integrator.reeval_fsal
    current_extrapolant!(u,t+dt,integrator)
  end # else uhold is previous value.

  rhs.t = t
  rhs.dt = dt
  nlres = alg.nlsolve(nl_rhs,u)

  copyto!(u,nlres)
  integrator.f.f2(rtmp1,nlres,integrator.p,t+dt)
  A = f.f1
  integrator.fsallast .= A*u .+ rtmp1
end


text:
# uses variables from setup_tests.jl

sym_implied_pop_m, sym_implied_pop_f = eqm_consistency(mgmkt.surplus, mgmkt.msingle, mgmkt.fsingle)
asym_implied_pop_m, asym_implied_pop_f = eqm_consistency(mgmkt2.surplus, mgmkt2.msingle, mgmkt2.fsingle)

# ensure convergence to valid equilibrium
@testset "validity of equilibrium population measures" begin
	match_tol = 1e-12
	@test mgmkt.matches ≈ match_matrix(mgmkt.surplus, mgmkt.msingle, mgmkt.fsingle) atol = match_tol
	@test sym_implied_pop_m ≈ mgmkt.mdist atol = match_tol
	@test sym_implied_pop_f ≈ mgmkt.fdist atol = match_tol

	@test mgmkt2.matches ≈ match_matrix(mgmkt2.surplus, mgmkt2.msingle, mgmkt2.fsingle) atol = match_tol
	@test asym_implied_pop_m ≈ mgmkt2.mdist atol = match_tol
	@test asym_implied_pop_f ≈ mgmkt2.fdist atol = match_tol
end

# check symmetry
@test mgmkt.msingle ≈ mgmkt.fsingle

# check that supermodular surplus results in positive assortative matching
@test mgmkt2.matches[1,1,1,1] + mgmkt2.matches[end,end,end,end] > mgmkt2.matches[1,1,end,end] + mgmkt2.matches[end,end,1,1]


text:
#script for solving a 1D BVP of the form
# -a0(x)*u''(x)+a1(x)*u(x)=f(x) for x ∈ (0,1)
# Here a0(x)=1, a1(x)=0, f(x)=4*π^2*sin(2*π*x)
# with Dirichlet boundary conditions u(0) = 0, u(1)=0
include("../src/Solver1D.jl")

a0(x) = 1
a1(x) = 0
f(x) = (4*π^2)*sin(2*π*x)
exact_sol(x) = sin(2*π*x)

ptLeft = 0.
ptRight = 1.
uPtLeft = 0.
uPtRight = 1.
nrb = nrbline(ptLeft, ptRight)

solPtLeft = exact_sol(ptLeft)
solPtRight = exact_sol(ptRight)

bound_left = Boundary1D("Dirichlet", ptLeft, uPtLeft, solPtLeft)
bound_right = Boundary1D("Dirichlet", ptRight, uPtRight, solPtRight)

numElem = 10
degP = 2
new_knots = collect(LinRange(ptLeft, ptRight, numElem+1)[2:end-1])
nrb = nrbdegelev(nrb, [degP-1])
nrb = nrbkntins(nrb, [new_knots])

#refine some more knots at the beginning
ptMid = (ptLeft+ptRight)/2
newer_knots = collect(LinRange(ptLeft, ptMid, numElem+1))[2:end-1]
newer_knots = setdifftol(newer_knots, new_knots)
#@show newer_knots
nrb = nrbkntins(nrb, [newer_knots])

prob_sin = Problem1D(f, a0, a1, [bound_left, bound_right], nrb)
IGAmesh = genMesh(prob_sin.domain)
plotBasisParam(IGAmesh)

gauss_rule = genGaussLegendre(degP+1)
stiff = assemble_stiff(IGAmesh, a0, gauss_rule)
mass = assemble_mass(IGAmesh, a1, gauss_rule)
rhs = assemble_rhs(IGAmesh, f, gauss_rule)
lhs, rhs = applyBCnurbs(prob_sin.boundary_cond, stiff, mass, rhs, nrb)
sol0 = lhs\rhs
plotSol(IGAmesh, real(sol0), "Poisson1D")
plotSolError(IGAmesh, real(sol0), exact_sol, "Poisson1D")


text:
#####
##### Types for describing an Instrument Operations (convolution, resampling)
#####

"""
    type AbstractInstrumentOperator
Abstract AbstractInstrumentOperator type
"""
abstract type AbstractInstrumentOperator end

"""
    type AbstractInstrument
Abstract Instrument type
"""
abstract type AbstractInstrument end

"""
    struct FixedKernelInstrument{FT}

A struct which provides all parameters for the convolution with a kernel, which is identical across the instrument spectral axis

# Fields
$(DocStringExtensions.FIELDS)
"""
struct FixedKernelInstrument{FT} <: AbstractInstrumentOperator
    "convolution Kernel"
    kernel::OffsetArray{FT,1}
    "Output spectral grid"
    ν_out::Array{FT,1}
end;

"""
    struct VariableInstrument{FT}

        A struct which provides all parameters for the convolution with a kernel, which varies across the instrument spectral axis

# Fields
$(DocStringExtensions.FIELDS)
"""
struct VariableKernelInstrument{FT,AA} <: AbstractInstrumentOperator
    "convolution Kernel"
    kernel::OffsetArray{FT,2,AA}
    "Output spectral grid"
    ν_out::Array{FT,1}
    "Output index grid"
    ind_out::Array{Int,1}
end;

struct FTSInstrument{FT} <: AbstractInstrument
    "Maximum Optical Path Difference (cm)"
    MOPD::FT
    "Field of View (rad)"
    FOV::FT
    "Assymmetry parameter"
    β::FT
end

"""
    type AbstractNoiseModel
Abstract AbstractNoiseModel type
"""
abstract type AbstractNoiseModel end

"""
    struct GratingNoiseModel{FT}

A struct which stores important variable to compute the noise of a grating spectrometer

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef  struct GratingNoiseModel <: AbstractNoiseModel
    "Integration time `[s]`"
    t_int::Unitful.Time
    "Detector pixel size (assuming quadratic) `[μm]`"
    detector_size::Unitful.Length
    "Detector quantum efficiency"
    Qₑ
    "Efficiency of the optical bench (Telescope, grating and other optical elements)"
    η
    "F-number"
    Fnumber
    "Spectral Sampling Interval SSI or Δλ `[μm]`"
    Δλ::Unitful.Length
    "Readout noise `[e⁻]`"
    σ_read
    "Dark current `[e⁻/s]`"
    dark_current::PerTime
    # "Slit width `[μm]`"
    # slit_width::Unitful.Length
end;

function createGratingNoiseModel(ET, DS, FPA_qe, effTransmission, fnumber, SSI, RN, DC)
    # @info "Creating instrument model"
    GratingNoiseModel(
       uconvert(u"s",ET),
       uconvert(u"μm",DS),
       FPA_qe, effTransmission,
       fnumber,
       uconvert(u"μm",SSI),
       RN,
       uconvert(u"s^-1",DC))
end

function Base.show(io::IO, m::GratingNoiseModel)
    compact = get(io, :compact, false)

    if !compact
        println("Instance of GratingNoiseModel:")
        println("Integration time         = ", m.t_int)
        println("Detector size            = ", m.detector_size)
        println("FPA Quantum efficiency   = ", m.Qₑ)
        println("Optical bench efficiency = ", m.η)
        println("F-number                 = ", m.Fnumber)
        println("Spectral Sampling        = ", m.Δλ)
        println("Readout noise            = ", m.σ_read)
        println("Dark_current             = ", m.dark_current)
    else
        show(io, "Instance of GratingNoiseModel")
    end
end

"""
    type AbstractL1File
Abstract AbstractL1File type
"""
abstract type AbstractL1File end
struct L1_OCO <: AbstractL1File
    geometry::Dict{String, Any}
    ils::Dict{String, Any}
    meteo::Dict{String, Any}
    measurement::Dict{String, Any}
end;

abstract type AbstractMeasurement end

struct MeasurementOCO{FT} <: AbstractMeasurement
    "Wavelength axis"
    SpectralGrid::Array{FT,1}
    "Band ID axis"
    BandID::Tuple
    "Spectral radiances (concatenated)"
    SpectralMeasurement::Array{FT,1}
    "Radiance Unit"
    MeasurementUnit::String
    "Wavelength unit"
    GridUnit::String
    "Latitude [°]"
    latitude::FT
    "Longitude [°]"
    longitude::FT
    "Viewing Zenith Angle [°]"
    vza::FT
    "Solar Zenith Angle [°]"
    sza::FT
    "Relative Azimuth Angle [°]"
    raa::FT
    "Gravity at surface"
    g₀::FT
    "Polarization Angle [°]"
    ϕ::FT
    "Mueller Coefficients"
    mueller::Array{FT,1}
    "Pressure profile half levels"
    p_half::Array{FT,1}
    "Pressure profile half levels"
    p_full::Array{FT,1}
    "Temperature profile full level"
    T::Array{FT,1}
    "Specific humidity full level"
    q::Array{FT,1}
    "Instrument Line Shape"
    ils::Tuple
end;


text:
#=

Denisity functions for Tapestree

Ignacio Quintero Mächler

t(-_-t)

September 23 2017

=#

"""
    logdexp(x::Float64, λ::Float64)

Compute the logarithmic transformation of the
**Exponential** density with mean `λ` for `x`.
"""
logdexp(x::Float64, λ::Float64) = (log(λ) - λ * x)::Float64

"""
    logdunifU(x::Float64)

Standard uniform distribution (`[0.0,1.0]`).
"""
logdunifU(x::Float64) = 0.0

"""
    logdunif(x::Float64, a::Float64, b::Float64)

Make function to compute the logarithmic transformation of the
**Uniform** density with lower bound `a` and upper bound `b` for `x`.
"""
function logdunif(x::Float64, a::Float64, b::Float64)
  if x < a
    return -Inf
  elseif x <= b
    return -log(b-a)::Float64
  else
    return -Inf
  end
end

"""
    llrdexp_x(xp::Float64, xc::Float64, λ::Float64)

Compute the loglik ratio of the
**Exponential** density for proposal
`xp` given current `xc` both with mean `λ`.
"""
llrdexp_x(xp::Float64, xc::Float64, λ::Float64) =
  λ * (xc - xp)

"""
    logdbeta(x::Float64, α::Float64, β::Float64)

Compute the logarithmic transformation of the
**Beta** density with shape `α` and shape `β` for `x`.
"""
logdbeta(x::Float64, α::Float64, β::Float64) =
  ((α-1.0) * log(x)                 +
  (β-1.0) * log(1.0 - x)           +
  log(gamma(α + β)/(gamma(α)*gamma(β))))

"""
    llrdbeta_x(xp::Float64, xc::Float64, α::Float64, β::Float64)

Compute the logarithmic ratio for the **Beta** density
with shape `α` and shape `β` between `xp` and `xc`.
"""
function llrdbeta_x(xp::Float64, xc::Float64, α::Float64, β::Float64)
  if !(0.0 < xp < 1.0)
    return -Inf
  else
    return ((α-1.0) * log(xp/xc) +
            (β-1.0) * log((1.0 - xp)/(1.0 - xc)))
  end
end

"""
    logdnorm(x::Float64, μ::Float64, σ²::Float64)

Compute the logarithmic transformation of the
**Normal** density with mean `μ` and variance `σ²` for `x`.
"""
logdnorm(x::Float64, μ::Float64, σ²::Float64) =
  -(0.5*log(2.0π) + 0.5*log(σ²) + (x - μ)^2/(2.0σ²))

"""
    logdnorm_tc(x::Float64, μ::Float64, σ²::Float64)

Compute the logarithmic transformation of the
**Normal** density with mean `μ` and variance `σ²` for `x`, up to a constant
"""
logdnorm_tc(x::Float64, μ::Float64, σ²::Float64) =
  -0.5*log(σ²) - (x - μ)^2/(2.0σ²)::Float64

"""
    llrdnorm_ωx(x::Float64, xi::Float64, μp::Float64, μc::Float64, σ²::Float64)

Compute the log-likelihood ratio for the **Normal** density
for `ωx` updates
"""
llrdnorm_ωx(x::Float64, xi::Float64, μp::Float64, μc::Float64, σ²::Float64) =
  (-(x - xi - μp)^2 + (x - xi - μc)^2)/(2.0σ²)

"""
    llrdnorm_σ²(x::Float64, μ::Float64, σ²p::Float64, σ²c::Float64)

Compute the log-likelihood ratio for the **Normal** density
for `σ²` updates
"""
llrdnorm_σ²(x::Float64, μ::Float64, σ²p::Float64, σ²c::Float64) =
  -0.5*(log(σ²p/σ²c) + (x - μ)^2*(1.0/σ²p - 1.0/σ²c))

"""
    llrdnorm_μ(x::Float64, μp::Float64, μc::Float64, σ²::Float64)

Compute the log-likelihood ratio for the **Normal** density
for `μ` updates
"""
llrdnorm_μ(x::Float64, μp::Float64, μc::Float64, σ²::Float64) =
  ((x - μc)^2 - (x - μp)^2)/(2.0σ²)

"""
    llrdnorm_x(xp::Float64, xc::Float64, μ::Float64, σ²::Float64)

Compute the log-likelihood ratio for the **Normal** density
for `x` updates
"""
llrdnorm_x(xp::Float64, xc::Float64, μ::Float64, σ²::Float64) =
  ((xc - μ)^2 - (xp - μ)^2)/(2.0σ²)

"""
    llrdnorm_xμ(xp::Float64, xc::Float64, μp::Float64, μc::Float64, σ²::Float64)

Compute the log-likelihood ratio for the **Normal** density
for `x` and `μ` updates
"""
llrdnorm_xμ(xp::Float64, xc::Float64, μp::Float64, μc::Float64, σ²::Float64) =
  ((xc - μc)^2 - (xp - μp)^2)/(2.0σ²)

"""
    logdtnorm(x::Float64, σ²::Float64)

Compute the log-likelihood density for the **Truncated Normal** density
with `a = -1` and `b = Inf`
"""
function logdtnorm(x::Float64, σ²::Float64)
  if x < -1.0
    return -Inf
  else
    return (-x^2/(2.0*σ²) - 0.5*log(2.0π) - log(0.5*sqrt(σ²)) -
            log(1.0 - erf_custom(-1.0/(sqrt(2.0σ²)))))
  end
end

"""
    llrtnorm_x(xp::Float64, xc::Float64, μ::Float64, σ²::Float64)

Compute the log-likelihood ratio for the **Truncated Normal** density
for `x` with `a = -1`, `b = Inf`
"""
function llrdtnorm_x(xp::Float64, xc::Float64, σ²::Float64)
  if xp < -1.0
    return -Inf
  else
    return (xc^2 - xp^2)/(2.0σ²)
  end
end

"""
    erf(x::Float64)

Compute the error function
"""
function erf_custom(x::Float64)
  z = abs(x)
  t = 1.0/(1.0+0.5*z)
  r = t*exp(-z*z-1.26551223 +
        t*(1.00002368 +
            t*(0.37409196 +
              t*(0.09678418 +
                t*(-0.18628806 +
                  t*(0.27886807 +
                    t*(-1.13520398 +
                      t*(1.48851587 +
                        t*(-0.82215223 +
                          t*0.17087277)))))))))

  return x >= 0.0 ? (1.0-r) : (r -1.0)
end

"""
    logdhcau(x::Float64, scl::Float64)

Compute the logarithmic transformation of the
**Half-Cauchy** density with scale `scl` for `x`.
"""
logdhcau(x::Float64, scl::Float64) =
  log(2.0 * scl/(π *(x * x + scl * scl)))

"""
    logdhcau1(x::Float64)

Compute the logarithmic transformation of the
**Half-Cauchy** density with scale of 1 for `x`.
"""
logdhcau1(x::Float64) =
  log(2.0/(π * (x * x + 1.)))

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
using Base: get_chunks_id
using PSSR
using Test

@testset "PSSR.jl" begin
    # Write your tests here.
    get_protein_data("C:/Users/honoa/Machine Learning/Protein")
end


text:
module PSID
using XLSX, DataDeps, DataFrames, CSV, LightXML, AbstractTrees, JSONTables, JSON3
using DataFramesMeta, SHA, DataStructures, Parsers

include("types.jl")
include("init.jl")
include("process_codebook.jl")
include("use_codebook.jl")
include("unzip_data.jl")
include("construct_alldata.jl")
export makePSID

end # module


text:
using ColorSchemes

# solve model for given rates and force interval
function ftsolve(M,rates,Fs,ts)
    sols = []
    for f in Fs
        # M.cons = f
        M.f = f
        sol = solve_prob(M,rates;tᵢ=ts)
        push!(sols,sol)
    end
    return sols
end

# take model solution and construct P and PDF surfaces for plotting
function construct_surf(sols,Fs,ts)
    p = []
    pdf = []
    for (i,f) in enumerate(Fs)

        pᵢ = sum(Array(sols[i](ts))',dims=2)
        pdfᵢ = p2pdf(pᵢ,sols[i])

        p = push!(p,pᵢ)
        pdf = push!(pdf,pdfᵢ)
    end
    p = reduce(hcat,p)
    pdf = reduce(hcat,pdf)
    return p, pdf
end

# compute pdf from p distribution
function p2pdf(p,sol)
    mean_lifetime(sol)
    pdf = p ./ mean_lifetime(sol)
end

function plot_ft_fit(M,rates,data;dt=0.1,dF=1.0)

    ts = 0:dt:maximum(data.t)
    Fs = 0:dF:maximum(data.F)
    # Fs = 0:dF:20
    sols = ftsolve(M,rates,Fs,ts)

    p, pdf = construct_surf(sols,Fs,ts)

    # plot contours

    p1 = plot(Fs,ts,log.(p),st=:contour,
                            fill=(true),#cgrad(:plasma)),
                            camera=(0,90),
                            clims=(-7,0),
                            ylabel="ln(p)",xlabel="F")
    scatter!(data.F,data.t,c=:white,markeralpha=0.7)
    display(p1)

    pex = surface(Fs,ts,(p),#,st=:contour,
                            fill=(true),#cgrad(:plasma)),
                            # camera=(0,90),
                            # clims=(-7,0),
                            zlabel="p",xlabel="F")
    display(pex)

    p2 = plot(Fs,ts,p,st=:contour,
                        fill=(true),#,cgrad(:plasma)),
                        camera=(0,90),
                        ylabel="p",xlabel="F")
    display(p2)

    p3 = plot(Fs,ts,log.(pdf),st=:contour,
                        fill=(true,cgrad(:plasma)),
                        camera=(0,90),
                        ylabel="pdf",xlabel="F")

    # plot lnp curves

    p4=plot()
    lnp = p2lnp.(p)
    for i in 1:1:size(lnp,2)
            plot!(ts,lnp[:,i],c = i,
                                xlims=(0,5),
                                ylims=(-7,0),
                                legend=false,
                                ylabel="ln(p)",xlabel="t")
    end

    # plot mean lifetime distribution
    t_means = mean_lifetime.(sols)
    p5 = plot(Fs,t_means,ylabel="<t>",xlabel="F")

    display(plot(p1,p2,p3,p4,p5))
    # display(p2)
    # display(p3)
    # display(p4)
    # display(p5)
    return sols
end

function plot_mean_ft_fit(M,k,data,binned;species=[],dt=0.1,dF=1.0)

    ts = 0:dt:maximum(data.t)
    Fs = 0:dF:maximum(data.F)

    M.solver = Tsit5  #change solver to Tsit5 (can handle ± std)
    sols = ftsolve(M,k,Fs,ts)
    mean_t = mean_lifetime.(sols)
    h = plot!(Fs,mean_lifetime.(sols),c=Color,
                    # ribbon=0.1,
                    markerstrokecolor=:auto,
                    linealpha = 0.3,
                    ylabel="<t>",
                    xlabel="F",
                    label=species)
    plot_F_vs_t(clean_dissociation_data(binned))
end

function extrapolate_fit(M,k,data,F_max;species=[],dt=0.1,dF=1.0)

    ts = 0:dt:maximum(data.t)
    Fs = 0:dF:maximum(F_max)

    # M.solver = Tsit5  #change solver to Tsit5 (can handle ± std)
    sols = ftsolve(M,k,Fs,ts)
    mean_t = mean_lifetime.(sols)
    h = plot!(Fs,log10.(mean_lifetime.(sols)),c=Color,
                    # ribbon=0.1,
                    markerstrokecolor=:auto,
                    linealpha = 1,
                    ylabel="<t>",
                    xlabel="F",
                    label=species)
end

function mean_ft_fit(M,k,data,F_max;species=[],dt=0.1,dF=1.0)

    ts = 0:dt:maximum(data.t)
    Fs = 0:dF:F_max

    M.solver = Tsit5  #change solver to Tsit5 (can handle ± std)
    sols = ftsolve(M,k,Fs,ts)
    mean_t = mean_lifetime.(sols)
end


text:
@testset "Hybrid DCPLossLess with BasicHybridDispatch formulation" begin
    device_model = DeviceModel(HybridSystem, BasicHybridDispatch)
    sys = PSB.build_system(PSITestSystems, "c_sys5_hybrid")

    # Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 672, 0, 384, 336, 192, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hybrid ACPPowerModel with BasicHybridDispatch formulation" begin
    device_model = DeviceModel(HybridSystem, BasicHybridDispatch)
    sys = PSB.build_system(PSITestSystems, "c_sys5_hybrid")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, sys)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 864, 0, 672, 624, 288, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hybrid DCPLossLess BasicHybridDispatch" begin
    device_model = DeviceModel(HybridSystem, StandardHybridDispatch)
    sys = PSB.build_system(PSITestSystems, "c_sys5_hybrid")

    # Parameters Testing
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 672, 0, 384, 336, 192, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Hybrid ACPPowerModel with StandardHybridDispatch formulation" begin
    device_model = DeviceModel(HybridSystem, StandardHybridDispatch)
    sys = PSB.build_system(PSITestSystems, "c_sys5_hybrid")

    # No Parameters Testing
    model = DecisionModel(MockOperationProblem, ACPPowerModel, sys)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 864, 0, 672, 624, 288, true)
    psi_checkobjfun_test(model, GAEVF)
end