text:
## --- Import the packages we'll be needing

    using StatGeochem, Statistics
    include("Utilities.jl")
    using Plots; gr(); default(fmt=:svg);
    using ProgressMeter: @showprogress

## --- Define physical parameters and set up model

    # Time step (Myr)
    dt = 0.1
    timevec = collect(0:dt:4000)

    # Calculate past mantle heat production
    p = Parameters()
    H = bse_heat_production(p,mcdbse,timevec) * p.Mm
    Qrm = H .* (p.Qrm_now / H[1])

## --- Run Newton's method

    nNewton = 20 # Number of Newton's method iterations
    dQrm_now = p.Qrm_now*1E-9 # Heat production step size
    p.Rc = 650E3 # Plate bending radius (km)

    Qrm_now_vec = Array{Float64}(undef, nNewton)
    Qrm_now_vec[1] = p.Qrm_now # Initial guess for modern heat production
    for i = 1:(nNewton-1)
        ll = coolingmodel_LL(p, Qrm * (Qrm_now_vec[i]/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma)
        llp = coolingmodel_LL(p, Qrm * ((Qrm_now_vec[i]+dQrm_now)/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma)
        dll_dQrm = (llp-ll) / dQrm_now
        Qrm_now_vec[i+1] = Qrm_now_vec[i] - ll/dll_dQrm/2
    end

    (Tm, Tc, Qm, Qc, dp) = coolingmodel(p, Qrm * (Qrm_now_vec[nNewton]/Qrm[1]), dt, timevec)
    h = plot(TrelObs_time, TrelObs, yerror=2*TrelObs_sigma, seriestype=:scatter, markersize=2, color=:darkblue, markerstrokecolor=:auto, label="Data")
    plot!(h, timevec/1000, Tm .- Tm[1], label="Model", xlabel="Age (Ma)", ylabel="\\Delta T (C)", framestyle=:box, legend=:topleft)
    xlims!(h,(0,4))
    display(h)
    print("Qrm convergence:\n")
    display(Qrm_now_vec)

## --- Explore different Ur at high d (8.0)

    # Get fresh parameters and set difficulty
    p = Parameters()
    p.Rc = (p.eta_L / (10.0^8))^(1/3)

    nSteps = 10
    Qrm_now_Min = 0.0
    Qrm_now_Max = 1.6E13
    Qrm_now_vec = reverse(range(Qrm_now_Min, Qrm_now_Max, length=nSteps))

    cmap = cgrad((plasma[40:212]))
    h = plot(TrelObs_time,TrelObs,yerror=2*TrelObs_sigma,seriestype=:scatter, markersize=2, color=:darkblue, markerstrokecolor=:auto, label="Data")
    plottime = downsample(timevec/1000,50)
    for i=1:nSteps
        (Tm, Tc, Qm, Qc, dp) = coolingmodel(p, Qrm * (Qrm_now_vec[i]/Qrm[1]), dt, timevec)
        plotTrel = downsample(Tm .- Tm[1], 50)
        if i==1
            plot!(h,plottime,plotTrel,line_z=Qrm_now_vec[i]*1E-12,line=(cmap),label="Model") #color=cmap[i]
        else
            plot!(h,plottime,plotTrel,line_z=Qrm_now_vec[i]*1E-12,line=(cmap),label="") #color=cmap[i]
        end
    end
    plot!(h,legend=:topleft,xlabel="Age (Ma)",ylabel="\\Delta T (C)",framestyle=:box,grid=:off,fg_color_legend=:white)
    ylims!(h,(-20,250))
    xlims!(h,(0,4))
    display(h)

    savefig(h,"d=8variableUr.pdf")

## --- Explore different Ur at low d (5.6)

    # Get fresh parameters and set difficulty
    p = Parameters()

    nSteps = 12
    Qrm_now_Min = 2.635E13
    Qrm_now_Max = 2.639E13
    p.Rc = (p.eta_L / (10.0^5.6))^(1/3)

    # Qrm_now_Min = 2.644E13
    # Qrm_now_Max = 2.6455E13
    # p.Rc = (p.eta_L / (10.0^5.5))^(1/3)

    Qrm_now_vec = reverse(range(Qrm_now_Min, Qrm_now_Max, length=nSteps))

    cmap = cgrad((plasma[40:212]))
    h = plot(TrelObs_time,TrelObs,yerror=2*TrelObs_sigma,seriestype=:scatter, markersize=2, color=:darkblue, markerstrokecolor=:auto, label="Data")
    plottime = downsample(timevec/1000,50)
    for i=1:nSteps
        (Tm, Tc, Qm, Qc, dp) = coolingmodel(p, Qrm * (Qrm_now_vec[i]/Qrm[1]), dt, timevec)
        plotTrel = downsample(Tm .- Tm[1], 50)
        if i==1
            plot!(h,plottime,plotTrel,line_z=Qrm_now_vec[i]*1E-12,line=(cmap),label="Model") #color=cmap[i]
        else
            plot!(h,plottime,plotTrel,line_z=Qrm_now_vec[i]*1E-12,line=(cmap),label="") #color=cmap[i]
        end
    end
    plot!(h,legend=:topleft,xlabel="Age (Ma)",ylabel="\\Delta T (C)",framestyle=:box,grid=:off,fg_color_legend=:white)
    ylims!(h,(-20,250))
    xlims!(h,(0,4))
    display(h)

    savefig(h,"d=5.6variableUr.pdf")

## --- Explore nl/Rc3 space - small range

    nNewton = 20
    n_difficulties = 100
    minDiff = 4.5
    maxDiff = 6.5
    difficulies = range(minDiff,maxDiff,length=n_difficulties)
    p = Parameters() # Default parameters
    descrip = "burnt-in"

    # p.Qc_now = 12E12 # Turn up present core heat flux from 6 TW
    # p.Tc_now = 6400+273.15 # Turn up present core temperature from 4400 K
    # p.Qm_now = 25E12 # Turn down present-day mantle heat flux from 35 TW.
    # descrip = "Qc12Qm25"

    # p.Qc_now = 20E12 # Turn up present core heat flux from 6 TW
    # p.Tc_now = 6400+273.15 # Turn up present core temperature from 4400 K
    # descrip = "Qc20"

    # p.Qm_now = 20E12 # Turn down present-day mantle heat flux from 35 TW.
    # descrip = "Qm20"

    Rc_t = (p.eta_L ./ (10.0.^difficulies)).^(1/3)
    Qrm_now_best = Array{Float64}(undef, size(Rc_t))
    residual_best = Array{Float64}(undef, size(Rc_t))
    residual_offset_mille = Array{Float64}(undef, size(Rc_t))
    residual_offset_cent = Array{Float64}(undef, size(Rc_t))

    cmap = cgrad(reverse(viridis[1:212]))
    h1 = plot(TrelObs_time,TrelObs,yerror=2*TrelObs_sigma,seriestype=:scatter, color=:darkblue, markersize=2, markerstrokecolor=:auto, label="Data")
    h2 = plot(xlabel="Age (Ma)",ylabel="Core temperature (K)",framestyle=:box)
    h3 = plot(xlabel="Age (Ma)",ylabel="Core heat flux (W)",framestyle=:box)
    h4 = plot(xlabel="Age (Ma)",ylabel="Mantle heat flux (W)",framestyle=:box)
    plottime = downsample(timevec/1000,50)
    @time for n=1:n_difficulties
        # Find best-fit Qrm_now for this Rc_t[n]
        p.Rc = Rc_t[n]
        Qrm_now_best[n] = coolingmodel_Newton_Qrm(p,nNewton,Qrm,dt,timevec,TrelObs_time*1000,TrelObs,TrelObs_sigma)

        # Keep track of results
        residual_best[n] = coolingmodel_resid(p, Qrm*(Qrm_now_best[n]/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma)
        residual_offset_mille[n] = coolingmodel_resid(p, Qrm*(Qrm_now_best[n]*1.001/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma) / 2
        residual_offset_mille[n] += coolingmodel_resid(p, Qrm*(Qrm_now_best[n]*0.999/Qrm[1]),dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma) / 2
        residual_offset_cent[n] = coolingmodel_resid(p, Qrm*(Qrm_now_best[n]*1.01/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma) / 2
        residual_offset_cent[n] += coolingmodel_resid(p, Qrm*(Qrm_now_best[n]*0.99/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma) / 2
        (Tm, Tc, Qm, Qc, dp) = coolingmodel(p, Qrm*(Qrm_now_best[n]/Qrm[1]), dt, timevec)

        # Plot results
        plotTrel = downsample(Tm .- Tm[1], 50)
        if n == n_difficulties
            plot!(h1,plottime,plotTrel,line_z=difficulies[n],line=(cmap),label="Model")
        else
            plot!(h1,plottime,plotTrel,line_z=difficulies[n],line=(cmap),label="")
        end
        plot!(h2,plottime,downsample(Tc,50),line_z=difficulies[n],line=(cmap),label="")
        plot!(h3,plottime,downsample(Qc,50),line_z=difficulies[n],line=(cmap),label="")
        plot!(h4,plottime,downsample(Qm,50),line_z=difficulies[n],line=(cmap),label="")
    end
    plot!(h1,legend=:topleft,xlabel="Age (Ma)",ylabel="\\Delta T (C)",framestyle=:box)
    xlims!(h1,(0,4))
    ylims!(h1,(-20,250))
    display(h1)
    savefig(h1,"$(descrip)_CoolingCurveVsDifficulty.pdf")
    display(h2)
    savefig(h2,"$(descrip)_CoolingCurveVsDifficulty_coretemp.pdf")
    display(h3)
    savefig(h3,"$(descrip)_CoolingCurveVsDifficulty_coreflux.pdf")
    display(h4)
    savefig(h4,"$(descrip)_CoolingCurveVsDifficulty_mantleflux.pdf")

    # Two-axis plot for best-fit d and Qrm_now
    h = plot(xlabel="log_{10}(\\eta_L/Rc^3)",ylabel="Sum squared residual",xlims=(minDiff,maxDiff),legend=:topright,fg_color_legend=:white)
    plot!(log10.(p.eta_L./Rc_t.^3),residual_best,label="Residual at best-fit Qrm",grid=:off)
    plot!(log10.(p.eta_L./Rc_t.^3),residual_offset_mille,label="with Qrm offset 1 per mille",ylims=(0,1.5*maximum(residual_best)),color=:blue)
    plot!(log10.(p.eta_L./Rc_t.^3),residual_offset_cent,label="with Qrm offset 1 percent",ylims=(0,1.5*maximum(residual_best)),color=:darkblue)
    plot!(twinx(),log10.(p.eta_L./Rc_t.^3),Qrm_now_best*1E-12,label="",grid=:off,color=:darkred,ylabel="Best-fit present Qrm (TW)",xlims=(minDiff,maxDiff))
    display(h)
    savefig(h,"$(descrip)_Qrm_now=f(Difficulty).pdf")

## --- Explore nl/Rc3 space -- Full range

    nNewton = 20
    n_difficulties = 100
    minDiff = 4
    maxDiff = 8
    difficulies = range(minDiff,maxDiff,length=n_difficulties)
    p = Parameters()
    descrip = "default"

    # p.Qc_now = 20E12 # Turn up present core heat flux from 6 TW
    # p.Tc_now = 6400+273.15 # Turn up present core temperature from 4400 K
    # descrip = "Qc20"

    p.Qm_now = 13E12 # Turn down present-day mantle heat flux from 35 TW.
    descrip = "Qm13"

    Rc_t = (p.eta_L ./ (10 .^ difficulies)).^(1/3)
    Qrm_now_best = Array{Float64}(undef, size(Rc_t))
    residual_best = Array{Float64}(undef, size(Rc_t))
    residual_offset_mille = Array{Float64}(undef, size(Rc_t))
    residual_offset_cent = Array{Float64}(undef, size(Rc_t))

    cmap = cgrad(reverse(viridis[1:212]))
    h = plot(TrelObs_time,TrelObs,yerror=2*TrelObs_sigma,seriestype=:scatter, color=:darkblue, markersize=2, markerstrokecolor=:auto, label="Data")
    plottime = downsample(timevec/1000,50)
    @time for n=1:n_difficulties
        # Find best-fit Qrm_now for this Rc_t[n]
        p.Rc = Rc_t[n]
        Qrm_now_best[n] = coolingmodel_Newton_Qrm(p,nNewton,Qrm,dt,timevec,TrelObs_time*1000,TrelObs,TrelObs_sigma)

        # Keep track of results
        residual_best[n] = coolingmodel_resid(p, Qrm*(Qrm_now_best[n]/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma)
        residual_offset_mille[n] = coolingmodel_resid(p, Qrm*(Qrm_now_best[n]*1.001/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma) / 2
        residual_offset_mille[n] += coolingmodel_resid(p, Qrm*(Qrm_now_best[n]*0.999/Qrm[1]),dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma) / 2
        residual_offset_cent[n] = coolingmodel_resid(p, Qrm*(Qrm_now_best[n]*1.01/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma) / 2
        residual_offset_cent[n] += coolingmodel_resid(p, Qrm*(Qrm_now_best[n]*0.99/Qrm[1]), dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma) / 2
        (Tm, Tc, Qm, Qc, dp) = coolingmodel(p, Qrm*(Qrm_now_best[n]/Qrm[1]), dt, timevec)

        # Plot results
        plotTrel = downsample(Tm .- Tm[1], 50)
        if n == n_difficulties
            plot!(h,plottime,plotTrel,line_z=difficulies[n],line=(cmap),label="Model")
        elseif mod(n,2)==0
            plot!(h,plottime,plotTrel,line_z=difficulies[n],line=(cmap),label="")
        end
    end
    plot!(h,legend=:topleft,xlabel="Age (Ma)",ylabel="\\Delta T (C)",framestyle=:box)
    xlims!(h,(0,4))
    ylims!(h,(-20,250))
    display(h)
    savefig(h,"$(descrip)_CoolingCurveVsDifficultyHighRes.pdf")

    # Two-axis plot for best-fit d and Qrm_now
    h = plot(xlabel="log_{10}(\\eta_L/Rc^3)",ylabel="Sum squared residual",xlims=(minDiff,maxDiff),legend=:topleft,fg_color_legend=:white)
    plot!(log10.(p.eta_L./Rc_t.^3),residual_best,label="Residual at best-fit Qrm",grid=:off)
    plot!(log10.(p.eta_L./Rc_t.^3),residual_offset_mille,label="with Qrm offset 1 per mille",ylims=(0,1.5*maximum(residual_best)),color=:blue)
    plot!(log10.(p.eta_L./Rc_t.^3),residual_offset_cent,label="with Qrm offset 1 percent",ylims=(0,1.5*maximum(residual_best)),color=:darkblue)
    plot!(twinx(),log10.(p.eta_L./Rc_t.^3),Qrm_now_best*1E-12,label="",grid=:off,color=:darkred,ylabel="Best-fit Qrm today",xlims=(minDiff,maxDiff))
    display(h)
    savefig(h,"$(descrip)_Ur=f(Difficulty)HighRes.pdf")

## --- MCMC-Newton hybrid inversion

    nSteps = 10^3
    burnin = 5*10^2
    d = 6 # Initial guess for difficulty
    d_step_sigma = 0.5 # Standard deviation of Gaussian proposal distribution by which we adjust difficulty

    p = Parameters() # Start with fresh parameters
    # p.Qm_now = 15E12 # Turn down present-day mantle heat flux from 35 TW.

    # Define a function to conduct a hybrid MCMC-Newton inversion,
    # adjusting d randomly, and finding the best-fit Qrm_now and Tm path
    function MCMC_Newton_SMC(nSteps, d, d_step_sigma, p)

        # Allocate arrays to record the stationary distributions
        d_dist = Array{Float64}(undef, nSteps)
        Qrm_now_dist = Array{Float64}(undef, nSteps)
        Trel_dist = Array{Float64}(undef, nSteps,length(TrelObs_time))
        ll_dist = Array{Float64}(undef, nSteps)
        acceptancedist = fill(false,nSteps)

        # First markov chain step
        p.Rc = (p.eta_L/(10.0^d))^(1/3)
        Qrm_now = coolingmodel_Newton_Qrm(p, nNewton, Qrm, dt, timevec, TrelObs_time*1000, TrelObs, TrelObs_sigma)
        Tm = coolingmodel(p, Qrm*(Qrm_now/Qrm[1]), dt, timevec)[1]
        Trel = linterp1(timevec,Tm .- Tm[1], TrelObs_time*1000)
        ll = sum( -(Trel - TrelObs).^2 ./ (2*TrelObs_sigma.^2) - log.(sqrt.(2*pi*TrelObs_sigma)))

        d_dist[1] = d
        Qrm_now_dist[1] = Qrm_now
        ll_dist[1] = ll

        @showprogress for n=1:nSteps
           d_prop = d + randn()*d_step_sigma
           p.Rc = (p.eta_L / (10.0^d_prop))^(1/3)
           Qrm_now_prop = coolingmodel_Newton_Qrm(p,nNewton,Qrm,dt,timevec,TrelObs_time*1000,TrelObs,TrelObs_sigma)
           Tm_prop = coolingmodel(p, Qrm*(Qrm_now_prop/Qrm[1]), dt, timevec)[1]
           Trel_prop = linterp1(timevec, Tm_prop .- Tm_prop[1], TrelObs_time*1000)
           ll_prop = sum( -(Trel_prop - TrelObs).^2 ./ (2*TrelObs_sigma.^2) - log.(sqrt.(2*pi*TrelObs_sigma)))

           # Accept or reject proposal based on likelihood
           if rand() < exp(ll_prop-ll)
               d_step_sigma = 2.9*abs(d_prop-d)
               Trel = Trel_prop
               d = d_prop
               Qrm_now = Qrm_now_prop
               ll = ll_prop
               acceptancedist[n] = true
           end

           # Record results
           Trel_dist[n,:]=Trel
           d_dist[n] = d
           Qrm_now_dist[n] = Qrm_now
           ll_dist[n] = ll
        end
        return (d_dist, Qrm_now_dist, Trel_dist, ll_dist, acceptancedist)
    end

    # Run the hybrid inversion
    (d_dist, Qrm_now_dist, Trel_dist, ll_dist, acceptancedist) = MCMC_Newton_SMC(nSteps, d, d_step_sigma, p)

    Qrm_now_mu = nanmean(Qrm_now_dist[burnin:end])
    Qrm_now_sigma = nanstd(Qrm_now_dist[burnin:end])

    d_mu = nanmean(d_dist[burnin:end])
    d_sigma = nanstd(d_dist[burnin:end])

    p.Rc = (p.eta_L/(10.0^d_mu))^(1/3)
    (Tm_mu, Tc_mu, Qm_mu, Qc_mu, dp_mu) = coolingmodel(p,Qrm*(Qrm_now_mu/Qrm[1]),dt,timevec)

    h = plot(TrelObs_time, TrelObs, yerror=2*TrelObs_sigma, seriestype=:scatter, color=:darkblue, markersize=2, markerstrokecolor=:auto, label="Data")
    plot!(h,timevec/1000, Tm_mu .- Tm_mu[1], label="Model: Qrm_now=$(round(Qrm_now_mu,sigdigits=4)),d=$(round(d_mu,sigdigits=4))", legend=:topleft)
    plot!(h,xlabel="Age (Ga)",ylabel="\\Delta T (C)")
    display(h)
    savefig(h,"MCMC-Newton Best Fit.pdf")

    print(" Qrm_now: $Qrm_now_mu +/- $Qrm_now_sigma\n d: $d_mu +/- $d_sigma\n\n")

## --- Full unconstrained MCMC inversion, backwards-time mode

    # While the results of forward and reverse-mode numerical models are fundamentally
    # identical due to time-reversal symmetry, they lead to different ways of
    # exploring the parameter space

    nSteps = 2*10^6
    burnin = 2*10^5
    jumpingsigmafactor = 2.718  # This can be adjusted to optimize acceptance likelihood

    # Start with fresh parameters
    p = Parameters()

    # Initial guesses
    # d = 4
    d = 8
    p.Rc = (p.eta_L/(10.0^d))^(1/3)

    p.Qm_initial = 43.5E12 # Initial mantle heat flux (W)
    p.Qc_initial = 25.8E12 # Initial core heat flux (W)
    p.Tm_initial = 1800.0 # Initial mantle temperature (K)
    p.Tc_initial = 8000.0 # Initial core temperature (K).
    p.Qm_initial_sigma = 20E12
    p.Qc_initial_sigma = 10E12
    p.Tm_initial_sigma = 200
    p.Tc_initial_sigma = 2000

    p.Qm_now = 35E12 # Present mantle heat flux (W)
    p.Qc_now = 10E12 # Present core heat flux (W)
    p.Tm_now = 1300+273.15 # Present mantle temperature (K)
    p.Tc_now = 6200+273.15 # Present core temperature (K). I/O boundary 6230 +/- 500 according to 10.1126/science.1233514
    p.Qm_now_sigma = 10E12
    p.Qc_now_sigma = 10E12
    p.Tm_now_sigma = 100.0
    p.Tc_now_sigma = 1000.0

    descrip = "Qm35-10_Qc10-10"

    # Define a function to conduct a full MCMC inversion
    function MCMC_SMC_bwd(nSteps, d, jumpingsigmafactor, p, Qrm)

        # Allocate arrays to record the stationary distributions
        acceptancedist = fill(false,nSteps)
        ll_dist = Array{Float64}(undef, nSteps)
        d_dist = Array{Float64}(undef, nSteps)
        Qrm_now_dist = Array{Float64}(undef, nSteps)
        Qm_initial_dist = Array{Float64}(undef, nSteps)
        Qc_initial_dist = Array{Float64}(undef, nSteps)
        Tc_initial_dist = Array{Float64}(undef, nSteps)
        Tm_initial_dist = Array{Float64}(undef, nSteps)
        Qm_now_dist = Array{Float64}(undef, nSteps)
        Qc_now_dist = Array{Float64}(undef, nSteps)
        Tc_now_dist = Array{Float64}(undef, nSteps)
        Tm_now_dist = Array{Float64}(undef, nSteps)
        Q_addtl_dist = Array{Float64}(undef, nSteps)
        Trel_dist = Array{Float64}(undef, nSteps,length(TrelObs_time))

        # Proposal distributions
        d_step_sigma = 0.1
        Qrm_now_step_sigma = 0.1E12
        Qm_step_sigma = 0.1E12
        Qc_step_sigma = 0.1E12
        Tc_step_sigma = 100
        Tm_step_sigma = 10
        Q_addtl_step_sigma = 0.1E12

        # First markov chain step / first proposal
        p_prop = deepcopy(p)
        d_prop = log10(p_prop.eta_L/p_prop.Rc^3)
        Qrm_now_prop = coolingmodel_Newton_Qrm(p_prop,20,Qrm,dt,timevec,TrelObs_time*1000,TrelObs,TrelObs_sigma)
        (Tm_prop, Tc_prop, Qm_prop, Qc_prop, dp_prop) = coolingmodel(p_prop, Qrm*(Qrm_now_prop/Qrm[1]), dt, timevec)
        Trel_prop = linterp1s(timevec, Tm_prop .- Tm_prop[1], TrelObs_time*1000)

        h = plot(TrelObs_time, TrelObs, yerror=2*TrelObs_sigma, seriestype=:scatter, color=:darkblue, markersize=2, markerstrokecolor=:auto, label="Data")
        plot!(h, timevec/1000, Tm_prop .- Tm_prop[1], label="Model: Qrm_now=$(round(Qrm_now_prop,sigdigits=3)),d=$(round(d,sigdigits=3))")
        plot!(h, xlabel="Age (Ga)", ylabel="\\Delta T (C)", grid=:off, legend=:topleft, fg_color_legend=:white)
        summarystring = " d:\t\t$(round(d,sigdigits=2))\n" *
            " Qrm_now: $(round(Qrm_now_prop,sigdigits=3))\n\n" *
            " Qc: $(round(p.Qc_now/1E12,sigdigits=3))\n" *
            " Qm: $(round(p.Qm_now/1E12,sigdigits=3))\n" *
            " Tc: $(round(Int,p.Tc_now))\n" *
            " Tm: $(round(Int,p.Tm_now))\n"
        plot!(h,1,1,seriestype=:scatter, ann=(0.16,130,text(summarystring,8,"Helvetica",:left)),label="")
        display(h)
        savefig(h,"MCMC_firstproposal.pdf")

        # Log likelihood of first proposal
        ll_prop = LL(Trel_prop,TrelObs,TrelObs_sigma*2) + LL(d_prop,5,5) +
           LL(Qrm_now_prop, p.Qrm_now, p.Qrm_now_sigma) +
           LL(p_prop.Qm_initial, p.Qm_initial, p.Qm_initial_sigma) +
           LL(p_prop.Qc_initial, p.Qc_initial, p.Qc_initial_sigma) +
           LL(p_prop.Tc_initial, p.Tc_initial, p.Tc_initial_sigma) +
           LL(p_prop.Tm_initial, p.Tm_initial, p.Tm_initial_sigma) +
           LL(Qm_prop[1],p.Qm_now,p.Qm_now_sigma) +
           LL(Qc_prop[1],p.Qc_now,p.Qc_now_sigma) +
           LL(Tc_prop[1],p.Tc_now,p.Tc_now_sigma) +
           LL(Tm_prop[1],p.Tm_now,p.Tm_now_sigma)

        # Accept and record first proposal
        ll = ll_prop
        d = d_prop
        Qrm_now = Qrm_now_prop
        Qm_now = p_prop.Qm_now
        Qc_now = p_prop.Qc_now
        Tc_now = p_prop.Tc_now
        Tm_now = p_prop.Tm_now
        Q_addtl = p_prop.Q_addtl
        Qm_initial = Qm_prop[end]
        Qc_initial = Qc_prop[end]
        Tc_initial = Tc_prop[end]
        Tm_initial = Tm_prop[end]
        Trel = copy(Trel_prop)

        ll_dist[1] = ll
        d_dist[1] = d
        Qrm_now_dist[1] = Qrm_now_prop
        Qm_initial_dist[1] = Qm_initial
        Qc_initial_dist[1] = Qc_initial
        Tc_initial_dist[1] = Tc_initial
        Tm_initial_dist[1] = Tm_initial
        Qm_now_dist[1] = Qm_now
        Qc_now_dist[1] = Qc_now
        Tc_now_dist[1] = Tc_now
        Tm_now_dist[1] = Tm_now
        Q_addtl_dist[1] = Q_addtl
        Trel_dist[1,:] = Trel

        # Run the Markov chain / Metropolis walker
        @showprogress for n=1:nSteps

            # Randomly choose which parameter to adjust. The main benefit of his
            # form of "Metropolis-in-Gibbs" / "Alternate conditional sampling"
            # is that it allows us to optimize the proposal (jumping) distributions
            # for each parameter individually
            case = ceil(rand()*6)

            # Propose new parameters
            Qrm_now_prop = case==1 ? abs(Qrm_now + randn()*Qrm_now_step_sigma) : Qrm_now
            d_prop  = case==2 ? abs(d + randn()*d_step_sigma - 1.0) + 1.0 : d
            p_prop.Rc = abs(p_prop.eta_L / (10.0^d_prop))^(1/3)
            p_prop.Qm_now = case==3 ? abs(Qm_now + randn()*Qm_step_sigma - 1E12) + 1E12 : Qm_now
            p_prop.Qc_now = case==4 ? abs(Qc_now + randn()*Qc_step_sigma - 1E12) + 1E12 : Qc_now
            p_prop.Tc_now = case==5 ? abs(Tc_now + randn()*Tc_step_sigma - 4000) + 4000 : Tc_now
            p_prop.Tm_now = case==6 ? abs(Tm_now + randn()*Tm_step_sigma - 1000) + 1000 : Tm_now
            p_prop.Q_addtl = case==7 ? abs(Q_addtl + randn()*Q_addtl_step_sigma) : Q_addtl

            # If proposal is so bad that we hit a math error, let ll be NaN to
            # auto-reject the proposal
            Trel_prop = NaN
            ll_prop = NaN
            try
                # Calculate proposed solution
                (Tm_prop, Tc_prop, Qm_prop, Qc_prop, dp_prop)  = coolingmodel(p_prop, Qrm*(Qrm_now_prop/Qrm[1]), dt, timevec)
                Trel_prop = linterp1s(timevec, Tm_prop .- Tm_prop[1], TrelObs_time*1000)

                # Calculate log likelihood of new proposal
                ll_prop = LL(Trel_prop,TrelObs,TrelObs_sigma*2) + LL(d_prop,5,5) +
                   LL(Qrm_now_prop, p.Qrm_now, p.Qrm_now_sigma) +
                   LL(p_prop.Qm_initial, p.Qm_initial, p.Qm_initial_sigma) +
                   LL(p_prop.Qc_initial, p.Qc_initial, p.Qc_initial_sigma) +
                   LL(p_prop.Tc_initial, p.Tc_initial, p.Tc_initial_sigma) +
                   LL(p_prop.Tm_initial, p.Tm_initial, p.Tm_initial_sigma) +
                   LL(Qm_prop[1], p.Qm_now, p.Qm_now_sigma) +
                   LL(Qc_prop[1], p.Qc_now, p.Qc_now_sigma) +
                   LL(Tc_prop[1], p.Tc_now, p.Tc_now_sigma) +
                   LL(Tm_prop[1], p.Tm_now, p.Tm_now_sigma)
            catch
                @info("Model failed to converge, trying again")
            end

            # Accept or reject proposal based on likelihood
            if rand() < exp(ll_prop-ll)
                # Update proposal (jumping) distributions
                Qrm_now_step_sigma = case==1 ? jumpingsigmafactor*abs(Qrm_now_prop - Qrm_now) : Qrm_now_step_sigma
                d_step_sigma  = case==2 ? jumpingsigmafactor*abs(d_prop-d) : d_step_sigma
                Qm_step_sigma = case==3 ? jumpingsigmafactor*abs(p_prop.Qm_now - Qm_now) : Qm_step_sigma
                Qc_step_sigma = case==4 ? jumpingsigmafactor*abs(p_prop.Qc_now - Qc_now) : Qc_step_sigma
                Tc_step_sigma = case==5 ? jumpingsigmafactor*abs(p_prop.Tc_now - Tc_now) : Tc_step_sigma
                Tm_step_sigma = case==6 ? jumpingsigmafactor*abs(p_prop.Tm_now - Tm_now) : Tm_step_sigma
                Q_addtl_step_sigma = case==7 ? jumpingsigmafactor*abs(p_prop.Q_addtl - Q_addtl) : Q_addtl_step_sigma

                # Accept results
                ll = ll_prop
                d = d_prop
                Qrm_now = Qrm_now_prop
                Qm_now = p_prop.Qm_now
                Qc_now = p_prop.Qc_now
                Tc_now = p_prop.Tc_now
                Tm_now = p_prop.Tm_now
                Q_addtl = p_prop.Q_addtl
                Trel = copy(Trel_prop)

                # Just to keep track
                p_prop.Qm_initial = Qm_prop[end]
                p_prop.Qc_initial = Qc_prop[end]
                p_prop.Tc_initial = Tc_prop[end]
                p_prop.Tm_initial = Tm_prop[end]

                # Record accptance
                acceptancedist[n] = true
            end

            # Record results
            ll_dist[n] = ll
            d_dist[n] = d
            Qrm_now_dist[n] = Qrm_now
            Qm_now_dist[n] = Qm_now
            Qc_now_dist[n] = Qc_now
            Tc_now_dist[n] = Tc_now
            Tm_now_dist[n] = Tm_now
            Qm_initial_dist[n] = p_prop.Qm_initial
            Qc_initial_dist[n] = p_prop.Qc_initial
            Tc_initial_dist[n] = p_prop.Tc_initial
            Tm_initial_dist[n] = p_prop.Tm_initial
            Q_addtl_dist[n] = Q_addtl
            Trel_dist[n,:] = Trel
        end

        return (acceptancedist, ll_dist, d_dist, Qrm_now_dist, Qm_initial_dist, Qc_initial_dist, Tc_initial_dist, Tm_initial_dist, Qm_now_dist, Qc_now_dist, Tc_now_dist, Tm_now_dist, Trel_dist)
    end

    # Run the full MCMC inversion
    (acceptancedist, ll_dist, d_dist, Qrm_now_dist, Qm_initial_dist, Qc_initial_dist,
    Tc_initial_dist, Tm_initial_dist, Qm_now_dist, Qc_now_dist, Tc_now_dist,
    Tm_now_dist, Trel_dist) =  MCMC_SMC_bwd(nSteps, d, jumpingsigmafactor, p, Qrm)

    # Summary statistics
    d_mu = nanmean(d_dist[burnin:end])
    d_sigma = nanstd(d_dist[burnin:end])
    Qrm_now_mu = nanmean(Qrm_now_dist[burnin:end])
    Qrm_now_sigma = nanstd(Qrm_now_dist[burnin:end])
    Qm_initial_mu = mean(Qm_initial_dist[burnin:end])
    Qm_initial_sigma = std(Qm_initial_dist[burnin:end])
    Qc_initial_mu = mean(Qc_initial_dist[burnin:end])
    Qc_initial_sigma = std(Qc_initial_dist[burnin:end])
    Tc_initial_mu = mean(Tc_initial_dist[burnin:end]) - 273.15
    Tc_initial_sigma = std(Tc_initial_dist[burnin:end])
    Tm_initial_mu = mean(Tm_initial_dist[burnin:end]) - 273.15
    Tm_initial_sigma = std(Tm_initial_dist[burnin:end])
    Qm_now_mu = mean(Qm_now_dist[burnin:end])
    Qm_now_sigma = std(Qm_now_dist[burnin:end])
    Qc_now_mu = mean(Qc_now_dist[burnin:end])
    Qc_now_sigma = std(Qc_now_dist[burnin:end])
    Tc_now_mu = mean(Tc_now_dist[burnin:end]) - 273.15
    Tc_now_sigma = std(Tc_now_dist[burnin:end])
    Tm_now_mu = mean(Tm_now_dist[burnin:end]) - 273.15
    Tm_now_sigma = std(Tm_now_dist[burnin:end])

    # Run model with average burnt-in parameters
    p.Rc = (p.eta_L/(10.0^d_mu))^(1/3)
    p.Qm_initial = Qm_initial_mu
    p.Qc_initial = Qc_initial_mu
    p.Tc_initial = Tc_initial_mu + 273.15
    p.Tm_initial = Tm_initial_mu + 273.15
    p.Qm_now = Qm_now_mu
    p.Qc_now = Qc_now_mu
    p.Tc_now = Tc_now_mu + 273.15
    p.Tm_now = Tm_now_mu + 273.15
    (Tm_mu, Tc_mu, Qm_mu, Qc_mu, dp_mu) = coolingmodel(p,Qrm*(Qrm_now_mu/Qrm[1]),dt,timevec)

    # Plot results
    lines_to_plot = 500
    cmap = cgrad(reverse(viridis[128:end]))
    dsfactor = round(Int, (nSteps-burnin)/lines_to_plot)
    zdata = downsample(d_dist[burnin:end], dsfactor)'
    Trel_dist_plot = downsample(Trel_dist[burnin:end,:], dsfactor,1)'
    h = plot(TrelObs_time, Trel_dist_plot, line_z=zdata, line=(cmap), linealpha=0.3, label="")
    plot!(h, TrelObs_time, TrelObs, yerror=2*TrelObs_sigma, seriestype=:scatter, color=:darkblue, markersize=2, markerstrokecolor=:auto, label="Data")
    plot!(h, timevec/1000, Tm_mu .- Tm_mu[1], color=:black, label="Model")
    plot!(h, xlabel="Age (Ma)", ylabel="\\Delta T (C)" ,xlims=(0,4), framestyle=:box, grid=:off, legend=:topleft, fg_color_legend=:white)
    summarystring = " d:\t\t$(round(d_mu,sigdigits=3))\n" *
        " Qrm0: $(round(Qrm_now_mu,sigdigits=3))\n\n" *
        " Qc: $(round(Qc_now_mu/1E12,sigdigits=3))\n" *
        " Qm: $(round(Qm_now_mu/1E12,sigdigits=3))\n" *
        " Tc: $(round(Int,Tc_now_mu))\n" *
        " Tm: $(round(Int,Tm_now_mu))\n"
    plot!(h,1,1,seriestype=:scatter, ann=(0.16,130,text(summarystring,8,"Helvetica",:left)),label="")
    display(h)
    savefig(h,"MCMC_bestfit_bwd_$descrip.pdf")

    # Print results
    resultstring = "Reverse-time MCMC inversion\n" *
     " d: $(round(d_mu,sigdigits=4)) +/- $(round(d_sigma,sigdigits=4))\n" *
     " Qrm_now: $(round(Qrm_now_mu,sigdigits=6)) +/- $(round(Qrm_now_sigma,sigdigits=6)) W\n" *
     " Qc_now: $(round(Qc_now_mu,sigdigits=4)) +/- $(round(Qc_now_sigma,sigdigits=4)) W\n" *
     " Qm_now: $(round(Qm_now_mu,sigdigits=4)) +/- $(round(Qm_now_sigma,sigdigits=4)) W\n" *
     " Tc_now: $(round(Tc_now_mu,sigdigits=4)) +/- $(round(Tc_now_sigma,sigdigits=4)) C\n" *
     " Tm_now: $(round(Tm_now_mu,sigdigits=4)) +/- $(round(Tm_now_sigma,sigdigits=4)) C\n" *
     " Qc_initial: $(round(Qc_initial_mu,sigdigits=4)) +/- $(round(Qc_initial_sigma,sigdigits=4)) W\n" *
     " Qm_initial: $(round(Qm_initial_mu,sigdigits=4)) +/- $(round(Qm_initial_sigma,sigdigits=4)) W\n" *
     " Tc_initial: $(round(Tc_initial_mu,sigdigits=4)) +/- $(round(Tc_initial_sigma,sigdigits=4)) C\n" *
     " Tm_initial: $(round(Tm_initial_mu,sigdigits=4)) +/- $(round(Tm_initial_sigma,sigdigits=4)) C\n"
    print(resultstring * " acceptance ratio: $(round(mean(acceptancedist[burnin:end]),sigdigits=2))\n")

## --- Full unconstrained MCMC inversion, forward-time mode

    # While the results of forward and reverse-mode numerical models are fundamentally
    # identical due to time-reversal symmetry, they lead to different ways of
    # exploring the parameter space

    nSteps = 1*10^6
    burnin = 2*10^5
    jumpingsigmafactor = 2.718  # This can be adjusted to optimize acceptance likelihood

    # Start with fresh parameters
    p = Parameters()

    # Initial guesses
    # d = 4
    d = 8
    p.Rc = (p.eta_L/(10.0^d))^(1/3)

    p.Qm_initial = 43.5E12 # Initial mantle heat flux (W)
    p.Qc_initial = 25.8E12 # Initial core heat flux (W)
    p.Tm_initial = 1800.0 # Initial mantle temperature (K)
    p.Tc_initial = 8000.0 # Initial core temperature (K).
    p.Qm_initial_sigma = 20E12
    p.Qc_initial_sigma = 10E12
    p.Tm_initial_sigma = 200
    p.Tc_initial_sigma = 2000

    p.Qm_now = 35E12 # Present mantle heat flux (W)
    p.Qc_now = 10E12 # Present core heat flux (W)
    p.Tm_now = 1300+273.15 # Present mantle temperature (K)
    p.Tc_now = 6200+273.15 # Present core temperature (K). I/O boundary 6230 +/- 500 according to 10.1126/science.1233514
    p.Qm_now_sigma = 10E12
    p.Qc_now_sigma = 10E12
    p.Tm_now_sigma = 100.0
    p.Tc_now_sigma = 1000.0

    descrip = "Qm35-10_Qc10-10"

    # Define a function to conduct a full MCMC inversion
    function MCMC_SMC_fwd(nSteps, d, jumpingsigmafactor, p, Qrm)

        # Allocate arrays to record the stationary distributions
        acceptancedist = fill(false,nSteps)
        ll_dist = Array{Float64}(undef, nSteps)
        d_dist = Array{Float64}(undef, nSteps)
        Qrm_now_dist = Array{Float64}(undef, nSteps)
        Qm_initial_dist = Array{Float64}(undef, nSteps)
        Qc_initial_dist = Array{Float64}(undef, nSteps)
        Tc_initial_dist = Array{Float64}(undef, nSteps)
        Tm_initial_dist = Array{Float64}(undef, nSteps)
        Qm_now_dist = Array{Float64}(undef, nSteps)
        Qc_now_dist = Array{Float64}(undef, nSteps)
        Tc_now_dist = Array{Float64}(undef, nSteps)
        Tm_now_dist = Array{Float64}(undef, nSteps)
        Q_addtl_dist = Array{Float64}(undef, nSteps)
        Trel_dist = Array{Float64}(undef, nSteps,length(TrelObs_time))

        # Proposal distributions
        d_step_sigma = 0.1
        Qrm_now_step_sigma = 0.1E12
        Qm_step_sigma = 0.1E12
        Qc_step_sigma = 0.1E12
        Tc_step_sigma = 100
        Tm_step_sigma = 10
        Q_addtl_step_sigma = 0.1E12

        # First markov chain step / first proposal
        p_prop = deepcopy(p)
        d_prop = log10(p_prop.eta_L/p_prop.Rc^3)
        Qrm_now_prop = coolingmodel_forward_Newton_Qrm(p_prop,20,Qrm,dt,timevec,TrelObs_time*1000,TrelObs,TrelObs_sigma)
        (Tm_prop, Tc_prop, Qm_prop, Qc_prop, dp_prop) = coolingmodel_forward(p_prop, Qrm*(Qrm_now_prop/Qrm[1]), dt, timevec)
        Trel_prop = linterp1s(timevec, Tm_prop .- Tm_prop[1], TrelObs_time*1000)

        h = plot(TrelObs_time, TrelObs, yerror=2*TrelObs_sigma, seriestype=:scatter, color=:darkblue, markersize=2, markerstrokecolor=:auto, label="Data")
        plot!(h, timevec/1000, Tm_prop .- Tm_prop[1], label="Model: Qrm_now=$(round(Qrm_now_prop,sigdigits=3)),d=$(round(d,sigdigits=3))")
        plot!(h, xlabel="Age (Ga)", ylabel="\\Delta T (C)", grid=:off, legend=:topleft, fg_color_legend=:white)
        summarystring = " d:\t\t$(round(d,sigdigits=2))\n" *
            " Qrm_now: $(round(Qrm_now_prop,sigdigits=3))\n\n" *
            " Qc_i: $(round(p.Qc_initial/1E12,sigdigits=1))\n" *
            " Qm_i: $(round(p.Qm_initial/1E12,sigdigits=1))\n" *
            " Tc_i: $(round(Int,p.Tc_initial))\n" *
            " Tm_i: $(round(Int,p.Tm_initial))\n"
        plot!(h,1,1,seriestype=:scatter, ann=(0.16,130,text(summarystring,8,"Helvetica",:left)),label="")
        display(h)
        savefig(h,"MCMC_firstproposal.pdf")

        # Log likelihood of first proposal
        ll_prop = LL(Trel_prop,TrelObs,TrelObs_sigma*2) + LL(d_prop,5,5) +
           LL(Qrm_now_prop, p.Qrm_now, p.Qrm_now_sigma) +
           LL(p_prop.Qm_initial, p.Qm_initial, p.Qm_initial_sigma) +
           LL(p_prop.Qc_initial, p.Qc_initial, p.Qc_initial_sigma) +
           LL(p_prop.Tc_initial, p.Tc_initial, p.Tc_initial_sigma) +
           LL(p_prop.Tm_initial, p.Tm_initial, p.Tm_initial_sigma) +
           LL(Qm_prop[1],p.Qm_now,p.Qm_now_sigma) +
           LL(Qc_prop[1],p.Qc_now,p.Qc_now_sigma) +
           LL(Tc_prop[1],p.Tc_now,p.Tc_now_sigma) +
           LL(Tm_prop[1],p.Tm_now,p.Tm_now_sigma)

        # Accept and record first proposal
        ll = ll_prop
        d = d_prop
        Qrm_now = Qrm_now_prop
        Qm_initial = p_prop.Qm_initial
        Qc_initial = p_prop.Qc_initial
        Tc_initial = p_prop.Tc_initial
        Tm_initial = p_prop.Tm_initial
        Q_addtl = p_prop.Q_addtl
        Qm_now = Qm_prop[1]
        Qc_now = Qc_prop[1]
        Tc_now = Tc_prop[1]
        Tm_now = Tm_prop[1]
        Trel = copy(Trel_prop)

        ll_dist[1] = ll
        d_dist[1] = d
        Qrm_now_dist[1] = Qrm_now_prop
        Qm_initial_dist[1] = Qm_initial
        Qc_initial_dist[1] = Qc_initial
        Tc_initial_dist[1] = Tc_initial
        Tm_initial_dist[1] = Tm_initial
        Qm_now_dist[1] = Qm_now
        Qc_now_dist[1] = Qc_now
        Tc_now_dist[1] = Tc_now
        Tm_now_dist[1] = Tm_now
        Q_addtl_dist[1] = Q_addtl
        Trel_dist[1,:] = Trel

        # Run the Markov chain / Metropolis walker
        @showprogress for n=1:nSteps

            # Randomly choose which parameter to adjust. The main benefit of his
            # form of "Metropolis-in-Gibbs" / "Alternate conditional sampling"
            # is that it allows us to optimize the proposal (jumping) distributions
            # for each parameter individually
            case = ceil(rand()*6)

            # Propose new parameters
            Qrm_now_prop = case==1 ? abs(Qrm_now + randn()*Qrm_now_step_sigma) : Qrm_now
            d_prop  = case==2 ? abs(d + randn()*d_step_sigma - 1.0) + 1.0 : d
            p_prop.Rc = abs(p_prop.eta_L / (10.0^d_prop))^(1/3)
            p_prop.Qm_initial = case==3 ? abs(Qm_initial + randn()*Qm_step_sigma - 1E12) + 1E12 : Qm_initial
            p_prop.Qc_initial = case==4 ? abs(Qc_initial + randn()*Qc_step_sigma - 1E12) + 1E12 : Qc_initial
            p_prop.Tc_initial = case==5 ? abs(Tc_initial + randn()*Tc_step_sigma - 4000) + 4000 : Tc_initial
            p_prop.Tm_initial = case==6 ? abs(Tm_initial + randn()*Tm_step_sigma - 1000) + 1000 : Tm_initial
            p_prop.Q_addtl = case==7 ? abs(Q_addtl + randn()*Q_addtl_step_sigma) : Q_addtl

            # If proposal is so bad that we hit a math error, let ll be NaN to
            # auto-reject the proposal
            Trel_prop = NaN
            ll_prop = NaN
            try
                # Calculate proposed solution
                (Tm_prop, Tc_prop, Qm_prop, Qc_prop, dp_prop)  = coolingmodel_forward(p_prop, Qrm*(Qrm_now_prop/Qrm[1]), dt, timevec)
                Trel_prop = linterp1s(timevec, Tm_prop .- Tm_prop[1], TrelObs_time*1000)

                # Calculate log likelihood of new proposal
                ll_prop = LL(Trel_prop,TrelObs,TrelObs_sigma*2) + LL(d_prop,5,5) +
                   LL(Qrm_now_prop, p.Qrm_now, p.Qrm_now_sigma) +
                   LL(p_prop.Qm_initial, p.Qm_initial, p.Qm_initial_sigma) +
                   LL(p_prop.Qc_initial, p.Qc_initial, p.Qc_initial_sigma) +
                   LL(p_prop.Tc_initial, p.Tc_initial, p.Tc_initial_sigma) +
                   LL(p_prop.Tm_initial, p.Tm_initial, p.Tm_initial_sigma) +
                   LL(Qm_prop[1], p.Qm_now, p.Qm_now_sigma) +
                   LL(Qc_prop[1], p.Qc_now, p.Qc_now_sigma) +
                   LL(Tc_prop[1], p.Tc_now, p.Tc_now_sigma) +
                   LL(Tm_prop[1], p.Tm_now, p.Tm_now_sigma)
            catch
                @info("Model failed to converge, trying again")
            end

            # Accept or reject proposal based on likelihood
            if rand() < exp(ll_prop-ll)
                # Update proposal (jumping) distributions
                Qrm_now_step_sigma = case==1 ? jumpingsigmafactor*abs(Qrm_now_prop - Qrm_now) : Qrm_now_step_sigma
                d_step_sigma  = case==2 ? jumpingsigmafactor*abs(d_prop-d) : d_step_sigma
                Qm_step_sigma = case==3 ? jumpingsigmafactor*abs(p_prop.Qm_initial - Qm_initial) : Qm_step_sigma
                Qc_step_sigma = case==4 ? jumpingsigmafactor*abs(p_prop.Qc_initial - Qc_initial) : Qc_step_sigma
                Tc_step_sigma = case==5 ? jumpingsigmafactor*abs(p_prop.Tc_initial - Tc_initial) : Tc_step_sigma
                Tm_step_sigma = case==6 ? jumpingsigmafactor*abs(p_prop.Tm_initial - Tm_initial) : Tm_step_sigma
                Q_addtl_step_sigma = case==7 ? jumpingsigmafactor*abs(p_prop.Q_addtl - Q_addtl) : Q_addtl_step_sigma

                # Accept results
                ll = ll_prop
                d = d_prop
                Qrm_now = Qrm_now_prop
                Qm_initial = p_prop.Qm_initial
                Qc_initial = p_prop.Qc_initial
                Tc_initial = p_prop.Tc_initial
                Tm_initial = p_prop.Tm_initial
                Q_addtl = p_prop.Q_addtl
                Trel = copy(Trel_prop)

                # These can be used to recalibrate Beta_c and Beta_d
                p_prop.Qm_now = Qm_prop[1]
                p_prop.Qc_now = Qc_prop[1]
                p_prop.Tc_now = Tc_prop[1]
                p_prop.Tm_now = Tm_prop[1]

                # Record accptance
                acceptancedist[n] = true
            end

            # Record results
            ll_dist[n] = ll
            d_dist[n] = d
            Qrm_now_dist[n] = Qrm_now
            Qm_initial_dist[n] = Qm_initial
            Qc_initial_dist[n] = Qc_initial
            Tc_initial_dist[n] = Tc_initial
            Tm_initial_dist[n] = Tm_initial
            Qm_now_dist[n] = p_prop.Qm_now
            Qc_now_dist[n] = p_prop.Qc_now
            Tc_now_dist[n] = p_prop.Tc_now
            Tm_now_dist[n] = p_prop.Tm_now
            Q_addtl_dist[n] = Q_addtl
            Trel_dist[n,:] = Trel
        end

        return (acceptancedist, ll_dist, d_dist, Qrm_now_dist, Qm_initial_dist, Qc_initial_dist, Tc_initial_dist, Tm_initial_dist, Qm_now_dist, Qc_now_dist, Tc_now_dist, Tm_now_dist, Trel_dist)
    end

    # Run the full MCMC inversion
    (acceptancedist, ll_dist, d_dist, Qrm_now_dist, Qm_initial_dist, Qc_initial_dist,
    Tc_initial_dist, Tm_initial_dist, Qm_now_dist, Qc_now_dist, Tc_now_dist,
    Tm_now_dist, Trel_dist) =  MCMC_SMC_fwd(nSteps, d, jumpingsigmafactor, p, Qrm)

    # Summary statistics
    d_mu = nanmean(d_dist[burnin:end])
    d_sigma = nanstd(d_dist[burnin:end])
    Qrm_now_mu = nanmean(Qrm_now_dist[burnin:end])
    Qrm_now_sigma = nanstd(Qrm_now_dist[burnin:end])
    Qm_initial_mu = mean(Qm_initial_dist[burnin:end])
    Qm_initial_sigma = std(Qm_initial_dist[burnin:end])
    Qc_initial_mu = mean(Qc_initial_dist[burnin:end])
    Qc_initial_sigma = std(Qc_initial_dist[burnin:end])
    Tc_initial_mu = mean(Tc_initial_dist[burnin:end]) - 273.15
    Tc_initial_sigma = std(Tc_initial_dist[burnin:end])
    Tm_initial_mu = mean(Tm_initial_dist[burnin:end]) - 273.15
    Tm_initial_sigma = std(Tm_initial_dist[burnin:end])
    Qm_now_mu = mean(Qm_now_dist[burnin:end])
    Qm_now_sigma = std(Qm_now_dist[burnin:end])
    Qc_now_mu = mean(Qc_now_dist[burnin:end])
    Qc_now_sigma = std(Qc_now_dist[burnin:end])
    Tc_now_mu = mean(Tc_now_dist[burnin:end]) - 273.15
    Tc_now_sigma = std(Tc_now_dist[burnin:end])
    Tm_now_mu = mean(Tm_now_dist[burnin:end]) - 273.15
    Tm_now_sigma = std(Tm_now_dist[burnin:end])

    # Run model with average burnt-in parameters
    p.Rc = (p.eta_L/(10.0^d_mu))^(1/3)
    p.Qm_initial = Qm_initial_mu
    p.Qc_initial = Qc_initial_mu
    p.Tc_initial = Tc_initial_mu + 273.15
    p.Tm_initial = Tm_initial_mu + 273.15
    p.Qm_now = Qm_now_mu
    p.Qc_now = Qc_now_mu
    p.Tc_now = Tc_now_mu + 273.15
    p.Tm_now = Tm_now_mu + 273.15
    (Tm_mu, Tc_mu, Qm_mu, Qc_mu, dp_mu) = coolingmodel_forward(p,Qrm*(Qrm_now_mu/Qrm[1]),dt,timevec)

    # Plot results
    lines_to_plot = 500
    cmap = cgrad(reverse(viridis[128:end]))
    dsfactor = round(Int, (nSteps-burnin)/lines_to_plot)
    zdata = downsample(d_dist[burnin:end], dsfactor)'
    Trel_dist_plot = downsample(Trel_dist[burnin:end,:], dsfactor,1)'
    h = plot(TrelObs_time, Trel_dist_plot, line_z=zdata, line=(cmap), linealpha=0.3, label="")
    plot!(h, TrelObs_time, TrelObs, yerror=2*TrelObs_sigma, seriestype=:scatter, color=:darkblue, markersize=2, markerstrokecolor=:auto, label="Data")
    plot!(h, timevec/1000, Tm_mu .- Tm_mu[1], color=:black, label="Model")
    plot!(h, xlabel="Age (Ma)", ylabel="\\Delta T (C)" ,xlims=(0,4), framestyle=:box, grid=:off, legend=:topleft, fg_color_legend=:white)
    summarystring = " d:\t\t$(round(d_mu,sigdigits=2))\n" *
        " Qrm_now: $(round(Qrm_now_mu,sigdigits=3))\n\n" *
        " Qc_i: $(round(Qc_initial_mu/1E12,sigdigits=1))\n" *
        " Qm_i: $(round(Qm_initial_mu/1E12,sigdigits=1))\n" *
        " Tc_i: $(round(Int,Tc_initial_mu))\n" *
        " Tm_i: $(round(Int,Tm_initial_mu))\n"
    plot!(h,1,1,seriestype=:scatter, ann=(0.16,130,text(summarystring,8,"Helvetica",:left)),label="")
    display(h)
    savefig(h,"MCMC_bestfit_fwd_$descrip.pdf")

    # Print results
    resultstring = "Forward-time MCMC inversion\n" *
     " d: $(round(d_mu,sigdigits=4)) +/- $(round(d_sigma,sigdigits=4))\n" *
     " Qrm_now: $(round(Qrm_now_mu,sigdigits=6)) +/- $(round(Qrm_now_sigma,sigdigits=6)) W\n" *
     " Qc_now: $(round(Qc_now_mu,sigdigits=4)) +/- $(round(Qc_now_sigma,sigdigits=4)) W\n" *
     " Qm_now: $(round(Qm_now_mu,sigdigits=4)) +/- $(round(Qm_now_sigma,sigdigits=4)) W\n" *
     " Tc_now: $(round(Tc_now_mu,sigdigits=4)) +/- $(round(Tc_now_sigma,sigdigits=4)) C\n" *
     " Tm_now: $(round(Tm_now_mu,sigdigits=4)) +/- $(round(Tm_now_sigma,sigdigits=4)) C\n" *
     " Qc_initial: $(round(Qc_initial_mu,sigdigits=4)) +/- $(round(Qc_initial_sigma,sigdigits=4)) W\n" *
     " Qm_initial: $(round(Qm_initial_mu,sigdigits=4)) +/- $(round(Qm_initial_sigma,sigdigits=4)) W\n" *
     " Tc_initial: $(round(Tc_initial_mu,sigdigits=4)) +/- $(round(Tc_initial_sigma,sigdigits=4)) C\n" *
     " Tm_initial: $(round(Tm_initial_mu,sigdigits=4)) +/- $(round(Tm_initial_sigma,sigdigits=4)) C\n"
    print(resultstring * " acceptance ratio: $(round(mean(acceptancedist[burnin:end]),sigdigits=2))\n")

## --- End of File

max_stars_repo_name: brenhinkeller/SecularCooling

max_issues_repo_name: brenhinkeller/SecularCooling

max_forks_repo_name: brenhinkeller/SecularCooling


text:
using DataFrames
using CSV
using StatsBase
using StableRNGs

# https://github.com/szilard/benchm-ml/blob/master/0-init/2-gendata.txt
# 2005 2006 2007

# 2005
url2005 = "https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/HG7NV7/JTFT25"

# 2006
url2006 = "https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/HG7NV7/EPIFFT"

# 2007
url2007 = "https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/HG7NV7/2BHLWK"

download(url2005, joinpath(@__DIR__, "data", "2005.csv.bz2"))
download(url2006, joinpath(@__DIR__, "data", "2006.csv.bz2"))
download(url2007, joinpath(@__DIR__, "data", "2007.csv.bz2"))

run(`bunzip2 ./data/2005.csv.bz2`)
run(`bunzip2 ./data/2006.csv.bz2`)
run(`bunzip2 ./data/2007.csv.bz2`)

df1 = CSV.file(joinpath(@__DIR__, "data", "2005.csv")) |> DataFrame;
df2 = CSV.file(joinpath(@__DIR__, "data", "2006.csv")) |> DataFrame;

df1 = df1[df1[!, :DepDelay] .!= "NA", :];
df2 = df2[df2[!, :DepDelay] .!= "NA", :];

df1[!, :dep_delayed_15min] = ifelse.(parse.(Int, df1[!, :DepDelay]) .>= 15, "Y", "N");
df2[!, :dep_delayed_15min] = ifelse.(parse.(Int, df2[!, :DepDelay]) .>= 15, "Y", "N");

cols = Symbol.(["Month", "DayofMonth", "DayOfWeek", "DepTime", "UniqueCarrier", "Origin", "Dest", "Distance","dep_delayed_15min"])

df1 = df1[!, cols];
df2 = df2[!, cols];

df = vcat(df1, df2);

for k in Symbol.(["Month", "DayofMonth", "DayOfWeek"])
    df[!, k] .= "c-" .* string.(df[!, k])
end

df[!, :DepTime] .= parse.(Int, df[!, :DepTime]);

CSV.write(joinpath(@__DIR__, "data", "air-process.csv"), df)

# Sampling

rng = StableRNG(2020)

ids1k = sample(rng, 1:size(df, 1), 1000, replace = false);
ids10k = sample(rng, 1:size(df, 1), 10000, replace = false);

df1k = df[ids1k, :];
df10k = df[ids10k, :];

CSV.write(joinpath(@__DIR__, "data", "air-1k.csv"), df1k);
CSV.write(joinpath(@__DIR__, "data", "air-10k.csv"), df10k);

###################
test data

df = urldownload("https://s3.amazonaws.com/benchm-ml--main/test.csv", true) |> DataFrame
CSV.write(joinpath(@__DIR__, "data", "air-test.csv"), df)

df = urldownload("https://s3.amazonaws.com/benchm-ml--main/train-0.1m.csv", true) |> DataFrame
CSV.write(joinpath(@__DIR__, "data", "air-0.1m.csv"), df)

df = urldownload("https://s3.amazonaws.com/benchm-ml--main/train-1m.csv", true) |> DataFrame
CSV.write(joinpath(@__DIR__, "data", "air-1m.csv"), df)

# optional
df = urldownload("https://s3.amazonaws.com/benchm-ml--main/train-10m.csv", true) |> DataFrame
CSV.write(joinpath(@__DIR__, "data", "air-10m.csv"), df)

ext: jl


text:
struct Propagator{T,D} <: AbstractParameterizedFunction{false}
    U::Matrix{Complex{T}}
    Δt::T
    dims::Dims{D}
end

dimsmatch(U::Propagator,A::QuObject) = U.dims==dims(A) || throw(DimensionMismatch("subspace dimensions must match"))

function (U::Propagator)(ψ::Ket)
    dimsmatch(U,ψ)
    return Ket(U.U*data(ψ),dims(ψ))
end

function (U::Propagator)(ψ::Ket, n::Integer)
    dimsmatch(U,ψ)
    return Ket(U.U^n*data(ψ),dims(ψ))
end

function (U::Propagator)(ρ::Operator)
    dimsmatch(U,ρ)
    vecρ = vec(data(ρ))
    d = prod(dims(ρ))
    return Operator(reshape(U.U*vecρ,(d,d)),dims(ρ))
end

function (U::Propagator)(ρ::Operator, n::Integer)
    dimsmatch(U,ρ)
    vecρ = vec(data(ρ))
    d = prod(dims(ρ))
    return Operator(reshape(U.U^n*vecρ,(d,d)),dims(ρ))
end

# Convert a Propagator to an Operator
Operator(U::Propagator) = Operator(U.U,U.dims)

# Convert a Liouvillian to a Propagator
function SchrodingerProp(L::Liouvillian,tspan,alg=Vern8();kwargs...)
    d = prod(dims(L))
    prob = ODEProblem(L,Matrix{ComplexF64}(I,d,d),tspan)
    sol  = solve(prob,alg;save_start=false,saveat=tspan[end],abstol=1E-8,reltol=1E-6,kwargs...)
    return Propagator(sol.u[end],float(tspan[2]-tspan[1]),dims(L))
end

function LindbladProp(L::Liouvillian,tspan,alg=Tsit5();kwargs...)
    d = prod(dims(L))^2
    prob = ODEProblem(L,Matrix{ComplexF64}(I,d,d),tspan)
    sol  = solve(prob,alg;save_start=false,saveat=tspan[end],abstol=1E-7,reltol=1E-5,kwargs...)
    return Propagator(sol.u[end],float(tspan[2]-tspan[1]),dims(L))
end


text:
module RRIFT

using Statistics

include("utils.jl")
export apply_mask, crop, negatives_to_zero!, ccc

using CancerImagingArchive: series, images
using DICOM: dcm_parse, lookup
using Perfusion: make_folder
include("download_gbm_data.jl")
export download_invivo_studies, download_invivo_masks, download_invivo_preprocessed
export get_mask

using MAT: matread, matwrite
using Perfusion: @extract, fit_relaxation, signal_to_concentration
include("preprocess.jl")
export compute_concentration, preprocess_dicom_to_mat, load_preprocessed_mat

using Statistics: mean, quantile, std
using Perfusion: fit_model
include("analysis.jl")
export process_patients

import Perfusion

fit_errm(; kwargs...) = Perfusion.fit_errm_lls(; kwargs...).estimates
fit_cerrm(; kwargs...) = Perfusion.fit_cerrm_lls(; kwargs...).estimates

function fit_cerrm_with_rrift(; kwargs...)
    return Perfusion.fit_rrift_with_cerrm(; kwargs...).estimates
end

positive_only_mask = Perfusion.positive_only_mask

export fit_errm, fit_cerrm, fit_cerrm_with_rrift
export positive_only_mask

end # module


text:
struct MaintenanceSchedulingIndices <: AbstractModelIndices
    hours::Vector{Int}
    plants::Vector{Plant}
    segments::Vector{Int}
    bids::Vector{Int}
end
plants(indices::MaintenanceSchedulingIndices) = indices.plants

function HydroModels.modelindices(data::MaintenanceSchedulingData, horizon::Horizon; areas::Vector{Area} = [0], rivers::Vector{River} = [:All])
    hours = collect(1:num_hours(horizon))
    plants = plants_in_areas_and_rivers(hydrodata(data), areas, rivers)
    if isempty(plants)
        error("No plants in given set of price areas and rivers")
    end
    segments = collect(1:2)
    bids = collect(1:length(data.bidlevels[1]))
    return MaintenanceSchedulingIndices(hours, plants, segments, bids)
end


text:
abstract type AbstractNLPEvaluator end

function optimize! end

max_stars_repo_name: exanauts/proxALM

max_issues_repo_name: exanauts/proxALM

max_forks_repo_name: exanauts/proxALM


text:
@testset "numerical checks" begin
    T = [-2.5 2 0.5; 1 -3 2; 1 2 -3]

    C = [0, 2.0, -6.0]
    S = PhaseSet(C)

    model = BoundedFluidQueue(T,S,nodes[end])
    am = augment_model(model)

    mtypes = (DGMesh,
        FVMesh,
        FRAPMesh)
    for mtype in mtypes
        @testset "mtype" begin
            nodes = collect(0:0.1:12)
            order = 5
            msh = mtype(nodes,order)

            ####
            # non-augmented model
            ####
            dq = DiscretisedFluidQueue(model,msh)
            generator = build_full_generator(dq)
            b = zeros(1,size(generator,1))
            b[1] = 1.0
            if mtype==FVMesh
                generator[:,1] .= [
                    ones(N₋(model.S));
                    repeat(Δ(msh),n_phases(model));
                    ones(N₊(model.S))]
            else
                generator[:,1] .= 1.0
            end
            stationary_coeffs = b/generator.B
            d = SFMDistribution(stationary_coeffs[:],dq)
            stationary_cdf_estimate = (mtype!=FRAPMesh) ? cdf(d) : cdf(d, normalised_closing_operator_cdf)
            (mtype==FRAPMesh) && (stationary_cdf_estimate_naive = cdf(d, naive_normalised_closing_operator_cdf))
            analytical_cdf = stationary_distribution_x(model)[3]
            x_vec = nodes[1]:0.23:nodes[end]
            @testset "numerics - cdf" begin
                for x in x_vec
                    @test isapprox(analytical_cdf(x),stationary_cdf_estimate.(x,1:3), rtol=1e-2)
                    (mtype==FRAPMesh) && (@test isapprox(analytical_cdf(x),stationary_cdf_estimate_naive.(x,1:3), rtol=1e-2))
                end
            end

            ####
            # augmented model
            ####
            dq_am = DiscretisedFluidQueue(am,msh)
            generator_am = build_full_generator(dq_am)
            b_am = zeros(1,size(generator_am,1))
            b_am[1] = 1.0
            if mtype==FVMesh
                generator_am[:,1] .= [
                    ones(N₋(am.S));
                    repeat(Δ(msh),n_phases(am));
                    ones(N₊(am.S))]
            else
                generator_am[:,1] .= 1.0
            end
            stationary_coeffs_am = b_am/generator_am.B
            d_am = SFMDistribution(stationary_coeffs_am[:],dq_am)
            stationary_cdf_estimate_am = (mtype!=FRAPMesh) ? cdf(d_am) : (stationary_cdf_estimate_am_naive = cdf(d_am, normalised_closing_operator_cdf))
            (mtype==FRAPMesh) && (stationary_cdf_estimate_am_naive = cdf(d_am, naive_normalised_closing_operator_cdf))
            analytical_cdf_am = stationary_distribution_x(am)[3]
            x_vec = nodes[1]:0.23:nodes[end]
            @testset "numerics - augmented model - cdf" begin
                for x in x_vec
                    @test isapprox(analytical_cdf_am(x),stationary_cdf_estimate_am.(x,1:4), rtol=1e-2)
                    (mtype==FRAPMesh) && (@test isapprox(analytical_cdf_am(x),stationary_cdf_estimate_am_naive.(x,1:4), rtol=1e-2))
                end
            end

            @testset "numerics - augmented model vs vanilla cdf" begin
                for x in x_vec
                    @test isapprox(stationary_cdf_estimate.(x,2:3),stationary_cdf_estimate_am.(x,3:4), rtol=1e-2 )
                    @test isapprox(stationary_cdf_estimate.(x,1),sum(stationary_cdf_estimate_am.(x,1:2)), rtol=1e-2 )
                    (mtype==FRAPMesh) && (@test isapprox(stationary_cdf_estimate_naive.(x,2:3),stationary_cdf_estimate_am_naive.(x,3:4), rtol=1e-2 ))
                    (mtype==FRAPMesh) && (@test isapprox(stationary_cdf_estimate_naive.(x,1),sum(stationary_cdf_estimate_am_naive.(x,1:2)), rtol=1e-2 ))
                end
            end
        end
    end
end


text:
using Catalyst, Test, OrdinaryDiffEq, LinearAlgebra

rn = @reaction_network begin
    (a1,a2), C <--> 0
    (k1,k2), A + B <--> C
    b1, 0 <-- B
end a1 a2 k1 k2 b1

rn2 = @reaction_network begin
    a1, C --> 0
    a2, 0 --> C
    k1, A + B --> C
    k2, C --> A + B
    b1, B --> 0
end a1 a2 k1 k2 b1

ps = ones(5)
u0 = [10.0,20.0,0.0]
tspan = (0,100.0)
oprob1 = ODEProblem(rn, u0, tspan, ps)
oprob2 = ODEProblem(rn2, u0, tspan, ps)
sol1 = solve(oprob1, Tsit5())
sol2 = solve(oprob2, Tsit5())

tv = range(tspan[1], tspan[2], length=100)
@test norm(sol1 - sol2, Inf) ≈ 0


text:
Xc1, yc1 = gen_classif(classes=["M", "F"])
Xc2, yc2 = gen_classif(classes=["A", "B", "C"])

@testset "LogRegClf" begin
    m, f = simple_test_classif_prob(LogisticClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:classes, :coef, :intercept)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
end

@testset "LogRegCVClf" begin
    m, f = simple_test_classif_prob(LogisticCVClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:classes, :coef, :intercept, :Cs, :l1_ratios, :coefs_paths, :scores, :C, :l1_ratio)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
end

@testset "PAClf" begin
    m, f = simple_test_classif(PassiveAggressiveClassifier(), Xc1, yc1; dummybinary=true)
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
end

@testset "PerceptronClf" begin
    m, f = simple_test_classif(PerceptronClassifier(), Xc2, yc2; dummybinary=true)
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
end

@testset "RidgeClf" begin
    m, f = simple_test_classif(RidgeClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
end

@testset  "RidgeCVClf" begin
    # NOTE: this is a hack to prevent getting a bunch of irrelevant sklearn warnings
    # which are due to version compat when using ridgecv (not relevant for us)
    old_stdout = stdout
    old_stderr = stderr
    _ = redirect_stdout()
    _ = redirect_stderr()
    m, f = simple_test_classif(RidgeCVClassifier(), Xc1, yc1)
    redirect_stdout(old_stdout)
    redirect_stderr(old_stderr)
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
end

# NOTE: SGD classifier with few points is tricky which is why we remove the dummy test
@testset "SGDClf" begin
    ((ScikitLearn.Skcore).pyimport("numpy.random")).seed(0)
    m, f = simple_test_classif(SGDClassifier(), Xc2, yc2; nodummy=true)
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
    ((ScikitLearn.Skcore).pyimport("numpy.random")).seed(0)
    m, f = simple_test_classif_prob(ProbabilisticSGDClassifier(), Xc2, yc2; nodummy=true)
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:target_scitype] == AbstractVector{<:Finite}
    @test !isempty(infos[:docstring])
end


text:
# BAM Reader
# ==========

"""
    BAMReader(input::IO; index=nothing)

Create a data reader of the BAM file format.

# Arguments
* `input`: data source
* `index=nothing`: filepath to a random access index (currently *bai* is Supported)
"""
type BAMReader{T} <: Bio.IO.AbstractReader
    stream::BGZFStreams.BGZFStream{T}
    header::SAMHeader
    start_offset::BGZFStreams.VirtualOffset
    refseqnames::Vector{String}
    refseqlens::Vector{Int}
    index::Nullable{BAI}
end

function BAMReader(input::IO; index=nothing)
    if isa(index, AbstractString)
        index = BAI(index)
    else
        if index != nothing
            error("unrecognizable index argument")
        end
    end
    reader = init_bam_reader(input)
    reader.index = index
    return reader
end

function Base.eltype{T}(::Type{BAMReader{T}})
    return BAMRecord
end

function Bio.IO.stream(reader::BAMReader)
    return reader.stream
end

function Base.show(io::IO, reader::BAMReader)
    println(io, summary(reader), ":")
      print(io, "  number of contigs: ", length(reader.refseqnames))
end

function header(reader::BAMReader; fillSQ::Bool=false)
    if fillSQ
        if !isempty(find(reader.header, "SQ"))
            throw(ArgumentError("SAM header already has SQ records"))
        end
        for (name, len) in zip(reader.refseqnames, reader.refseqlens)
            push!(reader.header, SAMMetaInfo("SQ", ["SN" => name, "LN" => len]))
        end
    end
    return reader.header
end

function Base.seek(reader::BAMReader, voffset::BGZFStreams.VirtualOffset)
    seek(reader.stream, voffset)
end

function Base.seekstart(reader::BAMReader)
    seek(reader.stream, reader.start_offset)
end

# Initialize a BAM reader by reading the header section.
function init_bam_reader(input::IO)
    stream = BGZFStreams.BGZFStream(input)

    # magic bytes
    B = read(stream, UInt8)
    A = read(stream, UInt8)
    M = read(stream, UInt8)
    x = read(stream, UInt8)
    if B != UInt8('B') || A != UInt8('A') || M != UInt8('M') || x != 0x01
        error("input was not a valid BAM file")
    end

    # SAM header
    textlen = read(stream, Int32)
    samreader = SAMReader(IOBuffer(read(stream, UInt8, textlen)))

    # reference sequences
    refseqnames = String[]
    refseqlens = Int[]
    n_refs = read(stream, Int32)
    for _ in 1:n_refs
        namelen = read(stream, Int32)
        data = read(stream, UInt8, namelen)
        seqname = unsafe_string(pointer(data))
        seqlen = read(stream, Int32)
        push!(refseqnames, seqname)
        push!(refseqlens, seqlen)
    end

    reader = BAMReader(
        stream,
        samreader.header,
        isa(input, Pipe) ? BGZFStreams.VirtualOffset(0, 0) : BGZFStreams.virtualoffset(stream),
        refseqnames,
        refseqlens,
        Nullable{BAI}())

    return reader
end

function Base.read!(reader::BAMReader, record::BAMRecord)
    unsafe_read(
        reader.stream,
        pointer_from_objref(record),
        BAM_FIXED_FIELDS_BYTES)
    dsize = data_size(record)
    if length(record.data) < dsize
        resize!(record.data, dsize)
    end
    unsafe_read(reader.stream, pointer(record.data), dsize)
    record.refseqnames = reader.refseqnames
    return record
end

function Base.start(reader::BAMReader)
    return BAMRecord()
end

function Base.done(reader::BAMReader, rec)
    return eof(reader)
end

function Base.next(reader::BAMReader, rec)
    read!(reader, rec)
    return copy(rec), rec
end


text:
using KeckResearch
using Test

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


text:
###Dispatch Formulations##

"""
This function creates the minimal themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{PowerSystems.ThermalGen}, category_formulation::Type{D}, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {D <: AbstractThermalDispatchForm, S <: PM.AbstractPowerFormulation}

    p_th = activepowervariables(m, sys.generators.thermal, sys.time_periods);

    varnetinjectiterate!(netinjection.var_active, p_th, sys.time_periods, sys.generators.thermal)

    activepower(m, sys.generators.thermal, category_formulation, system_formulation, sys.time_periods)

    cost = variablecost(m, sys.generators.thermal, category_formulation, system_formulation)

    add_to_cost!(m, cost)

end

function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{PowerSystems.ThermalGen}, category_formulation::Type{D}, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {D <: AbstractThermalDispatchForm, S <: AbstractACPowerModel}

    constructdevice!(m, netinjection, category, category_formulation, PM.AbstractPowerFormulation, sys; args...)

    q_th = reactivepowervariables(m, sys.generators.thermal, sys.time_periods);

    varnetinjectiterate!(netinjection.var_reactive, q_th, sys.time_periods, sys.generators.thermal)

    reactivepower(m, sys.generators.thermal, category_formulation, system_formulation, sys.time_periods)

end

function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{PowerSystems.ThermalGen}, category_formulation::Type{ThermalRampLimitDispatch}, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {S <: PM.AbstractPowerFormulation}

    constructdevice!(m, netinjection, category, AbstractThermalDispatchForm, PM.AbstractPowerFormulation, sys; args...)

    #rampargs = pairs((;(k=>v for (k,v) in pairs(args) if k in [:initalpower])...))

    rampconstraints(m, sys.generators.thermal, category_formulation, system_formulation, sys.time_periods; args...)

    return

end

###Commitment Formulations##

"""
This function creates the minimal the minimal thermal commitment formulation
"""
function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{PowerSystems.ThermalGen}, category_formulation::Type{D}, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {D <: AbstractThermalCommitmentForm, S <: AbstractDCPowerModel}

    p_th = activepowervariables(m, sys.generators.thermal, sys.time_periods);

    commitmentvariables(m, sys.generators.thermal, sys.time_periods)

    varnetinjectiterate!(netinjection.var_active, p_th, sys.time_periods, sys.generators.thermal)

    activepower(m, sys.generators.thermal, category_formulation, system_formulation, sys.time_periods)

    variable_cost = variablecost(m, sys.generators.thermal, AbstractThermalDispatchForm, system_formulation)

    commitment_cost = commitmentcost(m, sys.generators.thermal, category_formulation, system_formulation)

    add_to_cost!(m, variable_cost)

    add_to_cost!(m, commitment_cost)

end

"""
This function adds constraints to the minimal thermal commitment formulation
"""
function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{PowerSystems.ThermalGen}, category_formulation::Type{StandardThermalCommitment}, system_formulation::Type{S}, sys::PowerSystems.PowerSystem; args...) where {S <: AbstractDCPowerModel}

    constructdevice!(m, netinjection, category, AbstractThermalCommitmentForm, AbstractDCPowerModel, sys; args...)

    #commitargs = pairs((;(k=>v for (k,v) in pairs(args) if k in [:initalstatus,:initialonduration,:initialoffduration])...)) #this isn't strictly needed, could delete for cleanliness

    commitmentconstraints(m, sys.generators.thermal, category_formulation, system_formulation, sys.time_periods; args...)

    #rampargs = pairs((;(k=>v for (k,v) in pairs(args) if k in [:initalpower])...)) #this isn't strictly needed, could delete for cleanliness

    rampconstraints(m, sys.generators.thermal, category_formulation, system_formulation, sys.time_periods; args...)

    timeconstraints(m, sys.generators.thermal, category_formulation, system_formulation, sys.time_periods; args...)

end

"""
This function adds constraints to the minimal thermal commitment formulation
"""
function constructdevice!(m::JuMP.AbstractModel, netinjection::BalanceNamedTuple, category::Type{PowerSystems.ThermalGen}, category_formulation::Type{StandardThermalCommitment}, system_formulation::Type{CopperPlatePowerModel}, sys::PowerSystems.PowerSystem; args...)

    constructdevice!(m, netinjection, category, category_formulation, AbstractDCPowerModel, sys; args...)

end


text:
module IntroJuliaStatsDemo

# Actually using.
using Compat, Revise
# These are used in the examples, by having the module load (but not export) them, they will load faster later.
using DataFrames, Distributions, GLM, RCall, ForwardDiff

export  meaning_of_life

function __init__()
    Revise.track( @__FILE__ )
    ENV["EDITOR"] = "code"
end

# package code goes here
meaning_of_life() = 43

end # module


text:
"""
    ar_sum(x)
Reduce `k = length(x)` vectors through sum to a single vector.
"""
ar_sum(x) = sum(x)

"""
    ar_prod(x)
Reduce `k = length(x)` vectors through product to a single vector.
"""
ar_prod(x) = reduce((y, z) -> y .* z, x)

"""
    arithmetic_layer()
Generate the layer of arithmetic operations of the ICN. The operations are mutually exclusive, that is only one will be selected.
"""
function arithmetic_layer()
    arithmetics = LittleDict{Symbol, Function}(
        :sum => ar_sum,
        :prod => ar_prod,
    )

    return Layer(arithmetics, true)
end


text:
using WebIO, Widgets, JSExpr, Interact, AssetRegistry
using DigitalMusicology
using LightXML: XMLDocument

export pianorollwdg, veroviowdg, clear!, jumptonote!

notetodict(note::TimedNote, id=0) =
    Dict(:onset => onset(note),
         :offset => offset(note),
         :pitch => convert(Int,pitch(note)),
         :id => id)

dicttonote(dict) = TimedNote()

testnotes = [TimedNote(midip(0),0//1,1//1),
             TimedNote(midip(2),1//1,2//1),
             TimedNote(midip(4),2//1,3//1)]

jsdep(path...) = joinpath(@__DIR__, "..", "deps", path...)

"""
    pianorollwdg(notes; highlights=[])

Creates a widgets that shows `notes` in a pianoroll representation.
Optionally, groups of notes can be highlighted
by suppling an array of arrays of notes as `highlights`.
Notes and highlights can be changed and queried by accessing the
observables `:notes` and `:highlights` in the widget, e.g.

    pr[:notes][] = newnotes
"""
function pianorollwdg(notes; highlights=[], allowselect=false, zoomtohl=true)
    # https://d3js.org/d3.v4.min.js
    scp = Scope(imports=[joinpath(@__DIR__, "..", "deps", "d3.v4.min.js"),
                         joinpath(@__DIR__, "..", "deps", "pianoroll.js")])
    id = string("pianoroll", rand(UInt))

    onotes = Observable(scp, "onotes", notes)
    ohighlights = Observable{Vector{Vector}}(scp, "ohighlights", highlights)

    jnotes = map(onotes) do notes
        map(notetodict, notes, 1:length(notes))
    end
    setobservable!(scp, "jnotes", jnotes)
    jhighlights = map(ohighlights) do hls
        map(hl -> map(notetodict, hl), hls)
    end
    setobservable!(scp, "jhighlights", jhighlights)

    jselected = Observable(scp, "jselected", [])
    oselected = map(jselected) do sels
        map(sel -> notes[sel+1], sels)
    end

    jclear = Observable(scp, "jclear", nothing)

    importjs = @js function (d3)
        @var pr = document.getElementById($id)
        pianoroll($id, d3, $jnotes[], $jhighlights[], $allowselect)
        pr.updateSelectedOut = function (selected)
            $jselected[] = Array.from(selected)
        end
    end

    onimport(scp, importjs)

    onjs(scp["jnotes"], @js ns -> document.getElementById($id).updateNotes(ns))
    onjs(scp["jhighlights"],
         @js hls -> document.getElementById($id).updateHighlights(hls, $zoomtohl))
    onjs(scp["jclear"], @js clr -> document.getElementById($id).clearSel())

    lay(w) = node(:div,
        scope(w),
        id=id,
        style=Dict("background" => "white")
    )

    wdg = Widget{:pianoroll}([:notes => onotes, :highlights => ohighlights, :selected => oselected, :clear => jclear];
                             output=oselected, scope=scp, layout=lay)
end

clear!(pr::Widget{:pianoroll}) = pr[:clear][] = nothing

"""
    veroviowdg(input; highlights=[], allowselect=false, format="xml")

Takes as `input` a string in a format understood by verovio (e.g., MusicXML or MEI).
Returns a widget that renders the input using verovio and allows to
control the input as well as highlighted and selected notes using
the `:input`, `:highlights`, and `:selected` kewords.

Both highlights and selections use note ids,
so it makes sense to first set the `xml:id` attributes of notes in the input XML structure,
then feed the XML to this function and use the assigned note ids
to identify notes in the verovio SVG with notes used elsewhere.
While `:selection` is a list of note ids, `:highlights` is a list of note groups,
where each group is displayed in a different color.

The widget includes page controls, but the visible page cannot be directly controlled from Julia.
Instead, one can force the widget to show a specific note using `jumpto!(wdg, "note-id")`.

The selection can be cleared using `clear!(wdg)`.
"""
veroviowdg(xml::XMLDocument; args...) = veroviowdg(string(xml); args...)
function veroviowdg(input; highlights=[], allowselect=false, format="xml")
    # using the fix for emscripten-based dependencies:
    scp = Scope(imports=[jsdep("verovio-toolkit.js"),
                         jsdep("d3.v4.min.js"),
                         jsdep("veroviowdg.js")],
                systemjs_options=Dict(
                    "meta" =>
                    # don't hardcode hash, lookup in AssetRegistry if possible
                    Dict(AssetRegistry.getkey(jsdep("verovio-toolkit.js")) =>
                         Dict("format" => "global"))))
    id = string("verovioScore", rand(UInt))

    # observables
    oinput = Observable(scp, "input", input)
    ohighlights = Observable(scp, "highlights", highlights)
    oselected = Observable(scp, "selected", [])
    ojumpto = Observable(scp, "jumpto", "")

    # js setup
    onimport(scp, @js function(vrv, d3, wdg)
                 @var score = document.getElementById($id)
                 @var tk = @new vrv.verovio.toolkit()
                 console.log(wdg);
                 veroviowdg($id, tk, d3, $input, $highlights, $allowselect, $format)
                 #console.log("test2");
                 score.updateSelectedOut = function (selected)
                     $oselected[] = Array.from(selected)
                 end
                 score.updateHighlightsOut = function (highlights)
                     $ohighlights[] = Array.from(highlights);
                 end
             end)

    # handlers
    onjs(scp["input"], @js ns -> document.getElementById($id).updateInput(ns))
    onjs(scp["highlights"], @js hls -> document.getElementById($id).updateHighlightsIn(hls))
    onjs(scp["selected"], @js sel -> document.getElementById($id).updateSelectedIn(sel))
    onjs(scp["jumpto"], @js note -> document.getElementById($id).jumpToNote(note))

    # layout and widget
    lay(w) = node(:div, scope(w), id=id)
    wdg = Widget{:verovio}([:input => oinput,
                            :highlights => ohighlights,
                            :selected => oselected,
                            :jumpto => ojumpto];
                           scope=scp, layout=lay)
end

clear!(vrv::Widget{:verovio}) = vrv[:selected][] = []
jumptonote!(vrv::Widget{:verovio}, note) = vrv[:jumpto][] = note


text:
@with_kw mutable struct RMDP{S,A} <: MDP{S,A}
    amdp::MDP{S,A}
    π
    cost_fn = POMDPs.reward
    include_time_in_state = false
    dt = 0.1
    maxT = 100*dt
    disturbance_type=:arg #:arg if passed as argument, :noise if used as noise on the state
end

function POMDPs.initialstate(mdp::RMDP)
    s0 = initialstate(mdp.amdp)
    if mdp.include_time_in_state
        return ImplicitDistribution((rng) -> [0, rand(s0)...])
    else
        return s0
    end
end

get_s(mdp::RMDP, s) = mdp.include_time_in_state ? s[2:end] : s

POMDPs.actions(mdp::RMDP) = disturbances(mdp.amdp)
POMDPs.actionindex(mdp::RMDP, a) = disturbanceindex(mdp.amdp, a)

function POMDPs.isterminal(mdp::RMDP, s)
    isterm = isterminal(mdp.amdp, get_s(mdp, s))
    if mdp.include_time_in_state
        isterm = isterm || (s[1] > (mdp.maxT + mdp.dt/2))
    end
    isterm
end

function isfailure(mdp::RMDP, s)
    isfailure(mdp.amdp, get_s(mdp,s))
end

POMDPs.discount(mdp::RMDP) = 1f0

render(mdp::RMDP, s, a = nothing; kwargs...) = render(mdp.amdp, get_s(mdp,s), a; kwargs...)
render(mdp::RMDP; kwargs...) = render(mdp.amdp; kwargs...)

function POMDPs.gen(mdp::RMDP, s, x, rng::AbstractRNG = Random.GLOBAL_RNG; kwargs...)
    if mdp.disturbance_type == :arg
        sp, r = gen(mdp.amdp, get_s(mdp, s), action(mdp.π, get_s(mdp, s)), x, rng; kwargs...)
    elseif mdp.disturbance_type == :noise
        sp, r = gen(mdp.amdp, get_s(mdp, s), action(mdp.π, get_s(mdp, s) .+ x), rng; kwargs...)
    elseif mdp.disturbance_type == :both
        sp, r = gen(mdp.amdp, get_s(mdp, s), action(mdp.π, get_s(mdp, s) .+ x[2:end]), x[1], rng; kwargs...)
    else
        @error "Unrecognized disturbance type $(mdp.disturbance_type)"
    end

    if mdp.include_time_in_state
        t = s[1]
        sp = [t+mdp.dt, sp...]
    end
    (sp=sp, r=mdp.cost_fn(mdp, s, sp))
end

function POMDPs.transition(mdp::RMDP, s, x, rng::AbstractRNG = Random.GLOBAL_RNG; kwargs...)
    if mdp.disturbance_type == :arg
        t = transition(mdp.amdp, get_s(mdp, s), action(mdp.π, get_s(mdp, s)), x, rng; kwargs...)
    elseif mdp.disturbance_type == :noise
        t = transition(mdp.amdp, get_s(mdp, s), action(mdp.π, get_s(mdp, s) .+ x), rng; kwargs...)
    elseif mdp.disturbance_type == :both
        t = transition(mdp.amdp, get_s(mdp, s), action(mdp.π, get_s(mdp, s) .+ x[2:end]), x[1], rng; kwargs...)
    else
        @error "Unrecognized disturbance type $(mdp.disturbance_type)"
    end
    t
end

POMDPs.reward(mdp::RMDP, s, sp) = mdp.cost_fn(mdp.amdp, s, sp)


text:
module ComoniconOptions

export read_options, has_comonicon_toml, @asset_str

using TOML
using Configurations

include("types.jl")
include("utils.jl")

end


text:
@testset "action selectors" begin
    @testset "AlternateSelector" begin
        s = AlternateSelector(3)
        values = zeros(3)
        for _ in 1:3
            @test s(values) == 1
            @test s(values) == 2
            @test s(values) == 3
        end
    end

    @testset "EpsilonGreedySelector" begin
        @testset "test distribution" begin
            seed!(123)
            ϵ = 0.1
            s = EpsilonGreedySelector(ϵ)
            values = [1, 2, 0, -1]
            action_counts = countmap(s(values) for _ in 1:1000)
            @test isapprox(action_counts[1] / 1000, ϵ / length(values), atol=0.01)
            @test isapprox(action_counts[2] / 1000, 1 - ϵ + ϵ / length(values), atol=0.01)
            @test isapprox(action_counts[3] / 1000, ϵ / length(values), atol=0.01)
            @test isapprox(action_counts[4] / 1000, ϵ / length(values), atol=0.01)
        end

        # https://github.com/Ju-jl/ReinforcementLearningAnIntroduction.jl/issues/6
        @testset "test generator as input" begin
            seed!(123)
            ϵ = 0.1
            s = EpsilonGreedySelector(ϵ)
            values = (i for i in (1,2,0,-1))
            action_counts = countmap(s(values) for _ in 1:1000)
            @test isapprox(action_counts[1] / 1000, ϵ / length(values), atol=0.01)
            @test isapprox(action_counts[2] / 1000, 1 - ϵ + ϵ / length(values), atol=0.01)
            @test isapprox(action_counts[3] / 1000, ϵ / length(values), atol=0.01)
            @test isapprox(action_counts[4] / 1000, ϵ / length(values), atol=0.01)
        end
    end

    @testset "UpperConfidneceBound" begin
        s = UpperConfidenceBound(2)
        values = [1, 2]
        @test s(values) == 2
        @test s(values) == 1
        @test s(values) == 2
        @test s(values) == 2
        @test s(values) == 1
    end

    @testset "WeightedSample" begin
        seed!(123)
        s = WeightedSample()
        values = [0.2, 0.3, 0.5]
        action_counts = countmap(s(values) for _ in 1:1000)
        for i in 1:length(values)
            @test isapprox(action_counts[i]/1000, values[i], atol=0.1)
        end
    end
end


text:
include("../../LibSPP/librarySPP.jl")
include("./spp_dm2.jl")
include("../../../setSPP.jl")

function main()
    println("Etudiant(e)s : LE et OUSMAN")

    # Collecting the names of instances to solve located in the folder Data ----
    target = "../../Data"
    fnames = getfname(target)

    fres = splitdir(splitdir(pwd())[end-1])[end]
    io = open("../res/"*fres*".res", "w")
    for instance = 1:length(fnames)

        C, A = loadSPP(string(target,"/",fnames[instance]))

        resGRASP   = Array{Int64}(undef, 10)
        resReGRASP = Array{Int64}(undef, 10)

        for iter in 1:10
            t1 = @elapsed xGRASP, zGRASP = GRASP(C, A, 0.7, 10)
            t2 = @elapsed xReGRASP, zReGRASP = reactiveGRASP(C, A, [0.2,0.5,0.7,0.8,0.9,1], 120, 17)
            resGRASP[iter]   = zGRASP
            resReGRASP[iter] = zReGRASP
            #setSPP(C,A)
            println(io, fnames[instance], " | iter = ", iter, " | zGRASP = ", zGRASP, " | zReGRASP = ", zReGRASP, " | t1 = ", t1, " | t2 = ", t2)
        end
        println(io, "GRASP         : zMin = ", minimum(resGRASP), " | zMax = ", maximum(resGRASP), " | zMoy = ", sum(resGRASP[k] for k in 1:10)/10)
        println(io, "reactiveGRASP : zMin = ", minimum(resReGRASP), " | zMax = ", maximum(resReGRASP), " | zMoy = ", sum(resReGRASP[k] for k in 1:10)/10)
        println(io)
    end
    close(io)

end

main()

max_stars_repo_name: ducanh-le/metaSPPstu

max_issues_repo_name: ducanh-le/metaSPPstu

max_forks_repo_name: ducanh-le/metaSPPstu


text:
module VectorWrapper

using JavaCall

const Class = @jimport java.lang.Class
const JString = @jimport java.lang.String
const Object = @jimport java.lang.Object
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const StorelessUnivariateStatistic = @jimport org.hipparchus.stat.descriptive.StorelessUnivariateStatistic
const VectorialCovariance = @jimport org.hipparchus.stat.descriptive.vector.VectorialCovariance
const VectorialStorelessStatistic = @jimport org.hipparchus.stat.descriptive.vector.VectorialStorelessStatistic

export VectorialCovariance
export VectorialStorelessStatistic
export clear
export equals
export get_class
export get_dimension
export get_n
export get_result
export hash_code
export increment
export notify
export notify_all
export to_string
export wait

include("vectorial_covariance.jl")
include("vectorial_storeless_statistic.jl")

end


text:
Base.eltype(::D) where {T, D <: Union{DiscreteDomain{T}, ContinuousDomain{T}}} = T

function Base.convert(::Type{Intervals}, d::RangeDomain{T}) where {T <: Real}
    a, b = extrema(get_domain(d))
    return domain((Float64(a), true), (Float64(b), true))
end

function Base.convert(::Type{RangeDomain}, d::Intervals{T}) where {T <: Real}
    i = get_domain(d)[1]
    a = Int(value(i, :a) + (closed(i, :a) ? 0 : 1))
    b = Int(value(i, :b) - (closed(i, :b) ? 0 : 1))
    return domain(a:b)
end


text:
using DataFrames
using CSV
using Dates

# Compute the number of unique phenotypes in neighborhoods of increasing mutational size of a given genotype.
# Hypothesis:  Number unique phenotypes is postively related to complexity and
#   negatively related to robustness.

# Parallel map run_mutational_evolvability() over the goals of gl.  Output a dataframe csv file.
# Note that goals in gl may be replaced by random goals if no circuit is evolved to map to a goal.
function run_mutational_evolvability( gl::GoalList, p::Parameters, num_mutations::Integer, max_steps::Integer,  max_attempts::Integer;
    csvfile::String="" )
  gl_sav = deepcopy(gl)
  println("gl: ",gl)
  #complexity_robustness_geno_counts = map( g->mutational_evolvability( g, p, num_mutations, max_steps, max_attempts ), gl )
  complexity_robustness_geno_counts = pmap( g->mutational_evolvability( g, p, num_mutations, max_steps, max_attempts ), gl )
  println("length pairs: ",length(complexity_robustness_geno_counts))
  df = DataFrame()
  df.goal = map(x->x[1],complexity_robustness_geno_counts)
  df.complexity = map(x->x[2],complexity_robustness_geno_counts)
  df.robusness = map(x->x[3],complexity_robustness_geno_counts)
  for i = 1:num_mutations
    insertcols!(df,size(df)[2]+1,Symbol("geno_count$i")=>map(x->x[i+3],complexity_robustness_geno_counts))
  end
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      hostname = chomp(open("/etc/hostname") do f read(f,String) end)
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# num_mutations: ",num_mutations)
      println(f,"# max_steps: ",max_steps)
      println(f,"# max_attempts: ",max_attempts)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  df
end

# Run the below defined mutational_evolvability() for a single goal.
# Use neutral_evolution() to find the genotype that maps to the given goal (phenotype).
# If no circuit can be evolved that maps to g, then g is replace be a random goal (phenotype).
function mutational_evolvability( g::Goal, p::Parameters, num_mutations::Integer, max_steps::Integer,  max_attempts::Integer;
      max_goal_attempts::Int64=10)
  println("mut evolve g: ",g)
  nc = nothing
  goal_attempt = 0
  while goal_attempt < max_goal_attempts
    attempt = 0
    steps = max_steps
    while attempt < max_attempts && steps == max_steps
      c = random_chromosome(p,default_funcs(p.numinputs))
      attempt += 1
      (nc,steps) = neutral_evolution( c, g, max_steps )
      println("after neutral evolution: steps: ",steps)
    end
    if steps == max_steps
      g = randgoal(p.numinputs,p.numoutputs)
      println("new goal generated: ",g)
      goal_attempt += 1
    else
      break
    end
  end
  complexity = complexity5(nc)
  #println("output values(nc): ",output_values(nc),"  complexity: ",complexity)
  print_circuit(nc)
  (robustness,goal_sets) = mutational_evolvability(nc, num_mutations )
  Tuple(vcat([g,complexity,robustness],[length(goal_sets[end-j]) for j = (num_mutations-1):-1:0]))
end

function mutational_evolvability( c::Chromosome, num_mutations::Int64 )
  funcs = default_funcs( c.params.numinputs )
  robustness = -1.0   # establish scope
  g = output_values(c)
  goal_sets = Set{Goal}[]
  circuit_set = Set{Chromosome}[]
  goal_set = Set( Goal[] )
  for i = 1:num_mutations
    if i == 1
      goal_circuit_pairs = mutate_all( c, funcs, output_chromosomes=true )
      circuit_set = Set( map( x->x[2], goal_circuit_pairs ))
      goal_list =  map( x->x[1], goal_circuit_pairs )
      robustness = length(filter(x->x==g,goal_list))/length(goal_list)
      goal_set = Set( goal_list )
    else
      circuit_set_previous = deepcopy( circuit_set )
      circuit_set = Set{Chromosome}[]
      goal_set = Set{Goal}[]
      j = 0
      for c in circuit_set_previous
        j += 1
        goal_circuit_pairs = mutate_all( c, funcs, output_chromosomes=true )
        circuit_set = union( Set(map( x->x[2], goal_circuit_pairs )), circuit_set )
        goal_set = union( Set( map( x->x[1], goal_circuit_pairs ) ), goal_set )
        #println("j: ",j,"  lengths: cl: ",length(circuit_list),"  gs: ",length(goal_set))
      end
    end
    println("i: ",i,"  length(circuit_set): ",length(circuit_set))
    println("i: ",i,"  length(goal_set): ",length(goal_set))
    if i == 1
      push!(goal_sets,goal_set)
    else
      push!(goal_sets,union(goal_set,goal_sets[i-1]))
      println("i: ",i,"  goal: ",g,"  length(goal_sets[i]): ",length(goal_sets[i]))
    end
  end
  (robustness,goal_sets)
end


text:
struct VectorDiff{T1, T2} <: DeepDiff
    before::T1
    after::T2
    removed::Vector{Int}
    added::Vector{Int}
end

before(diff::VectorDiff) = diff.before
after(diff::VectorDiff) = diff.after
removed(diff::VectorDiff) = diff.removed
added(diff::VectorDiff) = diff.added
changed(diff::VectorDiff) = Int[]

Base.:(==)(d1::VectorDiff, d2::VectorDiff) = fieldequal(d1, d2)

# diffing an array is an application of the Longest Common Subsequence problem:
# https://en.wikipedia.org/wiki/Longest_common_subsequence_problem
function deepdiff(X::Vector, Y::Vector)
    # we're going to solve with dynamic programming, so let's first pre-allocate
    # our result array, which will store possible lengths of the common
    # substrings.

    lengths = zeros(Int, length(X)+1, length(Y)+1)

    for (j, v2) in enumerate(Y)
        for (i, v1) in enumerate(X)
            if isequal(v1, v2)
                lengths[i+1, j+1] = lengths[i, j] + 1
            else
                lengths[i+1, j+1] = max(lengths[i+1, j], lengths[i, j+1])
            end
        end
    end

    removed = Int[]
    added = Int[]
    backtrack(lengths, removed, added, X, Y, length(X), length(Y))

    VectorDiff(X, Y, removed, added)
end

# recursively trace back the longest common subsequence, adding items
# to the added and removed lists as we go
function backtrack(lengths, removed, added, X, Y, i, j)
    if i > 0 && j > 0 && isequal(X[i], Y[j])
        backtrack(lengths, removed, added, X, Y, i-1, j-1)
    elseif j > 0 && (i == 0 || lengths[i+1, j] ≥ lengths[i, j+1])
        backtrack(lengths, removed, added, X, Y, i, j-1)
        push!(added, j)
    elseif i > 0 && (j == 0 || lengths[i+1, j] < lengths[i, j+1])
        backtrack(lengths, removed, added, X, Y, i-1, j)
        push!(removed, i)
    end
end

# takes a function to be called for each item. The arguments given to the function
# are the items index, the state of the item (:removed, :added, :same) and a boolean
# for whether it's the last item. Indices are given for the `before` array when
# the state is :removed or :same, and for the `after` array when it's :added.
function visitall(f::Function, diff::VectorDiff)
    from = before(diff)
    to = after(diff)
    rem = removed(diff)
    add = added(diff)

    ifrom = 1
    ito = 1
    iremoved = 1
    iadded = 1

    while ifrom <= length(from) || ito <= length(to)
        if iremoved <= length(rem) && ifrom == rem[iremoved]
            ifrom += 1
            iremoved += 1
            f(ifrom-1, :removed, ifrom > length(from) && ito > length(to))
        elseif iadded <= length(add) && ito == add[iadded]
            ito += 1
            iadded += 1
            f(ito-1, :added, ifrom > length(from) && ito > length(to))
        else
            # not removed or added, must be in both
            ifrom += 1
            ito += 1
            f(ifrom-1, :same, ifrom > length(from) && ito > length(to))
        end
    end
end

function Base.show(io::IO, diff::VectorDiff)
    from = before(diff)
    to = after(diff)
    rem = removed(diff)
    add = added(diff)
    print(io, "[")

    visitall(diff) do idx, state, last
        if state == :removed
            printitem(io, from[idx], :red, "(-)")
            last || printstyled(io, ", ", color=:red)
        elseif state == :added
            printitem(io, to[idx], :green, "(+)")
            last || printstyled(io, ", ", color=:green)
        else
            printitem(io, from[idx])
            last || print(io, ", ")
        end
    end
    print(io, "]")
end

# prefix is printed if we're not using color
function printitem(io, v, color=:normal, prefix="")
    if hascolor(io)
        printstyled(io, v, color=color)
    else
        print(io, prefix, v)
    end
end

max_stars_repo_name: ssfrr/StructureDiffs

max_issues_repo_name: ssfrr/StructureDiffs

max_forks_repo_name: ssfrr/StructureDiffs


text:
module Granular

using JuLIP
using JuLIP.Potentials
using NeighbourLists
using StatsBase
using Hertz
import JuLIP: energy, forces, hessian_pos

export setup_cell

function setup_cell(N, ϕ, α; Lx = 1.0, Ly = 1.0, ratio=1.4, dim = 2)
    # create computational cell
    x = Lx * rand(N)'
    y = Ly * rand(N)'
    z = zeros(N)'
    X = [x; y; z]
    X = vecs(X)
    # create atoms
    # at = Atoms()
    # at = Atoms(X, zeros(X), ones(N), rand([1,2], N), _autocell(X),
    #         (false, false, false) )5 * eye(3), false)
    at = Atoms(:X, X)

    ρ = sqrt(2.0*ϕ*Lx*Ly/N/π/(1+ratio^2))

    idx = sample(1:N, Int(N/2); replace=false, ordered=true) # index for smaller discs
    at.Z = ones(N)
    at.Z[idx] = 2 # Here at.Z represents the tag of the grains,
                  # at.Z[i]=1, means the i-th grain is smaller
                  # at.Z[j]=2, means the j-th grain is bigger.

    M = ρ*(0.4*at.Z + 0.6); # Here, M does not mean masses but radius.
    at.M = M
    Rad = copy(M)
 # setup the basic information
    set_data!(at, :Num, N, Inf)
    set_data!(at, :ratio, ratio, Inf)
    set_data!(at, :Alpha, α, Inf)
    set_data!(at, :Phi, ϕ, Inf)
    set_data!(at, :Rad, Rad, Inf)

 # setup the shear type and initialize the strain
    StrainType = "SimpleShear"
    Strain = 0.0
    set_data!(at, :StrainType, StrainType, Inf)
    set_data!(at, :Strain, Strain, Inf)

 # initialize the pressure and energy
    set_data!(at, :Stress, zeros(dim, dim), Inf)
    set_data!(at, :Energy, 0.0, Inf)

 # setup the geometry information
    set_data!(at, :Length, Lx, Inf)
    set_data!(at, :Height, Ly, Inf)
    set_data!(at, :Dimension, dim, Inf)

    if StrainType == "SimpleShear"
        mycell = [Lx 0.0 0.0; Strain*Ly Ly 0.0; 0.0 0.0 1.0] # use mycell to differentiate this from JuLIP::cell
    elseif StrainType == "PureShear"
        mycell = [Lx 0.0 0.0; 0.0 Ly 0.0; 0.0 0.0 1.0]
    end
    set_cell!(at, mycell)

    # set_calculator!(at, V)
    if dim == 2
        set_constraint!(at, FixedCell2D(at))
    end
    set_pbc!(at, [true, true, true])

    return at
end

end

gran = Granular


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)

"""
    m = ElasticContactPairResponseMaterial(c_res, d_res, mu_k, mu_r, vsmall, wsmall)

Return an `ElasticContactPairResponseMaterial < AbstractContactPairMaterial` object.
"""
mutable struct ElasticContactPairResponseMaterial <: Modia3D.AbstractContactPairMaterial
    c_res::Float64
    c_geo::Float64
    n_geo::Float64
    d_res::Float64
    mu_k::Float64
    mu_r::Float64
    mu_r_geo::Float64
    vsmall::Float64
    wsmall::Float64
end

regularize(absv,vsmall) = absv >= vsmall ? absv : absv*(absv/vsmall)*(1.0 - (absv/vsmall)/3.0) + vsmall/3.0

resultantCoefficientOfRestitution(cor, abs_vreln, vsmall; cor_min=0.01) =
    cor + (cor_min - cor)*exp(log(0.01)*(abs_vreln/vsmall))

function resultantDampingCoefficient(cor, abs_vreln, vsmall; cor_min=0.01, d_max=1000.0)
    @assert(cor >= 0.0 && cor <= 1.0)
    @assert(abs_vreln >= 0.0)
    @assert(vsmall  > 0.0)
    @assert(cor_min > 0.0)
    @assert(d_max   > 0.0)

    cor_res = resultantCoefficientOfRestitution(cor,abs_vreln,vsmall;cor_min=cor_min)
    return min(d_max, 8.0*(1.0 - cor_res)/(5*cor_res*regularize(abs_vreln,vsmall)))
end

function elasticContactPairCoefficients(obj1::Object3D, obj2::Object3D)
  if typeof(obj1.feature.shape) == Modia3D.Shapes.Sphere && typeof(obj2.feature.shape) == Modia3D.Shapes.Sphere
    r1 = obj1.feature.shape.diameter*0.5
    r2 = obj2.feature.shape.diameter*0.5
    mu_r_geo = r1*r2/(r1 + r2)
    n_geo = 1.5
    c_geo = 4/3*sqrt(mu_r_geo)
  elseif typeof(obj1.feature.shape) == Modia3D.Shapes.Sphere && typeof(obj2.feature.shape) != Modia3D.Shapes.Sphere
    mu_r_geo = obj1.feature.shape.diameter*0.5
    n_geo = 1.5
    c_geo = 4/3*sqrt(mu_r_geo)
  elseif typeof(obj1.feature.shape) != Modia3D.Shapes.Sphere && typeof(obj2.feature.shape) == Modia3D.Shapes.Sphere
    mu_r_geo = obj2.feature.shape.diameter*0.5
    n_geo = 1.5
    c_geo = 4/3*sqrt(mu_r_geo)
  else
    mu_r_geo = 1.0
    n_geo = 1.0
    c_geo = 1.0
  end
  return (c_geo, n_geo, mu_r_geo)
end

"""
   resonseMaterial = contactStart(matPair::Modia3D.ElasticContactPairMaterial,
                                  obj1,obj2,rContact,contactNormal)

Return a `responseMaterial::ElasticContactPairResponseMaterial` object
at the start of a collision.
"""
function contactStart(matPair::Shapes.ElasticContactPairMaterial,
                      obj1::Object3D,
                      obj2::Object3D,
                      rContact::Frames.Vector3D,
                      contactNormal::Frames.Vector3D,
                      elasticContactReductionFactor::Float64)
    # Compute spring constant
    name1 = obj1.feature.contactMaterial
    name2 = obj2.feature.contactMaterial
    mat1  = Shapes.solidMaterialPalette[name1]
    mat2  = Shapes.solidMaterialPalette[name2]
    E1    = mat1.YoungsModulus
    E2    = mat2.YoungsModulus
    nu1   = mat1.PoissonsRatio
    nu2   = mat2.PoissonsRatio
    if E1 <= 0.0 || E2 <= 0.0 || nu1 <= 0.0 || nu1 >= 1.0 ||
       nu2 <= 0.0 || nu2 >= 1.0 || elasticContactReductionFactor <= 0.0
      responseMaterial = nothing
    else
      @assert(E1 > 0.0)
      @assert(E2 > 0.0)
      @assert(nu1 > 0.0 && nu1 < 1.0)
      @assert(nu2 > 0.0 && nu2 < 1.0)
      @assert(elasticContactReductionFactor > 0.0)
      c1 = E1/(1 - nu1^2)
      c2 = E2/(1 - nu2^2)
      c_res = elasticContactReductionFactor*c1*c2/(c1 + c2)

      # Compute damping constant
      delta_dot_start = normalRelativeVelocityAtContact(obj1, obj2, rContact, contactNormal)
      d_res = Modia3D.resultantDampingCoefficient(matPair.coefficientOfRestitution, abs(delta_dot_start), matPair.vsmall)

      # Determine other coefficients
      (c_geo, n_geo, mu_r_geo) = elasticContactPairCoefficients(obj1,obj2)
      responseMaterial = ElasticContactPairResponseMaterial(c_res, c_geo, n_geo, d_res,
                              matPair.slidingFrictionCoefficient,
                              matPair.rotationalResistanceCoefficient, mu_r_geo,
                              matPair.vsmall, matPair.wsmall)
    end
    return responseMaterial
end

contactEnd(mat::ElasticContactPairResponseMaterial,obj1,obj2)::Nothing = nothing
contactEnd(mat::Nothing,obj1,obj2)::Nothing = nothing

"""
    (f1,f2,t1,t2) = responseCalculation(material::ElasticContactPairResponseMaterial,
                            obj1, obj2, rContact, contactNormal, distanceWithHysteresis, time)

Compute contact forces `f1,f2` and contact torques `t1,t2` from the elastic contact pair
response material `material` the two penetrating objects `obj1, obj2`, the contact point `rContact`,
the contact normal `contactNormal` and the largest penetration depth `distanceWithHysteresis`
at time `time`.
"""
function responseCalculation(material::ElasticContactPairResponseMaterial, obj1::Object3D, obj2::Object3D,
                             rContact::Frames.Vector3D, e_n::Frames.Vector3D,
                             s::Float64, time::Float64, file, sim)::Tuple{Frames.Vector3D,Frames.Vector3D,Frames.Vector3D,Frames.Vector3D}
    # Material
    c_res    = material.c_res
    c_geo    = material.c_geo
    n_geo    = material.n_geo
    d_res    = material.d_res
    mu_k     = material.mu_k
    mu_r     = material.mu_r
    mu_r_geo = material.mu_r_geo
    vsmall   = material.vsmall
    wsmall   = material.wsmall

    ### signed velocity and relative velocity ####
    # Contact points and distances to local part frame (in world frame)
    r_rel1 = rContact - obj1.r_abs
    r_rel2 = rContact - obj2.r_abs

    # Velocities and angular velocities of contact frames in world frame
    w1 = obj1.R_abs'*obj1.w
    w2 = obj2.R_abs'*obj2.w
    v1 = obj1.v0 + cross(w1,r_rel1)
    v2 = obj2.v0 + cross(w2,r_rel2)

    # Velocities and angular velocities in normal and tangential direction
    w_rel   = w2 - w1
    e_w_reg = w_rel/Modia3D.regularize(norm(w_rel),wsmall)
    v_rel   = v2 - v1

    # delta_dot ... signed relative velocity in normal direction
    # v_t       ... relative velocity vector in tangential direction
    # delta     ... signed distance
    delta_dot = dot(v_rel,e_n)
    v_t       = v_rel - delta_dot*e_n
    e_t_reg   = v_t/Modia3D.regularize(norm(v_t),vsmall)
    delta     = -s

    if c_geo != 1.0 && n_geo != 1.0
      delta_comp = delta * sqrt(abs(delta))
    else
      delta_comp = delta
    end
    fn = -max(0.0, c_res * c_geo * delta_comp * (1 - d_res*delta_dot) )
    ft = -mu_k*fn*e_t_reg
    f1 = fn * e_n + ft
    f2 = -f1

    tau = -mu_r * mu_r_geo * fn * e_w_reg
    t1  = cross(r_rel1,f1) + tau
    t2  = cross(r_rel2,f2) - tau

    return (f1,f2,t1,t2)
end

responseCalculation(material::Nothing, obj1::Object3D, obj2::Object3D,
                    rContact::Frames.Vector3D, e_n::Frames.Vector3D,
                    s::Float64, time::Float64, file) =
                    (Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D)