text:
using Distributed

using Test

@everywhere begin
    using MLJBase
    using MLJTuning
    using ..Models
    import ComputationalResources: CPU1, CPUProcesses, CPUThreads
end

using Random
Random.seed!(1234*myid())
using .TestUtilities

N = 30
x1 = rand(N);
x2 = rand(N);
x3 = rand(N);
X = (; x1, x2, x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.4*rand(N);

m(K) = KNNRegressor(; K)
r = [m(K) for K in 13:-1:2]

# TODO: replace the above with the line below and post an issue on
# the failure (a bug in Distributed, I reckon):
# r = (m(K) for K in 13:-1:2)

@testset "constructor" begin
    @test_throws(MLJTuning.ERR_SPECIFY_RANGE,
                 TunedModel(model=first(r), tuning=Grid(), measure=rms))
    @test_throws(MLJTuning.ERR_SPECIFY_RANGE,
                 TunedModel(model=first(r), measure=rms))
    @test_throws(MLJTuning.ERR_BOTH_DISALLOWED,
                 TunedModel(model=first(r),
                            models=r, tuning=Explicit(), measure=rms))
    tm = @test_logs TunedModel(models=r, tuning=Explicit(), measure=rms)
    @test tm.tuning isa Explicit && tm.range ==r && tm.model == first(r)
    @test input_scitype(tm) == Unknown
    @test TunedModel(models=r, measure=rms) == tm
    @test_logs (:info, r"No measure") @test TunedModel(models=r) == tm

    @test_throws(MLJTuning.ERR_SPECIFY_MODEL,
                 TunedModel(range=r, measure=rms))
    @test_throws(MLJTuning.ERR_MODEL_TYPE,
                 TunedModel(model=42, tuning=Grid(),
                            range=r, measure=rms))
    @test_logs (:info, MLJTuning.INFO_MODEL_IGNORED) tm =
        TunedModel(model=42, tuning=Explicit(), range=r, measure=rms)
    @test_logs (:info, r"No measure") tm =
        TunedModel(model=first(r), range=r)
    @test_throws(MLJTuning.ERR_SPECIFY_RANGE_OR_MODELS,
                TunedModel(tuning=Explicit(), measure=rms))
    @test_throws(MLJTuning.ERR_NEED_EXPLICIT,
                 TunedModel(models=r, tuning=Grid()))
    @test_logs TunedModel(first(r), range=r, measure=rms)
    @test_logs(
        (:warn, MLJTuning.warn_double_spec(first(r), last(r))),
        TunedModel(first(r), model=last(r), range=r, measure=rms),
    )
    @test_throws(
        MLJTuning.ERR_TOO_MANY_ARGUMENTS,
        TunedModel(first(r), last(r), range=r, measure=rms),
    )
    tm = @test_logs TunedModel(model=first(r), range=r, measure=rms)
    @test tm.tuning isa RandomSearch
    @test input_scitype(tm) == Table(Continuous)

    @test_throws MLJTuning.ERR_UNINSTANTIATED_MODEL TunedModel(; model=KNNRegressor, range=r)
end

results = [(evaluate(model, X, y,
                     resampling=CV(nfolds=2),
                     measure=rms,
                     verbosity=0,)).measurement[1] for model in r]

@testset "measure compatibility check" begin
    tm = TunedModel(
        models=r,
        resampling=CV(nfolds=2),
        measures=cross_entropy
    )
    @test_logs((:error, r"Problem"),
               (:info, r""),
               (:info, r""),
               @test_throws ArgumentError fit(tm, 0, X, y))
end

@testset_accelerated "basic fit (CPU1)" accel begin
    printstyled("\n Testing progressmeter basic fit with $(accel) and CPU1 resampling \n", color=:bold)
    best_index = argmin(results)
    tm = TunedModel(
        models=r,
        resampling=CV(nfolds=2),
        measures=[rms, l1],
        acceleration=accel
    )
    verbosity = accel isa CPU1 ? 2 : 1
    fitresult, meta_state, _report = fit(tm, verbosity, X, y);
    history, _, state = meta_state;
    results2 = map(event -> event.measurement[1], history)
    @test results2 ≈ results
    @test fitresult.model == collect(r)[best_index]
    @test _report.best_model == collect(r)[best_index]
    @test _report.history[5] == MLJTuning.delete(history[5], :metadata)

    # training_losses:
    losses = training_losses(tm, _report)
    @test all(eachindex(losses)) do i
        minimum(results[1:i]) == losses[i]
    end
    @test MLJBase.iteration_parameter(tm) == :n

end

@static if VERSION >= v"1.3.0-DEV.573"
@testset_accelerated "Basic fit (CPUThreads)" accel begin
    printstyled("\n Testing progressmeter basic fit with $(accel) and CPUThreads resampling \n", color=:bold)
    tm = TunedModel(
        models=r,
        resampling=CV(nfolds=2),
        measures=[rms, l1],
        acceleration= CPUThreads(),
        acceleration_resampling=accel
    )
    fitresult, meta_state, report = fit(tm, 1, X, y);
    history, _, state = meta_state;
    results3 = map(event -> event.measurement[1], history)
    @test results3 ≈ results
end
end
@testset_accelerated "Basic fit (CPUProcesses)" accel begin
    printstyled("\n Testing progressmeter basic fit with $(accel) and CPUProcesses resampling \n", color=:bold)
    best_index = argmin(results)
    tm = TunedModel(
        models=r,
        resampling=CV(nfolds=2),
        measures=[rms, l1],
        acceleration=CPUProcesses(),
        acceleration_resampling=accel
    )
    fitresult, meta_state, report = fit(tm, 1, X, y);
    history, _, state = meta_state;
    results4 = map(event -> event.measurement[1], history)
    @test results4 ≈ results
end

@testset_accelerated(
    "under/over supply of models",
    accel,
    begin
    tm = TunedModel(
        models=r,
        measures=[rms, l1],
        acceleration=accel,
        resampling=CV(nfolds=2),
        n=4
    )
    mach = machine(tm, X, y)
    fit!(mach, verbosity=0)
    history = MLJBase.report(mach).history
    @test map(event -> event.measurement[1], history) ≈ results[1:4]

    tm.n += 2
    @test_logs((:info, r"^Updating"),
               (:info, r"^Attempting to add 2.*to 6"),
               fit!(mach, verbosity=1))
    history = MLJBase.report(mach).history
    @test map(event -> event.measurement[1], history) ≈ results[1:6]

    tm.n=100
    @test_logs (:info, r"Only 12") fit!(mach, verbosity=0)
    history = MLJBase.report(mach).history
    @test map(event -> event.measurement[1], history) ≈ results
    end
)

@everywhere begin

    # variation of the Explicit strategy that annotates the models
    # with metadata
    mutable struct MockExplicit <: MLJTuning.TuningStrategy end

    annotate(model) = (model, params(model)[1])

    _length(x) = length(x)
    _length(::Nothing) = 0
    MLJTuning.models(
        tuning::MockExplicit,
        model,
        history,
        state,
        n_remaining,
        verbosity
    ) = annotate.(state)[_length(history) + 1:end], state

    function default_n(tuning::Explicit, range)
        try
            length(range)
        catch MethodError
            DEFAULT_N
        end

    end
end

@testset_accelerated(
    "passing of model metadata",
    accel,
    begin
    tm = TunedModel(
        model=first(r),
        tuning=MockExplicit(),
        range=r,
        resampling=CV(nfolds=2),
        measures=[rms, l1],
        acceleration=accel
    )
    fitresult, meta_state, report = fit(tm, 0, X, y);
    history, _, state = meta_state;
    @test all(history) do event
    event.metadata == event.model.K
    end
    end
)

@testset "data caching" begin
    X = (x1= ones(5),);
    y = coerce(collect("abcaa"), Multiclass);

    m(b) = ConstantClassifier(testing=true, bogus=b)
    r = [m(b) for b in 1:5]

    tuned_model = TunedModel(
        models=r,
        resampling=Holdout(fraction_train=0.8),
        measure=log_loss,
        cache=true
    )

    # There is reformatting and resampling of X and y for training of
    # first model, and then resampling on the prediction side only
    # thereafter. Then, finally, for training on all supplied data, there
    # is reformatting and resampling again:
    @test_logs(
        (:info, "reformatting X, y"), # fit 1
        (:info, "resampling X, y"),   # fit 1
        (:info, "resampling X"),      # predict 1
        (:info, "resampling X"),      # predict 2
        (:info, "resampling X"),      # predict 3
        (:info, "resampling X"),      # predict 4
        (:info, "resampling X"),      # predict 5
        (:info, "reformatting X, y"), # fit on all data
        (:info, "resampling X, y"),   # fit one all data
        fit(tuned_model, 0, X, y)
    );

    # Otherwise, resampling and reformatting happen for every model
    # trained, and we get a reformat on the prediction side every
    # time:
    tuned_model.cache = false
    @test_logs(
        (:info, "reformatting X, y"), # fit 1
        (:info, "resampling X, y"),   # fit 1
        (:info, "reformatting X"),    # predict 1
        (:info, "reformatting X, y"), # fit 2
        (:info, "resampling X, y"),   # fit 2
        (:info, "reformatting X"),    # predict 2
        (:info, "reformatting X, y"), # fit 3
        (:info, "resampling X, y"),   # fit 3
        (:info, "reformatting X"),    # predict 3
        (:info, "reformatting X, y"), # fit 4
        (:info, "resampling X, y"),   # fit 4
        (:info, "reformatting X"),    # predict 4
        (:info, "reformatting X, y"), # fit 5
        (:info, "resampling X, y"),   # fit 5
        (:info, "reformatting X"),    # predict 5
        (:info, "reformatting X, y"), # fit on all data
        (:info, "resampling X, y"),   # fit on all data
        fit(tuned_model, 0, X, y)
    );
end

@testset "issue #128" begin
    X, y = make_regression(10, 2)
    dtc = DecisionTreeRegressor()
    r   = range(dtc, :max_depth, lower=1, upper=50);

    tmodel = TunedModel(model=dtc, ranges=[r, ],
                        tuning=Grid(resolution=50),
                        measure=mae,
                        n=48);
    mach = machine(tmodel, X, y)
    fit!(mach, verbosity=0);

    @test length(report(mach).history) == 48

    tmodel.n = 49
    fit!(mach, verbosity=0);

    @test length(report(mach).history) == 49
end

@testset_accelerated "Resampling reproducibility" accel begin
    X, y = make_regression(100, 2)
    dcr = DeterministicConstantRegressor()

    # Hold out reproducibility
    homodel = TunedModel(
        models=fill(dcr, 10),
        resampling=Holdout(rng=StableRNG(1234)),
        acceleration_resampling=accel,
        measure=mae
    )
    homach = machine(homodel, X, y)
    fit!(homach, verbosity=0);
    horep = report(homach)
    measurements = getproperty.(horep.history, :measurement)
    @test all(==(measurements[1]), measurements)

    # Cross-validation reproducibility
    cvmodel = TunedModel(
        models=fill(dcr, 10),
        resampling=CV(nfolds=5, rng=StableRNG(1234)),
        acceleration_resampling=accel,
        measure=mae
    )
    cvmach = machine(cvmodel, X, y)
    fit!(cvmach, verbosity=0);
    cvrep = report(cvmach)
    per_folds = getproperty.(cvrep.history, :per_fold)
    @test all(==(per_folds[1]), per_folds)
end

@testset "deterministic metrics for probabilistic models" begin

    # https://github.com/JuliaAI/MLJBase.jl/pull/599 allows mix of
    # deterministic and probabilistic metrics:
    X, y = MLJBase.make_blobs()
    model = DecisionTreeClassifier()
    range = MLJBase.range(model, :max_depth, values=[1,2])
    tmodel = TunedModel(model=model,
                        tuning=Grid(),
                        range=range,
                        measures=[MisclassificationRate(),
                                  LogLoss()])
    mach = machine(tmodel, X, y)
    @test_logs fit!(mach, verbosity=0)

end

@testset_accelerated "weights and class_weights are being passed" accel begin
    # we'll be tuning using 50/50 holdout
    X = (x=fill(1.0, 6),)
    y = coerce(["a", "a", "b", "a", "a", "b"], OrderedFactor)
    w = [1.0, 1.0, 100.0, 1.0, 1.0, 100.0]
    class_w = Dict("a" => 2.0, "b" => 100.0)

    model = DecisionTreeClassifier()

    # the first supports weights, the second class weights:
    ms=[MisclassificationRate(), MulticlassFScore()]

    resampling=Holdout(fraction_train=0.5)

    # without weights:
    tmodel = TunedModel(
        resampling=resampling,
        models=fill(model, 5),
        measures=ms,
        acceleration=accel
    )
    mach = machine(tmodel, X, y)
    fit!(mach, verbosity=0)
    measurement = report(mach).best_history_entry.measurement
    e = evaluate(model, X, y, measures=ms, resampling=resampling)
    @test measurement == e.measurement

    # with weights:
    tmodel.weights = w
    tmodel.class_weights = class_w
    fit!(mach, verbosity=0)
    measurement_weighted = report(mach).best_history_entry.measurement
    e_weighted = evaluate(model, X, y;
                          measures=ms,
                          resampling=resampling,
                          weights=w,
                          class_weights=class_w,
                          verbosity=-1)
    @test measurement_weighted == e_weighted.measurement

    # check both measures are different when they are weighted:
    @test !any(measurement .== measurement_weighted)
end

@testset "data caching at outer level suppressed" begin
    X, y = make_blobs()
    model = DecisionTreeClassifier()
    tmodel = TunedModel(models=[model,])
    mach = machine(tmodel, X, y)
    @test mach isa Machine{<:Any,false}
    fit!(mach, verbosity=-1)
    @test !isdefined(mach, :data)
    MLJBase.Tables.istable(mach.cache[end].fitresult.machine.data[1])
end

true

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
@testset "FakeAbsOrdQuoRing" begin
  G = small_group(8, 3)
  QG = QQ[G]
  ZG = Order(QG, basis(QG))

  for L in [ZG, maximal_order(ZG)]
    x = rand(L, -1:1)
    while iszero(det(representation_matrix(x)))
      x = rand(L, -1:1)
    end
    I = x * L
    R = @inferred Hecke.FakeAbsOrdQuoRing(L, I)
    @test (@inferred isone(@inferred one(R)))
    o = one(R)
    for i in 1:100
      a = @inferred rand(L, -1:1)
      aR = @inferred R(a)
      @test lift(aR) - a in I
      b = @inferred rand(L, -1:1)
      bR = @inferred R(b)
      @test lift(bR) - b in I
      @test o * aR == aR
      @test aR * o == aR
      @test aR * bR == R(a * b)
    end
  end
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 3b4135be-6682-11eb-1766-e772017ddc9b
begin
	using Plots
	using SpecialFunctions
	using LaTeXStrings
end

# ╔═╡ 88a17f66-667b-11eb-1d4d-57d100ff452e
md"
# Question:
Fluid bounded by two parallel plates extended to infinity. The lower wall is suddenly accelerated in the x-direction while the upper wall is stationary.

Governing Equation:

$\frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial y^2}$

Initial condition:\
t = 0; u = 0     for 0 < y ≤ h \
t = 0; u = $U_0$ for y = 0

Boundary condition: \
t ≥ 0; u = $U_0$ for y = 0 \
t ≥ 0; u = 0 for y = h

"

# ╔═╡ 291e4152-667e-11eb-3d13-031ad9572916
# Code for solving the problem
begin
	ν = 0.000217
	h = 0.04
	u_b = 40.0
	n = 41
    δx = h/(n-1)

	function makegrid(h,n,δx)
		for i in 1:n
			x[i] = 0.0 + (i-1)*δx
		end
		return x
	end

	function init(u,n)
		u      = zeros(n)
		u[1]   = u_b
		return u
	end

	function exact_sol(x,n,t_eval)
		@info(t_eval)
		uex = zeros(n)
		for i in 1:n
			η = x[i]/(2*sqrt(ν*t_eval))
			η1= h/(2*sqrt(ν*t_eval))
			uex[i] = u_b* (erfc(η) - erfc(2*η1 - η) + erfc(2*η1 + η)  - erfc(4*η1 - η) + erfc(4*η1 + η) )
		end
		return uex
	end

	#FTCS Time stepping
	function FTCS(δt,tmax,δx,u,n,ν,u_b)
		u_new = zeros(n)

	    #BC's
		u_new[1] = u_b
		u_new[n] = 0.0
		#Update
		t = 0.0
		nt= floor(Int,tmax/δt)
		for i in 1:nt
			for j in 2:n-1
				u_new[j] = u[j] + ν*(δt/(δx)^2) * (u[j+1] - 2*u[j] + u[j-1])
			end
			t = t + δt
			u = copy(u_new)
		end
		return u
	end

	#Dufort-Frankel Timestepping
	function DFF(δt,tmax,δx,u,n,ν,u_b)
		u_new = zeros(n)
		u_old = zeros(n)

	    #BC's
		u_new[1] = u_b
		u_new[n] = 0.0
		#Update
		t = 0.0
		nt= floor(Int,tmax/δt)
		d = ν*(δt/(δx)^2)

		u_old = copy(u)
		for i in 1:nt
			for j in 2:n-1
				u_new[j] = ((1.0-2.0*d)*u_old[j] + 2*d*(u[j+1] + u[j-1]))/(1.0+2.0*d)
			end
			t = t + δt
			u_old = copy(u)
			u = copy(u_new)
		end
		return u
	end

	function TriD(A,B,C,D,n)
		#/// n is only "n-2" in global
		a = copy(A)
		b = copy(B)
		c = copy(C)
		d = copy(D)
		u = zeros(n)
		for i in 2:n
			w = a[i]/b[i-1]
			b[i] = b[i] - w*c[i-1]
			d[i] = d[i] - w*d[i-1]
		end

		#Back solve
		u[n] = d[n]/b[n]
		for i in n-1:-1:1
			u[i] = (d[i] - c[i]*u[i+1])/b[i]
		end
		#@info(u)
		return u
	end

	#Laason implicit solve
	function LAimpl(δt,tmax,δx,u,n,ν,u_b)
		u_new = zeros(n)
		u_old = zeros(n)

	    #BC's
		u_new[1] = u_b
		u_new[n] = 0.0
		#Update
		t = 0.0
		nt= floor(Int,tmax/δt)
		d = ν*(δt/(δx)^2)

		A = zeros(n-2)
		B = zeros(n-2)
		C = zeros(n-2)
		D = zeros(n-2)
		for i in 1:nt
			A  .= -d
			B  .= 1+2*d
			C  .= -d
			D  .= copy(u[2:n-1])
			D[1]  = D[1]   + d*u_b
			D[n-2]= D[n-2] + d*0.0

			u_new[2:n-1] = TriD(A,B,C,D,n-2)

			u_new[1] = u_b
			u_new[n] = 0.0
			u = copy(u_new)
		end
		return u
	end

	#Crank-Nicolson implicit solve
	function CNimpl(δt,tmax,δx,u,n,ν,u_b)
		u_new = zeros(n)
		u_old = zeros(n)

	    #BC's
		u_new[1] = u_b
		u_new[n] = 0.0

		u[1] = u_b
		u[n] = 0.0
		#Update
		t = 0.0
		nt= floor(Int,tmax/δt)
		d = ν*(δt/(δx)^2)

		A = zeros(n-2)
		B = zeros(n-2)
		C = zeros(n-2)
		D = zeros(n-2)
		for i in 1:nt
			A  .= -0.5*d
			B  .= 1+1*d
			C  .= -0.5*d
			D  .= copy(u[2:n-1]) + 0.5*d*(u[3:n]-2.0*u[2:n-1]+u[1:n-2])

			D[1]  = D[1]   + 0.5*d*u_b
			D[n-2]= D[n-2] + 0.5*d*0.0

			u_new[2:n-1] = TriD(A,B,C,D,n-2)

			u_new[1] = u_b
			u_new[n] = 0.0
			u = copy(u_new)
		end
		return u
	end

	x = zeros(n)
	u = zeros(n)
	uex= zeros(n)

	#Run necessary stuff
	x = makegrid(h,n,δx)
	u = init(u,n)

end

# ╔═╡ 2e306e02-6719-11eb-2755-e364e8d201cc
md"# Velocity Profiles across schemes"

# ╔═╡ 270a7404-6713-11eb-2cf1-8b17120d4405
md"# Velocity Profiles for FTCS"

# ╔═╡ 1a86f16c-6713-11eb-2f1e-53e78d8e6a15
begin
	#Velocity Profiles for FTCS explicit method
	δt = 0.002
	tmax= [0.0,0.18,0.36,0.54,0.72,0.90,1.08]

	u_1 = FTCS(δt,tmax[1],δx,u,n,ν,u_b)
	u_2 = FTCS(δt,tmax[2],δx,u,n,ν,u_b)
	u_3 = FTCS(δt,tmax[3],δx,u,n,ν,u_b)
	u_4 = FTCS(δt,tmax[4],δx,u,n,ν,u_b)
	u_5 = FTCS(δt,tmax[5],δx,u,n,ν,u_b)
	u_6 = FTCS(δt,tmax[6],δx,u,n,ν,u_b)
	u_7 = FTCS(δt,tmax[7],δx,u,n,ν,u_b)

	gr()
	plot(title="Velocity Profile for dy = 0.001, dt = 0.002, d = 0.434")
	plot!(xlabel="u",ylabel="y")
	plot!(u_1,x,markershape = :ltriangle,label="t=0.0sec")
	plot!(u_2,x,markershape = :ltriangle,label="t=0.18sec")
	plot!(u_3,x,markershape = :ltriangle,label="t=0.36sec")
	plot!(u_4,x,markershape = :ltriangle,label="t=0.54sec")
	plot!(u_5,x,markershape = :ltriangle,label="t=0.72sec")
	plot!(u_6,x,markershape = :ltriangle,label="t=0.90sec")
	plot!(u_7,x,markershape = :ltriangle,label="t=1.08sec")

end

# ╔═╡ fb155b98-6715-11eb-0231-2bc9fde178fa
begin
	#Velocity Profiles for FTCS explicit method
	δt1 = 0.00232
	tmax1= [0.0,0.21,0.42,0.63,0.84,1.04,1.25]

	u1_1 = FTCS(δt1,tmax1[1],δx,u,n,ν,u_b)
	u1_2 = FTCS(δt1,tmax1[2],δx,u,n,ν,u_b)
	u1_3 = FTCS(δt1,tmax1[3],δx,u,n,ν,u_b)
	u1_4 = FTCS(δt1,tmax1[4],δx,u,n,ν,u_b)
	u1_5 = FTCS(δt1,tmax1[5],δx,u,n,ν,u_b)
	u1_6 = FTCS(δt1,tmax1[6],δx,u,n,ν,u_b)
	u1_7 = FTCS(δt1,tmax1[7],δx,u,n,ν,u_b)

	gr()
	plot(title="Velocity Profile for dy = 0.001, dt = 0.00232, d = 0.50344")
	plot!(xlabel="u",ylabel="y")
	plot!(u1_1,x,markershape = :ltriangle,label="t=0.0sec")
	plot!(u1_2,x,markershape = :ltriangle,label="t=0.21sec")
	plot!(u1_3,x,markershape = :ltriangle,label="t=0.42sec")
	plot!(u1_4,x,markershape = :ltriangle,label="t=0.63sec")
	plot!(u1_5,x,markershape = :ltriangle,label="t=0.84sec")
	plot!(u1_6,x,markershape = :ltriangle,label="t=1.04sec")
	plot!(u1_7,x,markershape = :ltriangle,label="t=1.25sec")

end

# ╔═╡ c3e4e836-6716-11eb-1fe8-4b01a8f19efe
md"# Velocity Profiles for Dufort-Frankel"

# ╔═╡ 800e1682-6716-11eb-2968-1f922c7f834c
begin
	#Velocity Profiles for DFF explicit method
	δt2 = 0.002
	tmax2= [0.0,0.18,0.36,0.54,0.72,0.90,1.08]

	u2_1 = DFF(δt2,tmax2[1],δx,u,n,ν,u_b)
	u2_2 = DFF(δt2,tmax2[2],δx,u,n,ν,u_b)
	u2_3 = DFF(δt2,tmax2[3],δx,u,n,ν,u_b)
	u2_4 = DFF(δt2,tmax2[4],δx,u,n,ν,u_b)
	u2_5 = DFF(δt2,tmax2[5],δx,u,n,ν,u_b)
	u2_6 = DFF(δt2,tmax2[6],δx,u,n,ν,u_b)
	u2_7 = DFF(δt2,tmax2[7],δx,u,n,ν,u_b)

	gr()
	plot(title="Velocity Profile for dy = 0.001, dt = 0.002, d = 0.434")
	plot!(xlabel="u",ylabel="y")
	plot!(u2_1,x,markershape = :ltriangle,label="t=0.0sec")
	plot!(u2_2,x,markershape = :ltriangle,label="t=0.18sec")
	plot!(u2_3,x,markershape = :ltriangle,label="t=0.36sec")
	plot!(u2_4,x,markershape = :ltriangle,label="t=0.54sec")
	plot!(u2_5,x,markershape = :ltriangle,label="t=0.72sec")
	plot!(u2_6,x,markershape = :ltriangle,label="t=0.90sec")
	plot!(u2_7,x,markershape = :ltriangle,label="t=1.08sec")

end

# ╔═╡ 39b664cc-6717-11eb-11b1-b74350c4341d
begin
	#Velocity Profiles for FTCS explicit method
	δt3 = 0.003
	tmax3= [0.0,0.21,0.42,0.63,0.84,1.04,1.25]

	u3_1 = DFF(δt3,tmax3[1],δx,u,n,ν,u_b)
	u3_2 = DFF(δt3,tmax3[2],δx,u,n,ν,u_b)
	u3_3 = DFF(δt3,tmax3[3],δx,u,n,ν,u_b)
	u3_4 = DFF(δt3,tmax3[4],δx,u,n,ν,u_b)
	u3_5 = DFF(δt3,tmax3[5],δx,u,n,ν,u_b)
	u3_6 = DFF(δt3,tmax3[6],δx,u,n,ν,u_b)
	u3_7 = DFF(δt3,tmax3[7],δx,u,n,ν,u_b)

	gr()
	plot(title="Velocity Profile for dy = 0.001, dt = 0.003, d = 0.651")
	plot!(xlabel="u",ylabel="y")
	plot!(u3_1,x,markershape = :ltriangle,label="t=0.0sec")
	plot!(u3_2,x,markershape = :ltriangle,label="t=0.21sec")
	plot!(u3_3,x,markershape = :ltriangle,label="t=0.42sec")
	plot!(u3_4,x,markershape = :ltriangle,label="t=0.63sec")
	plot!(u3_5,x,markershape = :ltriangle,label="t=0.84sec")
	plot!(u3_6,x,markershape = :ltriangle,label="t=1.04sec")
	plot!(u3_7,x,markershape = :ltriangle,label="t=1.25sec")

end

# ╔═╡ 99b7da0e-6717-11eb-061d-ef8e525fa7a6
md"# Velocity Profiles for Laason Implicit"

# ╔═╡ c2e5d336-6717-11eb-04f4-d9506627a208
begin
	#Velocity Profiles for LAimplicit method
	δt4 = 0.002
	tmax4= [0.0,0.18,0.36,0.54,0.72,0.90,1.08]

	u4_1 = LAimpl(δt4,tmax4[1],δx,u,n,ν,u_b)
	u4_2 = LAimpl(δt4,tmax4[2],δx,u,n,ν,u_b)
	u4_3 = LAimpl(δt4,tmax4[3],δx,u,n,ν,u_b)
	u4_4 = LAimpl(δt4,tmax4[4],δx,u,n,ν,u_b)
	u4_5 = LAimpl(δt4,tmax4[5],δx,u,n,ν,u_b)
	u4_6 = LAimpl(δt4,tmax4[6],δx,u,n,ν,u_b)
	u4_7 = LAimpl(δt4,tmax4[7],δx,u,n,ν,u_b)

	gr()
	plot(title="Velocity Profile for dy = 0.001, dt = 0.002, d = 0.434")
	plot!(xlabel="u",ylabel="y")
	plot!(u4_1,x,markershape = :ltriangle,label="t=0.0sec")
	plot!(u4_2,x,markershape = :ltriangle,label="t=0.18sec")
	plot!(u4_3,x,markershape = :ltriangle,label="t=0.36sec")
	plot!(u4_4,x,markershape = :ltriangle,label="t=0.54sec")
	plot!(u4_5,x,markershape = :ltriangle,label="t=0.72sec")
	plot!(u4_6,x,markershape = :ltriangle,label="t=0.90sec")
	plot!(u4_7,x,markershape = :ltriangle,label="t=1.08sec")
		plot!(exact_sol(x,n,0.18),x)
end

# ╔═╡ 0ce28b94-6718-11eb-1793-c5f06bfe5686
begin
	#Velocity Profiles for LaasonImplicit method
	δt5 = 0.01
	tmax5= [0.0,0.21,0.42,0.63,0.84,1.04,1.25]

	u5_1 = LAimpl(δt5,tmax5[1],δx,u,n,ν,u_b)
	u5_2 = LAimpl(δt5,tmax5[2],δx,u,n,ν,u_b)
	u5_3 = LAimpl(δt5,tmax5[3],δx,u,n,ν,u_b)
	u5_4 = LAimpl(δt5,tmax5[4],δx,u,n,ν,u_b)
	u5_5 = LAimpl(δt5,tmax5[5],δx,u,n,ν,u_b)
	u5_6 = LAimpl(δt5,tmax5[6],δx,u,n,ν,u_b)
	u5_7 = LAimpl(δt5,tmax5[7],δx,u,n,ν,u_b)

	gr()
	plot(title="Velocity Profile for dy = 0.001, dt = 0.01, d = 2.17")
	plot!(xlabel="u",ylabel="y")
	plot!(u3_1,x,markershape = :ltriangle,label="t=0.0sec")
	plot!(u3_2,x,markershape = :ltriangle,label="t=0.21sec")
	plot!(u3_3,x,markershape = :ltriangle,label="t=0.42sec")
	plot!(u3_4,x,markershape = :ltriangle,label="t=0.63sec")
	plot!(u3_5,x,markershape = :ltriangle,label="t=0.84sec")
	plot!(u3_6,x,markershape = :ltriangle,label="t=1.04sec")
	plot!(u3_7,x,markershape = :ltriangle,label="t=1.25sec")

end

# ╔═╡ 1eb83680-6719-11eb-16c6-89de2999c61c
md"# Velocity Profiles for Crank-Nicolson Implicit"

# ╔═╡ 1d27d028-6719-11eb-3ac6-7147ea8109ef
begin
	#Velocity Profiles for CNimplicit method
	δt6 = 0.002
	tmax6= [0.0,0.18,0.36,0.54,0.72,0.90,1.08]

	u6_1 = CNimpl(δt6,tmax6[1],δx,u,n,ν,u_b)
	u6_2 = CNimpl(δt6,tmax6[2],δx,u,n,ν,u_b)
	u6_3 = CNimpl(δt6,tmax6[3],δx,u,n,ν,u_b)
	u6_4 = CNimpl(δt6,tmax6[4],δx,u,n,ν,u_b)
	u6_5 = CNimpl(δt6,tmax6[5],δx,u,n,ν,u_b)
	u6_6 = CNimpl(δt6,tmax6[6],δx,u,n,ν,u_b)
	u6_7 = CNimpl(δt6,tmax6[7],δx,u,n,ν,u_b)

	gr()
	plot(title="Velocity Profile for dy = 0.001, dt = 0.002, d = 0.434")
	plot!(xlabel="u",ylabel="y")
	plot!(u6_1,x,markershape = :ltriangle,label="t=0.0sec")
	plot!(u6_2,x,markershape = :ltriangle,label="t=0.18sec")
	plot!(u6_3,x,markershape = :ltriangle,label="t=0.36sec")
	plot!(u6_4,x,markershape = :ltriangle,label="t=0.54sec")
	plot!(u6_5,x,markershape = :ltriangle,label="t=0.72sec")
	plot!(u6_6,x,markershape = :ltriangle,label="t=0.90sec")
	plot!(u6_7,x,markershape = :ltriangle,label="t=1.08sec")
    plot!(exact_sol(x,n,0.18),x)
end

# ╔═╡ 04589ef0-671a-11eb-3b04-4dc8b7a0fd43
begin
	#Velocity Profiles for CNImplicit method
	δt7 = 0.01
	tmax7= [0.0,0.21,0.42,0.63,0.84,1.04,1.25]

	u7_1 = CNimpl(δt7,tmax5[1],δx,u,n,ν,u_b)
	u7_2 = CNimpl(δt7,tmax5[2],δx,u,n,ν,u_b)
	u7_3 = CNimpl(δt7,tmax5[3],δx,u,n,ν,u_b)
	u7_4 = CNimpl(δt7,tmax5[4],δx,u,n,ν,u_b)
	u7_5 = CNimpl(δt7,tmax5[5],δx,u,n,ν,u_b)
	u7_6 = CNimpl(δt7,tmax5[6],δx,u,n,ν,u_b)
	u7_7 = CNimpl(δt7,tmax5[7],δx,u,n,ν,u_b)

	gr()
	plot(title="Velocity Profile for dy = 0.001, dt = 0.01, d = 2.17")
	plot!(xlabel="u",ylabel="y")
	plot!(u7_1,x,markershape = :ltriangle,label="t=0.0sec")
	plot!(u7_2,x,markershape = :ltriangle,label="t=0.21sec")
	plot!(u7_3,x,markershape = :ltriangle,label="t=0.42sec")
	plot!(u7_4,x,markershape = :ltriangle,label="t=0.63sec")
	plot!(u7_5,x,markershape = :ltriangle,label="t=0.84sec")
	plot!(u7_6,x,markershape = :ltriangle,label="t=1.04sec")
	plot!(u7_7,x,markershape = :ltriangle,label="t=1.25sec")
	plot!(exact_sol(x,n,0.21),x)

end

# ╔═╡ 08bc0ee8-671d-11eb-1e49-e392f826eba5
md"# Error Distribution for t = 0.18 sec, 1.08sec"

# ╔═╡ bb3fd91e-671d-11eb-1c84-d518b84da067
#Calculating Error
function err_calc(uexact,ucalc,n)
	err = zeros(n)
	for i in 1:n
		err[i] = (uexact[i] - ucalc[i])
	end
	return err
end

# ╔═╡ 24c371dc-671d-11eb-2611-c39d01c2ecf2
begin
	u_ex = exact_sol(x,n,0.18)
	u_ftcs = FTCS(0.002,0.18,δx,u,n,ν,u_b)
	u_dff  = DFF(0.002,0.18,δx,u,n,ν,u_b)
	u_la   = LAimpl(0.002,0.18,δx,u,n,ν,u_b)
	u_cn   = CNimpl(0.002,0.18,δx,u,n,ν,u_b)

	ftcs_err = err_calc(u_ex,u_ftcs,n)
	dff_err  = err_calc(u_ex,u_dff,n)
	la_err  = err_calc(u_ex,u_la,n)
	cn_err  = err_calc(u_ex,u_cn,n)

	plot(title="Error comparison for various schmes dt = 0.002 at t=0.18")
	plot!(xlabel="Error",ylabel="y",xlim=(-0.11,0.11))
	plot!(ftcs_err,x,markershape = :ltriangle,label="FTCS")
	plot!(dff_err,x,markershape = :ltriangle,label="Dufort-Frankel")
	plot!(la_err,x,markershape = :ltriangle,label="Laasonen-implicit")
	plot!(cn_err,x,markershape = :ltriangle,label="Crank-Nicolson-implicit")

end

# ╔═╡ 615146cc-672c-11eb-097d-41f289d3f896
begin
	u1_ex = exact_sol(x,n,1.08)
	u1_ftcs = FTCS(0.002,1.08,δx,u,n,ν,u_b)
	u1_dff  = DFF(0.002,1.08,δx,u,n,ν,u_b)
	u1_la   = LAimpl(0.002,1.08,δx,u,n,ν,u_b)
	u1_cn   = CNimpl(0.002,1.08,δx,u,n,ν,u_b)

	ftcs1_err = err_calc(u1_ex,u1_ftcs,n)
	dff1_err  = err_calc(u1_ex,u1_dff,n)
	la1_err  = err_calc(u1_ex,u1_la,n)
	cn1_err  = err_calc(u1_ex,u1_cn,n)

	plot(title="Error comparison for various schmes dt = 0.002 and t=1.08")
	plot!(xlabel="Error",ylabel="y",xlim=(-0.05,0.05))
	plot!(ftcs1_err,x,markershape = :ltriangle,label="FTCS")
	plot!(dff1_err,x,markershape = :ltriangle,label="Dufort-Frankel")
	plot!(la1_err,x,markershape = :ltriangle,label="Laasonen-implicit")
	plot!(cn1_err,x,markershape = :ltriangle,label="Crank-Nicolson-implicit")

end

# ╔═╡ dc70983e-672d-11eb-101b-47d6c05f02e2
md"# Error Distribution for different time-steps at t = 1.0sec for Lassonen"

# ╔═╡ 251d0ac8-672d-11eb-1d64-8daa66f59ffa
begin
	u2_ex = exact_sol(x,n,1.0)
	ula_1 = LAimpl(0.005,1.0,δx,u,n,ν,u_b)
	ula_2 = LAimpl(0.01,1.0,δx,u,n,ν,u_b)
	ula_3 = LAimpl(0.1,1.0,δx,u,n,ν,u_b)
	ula_4 = LAimpl(0.2,1.0,δx,u,n,ν,u_b)

	ula1_err  = err_calc(u2_ex,ula_1,n)
	ula2_err  = err_calc(u2_ex,ula_2,n)
	ula3_err  = err_calc(u2_ex,ula_3,n)
	ula4_err  = err_calc(u2_ex,ula_4,n)

	plot(title="Error comparison for different timesteps at t=1.0
		for Lassonen")
	plot!(xlabel="Error",ylabel="y")
	plot!(ula1_err,x,markershape = :ltriangle,label="dt=0.005")
	plot!(ula2_err,x,markershape = :ltriangle,label="dt=0.01")
	plot!(ula3_err,x,markershape = :ltriangle,label="dt=0.1")
	plot!(ula4_err,x,markershape = :ltriangle,label="dt=0.2")

end

# ╔═╡ Cell order:
# ╟─88a17f66-667b-11eb-1d4d-57d100ff452e
# ╠═3b4135be-6682-11eb-1766-e772017ddc9b
# ╠═291e4152-667e-11eb-3d13-031ad9572916
# ╟─2e306e02-6719-11eb-2755-e364e8d201cc
# ╟─270a7404-6713-11eb-2cf1-8b17120d4405
# ╟─1a86f16c-6713-11eb-2f1e-53e78d8e6a15
# ╟─fb155b98-6715-11eb-0231-2bc9fde178fa
# ╟─c3e4e836-6716-11eb-1fe8-4b01a8f19efe
# ╟─800e1682-6716-11eb-2968-1f922c7f834c
# ╟─39b664cc-6717-11eb-11b1-b74350c4341d
# ╟─99b7da0e-6717-11eb-061d-ef8e525fa7a6
# ╟─c2e5d336-6717-11eb-04f4-d9506627a208
# ╟─0ce28b94-6718-11eb-1793-c5f06bfe5686
# ╟─1eb83680-6719-11eb-16c6-89de2999c61c
# ╟─1d27d028-6719-11eb-3ac6-7147ea8109ef
# ╟─04589ef0-671a-11eb-3b04-4dc8b7a0fd43
# ╟─08bc0ee8-671d-11eb-1e49-e392f826eba5
# ╠═bb3fd91e-671d-11eb-1c84-d518b84da067
# ╟─24c371dc-671d-11eb-2611-c39d01c2ecf2
# ╟─615146cc-672c-11eb-097d-41f289d3f896
# ╟─dc70983e-672d-11eb-101b-47d6c05f02e2
# ╟─251d0ac8-672d-11eb-1d64-8daa66f59ffa


text:
@testset "Sample Parents" begin

    expObj = Exponential(0.5)
    kernelTest(x) = pdf.(expObj, x)
    bgFunction(x) = 0.5 * x

    @testset "Constant Background Rate" begin
        inferTests && @inferred HawkesProcesses.sample_parents([1,2,3], 0.5, 0.5, kernelTest)

        sp = HawkesProcesses.sample_parents([1,2,3], 0.5, 0.5, kernelTest)
        @test length(sp) == 4
        @test length(sp[1]) == 3
    end

    @testset "Function Background Rate" begin

    end

end


text:
#
# shm.jl --
#
# Management of shared memory for Julia.  Both POSIX and System V shared memory
# models are supported.
#
#------------------------------------------------------------------------------
#
# This file is part of IPC.jl released under the MIT "expat" license.
# Copyright (C) 2016-2018, Éric Thiébaut (https://github.com/emmt/IPC.jl).
#

"""
```julia
SharedMemory(id, len; perms=0o600, volatile=true)
```

yields a new shared memory object identified by `id` and whose size is `len`
bytes.  The identifier `id` can be a string starting by a `'/'` to create a
POSIX shared memory object or a System V IPC key to create a System V shared
memory segment.  In this latter case, the key can be `IPC.PRIVATE` to
automatically create a non-existing shared memory segment.

Keyword `perms` can be used to specify which access permissions are granted.
By default, only reading and writing by the user is granted.

Keyword `volatile` can be used to specify whether the shared memory is volatile
or not.  If non-volatile, the shared memory will remain accessible until
explicit destruction or system reboot.  By default, the shared memory is
destroyed when no longer in use.

To retrieve an existing shared memory object, call:

```julia
SharedMemory(id; readonly=false)
```

where `id` is the shared memory identifier (a string, an IPC key or a System V
IPC identifier of shared memory segment as returned by `ShmId`).  Keyword
`readonly` can be set true if only read access is needed.  Note that method
`shmid(obj)` may be called to retrieve the identifier of the shared memory
object `obj`.

Some methods are extended for shared memory objects.  Assuming `shm` is an
instance of `SharedMemory`, then:

```julia
pointer(shm)    # yields the base address of the shared memory
sizeof(shm)     # yields the number of bytes of the shared memory
shmid(shm)      # yields the identifier the shared memory
```

To ensure that shared memory object `shm` is eventually destroyed, call:

```julia
rm(shm)
```

See also [`shmid`](@ref), [`shmrm`](@ref).

"""
function SharedMemory(key::Key, len::Integer;
                      perms::Integer = S_IRUSR|S_IWUSR,
                      volatile::Bool = true)::SharedMemory{ShmId}
    # Create a new System V shared memory segment with given size and, at
    # least, read and write access for the caller.
    len ≥ 1 || throw_argument_error("bad number of bytes (", len, ")")
    flags = maskmode(perms) | (S_IRUSR|S_IWUSR|IPC_CREAT|IPC_EXCL)
    id = _shmget(key.value, len, flags)
    if id < 0
        throw_system_error("shmget")
    end

    # Attach shared memory segment to process address space.
    ptr = _shmat(id, C_NULL, 0)
    if ptr == Ptr{Cvoid}(-1)
        errno = Libc.errno()
        _shmctl(id, IPC_RMID, C_NULL)
        throw_system_error("shmat", errno)
    end
    if volatile && _shmctl(id, IPC_RMID, C_NULL) == -1
        errno = Libc.errno()
        _shmdt(ptr)
        throw_system_error("shmctl", errno)
    end

    # Instanciate Julia object.
    return SharedMemory{ShmId}(ptr, len, volatile, ShmId(id))
end

SharedMemory(key::Key; readonly::Bool = false, kwds...) =
    SharedMemory(ShmId(key, readonly); kwds...)

function SharedMemory(id::ShmId; readonly::Bool=false)::SharedMemory{ShmId}
    len = sizeof(id)
    ptr = shmat(id, readonly)
    return SharedMemory{ShmId}(ptr, len, false, id)
end

# Create a new POSIX shared memory object.
function SharedMemory(name::AbstractString,
                      len::Integer;
                      perms::Integer = S_IRUSR | S_IWUSR,
                      volatile::Bool = true) :: SharedMemory{String}
    # Make sure owner has read and write permissions (otherwise setting the
    # size will fail).
    mode = maskmode(perms) | (S_IRUSR | S_IWUSR)
    flags = O_CREAT | O_EXCL | O_RDWR
    return SharedMemory(name, flags, mode, len, volatile)
end

# Map an existing POSIX shared memory object.
function SharedMemory(name::AbstractString;
                      readonly::Bool=false) :: SharedMemory{String}
    flags = (readonly ? O_RDONLY : O_RDWR)
    mode = (readonly ? S_IRUSR : S_IRUSR|S_IWUSR)
    return SharedMemory(name, flags, mode, 0, false)
end

function SharedMemory(name::AbstractString,
                      flags::Integer,
                      mode::Integer,
                      len::Integer = 0,
                      volatile::Bool = false)

    # Create a new POSIX shared memory object?
    create = ((flags & O_CREAT) != 0)
    if create
        len ≥ 1 || throw_argument_error("bad number of bytes (", len, ")")
    end

    # Open shared memory and set or get its size.
    fd = _shm_open(name, flags, mode)
    if fd == -1
        throw_system_error("shm_open")
    end
    local nbytes::Int = 0
    if create
        # Set the size of the new shared memory object.
        if _ftruncate(fd, len) != SUCCESS
            errno = Libc.errno()
            _close(fd)
            _shm_unlink(name)
            throw_system_error("ftruncate", errno)
        end
        nbytes = Int(len)
    else
        # Get the size of the existing shared memory object.
        try
            nbytes = Int(filesize(fd))
        catch err
            _close(fd)
            rethrow(err)
        end
    end

    # Map the shared memory.  Note that `prot = PROT_NONE` should never occur.
    access = flags & (O_RDONLY|O_WRONLY|O_RDWR)
    prot = (access == O_RDONLY ? PROT_READ :
            access == O_WRONLY ? PROT_WRITE :
            access == O_RDWR   ? PROT_READ|PROT_WRITE : PROT_NONE)
    ptr = _mmap(C_NULL, nbytes, prot, MAP_SHARED, fd, 0)
    if ptr == MAP_FAILED
        errno = Libc.errno()
        _close(fd)
        if create
            _shm_unlink(name)
        end
        throw_system_error("mmap", errno)
    end

    # File descriptor can be closed.
    if _close(fd) != SUCCESS
        errno = Libc.errno()
        if create
            _shm_unlink(name)
        end
        _munmap(ptr, len)
        throw_system_error("close", errno)
    end

    # Return the shared memory object.
    return SharedMemory{String}(ptr, nbytes, volatile, String(name))
end

function _destroy(obj::SharedMemory{String})
    if obj.volatile
        _shm_unlink(obj.id)
    end
    _munmap(obj.ptr, obj.len)
    if PARANOID
        obj.ptr = C_NULL
        obj.len = 0
        obj.volatile = false
        obj.id = ""
    end
end

function _destroy(obj::SharedMemory{ShmId})
    _shmdt(obj.ptr)
    if PARANOID
        obj.ptr = C_NULL
        obj.len = 0
        obj.volatile = false
        obj.id = ShmId(-1)
    end
end

Base.sizeof(obj::SharedMemory) = obj.len
Base.pointer(obj::SharedMemory) = obj.ptr

# The short version of `show` if also used for string interpolation in
# scripts so that it is not necessary to extend methods:
#     Base.convert(::Type{String}, obj::T)
#     Base.string(obj::T)

Base.show(io::IO, obj::SharedMemory{String}) =
    print(io, "SharedMemory(\"", obj.id, "\"; len=", obj.len,
          ", ptr=Ptr{Cvoid}(0x", string(convert(Int, obj.ptr), base=16),"))")

Base.show(io::IO, obj::SharedMemory{ShmId}) =
    print(io, "SharedMemory(", obj.id, "; len=", obj.len,
          ", ptr=Ptr{Cvoid}(0x", string(convert(Int, obj.ptr), base=16),"))")

Base.show(io::IO, ::MIME"text/plain", obj::SharedMemory) = show(io, obj)

"""
```julia
shmid(arg)
```

yield the identifier of an existing POSIX shared memory object or Sytem V
shared memory segment identifed by `arg` or associated with `arg`.  Argument
can be:

* An instance of `SharedMemory`.

* An instance of `WrappedArray` whose contents is stored into shared memory.

* A string starting with a `'/'` (and no other `'/'`) to identify a POSIX
  shared memory object.

* An instance of `IPC.Key` to specify a System V IPC key associated with a
  shared memory segment.  In that case, an optional second argument `readonly`
  can be set `true` to only request read-only access; otherwise read-write
  access is requested.

* An instance of `ShmId` to specify a System V shared memory segment.

See also: [`SharedMemory`](@ref), [`shmrm`](@ref).

"""
shmid(shm::SharedMemory) = shm.id
shmid(arr::WrappedArray{T,N,<:SharedMemory}) where {T,N} = shmid(arr.mem)
shmid(id::ShmId) = id
shmid(key::Key, args...) = ShmId(key, args...)

Base.rm(shm::SharedMemory) = shmrm(shm)
Base.rm(id::ShmId) = shmrm(id)
Base.rm(::Type{SharedMemory}, key::Key) = shmrm(key)
Base.rm(::Type{SharedMemory}, name::AbstractString) = shmrm(name)

"""

```julia
shmrm(arg)
```

removes the shared memory associated with `arg`.  If `arg` is a name, the
corresponding POSIX named shared memory is unlinked.  If `arg` is a key or
identifier of a BSD shared memory segment, the segment is marked to be
eventually destroyed.  Argument `arg` can also be a `SharedMemory` object.

The `rm` method may also be called to remove an existing shared memory segment
or object.  There are several possibilities:

```julia
rm(SharedMemory, name)  # `name` identifies a POSIX shared memory object
rm(SharedMemory, key)   # `key` is associated with a BSD shared memory segment
rm(id)                  # `id` is the identifier of a BSD shared memory segment
rm(shm)                 # `shm` is an instance of `SharedMemory`
```

See also: [`SharedMemory`](@ref), [`shmid`](@ref), [`shmat`](@ref).

"""
shmrm(shm::SharedMemory) = shmrm(shmid(shm))
shmrm(key::Key) = shmrm(shmid(key))

shmrm(name::AbstractString) = begin
    if _shm_unlink(name) != SUCCESS
        errno = Libc.errno()
        if errno != Libc.ENOENT
            throw_system_error("shm_unlink", errno)
        end
    end
end

shmrm(id::ShmId) = begin
    if _shmrm(name) != SUCCESS
        # Only throw an error if not an already removed shared memory segment.
        errno = Libc.errno()
        if errno != Libc.EIDRM
            throw_system_error("shmctl", errno)
        end
    end
end

#------------------------------------------------------------------------------
# System V Shared Memory
#
# Methods which directly call C functions are prefixed by an underscore and do
# not perform any checking of their arguments (except that type of arguments
# must be correct).
#
# Higher level versions are not prefixed by an underscore and throw a
# `SystemError` if the value returned by the C function indicates an error.
# Their arguments may be slightly different.
#

Base.convert(::Type{Cint}, id::ShmId) = id.value
Base.convert(::Type{T}, id::ShmId) where {T<:Integer} = convert(T, id.value)

Base.show(io::IO, id::ShmId) = print(io, "ShmId(", id.value, ")")
Base.show(io::IO, ::MIME"text/plain", id::ShmId) = show(io, id)

function Base.sizeof(id::ShmId)
    # Note it is faster to create a Julia array of small size than calling
    # malloc/free or using a tuple unless the number of elements of the tuple
    # is small.
    buf = _workspace(_sizeof_struct_shmid_ds)
    shmctl(id, IPC_STAT, buf)
    return _peek(_typeof_shm_segsz, buf, _offsetof_shm_segsz)
end

"`_workspace(size)` yields an array of `size` uninitialized bytes."
@inline _workspace(size::Integer) = Array{UInt8}(undef, size)

"""
# Get the identifier of an existing System V shared memory segment

The following calls:

```julia
ShmId(id)                  -> id
ShmId(arr)                 -> id
ShmId(key, readlony=false) -> id
```

yield the the identifier of the existing System V shared memory segment
associated with the value of the first argument.  `id` is the identifier of the
shared memory segment, `arr` is an array attached to a System V shared memory
segment and `key` is the key associated with the shared memory segment.  In
that latter case, `readlony` can be set `true` to only request read-only
access; otherwise read-write access is requested.

See also: [`shmid`](@ref), [`shmget`](@ref).

"""
ShmId(id::ShmId) = id
ShmId(shm::SharedMemory{ShmId}) = shmid(shm)
ShmId(arr::WrappedArray{T,N,SharedMemory{ShmId}}) where {T,N} = shmid(arr)
ShmId(key::Key, readonly::Bool = false) =
    shmget(key, 0, (readonly ? S_IRUSR : (S_IRUSR|S_IWUSR)))

"""
# Get or create a System V shared memory segment

The call:

```julia
shmget(key, siz, flg) -> id
```

yields the identifier of the shared memory segment associated with the value of
the argument `key`.  A new shared memory segment, with size equal to the value
of `siz` (possibly rounded up to a multiple of the memory page size
`IPC.PAGE_SIZE`), is created if `key` has the value `IPC.PRIVATE` or `key`
isn't `IPC.PRIVATE`, no shared memory segment corresponding to `key` exists,
and `IPC_CREAT` is specified in argument `flg`.

Arguments are:

* `key` is the System V IPC key associated with the shared memory segment.

* `siz` specifies the size (in bytes) of the shared memory segment (may be
  rounded up to multiple of the memory page size).

* `flg` is a bitwise combination of flags.  The least significant 9 bits
  specify the permissions granted to the owner, group, and others.  These bits
  have the same format, and the same meaning, as the mode argument of `chmod`.
  Bit `IPC_CREAT` can be set to create a new segment.  If this flag is not
  used, then `shmget` will find the segment associated with `key` and check to
  see if the user has permission to access the segment.  Bit `IPC_EXCL` can be
  set in addition to `IPC_CREAT` to ensure that this call creates the segment.
  If `IPC_EXCL` and `IPC_CREAT` are both set, the call will fail if the segment
  already exists.

"""
function shmget(key::Key, siz::Integer, flg::Integer)
    id = _shmget(key.value, siz, flg)
    systemerror("shmget", id < 0)
    return ShmId(id)
end

_shmget(key::Integer, siz::Integer, flg::Integer) =
    ccall(:shmget, Cint, (_typeof_key_t, Csize_t, Cint), key, siz, flg)

"""
```julia
shmat(id, readonly) -> ptr
```

attaches a shared memory segment to the address space of the caller.  Argument
`id` is the identifier of the shared memory segment.  Boolean argument
`readonly` specifies whether to attach the segment for read-only access;
otherwise, the segment is attached for read and write accesses and the process
must have read and write permissions for the segment.  The returned value is
the pointer to access the shared memory segment.

See also: [`shmdt`](@ref), [`shmrm`](@ref).

"""
function shmat(id::ShmId, readonly::Bool)
    flg = (readonly ? SHM_RDONLY : zero(SHM_RDONLY))
    ptr = _shmat(id.value, C_NULL, flg)
    systemerror("shmat", ptr == Ptr{Cvoid}(-1))
    return ptr
end

_shmat(id::Integer, ptr::Ptr{Cvoid}, flg::Integer) =
    ccall(:shmat, Ptr{Cvoid}, (Cint, Ptr{Cvoid}, Cint), id, ptr, flg)

"""
```julia
shmdt(ptr)
```

detaches a System V shared memory segment from the address space of the caller.
Argument `ptr` is the pointer returned by a previous `shmat()` call.

See also: [`shmdt`](@ref), [`shmget`](@ref).

"""
shmdt(ptr::Ptr{Cvoid}) = systemerror("shmdt", _shmdt(ptr) != SUCCESS)

_shmdt(ptr::Ptr{Cvoid}) =
    ccall(:shmdt, Cint, (Ptr{Cvoid},), ptr)

"""

```julia
shmctl(id, cmd, buf)
```

performs the control operation specified by `cmd` on the System V shared memory
segment whose identifier is given in `id`.  The `buf` argument is a pointer to
a `shmid_ds` C structure.

See also [`shminfo`](@ref), [`shmcfg`](@ref) and [`shmrm`](@ref).

"""
shmctl(id::ShmId, cmd::Integer, buf::Union{DenseVector,Ptr}) =
    systemerror("shmctl", _shmctl(id.value, cmd, buf) == -1)

_shmctl(id::Integer, cmd::Integer, buf::Union{DenseVector,Ptr}) =
    ccall(:shmctl, Cint, (Cint, Cint, Ptr{Cvoid}), id, cmd, buf)

_shmrm(id::Integer) = _shmctl(id, IPC_RMID, C_NULL)

"""

```julia
shmcfg(arg, perms) -> id
```

changes the access permissions of a System V IPC shared memory segment.
Argument `perms` specifies bitwise flags with the new permissions.  The first
argument `arg` can be the identifier of the shared memory segment, a shared
array attached to the shared memory segment or the System V IPC key associated
with the shared memory segment.  In all cases, the identifier of the shared
memory segment is returned.

See also [`ShmId`](@ref), [`shmget`](@ref), [`shmctl`](@ref) and
[`SharedMemory`](@ref).

"""
function shmcfg(id::ShmId, perms::_typeof_shm_perm_mode)
    mask = convert(_typeof_shm_perm_mode, MASKMODE)
    buf = _workspace(_sizeof_struct_shmid_ds)
    shmctl(id, IPC_STAT, buf)
    mode = _peek(_typeof_shm_perm_mode, buf, _offsetof_shm_perm_mode)
    if (mode & mask) != (perms & mask)
        _poke!(_typeof_shm_perm_mode, buf, _offsetof_shm_perm_mode,
               (mode & ~mask) | (perms & mask))
        shmctl(id, IPC_SET, buf)
    end
    return id
end

shmcfg(id::ShmId, perms::Integer) =
    shmcfg(id, convert(_typeof_shm_perm_mode, perms))

shmcfg(arg::Union{ShmArray,Key}, perms::Integer) =
    shmcfg(ShmId(arg), perms)

"""

`ShmInfo` is the structure used to store information about a System V shared memory segment.
The call `ShmInfo(arg)` is equivalent to [`shminfo(arg)`](@ref shminfo).

"""
ShmInfo(arg::Union{ShmId,ShmArray,Key}) = shminfo(arg)

"""

```julia
shminfo(arg) -> info
```

yields information about the System V shared memory segment identified or
associated with `arg` which can be the identifier of the shared memory segment,
a shared array attached to the shared memory segment or the System V IPC key
associated with the shared memory segment.

See also [`ShmInfo`](@ref), [`ShmId`](@ref), [`shmget`](@ref), [`shmat`](@ref),
[`SharedMemory`](@ref).

"""
shminfo(arg::Union{ShmId,ShmArray,Key}) = shminfo!(arg, ShmInfo())

"""

```julia
shminfo!(arg, info) -> info
```

overwrites `info` (an instance of [`ShmInfo`](@ref)) with the information about
the System V shared memory segment identified or associated with `arg`.  See
[`shminfo`](@ref) for more details.

"""
function shminfo!(id::ShmId, info::ShmInfo)
    buf = _workspace(_sizeof_struct_shmid_ds)
    shmctl(id, IPC_STAT, buf)
    info.atime  = _peek(_typeof_time_t,        buf, _offsetof_shm_atime)
    info.dtime  = _peek(_typeof_time_t,        buf, _offsetof_shm_dtime)
    info.ctime  = _peek(_typeof_time_t,        buf, _offsetof_shm_ctime)
    info.segsz  = _peek(_typeof_shm_segsz,     buf, _offsetof_shm_segsz)
    info.id     = id.value
    info.cpid   = _peek(_typeof_pid_t,         buf, _offsetof_shm_cpid)
    info.lpid   = _peek(_typeof_pid_t,         buf, _offsetof_shm_lpid)
    info.nattch = _peek(_typeof_shmatt_t,      buf, _offsetof_shm_nattch)
    info.mode   = _peek(_typeof_shm_perm_mode, buf, _offsetof_shm_perm_mode)
    info.uid    = _peek(_typeof_uid_t,         buf, _offsetof_shm_perm_uid)
    info.gid    = _peek(_typeof_gid_t,         buf, _offsetof_shm_perm_gid)
    info.cuid   = _peek(_typeof_uid_t,         buf, _offsetof_shm_perm_cuid)
    info.cgid   = _peek(_typeof_gid_t,         buf, _offsetof_shm_perm_cgid)
    return info
end

shminfo!(arg, info::ShmInfo) = shminfo!(ShmId(arg), info)

shminfo!(key::Key, info::ShmInfo) = shminfo!(ShmId(key, true), info)


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2019 Varnerlab
# Robert Frederick School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

function contains(string,token)
    return occursin(token,string)
end

function eye(number_of_rows,number_of_cols)
    return Matrix{Float64}(I,number_of_rows, number_of_cols)
end

# load the files -
searchdir(path,key) = filter(x->contains(x,key),readdir(path))

# Function to set some values for the weights, gene lengths etc -
function customize_data_dictionary(default_data_dictionary::Dict{String,Any}, host_type::Symbol)::Dict{String,Any}

  # make a deepcopy -
  customized_data_dictionary = deepcopy(default_data_dictionary)

  # set the gene lengths -
  gene_length_array = customized_data_dictionary["gene_coding_length_array"]
  gene_length_array[1] = 744.0    # clssrA
  gene_length_array[2] = 711.0    # deGFPssrA
  gene_length_array[3] = 720.0    # S28
  gene_length_array[4] = 720.0    # S70 => doesn't matter, gene expression = 0

  # set the protein lengths -
  protein_coding_length_array = customized_data_dictionary["protein_coding_length_array"]
  protein_coding_length_array[1] = 248.0  # clssrA
  protein_coding_length_array[2] = 237.0  # deGFPssrA
  protein_coding_length_array[3] = 240.0  # S28
  protein_coding_length_array[4] = 240.0  # S70 => doesn't matter, translation = 0

  # set ICs -
  initial_condition_array = customized_data_dictionary["initial_condition_array"]
  initial_condition_array[1] = 0.001  # clssrA gene
  initial_condition_array[2] = 0.008  # deGFP gene
  initial_condition_array[3] = 0.0015 # S28 gene
  initial_condition_array[4] = 0.0    # no S70 gene, just S70 protein in the extract
  #initial_condition_array[11] = 0.020 # sigma28 protein in muM
  initial_condition_array[12] = 0.035 # S70 protein in muM

  # setup the W's -
  control_parameter_dictionary = customized_data_dictionary["control_parameter_dictionary"]
  control_parameter_dictionary["W_clssrA_RNAP"] = 0.000014
  control_parameter_dictionary["W_clssrA_sigma_28"] = 1.0
  control_parameter_dictionary["W_deGFPssrA_RNAP"] = 0.000014
  control_parameter_dictionary["W_deGFPssrA_sigma_70"] = 10.0
  control_parameter_dictionary["W_deGFPssrA_clssrA"] = 100.0
  control_parameter_dictionary["W_sigma_28_RNAP"] = 0.00014
  control_parameter_dictionary["W_sigma_28_sigma_70"] = 100.0
  control_parameter_dictionary["W_sigma_28_clssrA"] = 0.1
  control_parameter_dictionary["W_sigma_70_RNAP"] = 0.0

  # setup the binding parameters -
  bP = readdlm("bP-Opt.dat")
  binding_parameter_dictionary = customized_data_dictionary["binding_parameter_dictionary"]
  bp_symbol_name_array = String[]
  push!(bp_symbol_name_array, "n_clssrA_sigma_28")
  push!(bp_symbol_name_array, "K_clssrA_sigma_28")
  push!(bp_symbol_name_array, "n_deGFPssrA_sigma_70")
  push!(bp_symbol_name_array, "K_deGFPssrA_sigma_70")
  push!(bp_symbol_name_array, "n_deGFPssrA_clssrA")
  push!(bp_symbol_name_array, "K_deGFPssrA_clssrA")
  push!(bp_symbol_name_array, "n_sigma_28_sigma_70")
  push!(bp_symbol_name_array, "K_sigma_28_sigma_70")
  push!(bp_symbol_name_array, "n_sigma_28_clssrA")
  push!(bp_symbol_name_array, "K_sigma_28_clssrA")
  for (index,key_text) in enumerate(bp_symbol_name_array)

      # get value -
      value = bP[index]

      # cache -
      binding_parameter_dictionary[key_text] = value
  end

  # setup degradation_modifier_array -
  degradation_modifier_array = [
    0.0	    ;	# 1	clssrA
    0.0	    ;	# 2	deGFPssrA
    0.0	    ;	# 3	sigma_28
    0.0	    ;	# 4	sigma_70
    0.05    ;	# 5	mRNA_clssrA 0.0085
    0.05	;	# 6	mRNA_deGFPssrA 0.05
    0.05	;	# 7	mRNA_sigma_28
    1.0	    ;	# 8	mRNA_sigma_70
    30.0	;	# 9	protein_clssrA
    28.0	;	# 10	protein_deGFPssrA
    1.5	    ;	# 11	protein_sigma_28
    1.0	    ;	# 12	protein_sigma_70
  ]
  customized_data_dictionary["degradation_modifier_array"] = degradation_modifier_array
  biophysical_constants_dictionary = customized_data_dictionary["biophysical_constants_dictionary"]
  species_symbol_type_array = customized_data_dictionary["species_symbol_type_array"]
  protein_coding_length_array = customized_data_dictionary["protein_coding_length_array"]
  gene_coding_length_array = customized_data_dictionary["gene_coding_length_array"]

  # get gene IC -
  idx_gene = findall(x->x==:gene,species_symbol_type_array)
  gene_abundance_array = initial_condition_array[idx_gene]

  # Dilution degrdation matrix -
  dilution_degradation_matrix = build_dilution_degradation_matrix(biophysical_constants_dictionary,species_symbol_type_array,degradation_modifier_array)
  customized_data_dictionary["dilution_degradation_matrix"] = dilution_degradation_matrix

  # time constant modifier -
  time_constant_modifier_array = [
    0.0	                    ;	  # 1	clssrA
    0.0	                    ;	  # 2	deGFPssrA
    0.0	                    ;	  # 3	sigma_28
    0.0                     ;	  # 4	sigma_70
    0.5549555127597618      ;	  # 5	mRNA_clssrA
    0.9961201331095983      ;	  # 6	mRNA_deGFPssrA
    0.9573079006791378      ;	  # 7	mRNA_sigma_28
    2.47595196225463	    ;	  # 8	mRNA_sigma_70
    1.2340314042937286	    ;	  # 9	protein_clssrA
    2.0833904410438735 	    ;	  # 10	protein_deGFPssrA
    4.99999996459193	    ;	  # 11	protein_sigma_28
    2.47595196225463	    ;	  # 12	protein_sigma_70
  ]
  customized_data_dictionary["time_constant_modifier_array"] = time_constant_modifier_array

  # Precompute the translation parameters -
  translation_parameter_array = precompute_translation_parameter_array(biophysical_constants_dictionary, protein_coding_length_array, time_constant_modifier_array,host_type)
  customized_data_dictionary["translation_parameter_array"] = translation_parameter_array

  # Precompute the kinetic limit of transcription -
  transcription_kinetic_limit_array = precompute_transcription_kinetic_limit_array(biophysical_constants_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array, host_type)
  customized_data_dictionary["transcription_kinetic_limit_array"] = transcription_kinetic_limit_array

  # return -
  return customized_data_dictionary
end

function load_experimental_data_dictionary(base_path::String)

  # create experimental data dictionary -
  exp_data_dictionary = Dict{String,DataFrame}()

  # load data -
  # deGFP mRNA -
  deGFP_mRNA_path = "$(base_path)/deGFP-ssrA-mRNA-pooled.csv"
  deGFP_mRNA_data_table = CSV.read(deGFP_mRNA_path)

  # sigma28 mRNA -
  sigma28_mRNA_path = "$(base_path)/sigma28-mRNA-pooled.csv"
  sigma28_mRNA_data_table = CSV.read(sigma28_mRNA_path)

  # c1-ssrA mRNA -
  c1_ssrA_mRNA_path = "$(base_path)/c1-ssrA-mRNA-pooled.csv"
  c1_ssrA_mRNA_data_table = CSV.read(c1_ssrA_mRNA_path)

  # deGFP protein -
  deGFP_protein_path = "$(base_path)/deGFP-protein-pooled.csv"
  deGFP_protein_data_table = CSV.read(deGFP_protein_path)

  # package -------------------------------------------------- #
  exp_data_dictionary["deGFP_mRNA_data_table"] = deGFP_mRNA_data_table
  exp_data_dictionary["sigma28_mRNA_data_table"] = sigma28_mRNA_data_table
  exp_data_dictionary["c1_ssrA_mRNA_data_table"] = c1_ssrA_mRNA_data_table
  exp_data_dictionary["deGFP_protein_data_table"] = deGFP_protein_data_table
  return exp_data_dictionary
  # ---------------------------------------------------------- #
end

function sort_parameter_ensemble_array(parameter_ensemble_array::Array{Float64,2})

  # The last row is the error -
  error_array = parameter_ensemble_array[end,:]

  # what is the perm index?
  idx_sort = sortperm(error_array) # default: small to big

  # shuffle cols -
  sorted_parameter_array = parameter_ensemble_array[:,idx_sort]

  # return -
  return sorted_parameter_array
end

function calculate_fisher_information_matrix(scaled_sensitivity_array,measurement_weight_array,id_parameter_index_array)

  # Get the size of the system -
  (number_of_rows,number_of_cols) = size(scaled_sensitivity_array);

  # grab the id cols -
  id_sensitivity_array = scaled_sensitivity_array[:,id_parameter_index_array]

  # Create the FIM (fisher_information_matrix) -
  FIM = transpose(id_sensitivity_array)*measurement_weight_array*id_sensitivity_array

  # return the FIM to the caller -
  return FIM
end

function estimate_identifiable_parameters(scaled_sensitivity_array,epsilon)

  # Get the size of the system -
  (number_of_rows,number_of_cols) = size(scaled_sensitivity_array);
  X = zeros(number_of_rows,1)
  pset = Int64[]

  R = scaled_sensitivity_array

  for col_index = 1:number_of_cols

    # get mag of col -
    local_m_array = zeros(number_of_cols)
    for inner_col_index = 1:number_of_cols

      value = R[:,inner_col_index]'*R[:,inner_col_index]
      local_m_array[inner_col_index] = value[1]
    end

    # what is the maximum element?
    max_value = maximum(local_m_array)
    max_index = indmax(local_m_array)

    # check tolerance -
    if max_value>epsilon

      # Grab this parameter -
      push!(pset,max_index)

      # Grab this col -
      X = [X scaled_sensitivity_array[:,max_index]]

      # remove leading col if first time through -
      if (col_index == 1)
        X = X[:,2:end]
      end

      # Create local array -
      Shat=X*inv(X'*X)*X'*scaled_sensitivity_array

      # Update R -
      R = scaled_sensitivity_array - Shat
    end
  end

  return sort!(pset)
end

function calculate_sensitivity_array(path_to_senstivity_files,file_pattern,time_skip,data_dictionary)

  # what is my system dimension?
  number_of_states = data_dictionary["number_of_states"]

  # load the files -
  searchdir(path,key) = filter(x->contains(x,key),readdir(path))

  # block_dictionary -
  block_dictionary = Dict()

  # build file list -
  list_of_files = searchdir(path_to_senstivity_files,file_pattern)
  number_of_files = length(list_of_files)
  time_array = []
  for file_index = 1:number_of_files

    # Build path -
    path_to_data_file = path_to_senstivity_files*"/"*file_pattern*string(file_index)*".dat"

    # Load file -
    local_data_array = readdlm(path_to_data_file)

    # split -
    time_array = local_data_array[:,1]
    X = local_data_array[:,2:end]
    state_array = X[:,1:number_of_states]
    sensitivity_array = X[:,(number_of_states+1):end]
    scaled_sensitivity_block = scale_sensitivity_array(time_array,state_array,sensitivity_array,file_index,data_dictionary)

    # store the transpose -
    key_symbol = file_pattern*string(file_index)*".dat"
    block_dictionary[key_symbol] = transpose(scaled_sensitivity_block)
  end

  # what is my system dimension?
  number_of_timesteps = length(time_array)
  number_of_parameters = number_of_files

  # initialize -
  sensitivity_array = zeros(number_of_states,number_of_parameters)
  sample_time_array = Float64[]
  for time_step_index = 1:time_skip:number_of_timesteps

    time_value = time_array[time_step_index]
    push!(sample_time_array,time_value)

    local_sens_block = zeros(number_of_states,number_of_parameters)
    for parameter_index = 1:number_of_parameters

      # get the block from the dictionary -
      key_symbol = file_pattern*string(parameter_index)*".dat"
      block = block_dictionary[key_symbol]

      # grab the col -
      block_col = block[:,time_step_index]

      for state_index = 1:number_of_states
        local_sens_block[state_index,parameter_index] = block_col[state_index]
      end
    end

    # add block to s array -
    sensitivity_array = [sensitivity_array ; local_sens_block]
  end

  # cutoff leading block -
  sensitivity_array = sensitivity_array[(number_of_states+1):end,:]
  return (sample_time_array,sensitivity_array)
end

function calculate_average_scaled_sensitivity_array(path_to_senstivity_files,file_pattern,data_dictionary)

  # what is my system dimension?
  number_of_states = data_dictionary["number_of_states"]

  # initialize -
  average_scaled_sensitivity_array = zeros(number_of_states,1)

  # load the files -
  searchdir(path,key) = filter(x->contains(x,key),readdir(path))

  # build file list -
  list_of_files = searchdir(path_to_senstivity_files,file_pattern)
  number_of_files = length(list_of_files)
  for file_index = 1:number_of_files

    # Build path -
    path_to_data_file = path_to_senstivity_files*"/"*file_pattern*string(file_index)*".dat"

    # Load file -
    local_data_array = readdlm(path_to_data_file)

    # split -
    time_array = local_data_array[:,1]
    X = local_data_array[:,2:end]
    state_array = X[:,1:number_of_states]
    sensitivity_array = X[:,(number_of_states+1):end]
    scaled_sensitivity_array = scale_sensitivity_array(time_array,state_array,sensitivity_array,file_index,data_dictionary)

    # time average -
    average_sensitivity_col = time_average_array(time_array,scaled_sensitivity_array)

    # grab -
    average_scaled_sensitivity_array = [average_scaled_sensitivity_array average_sensitivity_col]
  end

  # trim leading col -
  average_scaled_sensitivity_array = average_scaled_sensitivity_array[:,2:end]
  return average_scaled_sensitivity_array
end

function time_average_array(time_array,data_array)

  # what is the delta T?
  delta_time = (time_array[end] - time_array[1])

  # initialize -
  average_array = Float64[]

  # what is the size of the array?
  (number_of_timesteps,number_of_states) = size(data_array)
  for state_index = 1:number_of_states

    # grab the data column -
    data_col = data_array[:,state_index]

    # average -
    average_value = (1/delta_time)*trapz(time_array,data_col)

    # push -
    push!(average_array,average_value)
  end

  return average_array

end

function scale_sensitivity_array(time_array,state_array,sensitivity_array,parameter_index,data_dictionary)

  # what is small?
  epsilon = 1e-6

  # initialize -
  (number_of_timesteps,number_of_states) = size(state_array)
  scaled_sensitivity_array = zeros(number_of_timesteps,number_of_states)

  # What is the nominal parameter value?
  parameter_name_mapping_array = data_dictionary["parameter_name_mapping_array"]
  parameter_name = parameter_name_mapping_array[parameter_index]

  # build the total parameter dictionary -
  binding_parameter_dictionary = data_dictionary["binding_parameter_dictionary"]
  control_parameter_dictionary = data_dictionary["control_parameter_dictionary"]
  txtl_parameter_dictionary = data_dictionary["biophysical_constants_dictionary"]
  total_parameter_dictionary = merge(binding_parameter_dictionary,control_parameter_dictionary,txtl_parameter_dictionary)

  # Grab the default value -
  default_parameter_value = total_parameter_dictionary[parameter_name]

  # main loop -
  for (time_index,time_value) in enumerate(time_array)

    for state_index = 1:number_of_states

      state_value = state_array[time_index,state_index]
      if (state_value<epsilon)
        state_value = epsilon
      end

      old_value = sensitivity_array[time_index,state_index]
      new_value = old_value*(default_parameter_value/state_value)
      scaled_sensitivity_array[time_index,state_index] = new_value
    end
  end

  return scaled_sensitivity_array
end

function calculate_jacobian(time,state_array,data_dictionary)

  # what is the size of the system?
  number_of_states = length(state_array)

  # calculate each row of the jacobian -
  jacobian_array = zeros(1,number_of_states)
  for (state_index,state_value) in enumerate(state_array)

    jacobian_row = calculate_jacobian_row(time,state_array,state_index,data_dictionary)
    jacobian_array = [jacobian_array  ; transpose(jacobian_row)]
  end

  jacobian_array = jacobian_array[2:end,:]
  return jacobian_array
end

function calculate_bmatrix(time,state_array,data_dictionary)

  # what is the size of the system?
  parameter_name_mapping_array = data_dictionary["parameter_name_mapping_array"]
  number_of_parameters = length(parameter_name_mapping_array)
  number_of_states = length(state_array)

  # calculate each row of the jacobian -
  b_array = zeros(number_of_parameters,1)
  for (state_index,state_value) in enumerate(state_array)

    bmatrtix_row = calculate_bmatrix_row(time,state_array,state_index,data_dictionary)
    b_array = [b_array bmatrtix_row]
  end

  b_array = b_array[:,2:end]
  return transpose(b_array)
end

function calculate_bmatrix_row(time,state_array,balance_index,data_dictionary)

  # define some constants -
  epsilon = 1e-6
  delta = 0.01

  # parameter name dictionary -
  parameter_name_mapping_array = data_dictionary["parameter_name_mapping_array"]
  binding_parameter_dictionary = data_dictionary["binding_parameter_dictionary"]
  control_parameter_dictionary = data_dictionary["control_parameter_dictionary"]
  txtl_parameter_dictionary = data_dictionary["biophysical_constants_dictionary"]
  number_of_parameters = length(parameter_name_mapping_array)

  # how many binding parameters do we have?
  number_of_binding_parameters = length(binding_parameter_dictionary)
  number_of_control_parameters = length(control_parameter_dictionary)

  # create a mega dictionary -
  total_parameter_dictionary = merge(binding_parameter_dictionary,control_parameter_dictionary,txtl_parameter_dictionary)

  # create delta parameter array -
  parameter_delta_array = Float64[]
  for (parameter_index,parameter_name) in enumerate(parameter_name_mapping_array)

    # Grab the default value -
    default_parameter_value = total_parameter_dictionary[parameter_name]

    # state perturbation -
    peturbed_parameter_value = default_parameter_value*(delta);

    #@show (peturbed_parameter_value,default_parameter_value)

    # check -
    if (peturbed_parameter_value<epsilon)
      peturbed_parameter_value = epsilon
    end

    # capture -
    push!(parameter_delta_array,(peturbed_parameter_value))
  end

  # Create the diag array -
  diag_delta_array = Diagonal(vec(parameter_delta_array))

  # Create bVec -
  f_nominal = calculate_balances(time,vec(state_array),data_dictionary)

  # estimate the perturbed balances -
  rhs_delta_array = Float64[]
  for (parameter_index,parameter_name) in enumerate(parameter_name_mapping_array)

    # copy -
    local_data_dictionary = deepcopy(data_dictionary)

    # Grab the default value -
    default_parameter_value = total_parameter_dictionary[parameter_name]

    # update the state -
    perturbed_parameter_array = zeros(number_of_parameters)
    for local_index = 1:number_of_parameters
      if (parameter_index == local_index)
        perturbed_parameter_array[parameter_index] = default_parameter_value*(1+delta);
      else
        local_parameter_name = parameter_name_mapping_array[local_index]
        perturbed_parameter_array[local_index] = total_parameter_dictionary[local_parameter_name]
      end
    end

    if (parameter_index<=number_of_binding_parameters)

      # we are in the binding section -
      local_data_dictionary["binding_parameter_dictionary"][parameter_name] = perturbed_parameter_array[parameter_index]
    elseif (parameter_index>number_of_binding_parameters && parameter_index<=(number_of_binding_parameters+number_of_control_parameters))
      # we are in the control section -
      local_data_dictionary["control_parameter_dictionary"][parameter_name] = perturbed_parameter_array[parameter_index]
    else
      # we are in the bar parameter section -
      local_data_dictionary[parameter_name] = perturbed_parameter_array[parameter_index]
    end

    # calculate the perturbed balances -
    f_perturbed = calculate_balances(time,vec(state_array),local_data_dictionary)
    f_perturbed = f_perturbed[balance_index] - f_nominal[balance_index]

    # capture -
    push!(rhs_delta_array,f_perturbed)
  end

  # calculate the bmatrix row -
  bmatrix_row = diag_delta_array\rhs_delta_array

  # return -
  return bmatrix_row
end

function finite_diff_jacobian(time,state_array,data_dictionary)

  # define some constants -
  epsilon = 1e-6
  delta = 0.05
  number_of_states = length(state_array)

  # initialize -
  jacobian_array = zeros(number_of_states,number_of_states)

  # nominal -
  f_nominal = calculate_balances(time,vec(state_array),data_dictionary)

  #@show state_array

  for row_index = 1:number_of_states
    for col_index = 1:number_of_states

      perturbed_state_array = zeros(number_of_states)
      for perturbation_index = 1:number_of_states

        if (col_index == perturbation_index)
          perturbed_state_array[col_index] = state_array[col_index]*(1+delta)
        else
          perturbed_state_array[perturbation_index] = state_array[perturbation_index]
        end
      end

      #@show perturbed_state_array

      # calculate the balances -
      f_perturbed = calculate_balances(time,vec(perturbed_state_array),data_dictionary)

      # calculate the entry -
      perturbation_size = state_array[col_index]*delta
      if (perturbation_size<epsilon)
        perturbation_size = epsilon
      end
      jacobian_array[row_index,col_index] = (f_perturbed[row_index] - f_nominal[row_index])/(perturbation_size)
    end
  end

  return jacobian_array
end

function calculate_jacobian_row(time,state_array,balance_index,data_dictionary)

  # define some constants -
  epsilon = 1e-6
  delta = 0.001
  number_of_states = length(state_array)

  # Create the delta state array -
  state_delta_array = Float64[]
  for (state_index,state_value) in enumerate(state_array)

    # state perturbation -
    peturbed_state = state_value*(delta);

    # check -
    if (peturbed_state<epsilon)
      peturbed_state = epsilon
    end

    # capture -
    push!(state_delta_array,peturbed_state)
  end

  # Create the diag array -
  diag_delta_array = Diagonal(vec(state_delta_array))

  # Create bVec -
  f_nominal = calculate_balances(time,vec(state_array),data_dictionary)

  # estimate the perturbed balances -
  rhs_delta_array = Float64[]
  for (state_index,state_value) in enumerate(state_array)

    # update the state -
    perturbed_state_array = zeros(number_of_states)

    for local_index = 1:number_of_states
      if (state_index == local_index)
        perturbed_state_array[state_index] = state_value*(1+delta);
      else
        perturbed_state_array[local_index] = state_array[local_index]
      end
    end

    # calculate the perturbed balances -
    f_perturbed = calculate_balances(time,vec(perturbed_state_array),data_dictionary)
    f_perturbed_delta = f_perturbed[balance_index] - f_nominal[balance_index]

    # capture -
    push!(rhs_delta_array,f_perturbed_delta)
  end

  # calculate the jacobian row -
  jacobian_row = diag_delta_array\rhs_delta_array

  # return -
  return jacobian_row
end

function build_biophysical_dictionary(path_to_biophysical_constants_file::String, host_type::Symbol)::Dict{String,Any}

    # check the path - is it legit?
    check_file_existence(path_to_biophysical_constants_file)

    # Load the model dictionary -
    model_dictionary = JSON.parsefile(path_to_biophysical_constants_file)

    # check the host type in the biophysical file - need to match
    file_host_type = Symbol(model_dictionary["host_type"])
    if (file_host_type != host_type)
        throw(error("inconsistent host_type. Expected $(host_type) found $(file_host_type)"))
    end

    # initialize -
    biophysical_dictionary = Dict{String,Any}()

    # convert -
    list_of_constant_dictionaries = model_dictionary["biophysical_constants"]
    for (key_string, local_dictionary) in list_of_constant_dictionaries

        # get the value -
        value = parse(Float64, local_dictionary["value"])

        # cache -
        biophysical_dictionary[key_string] = value
    end

    # return -
    return biophysical_dictionary
end

function trapz(x::Array{Float64}, y::Array{Float64})

    # Trapezoidal integration rule
    local n = length(x)
    if (length(y) != n)
        error("Vectors 'x', 'y' must be of same length")
    end

    r = 0
    if n == 1; return r; end
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    #= correction -h^2/12 * (f'(b) - f'(a))
    ha = x[2] - x[1]
    he = x[end] - x[end-1]
    ra = (y[2] - y[1]) / ha
    re = (y[end] - y[end-1]) / he
    r/2 - ha*he/12 * (re - ra)
    =#
    return r/2
end

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
using SpectralMethodFFT

struct TimeIntegTools

    klm_A::KFunc
    klm_B::KFunc
    klm_C::KFunc

    function TimeIntegTools(klm_A, klm_B, klm_C)
        new(klm_A, klm_B, klm_C)
    end

    function TimeIntegTools(
        Re⁻¹::Real, dt::Real, c::ConfigFFT)

        xlens = xlens_xranges(c.xranges)

        klm_K = klm_Kgen(c)
        klm_L = klm_Lgen(c)
        klm_M = klm_Mgen(c)

        klm_A = exp(
            - Re⁻¹ * dt / 2.0
            * sum(
                (2π .* [klm_K, klm_L, klm_M] ./ xlens).^2
            )
        )

        klm_B = exp(
            - Re⁻¹ * dt / 2.0
            * (2π * klm_L / xlens[2])^2
        )

        function f(x)
            if abs(x) < eps()
                dt/2.0 - x*dt^2/8.0
            else
                (1 - exp(-x*dt/2.0)) / x
            end
        end

        xyz_Y = xyz_Ygen(c)
        klm_F = klm_xyz(sin(π/2 * xyz_Y))

        klm_C = KFunc(f.(Re⁻¹ * klm_L.vals), c) * klm_F

        TimeIntegTools(klm_A, klm_B, klm_C)
    end

end

function K_K_K_K_lin_K_K_K_K(
        klm_Ψ::KFunc{T,3},
        klm_Φ::KFunc{T,3},
        klm_Ux::KFunc{T,3},
        klm_Uz::KFunc{T,3},
        tit::TimeIntegTools
    ) where T

    return (
        tit.klm_A * klm_Ψ,
        tit.klm_A * klm_Φ,
        tit.klm_B * klm_Ux + tit.klm_C,
        tit.klm_B * klm_Uz
    )

end

function klm_Δxz_klm(f::KFunc{T,3}) where T

    return (
        (klm_∂x_klm ∘ klm_∂x_klm)(f)
        +
        (klm_∂z_klm ∘ klm_∂z_klm)(f)
    )

end

function klm_Δxz⁻¹_klm(f::KFunc{T,3}) where T

    c = f.config
    xlens = xlens_xranges(c.xranges)

    lap = - 4π^2 * (
        ( c.Kcoords[1] / xlens[1] ) .^ 2
        +
        ( c.Kcoords[3] / xlens[3] ) .^ 2
    )
    gvals = f.vals ./ lap

    zero_indiceses = [(
        map(c.ngrids) do x
            if x % 2 == 0
                return (1, x÷2 + 1)
            else
                return (1,)
            end
        end
    )...]
    zero_indiceses = [
        dim==2 ? tuple(1:c.ngrids[2]...) : zero_indiceses[dim]
        for dim = 1:3
    ]

    for indices in Iterators.product(zero_indiceses...)
        gvals[indices...] = 0.
    end

    return KFunc(gvals, c)

end

function K3_K3_uω_ΨΦUxUz_K_K_K_K(
        klm_Ψ::KFunc{T,3},
        klm_Φ::KFunc{T,3},
        klm_Ux::KFunc{T,3},
        klm_Uz::KFunc{T,3}
    ) where T

    klm_ψ = klm_Δxz⁻¹_klm(klm_Ψ)
    klm_uy = - klm_Δ⁻¹_klm(klm_Φ)
    klm_ϕ = - klm_Δxz⁻¹_klm(klm_uy)

    # --- u ---
    klm_ux = (
        - klm_∂z_klm(klm_ψ)
        + klm_∂x_klm(klm_∂y_klm(klm_ϕ))
        + klm_Ux
    )
    klm_uz = (
        + klm_∂x_klm(klm_ψ)
        + klm_∂z_klm(klm_∂y_klm(klm_ϕ))
        + klm_Uz
    )

    klm3_u = [
        klm_ux,
        klm_uy,
        klm_uz
    ]

    # --- ω ---
    klm3_ω = klm3_rot_klm3(klm3_u)

    return klm3_u, klm3_ω

end

function K_K_K_K_nonlin_K_K_K_K(
        klm_Ψ::KFunc{T,3},
        klm_Φ::KFunc{T,3},
        klm_Ux::KFunc{T,3},
        klm_Uz::KFunc{T,3}
    ) where T

    klm3_u, klm3_ω = (
        K3_K3_uω_ΨΦUxUz_K_K_K_K(
            klm_Ψ, klm_Φ, klm_Ux, klm_Uz
        )
    )

    klm_ux, klm_uy, klm_uz = klm3_u

    klm3_tmp1 = klm3_rot_klm3(klm3_ω × klm3_u)
    klm3_tmp2 = - klm3_rot_klm3(klm3_tmp1)

    klm_nonlin3 = - klm_∂y_klm(klm_mean_xz_klm(
        klm_uy ∗ klm_ux
    ))
    klm_nonlin4 = - klm_∂y_klm(klm_mean_xz_klm(
        klm_uy ∗ klm_uz
    ))

    return (
        klm3_tmp1[2],
        klm3_tmp2[2],
        klm_nonlin3,
        klm_nonlin4
    )

end

function klm_mean_xz_klm(f::KFunc{T,3}) where T

    ngrids = f.config.ngrids
    max_nwaves = (
        0, ngrids[2] ÷ 2 + 1, 0
    )
    return klm_lowpass_klm(f, max_nwaves)

end

function K_K_K_K_RK4_K_K_K_K(
        klm_Ψ::KFunc{T,3},
        klm_Φ::KFunc{T,3},
        klm_Ux::KFunc{T,3},
        klm_Uz::KFunc{T,3},
        tit::TimeIntegTools,
        dt::Real
    ) where T

    klm4_ΨΦUxUz = (
        klm_Ψ, klm_Φ, klm_Ux, klm_Uz
    )

    klm4_ΨΦUxUz1 = (
        dt .* K_K_K_K_nonlin_K_K_K_K(
            klm4_ΨΦUxUz...
        )
    )

    klm4_ΨΦUxUz2 = (
        dt .* K_K_K_K_nonlin_K_K_K_K(
            K_K_K_K_lin_K_K_K_K(
                (
                    klm4_ΨΦUxUz .+ klm4_ΨΦUxUz1 ./ 2
                )...,
                tit
            )...
        )
    )

    klm4_ΨΦUxUz3 = (
        dt .* K_K_K_K_nonlin_K_K_K_K(
            (
                K_K_K_K_lin_K_K_K_K(
                    klm4_ΨΦUxUz..., tit
                )
                .+ klm4_ΨΦUxUz2 ./ 2
            )...
        )
    )

    klm4_ΨΦUxUz4 = (
        dt .* K_K_K_K_nonlin_K_K_K_K(
            K_K_K_K_lin_K_K_K_K(
                (
                    K_K_K_K_lin_K_K_K_K(
                        klm4_ΨΦUxUz..., tit
                    )
                    .+ klm4_ΨΦUxUz3
                )...,
                tit
            )...
        )
    )

    return (
        K_K_K_K_lin_K_K_K_K(
            (
                K_K_K_K_lin_K_K_K_K(
                    (
                        klm4_ΨΦUxUz .+ klm4_ΨΦUxUz1 ./ 6
                    )...,
                    tit
                )
                .+
                (
                    klm4_ΨΦUxUz2 .+ klm4_ΨΦUxUz3
                ) ./ 3
            )...,
            tit
        )
        .+ klm4_ΨΦUxUz4 ./ 6
    )

end

function main()

    # *************************
    #   Setting
    # *************************
    # --- Space ---
    xranges = ((0., 16.), (-1., 3.), (0., 16.))
    ngrids = (64, 16, 64)

    # --- Time ---
    t_st = 0.0
    t_ed = 1.0
    nt = 1000
    dt = (t_ed - t_st) / nt
    Re⁻¹ = 1 / (2.56 * 180)

    # *************************
    #   Initialization
    # *************************
    config = ConfigFFT(ngrids, xranges)
    tit = TimeIntegTools(Re⁻¹, dt, config)

    xyz_X = xyz_Xgen(config)
    xyz_Y = xyz_Ygen(config)
    xyz_Z = xyz_Zgen(config)

    xyz_ψ = sin(π/2 * xyz_Y)*sin(exp(π*xyz_X))*cos(exp(π*xyz_Z))
    xyz_ϕ = sin(π/2 * xyz_Y)*cos(exp(π*xyz_X))*sin(exp(π*xyz_Z))
    xyz_Ux = XFunc(zeros(config.ngrids), config)
    xyz_Uz = XFunc(zeros(config.ngrids), config)

    klm_ψ = klm_xyz(xyz_ψ)
    klm_ϕ = klm_xyz(xyz_ϕ)
    klm_Ux = klm_xyz(xyz_Ux)
    klm_Uz = klm_xyz(xyz_Uz)

    klm_Ψ = klm_Δxz_klm(klm_ψ)
    klm_Φ = klm_Δ_klm(klm_Δxz_klm(klm_ϕ))

    for it = 0:nt

        # output here
        println("it=", it)
        klm3_u, klm3_ω = K3_K3_uω_ΨΦUxUz_K_K_K_K(
            klm_Ψ, klm_Φ, klm_Ux, klm_Uz
        )
        xyz3_u = xyz_klm.(klm3_u)
        println(xyz3_u[2][1, 3, 1])

        if it < nt
            klm_Ψ, klm_Φ, klm_Ux, klm_Uz = (
                K_K_K_K_RK4_K_K_K_K(
                    klm_Ψ, klm_Φ, klm_Ux, klm_Uz,
                    tit, dt
                )
            )
        end
    end

end

main()

max_stars_repo_name: mat-der-D/SpectralMethodFFTW

max_issues_repo_name: mat-der-D/SpectralMethodFFTW

max_forks_repo_name: mat-der-D/SpectralMethodFFTW


text:
export git_otype, message, tree, tree_id,
       author, committer, parent_id, parent_count, parents

function message(c::GitCommit, raw::Bool=false)
    local msg_ptr::Ptr{Uint8}
    msg_ptr = raw? ccall((:git_commit_message_raw, libgit2), Ptr{Uint8}, (Ptr{Void},), c) :
                   ccall((:git_commit_message, libgit2), Ptr{Uint8}, (Ptr{Void},), c)
    if msg_ptr == C_NULL
        return nothing
    end
    return utf8(bytestring(msg_ptr))
end

function GitTree(c::GitCommit)
    tree_ptr = Ptr{Void}[0]
    @check ccall((:git_commit_tree, libgit2), Cint, (Ptr{Ptr{Void}}, Ptr{Void}), tree_ptr, c.ptr)
    return GitTree(tree_ptr[1])
end

git_tree_id(c::GitCommit) = Oid(ccall((:git_commit_tree_id, libgit2), Ptr{Uint8}, (Ptr{Void},), c))

function author(c::GitCommit)
    ptr = ccall((:git_commit_author, libgit2), Ptr{SignatureStruct}, (Ptr{Void},), c)
    @assert ptr != C_NULL
    return Signature(ptr)
end

function committer(c::GitCommit)
    ptr = ccall((:git_commit_committer, libgit2), Ptr{SignatureStruct}, (Ptr{Void},), c)
    @assert ptr != C_NULL
    return Signature(ptr)
end

Base.parent(c::GitCommit, n::Integer) = begin
    n >= 1 || throw(ArgumentError("n must be greater than or equal to 1"))
    commit_ptr = Ptr{Void}[0]
    @check ccall((:git_commit_parent, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Cuint), commit_ptr, c, n-1)
    return GitCommit(commit_ptr[1])
end

parent_count(c::GitCommit) = int(ccall((:git_commit_parentcount, libgit2), Csize_t, (Ptr{Void},), c))
parents(c::GitCommit) = (n = parent_count(c); GitCommit[parent(c, i) for i=1:n])

function parent_id(c::GitCommit, n::Integer)
    n >= 1 || throw(ArgumentError("n must be greater than or equal to 1"))
    oid_ptr = ccall((:git_commit_parent_id, libgit2), Ptr{Uint8}, (Ptr{Void}, Cuint), c, n-1)
    return Oid(oid_ptr)
end


text:
# File: fMean.jl
# Author:
# Date: 29-June-2020

# MEAN
function fMean(x,f)
# -------------------------------------------------------------------
# x:     returns
# f:     reporting frequency
# m1:    mean
# -------------------------------------------------------------------

    m1 = mean(x).*f
    return m1
end


text:
module utils

export @swap!
export dec2bin

macro swap!(a::Symbol,b::Symbol)
	blk = quote
		if typeof($(esc(a))) != typeof($(esc(b)))
			throw(ArgumentError("Arrays of different type"))
		else
			c = $(esc(a))
			$(esc(a)) = $(esc(b))
			$(esc(b)) = c
		end
	end
	return blk
end

function dec2bin(number::Int64, size::Int64)
	binary = Array{Int64, 1}(undef, size)
	quotient = number
	remainder = 0
	for i=size:-1:1
		remainder = mod(quotient, 2)
		quotient = div(quotient, 2)
		binary[i] = remainder
	end
	return binary
end

end


text:
# Autogenerated wrapper script for ldid_jll for armv6l-linux-gnueabihf-cxx03
export ldid

using libplist_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("ldid")
JLLWrappers.@declare_executable_product(ldid)
function __init__()
    JLLWrappers.@generate_init_header(libplist_jll, OpenSSL_jll)
    JLLWrappers.@init_executable_product(
        ldid,
        "bin/ldid",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ArgParse
using Formatting
using JSON,TOML
using DataStructures

function runOneCmd(cmd)
    println(">> ", string(cmd))
    run(cmd)
end

function pleaseRun(cmd)
    if isa(cmd, Array)
        for i = 1:length(cmd)
            runOneCmd(cmd[i])
        end
    else
        runOneCmd(cmd)
    end
end

function getCESMConfig(path, ids)

    wdir = pwd()

    cd(path)

    d = OrderedDict()
    for id in ids
        d[id] = readchomp(`./xmlquery $(id) -silent -valonly`)
    end

    cd(wdir)

    return d
end

function setCESMConfig(path, filename, kv_list)

    wdir = pwd()
    cd(path)
    for (key, val) in kv_list
        pleaseRun(`
            ./xmlchange -f $(filename) -id $(key) -val $(val)
        `)
    end
    cd(wdir)
end

s = ArgParseSettings()
@add_arg_table s begin

    "--project"
        help = "Project code sent to PBS scheduling system"
        arg_type = String
        required = true

    "--casename"
        help = "Casename"
        arg_type = String
        required = true

    "--root"
        help = "The folder where case folders are contained"
        arg_type = String
        required = true

    "--ncpu"
        help = "The number of CPUs that IOM can use"
        arg_type = Int64
        required = true

    "--walltime"
        help = "Walltime sent to PBS scheduling system"
        arg_type = String
        required = true

    "--queue"
        help = "Queue sent to PBS scheduling system"
        arg_type = String
        default = "economy"

    "--resolution"
        help = "Casename"
        arg_type = String
        required = true

    "--compset"
        help = "Casename"
        arg_type = String
        required = true

    "--machine"
        help = "Casename"
        arg_type = String
        required = true

    "--cesm-root"
        help = "Casename"
        arg_type = String
        required = true

    "--user-namelist-dir"
        help = "Casename"
        arg_type = String
        default = ""

    "--cesm-env"
        help = "The TOML file used to set env_mach_pes.xml, env_run.xml"
        arg_type = String
        default = ""

    "--build"
        help = "If set then will try to build the case."
        action = :store_true

end

parsed = DataStructures.OrderedDict(parse_args(ARGS, s))

JSON.print(parsed, 4)

if parsed["cesm-env"] != ""
    overwrite_cesm_env = TOML.parsefile(parsed["cesm-env"]) |> DataStructures.OrderedDict
    println("Loaded env file ", parsed["cesm-env"])
    JSON.print(overwrite_cesm_env, 4)
else
    overwrite_cesm_env = Dict()
end

mkpath(parsed["root"])
cd(parsed["root"])

if isdir(parsed["casename"])
    throw(ErrorException(format("Error: `{:s}` already exists.", parsed["casename"])))
end

pleaseRun(`
$(parsed["cesm-root"])/scripts/create_newcase
    -case    $(parsed["casename"])
    -compset $(parsed["compset"])
    -res     $(parsed["resolution"])
    -mach    $(parsed["machine"])
`)

if !isdir(parsed["casename"])
    throw(ErrorException(format("Error: `{:s}` is not created")))
end

println("Entering case folder $(parsed["casename"])")
cd(parsed["casename"])

if parsed["user-namelist-dir"] != ""
    println("copy user namelist")
    pleaseRun(`
        cp $(parsed["user-namelist-dir"])/user_nl_\* .
    `)
end

for k in ["env_run", "env_mach_pes"]
    if haskey(overwrite_cesm_env, k)
        setCESMConfig(pwd(), "$(k).xml", overwrite_cesm_env[k])
    end
end

pleaseRun(`./cesm_setup`)
pleaseRun(`./cesm_setup -clean`)
pleaseRun(`./cesm_setup`)

cesm_env = getCESMConfig(pwd(), [
    "CASEROOT",
    "RUNDIR",
    "DIN_LOC_ROOT",
    "DOUT_S_ROOT",
    "OCN_DOMAIN_FILE",
    "OCN_DOMAIN_PATH",
    "TOTALPES",
    "MAX_TASKS_PER_NODE",
])

println("env data: ")
JSON.print(cesm_env, 4)

# The $casename.run file is created by cesm_setup and mk_batch files in CESM1 util codes
cp("$(parsed["casename"]).run", "$(parsed["casename"]).cesm.run", force=true)

cesm_nodes = ceil(Int64, parse(Float64, cesm_env["TOTALPES"]) / parse(Float64, cesm_env["MAX_TASKS_PER_NODE"]))

open(joinpath(cesm_env["CASEROOT"], "$(parsed["casename"]).run"), "w") do io
    write(io, """
#PBS -A $(parsed["project"])
#PBS -N $(parsed["casename"])
#PBS -q $(parsed["queue"])
#PBS -l select=$(cesm_nodes+1):ncpus=$(cesm_env["MAX_TASKS_PER_NODE"]):mpiprocs=$(cesm_env["MAX_TASKS_PER_NODE"]):ompthreads=1
#PBS -l walltime="$(parsed["walltime"])"
#PBS -j oe
#PBS -S /bin/bash

#!/bin/bash

bash $(cesm_env["CASEROOT"])/$(parsed["casename"]).destroy_tunnel
/bin/csh $(cesm_env["CASEROOT"])/$(parsed["casename"]).cesm.run &
$(cesm_env["CASEROOT"])/$(parsed["casename"]).ocn.run &
wait

""")
end

open(joinpath(cesm_env["CASEROOT"],"$(parsed["casename"]).ocn.run"), "w") do io
    write(io, """
#!/bin/bash

ml av
ml load openmpi/4.0.3
ml load julia/1.6.0
julia --project -e 'ENV["JULIA_MPI_BINARY"]="system"; using Pkg; Pkg.build("MPI"; verbose=true)'

LID=\$( date +%y%m%d-%H%M%S )

caseroot="$(cesm_env["CASEROOT"])"
caserun="$(cesm_env["RUNDIR"])"
archive_log_dir="$(cesm_env["DOUT_S_ROOT"])/ocn/logs"
logfile="iom.log.\${LID}"

mpiexec -n $(parsed["ncpu"]) julia --project \${caseroot}/IOM/src/CESM_driver/main.jl --config-file=\${caseroot}/config.toml &> \${caserun}/\${logfile}

ret_code=\$?
echo "Return code: \$ret_code"

if [ \$ret_code -eq 0 ]; then
    echo "Program ends succefully. Move log file to arcihve."
    mv "\${caserun}/\${logfile}" "\${archive_log_dir}/"
else
    echo "Program ends abnomally. Please check."
fi

""")
end

# Create sh that removes x_tmp. This is useful when run on local machine for developing
open(joinpath(cesm_env["CASEROOT"], "$(parsed["casename"]).destroy_tunnel"), "w") do io
    write(io, """
#!/bin/bash
rm -rf $(cesm_env["RUNDIR"])/x_tmp/*
""", )
end

# If user clean_build and rerun cesm_setup, the run file will be overwritten.
cp("$(parsed["casename"]).run", "backup_$(parsed["casename"]).run", force=true)
open(joinpath(cesm_env["CASEROOT"],"$(parsed["casename"]).recover_run_file"), "w") do io
    write(io, """
#!/bin/bash
cp "$(parsed["casename"]).run" "$(parsed["casename"]).cesm.run"
cp "backup_$(parsed["casename"]).run" "$(parsed["casename"]).run"

""")
end

pleaseRun(`chmod +x $(parsed["casename"]).run`)
pleaseRun(`chmod +x $(parsed["casename"]).ocn.run`)
pleaseRun(`chmod +x $(parsed["casename"]).destroy_tunnel`)
pleaseRun(`chmod +x $(parsed["casename"]).recover_run_file`)

open(joinpath(cesm_env["CASEROOT"],"user_nl_cam"), "w") do io
    write(io, """
nhtfrq=0
empty_htapes=.true.
fincl1='SOLIN','FSUTOA','FLUT','FLUTC','CLDHGH','CLDLOW','CLDMED','CLDTOT','FLDS','FLNS','FSDS','FSDSC','FSNS','FSNSC','FLNT','FSNT','ICEFRAC','LHFLX','LWCF','PHIS','PRECC','PRECL','PRECSC','PRECSL','PS','PSL','QFLX','QREFHT','SHFLX','SNOWHICE','SNOWHLND','SWCF','TAUX','TAUY','TMQ','TREFHT','TS','U10','CLOUD','OMEGA','Q','T','U','V','VT','VU','VQ','Z3','SST'
""")
end

pleaseRun(`git clone --branch "dev/cesm-coupling" https://github.com/meteorologytoday/IOM.git`)

cd(joinpath("SourceMods", "src.docn"))
pleaseRun(`ln -s ../../IOM/src/CESM_driver/cesm1_tb_docn_comp_mod.F90 ./docn_comp_mod.F90`)
pleaseRun(`ln -s ../../IOM/src/CESM_driver/ProgramTunnel .`)

cd(joinpath("..", ".."))

if parsed["build"]
    pleaseRun(`./$(parsed["casename"]).build`)
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM