text:
using TimeSeriesCausality
using Test
using Statistics: mean
using Distributions: MvNormal

@time @testset "granger_estimation.jl" begin
    time_steps = 256 * 1024
    segment_length = 1024
    noise_cov = MvNormal([0.25 0.0; 0.0 0.64])  # uncorrelated cov matrix of noise
    noise = rand(noise_cov, time_steps)'

    # initial values
    signal = zeros(time_steps, 2)
    signal[1:3, :] = rand(3, 2) + noise[1:3, :]

    # simulation
    designer = [0.4 -0.6  0.8  0.0  0.0  0.0;
                0.5  0.9  0.0  0.0  0.0  0.7]

    for t in 4:time_steps
        signal[t, :] = designer * reshape(signal[t-3:t-1, :], :, 1) + noise[t, :]
    end

    covariances = TimeSeriesCausality.est_sig_covs(signal, 3)
    ar_factors, p_error = TimeSeriesCausality.mvar_est(covariances)

    begin
        @test all(isapprox.(ar_factors[:, 1], designer[:, 3], atol = 0.01))
        @test all(isapprox.(ar_factors[:, 2], designer[:, 6], atol = 0.01))
        @test all(isapprox.(ar_factors[:, 3], designer[:, 2], atol = 0.01))
        @test all(isapprox.(ar_factors[:, 4], designer[:, 5], atol = 0.01))
        @test all(isapprox.(ar_factors[:, 5], designer[:, 1], atol = 0.01))
        @test all(isapprox.(ar_factors[:, 6], designer[:, 4], atol = 0.01))
    end

    granger_idx, std_error = granger_est(signal, segment_length; order = 3, method = "jackknife")
    @test granger_idx ≈ 0.59 atol = 0.02
    @test std_error ≈ 0.0035 atol = 0.002

    order_range = 1:7
    best_aic = argmin(granger_aic(signal, segment_length, order_range))
    @test best_aic == 3
    best_bic = argmin(granger_bic(signal, segment_length, order_range))
    @test best_bic == 3
end

@time @testset "utils.jl" begin

    # tests of int ####################################################
    @test TimeSeriesCausality.int(3.14) == 3
    @test TimeSeriesCausality.int(-2.72) == -2
    @test TimeSeriesCausality.int(0.0) == 0

    # tests of dropmean ###############################################
    test_a = rand(13)
    @test TimeSeriesCausality.dropmean(test_a, 1) == mean(test_a; dims = 1)
    test_a = rand(3, 5)
    @test TimeSeriesCausality.dropmean(test_a, 1) == mean(test_a; dims = 1)[1, :]
    test_a = rand(7, 11)
    @test TimeSeriesCausality.dropmean(test_a, 2) == mean(test_a; dims = 2)[:, 1]

    # tests of squeeze ################################################
    test_a = rand(17)
    @test TimeSeriesCausality.squeeze(test_a) == test_a
    test_a = rand(3, 5, 7)
    @test TimeSeriesCausality.squeeze(test_a) == test_a
    test_a = rand(3, 1, 7)
    @test TimeSeriesCausality.squeeze(test_a) == test_a[:, 1, :]

    # tests of detrend ################################################
    # testing the inplace
    test_x = [0:0.01:2.0;]
    zeros_x = zeros(size(test_x))
    @test all(isapprox.(TimeSeriesCausality.detrend!(test_x, 1), zeros_x, atol = 0.001))
    @test all(isapprox.(test_x, zeros_x, atol = 0.001))

    test_x = [0:0.01:2.0;]
    @test all(isapprox.(TimeSeriesCausality.detrend!(test_x, 0), [-1.0:0.01:1.0;], atol = 0.001))
    test_x = [0:0.01:2.0;]
    @test all(isapprox.(TimeSeriesCausality.detrend!(test_x, 1), zeros_x, atol = 0.001))

    test_x = [0:0.01:(2pi);]
    test_y = sin.(test_x)
    @test all(isapprox.(TimeSeriesCausality.detrend!(test_y, 0), test_y, atol = 0.001))
    test_x = [0:0.01:(2pi);]
    test_y = sin.(test_x)
    @test all(isapprox.(TimeSeriesCausality.detrend!(test_y, 1), test_y, atol = 0.001))

    test_x = [0:0.01:(2pi);]
    test_y = cos.(test_x) .+ test_x
    aux_x = [0:0.01:(2pi);] .- pi
    aux_y = cos.(test_x) .+ aux_x
    @test all(isapprox.(TimeSeriesCausality.detrend!(test_y, 0), aux_y, atol = 0.001))
    test_x = [0:0.01:(2pi);]
    test_y = cos.(test_x) .+ test_x
    @test all(isapprox.(TimeSeriesCausality.detrend!(test_y, 1), cos.(test_x), atol = 0.01))

    # tests of hanning_fun ############################################
    # the arrays are from MATLAB `hanning` function.
    hann_12 = [
        0.057271987173395
        0.215967626634422
        0.439731659872338
        0.677302443521268
        0.874255374085551
        0.985470908713026
        0.985470908713026
        0.874255374085551
        0.677302443521268
        0.439731659872338
        0.215967626634422
        0.057271987173395
    ]
    hann_13 = [
        0.049515566048790
        0.188255099070633
        0.388739533021843
        0.611260466978157
        0.811744900929367
        0.950484433951210
        1.000000000000000
        0.950484433951210
        0.811744900929367
        0.611260466978157
        0.388739533021843
        0.188255099070633
        0.049515566048790
    ]
    hanning_window_12 = TimeSeriesCausality.hanning_fun(12)
    hanning_window_13 = TimeSeriesCausality.hanning_fun(13)
    @test all(hanning_window_12 .≈ hann_12)
    @test all(hanning_window_13 .≈ hann_13)
    @test hanning_window_12[1] == hanning_window_12[end]
    @test hanning_window_13[1] == hanning_window_13[end]
    @test hanning_window_12[1] > 0.0
    @test hanning_window_13[1] > 0.0

end
@time @testset "phase_slope_index.jl" begin
    # tests of psi_est ###############################################
    # two random signals
    signal = [[randn(1000000);] [randn(1000000);]]
    psi, _ = psi_est(signal, 100; subave = true, segave = true)
    @test all(isapprox(psi, zeros(2, 2); atol = 0.01))

    # induced causality should always be inferred independent of optional arguments
    for eplen_ in [0, 200], detrend_ in [true, false]
        for segave_ in [true, false], subave_ in [true, false]
            for method_ in ["jackknife", "bootstrap"]
                # channel 2 being exactly channel 1
                ch1_ = randn(1000)
                signal = [[ch1_;] [ch1_;]]
                psi, _ = psi_est(
                    signal,
                    100;
                    method = method_,
                    eplen = eplen_,
                    detrend = detrend_,
                    segave = segave_,
                    subave = subave_,
                )
                @test psi == zeros(2, 2)

                # channel 1 leading channel 2
                ch1_ = randn(100000)
                signal = [[ch1_[2:end];] [ch1_[1:(end-1)];]]
                psi, _ = psi_est(
                    signal,
                    100;
                    method = method_,
                    eplen = eplen_,
                    detrend = detrend_,
                    segave = segave_,
                    subave = subave_,
                )
                @test psi[1, 1] == 0.0 && psi[2, 2] == 0.0
                @test psi[1, 2] > 1.0 && psi[2, 1] + psi[1, 2] == 0.0

                # channel 2 leading channel 1
                ch2_ = randn(100000)
                signal = [[ch2_[1:(end-1)];] [ch2_[2:end];]]
                psi, _ = psi_est(
                    signal,
                    100;
                    method = method_,
                    eplen = eplen_,
                    detrend = detrend_,
                    segave = segave_,
                    subave = subave_,
                )
                @test psi[1, 1] == 0.0 && psi[2, 2] == 0.0
                @test psi[1, 2] < -1.0 && psi[2, 1] + psi[1, 2] == 0.0
            end
        end
    end

    # test freqlist ###################################################
    ch1_ = randn(100000)
    signal = [[ch1_[2:end];] [ch1_[1:(end-1)];]]
    psi_range, _ = psi_est(signal, 100; freqlist = 1:49)
    psi_default, _ = psi_est(signal, 100)  # default is based on seglen
    @test all(psi_range == psi_default)

    # auxiliary function for testing
    """
        sin_sum = v_sin(X, A, Ω, Φ)

    returns an array summed over series of Sin waves with A-amplitudes,
        Ω-frequencies and Φ-phase delay
    """
    v_sin(X, A, Ω, Φ) = sum([a .* sin.(ω .* X .+ ϕ) for (a, ω, ϕ) in zip(A, Ω, Φ)])

    nsamples_ = 200000
    fs_ = 100  # sampling frequency (Hz)
    f_range_ = 12:30  # beta frequency band [12-30] Hz
    Ω_ = f_range_ * 2pi
    X_ = range(0.0, nsamples_ / fs_; length = nsamples_)
    A_ = rand(0.1:0.01:1.0, size(f_range_, 1))
    Φ_ = rand(size(f_range_, 1))
    Y_ = v_sin(X_, A_, Ω_, Φ_)
    Y_ch1, Y_ch2 = Y_[5:end], Y_[1:(end-4)]
    noise_ = randn(size(Y_ch1))
    signal = [Y_ch1 Y_ch2] .+ noise_
    seglen_ = fs_
    psi_full_range, _ = psi_est(signal, seglen_)
    psi_low, _ = psi_est(signal, seglen_; freqlist = 1:11)
    psi_freq, _ = psi_est(signal, seglen_; freqlist = f_range_)
    psi_high, _ = psi_est(signal, seglen_; freqlist = 31:49)
    @test psi_full_range[1, 2] > 1.0  # ch1 leading ch2
    @test psi_freq[1, 2] > psi_full_range[1, 2]  # PSI higher in target frequency range
    @test psi_freq[1, 2] > psi_low[1, 2]  # PSI higher in target frequency range
    @test psi_freq[1, 2] > psi_high[1, 2]  # PSI higher in target frequency range

    # ndims(data) should be 2
    signal = randn(100, 3, 2)
    @test_throws DimensionMismatch psi_est(signal, 100)

    # squeeze
    signal = randn(100, 1, 4)
    @test psi_est(signal, 100)[1] == psi_est(signal[:, 1, :], 100)[1]

    # size(data, 1) >!> seglen
    signal = randn(50, 3)
    @test_throws DimensionMismatch psi_est(signal, 100)

    # if eplen == 0 then no std estimation
    signal = [[randn(1000000);] [randn(1000000);]]
    _, psi_std = psi_est(signal, 100)
    @test all(isnan.(psi_std))

    # if continuous data, subave shall not be true otherwise NaN is returned
    signal = [[randn(1000000);] [randn(1000000);]]
    psi, _ = psi_est(signal, 100; subave = true)
    @test !all(isnan.(psi))

end


text:
if !isdefined(@__MODULE__, :IoK8sApimachineryPkgUtilIntstrIntOrString)
    const IoK8sApimachineryPkgUtilIntstrIntOrString = String
else
    @warn("Skipping redefinition of IoK8sApimachineryPkgUtilIntstrIntOrString to String")
end

convert(::Type{IoK8sApimachineryPkgUtilIntstrIntOrString}, v::T) where {T<:Integer} = string(v)

ext: jl


text:
################################################################################
 # Impulses
################################################################################
# see joints/joint.jl method

################################################################################
 # Derivatives
################################################################################
function impulse_transform_jacobian(relative::Symbol, jacobian::Symbol,
    joint::Rotational{T,Nλ},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion, p; attjac=true) where {T,Nλ}

    if relative == :parent
        if jacobian == :parent
            # ∂(Ja'*p)/∂(xa,qa)
            Z3 = szeros(T,3,3)
            ∂Q∂qa = ∂VLᵀmat∂q(Tmat() * Rᵀmat(qb) * LVᵀmat(joint.orientation_offset) * p) * LVᵀmat(qa)
            return cat(I(3), 0.5 * I(3), dims=(1,2)) * [Z3 Z3; Z3 ∂Q∂qa]
        elseif jacobian == :child
            # ∂(Ja'*p)/∂(xb,qb)
            Z3 = szeros(T,3,3)
            ∇Qqb = VLᵀmat(qa) * Tmat(T) * ∂Rᵀmat∂q(LVᵀmat(joint.orientation_offset) * p) * LVᵀmat(qb)
            return cat(I(3), 0.5 * I(3), dims=(1,2)) * [Z3 Z3; Z3 ∇Qqb]
        end
    elseif relative == :child
        if jacobian == :parent
            # ∂(Jb'*p)/∂(xa,qa)
            Z3 = szeros(T,3,3)
            ∇Qqa = VLᵀmat(qb) * ∂Lmat∂q(LVᵀmat(joint.orientation_offset) * p) * LVᵀmat(qa)
            return cat(I(3), 0.5 * I(3), dims=(1,2)) * [Z3 Z3; Z3 ∇Qqa]
        elseif jacobian == :child
        # ∂(Jb'*p)/∂(xb,qb)
            Z3 = szeros(T,3,3)
            ∇Qqb = ∂VLᵀmat∂q(Lmat(qa) * LVᵀmat(joint.orientation_offset) * p) * LVᵀmat(qb)
            return cat(I(3), 0.5 * I(3), dims=(1,2)) * [Z3 Z3; Z3 ∇Qqb]
        end
    end
end


text:
using LinearAlgebra, Statistics, Printf
# commandline arguments.  In reality, use ArgParse equivalents
#
num = parse(Int64,ARGS[2])
ϵ = randn(num)
m = mean(ϵ)
s = std(ϵ)
println("mean = $m, std = $s of $num values")

open("myout.txt","w") do f
    write(f, @sprintf("mean = %s, std = %s of %s\n", m, s, num))
end

i = 0
while i <= 60
        println(i)
        sleep(1)
        global i += 1
    end


text:
"""
    f = get_inb(pedlist [,check, sorthere, method])
    f = get_inb(s, d [,check, sorthere, method])
    f = get_inb(Tv, pedlist [,check, sorthere, method])

Computes inbreeding coefficients for animals in a pedigree list `pedlist`.
The list is an integer matrix (2 x *n*, where *n* is the number of animals).
Or, you can provide the lists of sires and dams (`s` and `d`, respectively).
The results are in the vector `f` with the type `Tv` (Float64 by default).
Unknown parent groups (UPGs) will be replaced with a missing code (`0`).

The pedigree should be coded as parents precede their progeny.
All descendent mush have a larger code than their ancestors.
By default, the function checks if the pedigree is "permutaed", and stops if it is wrong.
You can omit the check by `check=false`.
If you want to permute the pedigree chronologically within this function, put `sorthere=true`.

The default method is Meuwissen and Luo (1992).
This is the only method to be supported in this version.

### Examples
```jldoctest
julia> using PedigreeTools;

julia> pedlist,idtable = read_ped("pedigree.txt");

# assuming animals are chronologically sorted.
julia> f = get_inb(pedlist);

# sort animals here.
julia> f = get_inb(pedlist, sorthere=true);
```
"""
function get_inb(Tv::DataType, pedlist::Matrix{Ti};
                 check=true, sorthere=false, method::String="MeuwissenAndLuo") where Ti<:Integer
   # make a copy of pedigree
   ped = normalize_pedigree(pedlist)

   # pedigree checks
   if check
      if sorthere
         # being sorted so no checks for `parentsfirst`
         if !check_ped(ped)
            error("error in pedigree")
         end
      else
         # assuming to be already sorted
         if !check_ped(ped,parentsfirst=true)
            error("error in pedigree, or not sorted")
         end
      end
   end
   # computations
   if method=="MeuwissenAndLuo"
      # This algorithm rewrites the pedigree list so
      # we use the copy of it to keep the original unchanged.
      if sorthere
         perm,invp = find_ped_order(ped)
         permute_ped!(invp,ped)
         f = kernel_meuwissen_and_luo!(Tv,ped)
         f .= f[invp]
         return f
      else
         return kernel_meuwissen_and_luo!(Tv,ped)
      end
      #return get_inb_ml(Tv,pedlist)
   else
      throw(ArgumentError("unsupported method for inbreeding computations"))
   end
end

# default = Float64
function get_inb(pedlist::Matrix{Ti}, Tv::DataType=Float64;
                 check=false, sorthere=false, method::String="MeuwissenAndLuo") where Ti<:Integer
   return get_inb(Tv, pedlist, check=check, sorthere=sorthere, method=method)
end

# pedigree as vectors
function get_inb(s::Vector{Ti}, d::Vector{Ti}, Tv::DataType=Float64;
                 check=false, sorthere=false, method::String="MeuwissenAndLuo") where Ti<:Integer
   pedlist = [s'; d']
   return get_inb(Tv, pedlist, check=check, sorthere=sorthere, method=method)
end

# normalize pedigree
# replace a negative/large ID with 0
function normalize_pedigree(pedlist::Matrix{Ti}) where {Ti<:Integer}
   n::Ti = size(pedlist,2)
   ped = zeros(Ti,2,n)
   sid::Ti = 0
   did::Ti = 0
   # normalize pedigree; replace negative/large id with 0
   for i=1:n
      sid = pedlist[1,i]
      did = pedlist[2,i]
      if sid<1 || sid>n
         sid = 0
      end
      if did<1 || did>n
         did = 0
      end
      ped[1,i] = sid
      ped[2,i] = did
   end
   return ped
end

# kernel of Meuwissen and Luo (1992)
# the argument "ped" being rewritten
function kernel_meuwissen_and_luo!(Tv::DataType, ped::Matrix{Ti}) where Ti<:Integer
   n::Ti = size(ped,2)
   #f = OffsetArray{Tv}(0:n)
   f = OffsetArray{Tv}(undef,0:n)
   i::Ti = 0
   j::Ti = 0
   k::Ti = 0
   s0::Ti = 0
   d0::Ti = 0
   ks::Ti = 0
   kd::Ti = 0
   fi::Tv = 0.0
   r::Tv = 0.0
   point = zeros(Ti,n)
   T = zeros(Tv,n)
   B = zeros(Tv,n)

   # start
   f[0] = -1
   for i=1:n
      s0 = ped[1,i]
      d0 = ped[2,i]
      ped[1,i] = max(s0,d0)
      ped[2,i] = min(s0,d0)
      B[i] = 0.5 - 0.25*(f[s0]+f[d0])

      if s0==0 || d0==0
         f[i] = 0.0
      else
         fi = -1.0
         T[i] = 1.0
         j = i
         while(j!=0)
            k = j
            r = 0.5 * T[k]
            ks = ped[1,k]
            kd = ped[2,k]
            if ks != 0
               while(point[k]>ks)
                  k = point[k]
               end
               T[ks] = T[ks] + r
               if ks != point[k]
                  point[ks] = point[k]
                  point[k] = ks
               end
               if kd != 0
                  while(point[k]>kd)
                     k = point[k]
                  end
                  T[kd] = T[kd] + r
                  if kd != point[k]
                     point[kd] = point[k]
                     point[k] = kd
                  end
               end
            end
            fi = fi + T[j]*T[j]*B[j]
            T[j] = 0.0
            k = j
            j = point[j]
            point[k] = 0
         end
         f[i] = fi
      end
   end
   return f[1:n]
end

"""
    c = get_inbupg_code(s, d, f)
    c = get_inbupg_code(pedlist, f)

Returns the "inb/upg code" out of sire ID `s`, dam ID `d`, and a list of inbreeding coefficient `f`.
Instead of single ID, you can provide lists of sires and dams (unified list `pedlist` or separate lists `s` and `d`) to obtain the code for all individuals.
The code is integer, and it is useful only for the BLUPF90 programs.
"""
function get_inbupg_code(s::Int,d::Int,f::Vector{Float64})
   if s>0
      ms = 0.0
      fs = f[s]
   else
      ms = 1.0
      fs = 0.0
   end
   if d>0
      md = 0.0
      fd = f[d]
   else
      md = 1.0
      fd = 0.0
   end
   return Int(round( 4000/((1+ms)*(1-fs) + (1+md)*(1-fd)) ))
end

function get_inbupg_code(pedlist::Matrix{Int},f::Vector{Float64})
   if size(pedlist,2)<size(f,1)
      throw(DimensionMismatch("pedlist and inbreeding lists"))
   end
   n = size(pedlist,2)
   code = zeros(Int,n)
   for j=1:n
      code[j] = get_inbupg_code(pedlist[1,j],pedlist[2,j],f)
   end
   return code
end

function get_inbupg_code(s::Vector{Int},d::Vector{Int},f::Vector{Float64})
   if size(s,1)!=size(d,1) || size(s,1)<size(f,1)
      throw(DimensionMismatch("sire, dam, and inbreeding lists"))
   end
   pedlist = [s'; d']
   return get_inbupg_code(pedlist,f)
end


text:
abstract MarkovNetSampler

"""
The MarkovNetGibbsSampler type houses the parameters of the Gibbs sampling algorithm.  The parameters are defined below:

burn_in:  The first burn_in samples will be discarded.  They will not be returned.
The thinning parameter does not affect the burn in period.
This is used to ensure that the Gibbs sampler converges to the target stationary distribution before actual samples are drawn.

thinning: For every thinning + 1 number of samples drawn, only the last is kept.
Thinning is used to reduce autocorrelation between samples.
Thinning is not used during the burn in period.
e.g. If thinning is 1, samples will be drawn in groups of two and only the second sample will be in the output.

evidence: the assignment that all samples must be consistent with (ie, Assignment(:A=>1) means all samples must have :A=1).
Use to sample conditional distributions.

initial_sample:  The inital assignment to variables to use.  If null, the initial sample is chosen at random
"""
type MNGibbsSampler <: MarkovNetSampler

    evidence::Assignment
    burn_in::Int
    thinning::Int
    initial_sample::Nullable{Assignment}

    function MNGibbsSampler(evidence::Assignment=Assignment();
                             burn_in::Int=100,
                             thinning::Int=0,
                             initial_sample::Nullable{Assignment}=Nullable{Assignment}()
                             )
        new(evidence, burn_in, thinning, initial_sample)
    end
end

"""
Implements Gibbs sampling for MarkovNets.
This Gibbs sample only supports discrete MarkovNets, and samples are
drawn following a Categorical Distribution with probabilities
equal to the normalized potentials

Sampling requires an MarkovNetGibbsSampler object which contains the parameters
"""
function Base.rand(mn::MarkovNet, sampler::MNGibbsSampler, nsamples::Integer)

    return gibbs_sample(mn, nsamples, sampler.burn_in, thinning=sampler.thinning, evidence=sampler.evidence, initial_sample=sampler.initial_sample)
end

function gibbs_sample(mn::MarkovNet, nsamples::Integer, burn_in::Integer;
                      thinning::Integer=0,
                      evidence::Assignment=Assignment(),
                      initial_sample::Nullable{Assignment}=Nullable{Assignment}()
                      )
    # Check parameters for correctness
    nsamples > 0 || throw(ArgumentError("nsamples parameter less than 1"))
    burn_in >= 0 || throw(ArgumentError("Negative burn_in parameter"))
    if ~ isnull(initial_sample)
        init_sample = get(initial_sample)
        for vertex in vertices(mn.ug)
            haskey(init_sample, Symbol(vertex)) || throw(ArgumentError("Gibbs sample initial_sample must be an assignment with all variables in the Bayes Net"))
        end
        for (vertex, value) in evidence
            init_sample[vertex] == value || throw(ArgumentError("Gibbs sample initial_sample was inconsistent with evidence"))
        end
    end

    # reduce the MarkovNet according to the evidence:
    mn_reduced = evidence_reduce(mn, evidence)

    if isnull(initial_sample)
        # Hacky
        initial_sample = Assignment()
        for factor in mn_reduced.factors
            for (i, vertex) in enumerate(factor.dimensions)
                n_categories = length(factor.potential[:, i])
                initial_sample[vertex.name] = rand(Categorical(n_categories))
            end
        end
    end

    # create the data frame
    t = Dict{Symbol, Vector{Any}}()
    for name in keys(initial_sample)
        t[name] = Any[]
    end

    # initialize the sample to our initial sample
    current_sample = initial_sample

    # burn in, if present
    if burn_in != 0
        for burn_in_sample in 1:burn_in
            current_sample = gibbs_sample_loop(mn_reduced, current_sample)
        end
    end

    # main loop
    for sample_iter in 1:nsamples

        # first skip over the thinning
        for skip_iter in 1:thinning
            current_sample = gibbs_sample_loop(mn_reduced, current_sample)
        end

        # real loop, we store the values in the dict
        current_sample = gibbs_sample_loop(mn_reduced, current_sample, t)

    end

    return DataFrame(t)
end

function gibbs_sample_loop(mn::MarkovNet, current_sample::Assignment,
                           results::Dict{Symbol, Vector{Any}} = Dict{Symbol, Vector{Any}}())::Assignment
    sample::Assignment = Assignment()
    other_assg::Assignment = Assignment()

    for v in vertices(mn.ug)

        other_assg = copy(current_sample) # this is faster than filter
        delete!(other_assg, Symbol(v))
        new_f = reduce(*, [mn.factors[factor_index][other_assg] for factor_index in mn.name_to_factor_indices[Symbol(v)]])
        proba = new_f.potential[:]
        proba = proba/sum(proba)
        sample_index = rand(Categorical(proba))
        current_sample[Symbol(v)] = sample_index
        if ~ isempty(results)
            push!(results[Symbol(v)], current_sample[Symbol(v)])
        end
    end
    return current_sample
end


text:
module TestFoldsDagger

using Test

for file in readdir(@__DIR__)
    if match(r"^test_.*\.jl$", file) !== nothing
        include(joinpath(@__DIR__, file))
    end
end

end # module


text:
using ..Serde: QuantizedPayloadSerde

"""
    QuantizedConfig{T <: Unsigned}(base::BaseConfig) <: Configuration

Configuration for compression through quantization. The quantization type needs
to be speficied when the configuration is instanciated. This configuration
instanciates a  [QuantizedPayloadSerde](@ref).
"""
struct QuantizedConfig{T <: Unsigned} <: Configuration
    base::BaseConfig
    payload_serde::QuantizedPayloadSerde{T}

    QuantizedConfig{T}(base::BaseConfig) where T <: Unsigned =
        new(base, QuantizedPayloadSerde{T}())
end


text:
using LinearAlgebra
using Makie

μ1 = μ2 = 0.0

𝛍 = [
    μ1
    μ2
]

σ11 = σ22 = 1.0
σ12 = σ21 = 0.0

Σ = [
    σ11 σ12
    σ21 σ22
]

function 𝒩(x1, x2; 𝛍 = 𝛍, Σ = Σ)
    𝐱 = [x1, x2]
    Λ = inv(Σ)
    D = length(𝛍)
    return √(det(Λ)) * exp(-0.5 * (𝐱 .- 𝛍)' * Λ * (𝐱 .- 𝛍)) / √(2π)^D
end

x1min, x1max = -1.0, 1.0
x2min, x2max = -1.0, 1.0
x1 = range(x1min, x1max, length = 50)
x2 = range(x2min, x2max, length = 50)
z = 𝒩.(x1', x2)
wireframe(x1, x2, z, colormap = :Spectral)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# Test heat equation agrees with analytic solution to problem 55 on page 28 in https://ocw.mit.edu/courses/mathematics/18-303-linear-partial-differential-equations-fall-2006/lecture-notes/heateqni.pdf
using MPI
using OrderedCollections
using StaticArrays
using Statistics
using Dierckx
using Test

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations
using ClimateMachine.Land.SoilHeatParameterizations
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.DGMethods: BalanceLaw, LocalGeometry
using ClimateMachine.MPIStateArrays
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.BalanceLaws

function init_soil!(land, state, aux, localgeo, time)
    FT = eltype(state)
    ϑ_l, θ_i = get_water_content(land.soil.water, aux, state, time)
    θ_l = volumetric_liquid_fraction(ϑ_l, land.soil.param_functions.porosity)
    param_set = parameter_set(land)
    ρc_s = volumetric_heat_capacity(
        θ_l,
        θ_i,
        land.soil.param_functions.ρc_ds,
        param_set,
    )

    state.soil.heat.ρe_int = FT(volumetric_internal_energy(
        θ_i,
        ρc_s,
        land.soil.heat.initialT(aux),
        param_set,
    ))
end

const FT = Float64

@testset "Heat analytic unit test" begin
    ClimateMachine.init()

    soil_param_functions = SoilParamFunctions(
        FT;
        porosity = 0.495,
        ν_ss_gravel = 0.1,
        ν_ss_om = 0.1,
        ν_ss_quartz = 0.1,
        ρc_ds = 0.43314518988433487,
        κ_solid = 8.0,
        ρp = 2700.0,
        κ_sat_unfrozen = 0.57,
        κ_sat_frozen = 2.29,
    )

    heat_surface_state = (aux, t) -> eltype(aux)(0.0)

    tau = FT(1) # period (sec)
    A = FT(5) # amplitude (K)
    ω = FT(2 * pi / tau)
    heat_bottom_state = (aux, t) -> A * cos(ω * t)
    T_init = (aux) -> eltype(aux)(0.0)

    soil_water_model = PrescribedWaterModel(
        (aux, t) -> eltype(aux)(0.0),
        (aux, t) -> eltype(aux)(0.0),
    )

    bc = LandDomainBC(
        bottom_bc = LandComponentBC(soil_heat = Dirichlet(heat_bottom_state)),
        surface_bc = LandComponentBC(soil_heat = Dirichlet(heat_surface_state)),
    )
    soil_heat_model = SoilHeatModel(FT; initialT = T_init)

    m_soil = SoilModel(soil_param_functions, soil_water_model, soil_heat_model)
    sources = ()
    m = LandModel(
        param_set,
        m_soil;
        boundary_conditions = bc,
        source = sources,
        init_state_prognostic = init_soil!,
    )

    N_poly = 5
    nelem_vert = 10

    # Specify the domain boundaries
    zmax = FT(1)
    zmin = FT(0)

    driver_config = ClimateMachine.SingleStackConfiguration(
        "LandModel",
        N_poly,
        nelem_vert,
        zmax,
        param_set,
        m;
        zmin = zmin,
        numerical_flux_first_order = CentralNumericalFluxFirstOrder(),
    )

    t0 = FT(0)
    timeend = FT(2)
    dt = FT(1e-4)
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )
    mygrid = solver_config.dg.grid
    aux = solver_config.dg.state_auxiliary
    ClimateMachine.invoke!(solver_config)
    t = ODESolvers.gettime(solver_config.solver)

    z_ind = varsindex(vars_state(m, Auxiliary(), FT), :z)
    z = Array(aux[:, z_ind, :][:])

    T_ind = varsindex(vars_state(m, Auxiliary(), FT), :soil, :heat, :T)
    T = Array(aux[:, T_ind, :][:])

    num =
        exp.(sqrt(ω / 2) * (1 + im) * (1 .- z)) .-
        exp.(-sqrt(ω / 2) * (1 + im) * (1 .- z))
    denom = exp(sqrt(ω / 2) * (1 + im)) - exp.(-sqrt(ω / 2) * (1 + im))
    analytic_soln = real(num .* A * exp(im * ω * timeend) / denom)
    MSE = mean((analytic_soln .- T) .^ 2.0)
    @test eltype(aux) == FT
    @test MSE < 1e-5
end


text:
# Key methods to add for a new dimensional data type

"""
    dims(x)

Return a tuple of the dimensions for a dataset. These can
contain the coordinate ranges if `bounds()` and `select()` are to be used,
or you want them to be shown on plots in place of the array indices.

They can also contain a units string or unitful unit to use and plot
dimension units.

This is the only method required for this package to work. It probably
requires defining a dims field on your object to store dims in.
"""
function dims end
dims(x::T) where T = error("`dims` not defined for type $T")
dims(x::Nothing) = nothing

"""
    refdims(x)

Reference dimensions for an array that is a slice or view of another
array with more dimensions.

`slicedims(a, dims)` returns a tuple containing the current new dimensions
and the new reference dimensions. Refdims can be stored in a field or disgarded,
as it is mostly to give context to plots. Ignoring refdims will simply leave some captions empty.  """
function refdims end
refdims(x) = ()
"""
    rebuild(x::AbstractDimensionalArray, data, [dims], [refdims])
    rebuild(x::AbstractDimension, val, [grid], [metadata])
    rebuild(x; kwargs...)

Rebuild an object struct with updated values.
"""
function rebuild end
rebuild(x; kwargs...) = ConstructionBase.setproperties(x, (;kwargs...))

"""
    val(x)

Return the contained value of a wrapper object, otherwise just returns the object.
"""
function val end
val(x) = x

"""
    metadata(x)

Return the metadata of a dimension or data object.
"""
function metadata end

"""
    units(x)

Return the units of a dimension. This could be a string, a unitful unit, or nothing.
"""
function units end
units(x) = nothing
units(xs::Tuple) = map(units, xs)

"""
    name(x)

Get the name of data or a dimension.
"""
function name end
name(x) = name(typeof(x))
name(x::Type) = ""
name(xs::Tuple) = map(name, xs)

"""
    shortname(x)

Get the short name of array data or a dimension.
"""
function shortname end
shortname(x) = shortname(typeof(x))
shortname(xs::Tuple) = map(shortname, xs)
shortname(x::Type) = ""

"""
    label(x)

Get a plot label for data or a dimension. This will include the name and units
if they exist, and anything else that should be shown on a plot.
"""
function label end
label(x) = string(name(x), (units(x) === nothing ? "" : string(" ", units(x))))
label(xs::Tuple) = join(map(label, xs), ", ")


text:
"""
    plot_field(femesh, ξ)

Plot field `ξ` on the finite element mesh.
"""
function plot_field(femesh::FEMesh, ξ)
    ncompartment, nboundary = size(femesh.facets)
    ξ_cmpts = split_field(femesh, ξ)
    scene = nothing
    first = true
    for icmpt = 1:ncompartment, iboundary = 1:nboundary
        facets = femesh.facets[icmpt, iboundary]
        points = femesh.points[icmpt]
        colors = abs.(ξ_cmpts[icmpt])
        # colors = fill(1.0 * iboundary, length(colors))
        # colors = iboundary / nboundary
        if first
            scene = mesh(points', facets', color = colors, shading = false)
            first = false
        else
            mesh!(points', facets', color = colors, shading = false)
        end
    end
    scene
end


text:
using Test
import ArchGDAL; const AG = ArchGDAL
import DataStreams

@testset "DataStream Support" begin
    df = AG.read("data/point.geojson") do dataset
        DataStreams.Data.close!(DataStreams.Data.stream!(
            AG.Source(AG.getlayer(dataset,0)), DataStreams.Data.Table
        ))
    end
    @test df.FID == [2.0, 3.0, 0.0, 3.0]
    @test df.pointname == ["point-a", "point-b", "a", "b"]
    @test AG.toWKT.(df.geometry0) == [
        "POINT (100 0)",
        "POINT (100.2785 0.0893)",
        "POINT (100 0)",
        "POINT (100.2785 0.0893)"
    ]
end


text:
using Bolt
𝕡 = CosmoParams()
bg = Background(𝕡)
𝕣 = Bolt.RECFAST(bg=bg, Yp=𝕡.Y_p, OmegaB=𝕡.Ω_b)  #  𝕣 = Bolt.Peebles()
ih = IonizationHistory(𝕣, 𝕡, bg)

Nz = 1000
dz = 10000/float(Nz)
z = (10000 - dz):(-dz):0.0
##

using PyPlot
clf()
plot(z, ih.Tmat.(z2x.(z)), "-", label=raw"$T_{\mathrm{mat}}$")
plot(z, ih.Trad.(z2x.(z)), "--", label=raw"$T_{\mathrm{rad}}$")

yscale("log")
xscale("log")
legend()
ylim(1, 2e4)
xlim(10, 10000)
xlabel("redshift")
ylabel("temperature [K]")
gcf()


text:
addprocs(5)
################
#My local machine has 6 processors
#CYCLOPS is written to use them
#################
@everywhere basedir=homedir();
@everywhere netdir=string(basedir,"/Desktop/CYCLOPS_UPDATE_SHARE");
@everywhere cd(netdir);

using StatsBase
using MultivariateStats
using Distributions

@everywhere include("CYCLOPS_v6_2a_AutoEncoderModule_multi.jl")
@everywhere include("CYCLOPS_v6_2a_seed.jl")
@everywhere include("CYCLOPS_v6_2a_PreNPostprocessModule.jl")
@everywhere include("CYCLOPS_v6_2a_CircularStats_U.jl")
@everywhere include("CYCLOPS_v6_2a_MultiCoreModule_Smooth.jl")

using CYCLOPS_v6_2a_AutoEncoderModule_multi
using CYCLOPS_v6_2a_Seed
using CYCLOPS_v6_2a_PreNPostprocessModule
using CYCLOPS_v6_2a_CircularStats_U
using CYCLOPS_v6_2a_MultiCoreModule_Smooth

using PyPlot
#indirlaval=string(basedir,"/Documents/LungWork_March2015/LAVAL/UsableData");
#indirgrng=string(basedir,"/Documents/LungWork_March2015/GRNG/UsableData");
#homologuedir=string(basedir,"/Google Drive/BHTC_Homologues")

indirlaval=netdir;
indirgrng=netdir;
homologuedir=netdir;

############################################################
Frac_Var=0.85 # Set Number of Dimensions of SVD to maintain this fraction of variance
DFrac_Var=0.03 # Set Number of Dimensions of SVD to so that incremetal fraction of variance of var is at least this much
N_best =10 # Number of random initial conditions to try for each optimization

total_background_num=10; # Number of background runs for global background refrence distribution (bootsrap) for reall runs should be much higher
n_cores=5; # Number of machine cores
############################################################

cd(homologuedir)
seed_homologues1=readcsv("Human_UbiquityCyclers.csv");
homologue_symbol_list1=seed_homologues1[2:end,2];
###########################################################
Seed_MinCV 				= 0.14
Seed_MaxCV 				= .7
Seed_Blunt				=.975
MaxSeeds                =10000
###########################################################

srand(12345);

cd(indir);
fullnonseed_data=readcsv("Annotated_Unlogged_BA11Data.csv");

alldata_probes=fullnonseed_data[4:end ,1];
alldata_symbols=fullnonseed_data[4:end ,2];

alldata_times=fullnonseed_data[3,4:end];
alldata_subjects=fullnonseed_data[2,4:end];
alldata_samples=fullnonseed_data[1,4:end];

alldata_data=fullnonseed_data[4:end , 4:end];
alldata_data=Array{Float64}(alldata_data);

n_samples=length(alldata_times)
n_probes=length(alldata_probes)
timestamped_samples=setdiff(1:n_samples,findin(alldata_times,["NA"]))
cutrank=n_probes-MaxSeeds

Seed_MinMean										= (sort(vec(mean(alldata_data,2))))[cutrank]

##################################
# This extracts the genes from the dataset that were felt to have a high likelyhood to be cycling - and also had a reasonable coefficient of variation in this data sets
##################################
seed_symbols1, seed_data1 							= getseed_homologuesymbol_brain(fullnonseed_data,homologue_symbol_list1,Seed_MaxCV,Seed_MinCV,Seed_MinMean,Seed_Blunt);
seed_data1 											= dispersion!(seed_data1)
outs1, norm_seed_data1								= GetEigenGenes(seed_data1,Frac_Var,DFrac_Var,30)

##################################
# This example creates a "balanced autoencoder" where the eigengenes ~ principle components are encoded by a single phase angle
# This is the version of CYCLOPS published in 2017
##################################
estimated_phaselist1,bestnet1,global_var_metrics1 	= CYCLOPS_Order(outs1,norm_seed_data1,N_best)
estimated_phaselist1								= mod.(estimated_phaselist1 + 2*pi,2*pi)
estimated_magnitudes1								= ExtractMagnitude(norm_seed_data1,bestnet1)
###################################

###################################
#this is the command to evaluat the fit of the data (with a balanced autoencoder and only a single phase angle) to that we observe with shuffled data that does not (by design)
#have a circualar structure
#I suspect that some cross validation statistics are better suited - (especially if more complex network architectures are used)
#This is also very time/resource intensive

#I also think a cross validation error is the right way to choose the proper network complexity
#We can discuss this
###################################
#pvals=multicore_backgroundstatistics_global_eigen(seed_data1,outs1,N_best,total_background_num,global_metrics1)
####################################

truetimes=mod.(Array{Float64}(alldata_times[timestamped_samples]),24)

estimated_phaselist1aa=estimated_phaselist1[timestamped_samples]

shiftephaselist1=best_shift_cos(estimated_phaselist1aa,truetimes,"hours")

#shiftephaselist3=best_shift_cos(estimated_phaselist3a,truetimes,"hours")

close()
scatter(truetimes,shiftephaselist1,alpha=.75,s=14)
title("Eigengenes Encoded by Single Phase")
ylabp=[0,pi/2,pi,3*pi/2,2*pi]
ylabs=[0, "","π","","2π"]
xlabp=[0,6,12,18,24]
xlabs=["0", "6","12","18","24"]
ylabel("CYCLOPS Phase (radians)", fontsize=14)
xlabel("Hour of Death", fontsize=14)
xticks(xlabp, xlabs)
yticks(ylabp, ylabs)

suptitle("CYCLOPS Phase Prediction: Human Frontal Cortex", fontsize=18)

###################
#Asses Error
#####################

errors=Circular_Error_List(2*pi*truetimes/24,shiftephaselist1)
hrerrors=(12/pi)*abs.(errors)
mean(hrerrors)
median(hrerrors)
sqrt(var(hrerrors))
sort(hrerrors)[Integer(round(.75*length(hrerrors)))]

correlations=Jammalamadka_Circular_CorrelationMeasures(2*pi*truetimes/24,shiftephaselist1)

###################################################
clockgenes=["C1orf51","ARNTL","NR1D1","CLOCK","PER3"]
clockrows=findin(alldata_symbols,clockgenes);
clockrows=clockrows[[2,4,1,3,5]] # reorder for plotting

clockdata=alldata_data[clockrows,:];
clockannotations=alldata_symbols[clockrows,:];
clockannotations[1]="CHRONO"  # C1orf51 renamed CHRONO

r1=reshape(alldata_samples[timestamped_samples],1,length(timestamped_samples))
r3=alldata_data[:,timestamped_samples]

fulltimed_data1=vcat(r1,r3)
c1=vcat(["Probe"],alldata_probes)
c2=vcat(["Gene Symbol"],alldata_symbols)
c3=vcat(["Gene Symbol"],alldata_symbols)

fulltimed_data=hcat(c1,c2,c3,fulltimed_data1)

MultiCore_Cosinor_Statistics(r3,shiftephaselist1,24)
cosinor=Compile_MultiCore_Cosinor_Statistics(fulltimed_data,shiftephaselist1,4,24)
cosinor_clock=cosinor[clockrows+1,:];  ## the +1 is because the first row of this array is the headings

for n in 1:4
	m1=2*n-1
	m2=2n
	subplot(5,2,m1)
		scatter(truetimes,clockdata[n,timestamped_samples],alpha=.75,s=7,color="DarkRed")
		title(clockannotations[n])
		xlabp=[0,24]
		xlabs=["0","24"]
		ylabp=[]
		ylabs=[]
		xticks(xlabp, xlabs)
		yticks(ylabp, ylabs)

	subplot(5,2,m2)
        scatter(shiftephaselist1,clockdata[n,timestamped_samples],alpha=.75,s=7)
        title(clockannotations[n])
		xlabp=[0,2*pi]
		xlabs=[0,"2π"]
		ylabp=[]
		ylabs=[]
		xticks(xlabp, xlabs)
		yticks(ylabp, ylabs)
	end
n=5
m1=2*n-1
m2=2n
subplot(5,2,m1)
	scatter(truetimes,clockdata[n,timestamped_samples],alpha=.75,s=7,color="DarkRed")
	title(clockannotations[n])
	xlabp=[0,24]
	xlabs=["0","24"]
	ylabp=[]
	ylabs=[]
	xticks(xlabp, xlabs)
	yticks(ylabp, ylabs)
	xlabel("Hour of Death", fontsize=18,color="DarkRed")

subplot(5,2,m2)
    scatter(shiftephaselist1,clockdata[n,timestamped_samples],alpha=.75,s=7)
    title(clockannotations[n])
	xlabp=[0,2*pi]
	xlabs=[0,"2π"]
	ylabp=[]
	ylabs=[]
	xticks(xlabp, xlabs)
	yticks(ylabp, ylabs)
	xlabel("CYCLOPS Phase", fontsize=18,color="DarkBlue")

text(-4*pi,600,"Relative Expression",rotation="vertical",fontsize=19)

#suptitle("Gene Expression As A Function of TOD and CYCLOPS Phase ", fontsize=14, fontweight="bold")

############################

for n in 1:15
	m1=2*n-1
	m2=2n
	subplot(5,8,m1)
		scatter(truetimes,seed_data1[n,timestamped_samples],alpha=.75,s=7,color="DarKRed")
		title(seed_symbols1[n])
		xlabp=[0,24]
		xlabs=["0","24"]
		xticks(xlabp, xlabs)

	subplot(5,8,m2)
        scatter(shiftephaselist1,seed_data1[n,timestamped_samples],alpha=.75,s=7)
        title(seed_symbols1[n])
		xlabp=[0,2*pi]
		xlabs=[0,"2π"]
		yticks=[]
		ylabs=[]
		xticks(xlabp, xlabs)

	end
####################

title("Selected Clock Genes")

max_stars_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_issues_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_forks_repo_name: AdithyaVenkateshMohan/CYCLOPS


text:
module CARMAKepler

using ..Celerite
using Distributions
using Ensemble
using HDF5
using ..Kepler

import Base:
    write, read

struct MultiEpochPosterior
    ts::Array{Array{Float64, 1}, 1}
    ys::Array{Array{Float64, 1}, 1}
    dys::Array{Array{Float64, 1}, 1}

    inds::Array{Array{Int, 1}, 1}

    allts::Array{Float64, 1}
    allys::Array{Float64, 1}
    alldys::Array{Float64, 1}

    P_min::Float64
    P_max::Float64

    K_min::Float64
    K_max::Float64

    ndrw::Int
    nosc::Int

    rate_min::Float64
    rate_max::Float64

    f_min::Float64
    f_max::Float64

    rms_min::Float64
    rms_max::Float64

    Q_min::Float64
    Q_max::Float64
end

function MultiEpochPosterior(ts, ys, dys, per_min, per_max, ndrw, nosc, rate_min, rate_max, f_min, f_max, Q_max)
    allts = vcat(ts...)
    allys = vcat(ys...)
    alldys = vcat(dys...)

    allinds = sortperm(allts)

    iallinds = invperm(allinds)

    inds = []
    i = 1
    for t in ts
        push!(inds, iallinds[i:i+size(t,1)-1])
        i = i+size(t,1)
    end

    allts = allts[allinds]
    allys = allys[allinds]
    alldys = alldys[allinds]

    rms_max = maximum([std(y) for y in ys])
    rms_min = minimum([std(y) for y in ys])

    T = allts[end]-allts[1]
    dtmin = minimum(diff(allts))

    MultiEpochPosterior(ts, ys, dys, inds, allts, allys, alldys, per_min, per_max, rms_min/100.0, 10.0*rms_max, ndrw, nosc, rate_min, rate_max, f_min, f_max, rms_min/100.0, rms_max*10.0, 1.0, Q_max)
end

struct MultiEpochParams
    mu::Array{Float64, 1}
    nu::Array{Float64, 1}

    K::Float64
    P::Float64
    e::Float64
    omega::Float64
    chi::Float64

    drw_rms::Array{Float64,1}
    drw_rate::Array{Float64,1}

    osc_rms::Array{Float64, 1}
    osc_freq::Array{Float64, 1}
    osc_Q::Array{Float64, 1}
end

function nparams(post::MultiEpochPosterior)
    2*size(post.ts, 1) + 5 + 2*post.ndrw + 3*post.nosc
end

function to_params(post::MultiEpochPosterior, x::Array{Float64, 1})
    i = 1
    mu = zeros(size(post.ts, 1))
    nu = zeros(size(post.ts, 1))
    for j in 1:size(post.ts, 1)
        mu[j] = x[i]
        nu[j] = Parameterizations.bounded_value(x[i+1], 0.5, 2.0)
        i = i + 2
    end

    K = Parameterizations.bounded_value(x[i], post.K_min, post.K_max)
    P = Parameterizations.bounded_value(x[i+1], post.P_min, post.P_max)
    ecosw, esinw = Parameterizations.unit_disk_value(x[i+2:i+3])
    e = sqrt(ecosw*ecosw + esinw*esinw)
    omega = atan(esinw, ecosw)
    chi = Parameterizations.bounded_value(x[i+4], 0.0, 1.0)
    i = i + 5

    drw_rms = zeros(post.ndrw)
    drw_rate = zeros(post.ndrw)

    rmin = post.rate_min
    for j in 1:post.ndrw
        drw_rms[j] = Parameterizations.bounded_value(x[i], post.rms_min, post.rms_max)
        drw_rate[j] = Parameterizations.bounded_value(x[i+1], rmin, post.rate_max)
        i = i + 2
        rmin = drw_rate[j]
    end

    osc_rms = zeros(post.nosc)
    osc_freq = zeros(post.nosc)
    osc_Q = zeros(post.nosc)
    freq_min = post.f_min
    for j in 1:post.nosc
        osc_rms[j] = Parameterizations.bounded_value(x[i], post.rms_min, post.rms_max)
        osc_freq[j] = Parameterizations.bounded_value(x[i+1], freq_min, post.f_max)
        osc_Q[j] = Parameterizations.bounded_value(x[i+2], post.Q_min, post.Q_max)
        i = i + 3
        freq_min = osc_freq[j]
    end

    MultiEpochParams(mu, nu, K, P, e, omega, chi, drw_rms, drw_rate, osc_rms, osc_freq, osc_Q)
end

function to_array(post::MultiEpochPosterior, p::MultiEpochParams)
    x = zeros(nparams(post))

    i = 1
    for j in 1:size(post.ts,1)
        x[i] = p.mu[j]
        x[i+1] = Parameterizations.bounded_param(p.nu[j], 0.5, 2.0)
        i = i + 2
    end

    x[i] = Parameterizations.bounded_param(p.K, post.K_min, post.K_max)
    x[i+1] = Parameterizations.bounded_param(p.P, post.P_min, post.P_max)
    x[i+2:i+3] = Parameterizations.unit_disk_param([p.e*cos(p.omega), p.e*sin(p.omega)])
    x[i+4] = Parameterizations.bounded_param(p.chi, 0.0, 1.0)
    i = i + 5

    rmin = post.rate_min
    for j in 1:post.ndrw
        x[i] = Parameterizations.bounded_param(p.drw_rms[j], post.rms_min, post.rms_max)
        x[i+1] = Parameterizations.bounded_param(p.drw_rate[j], rmin, post.rate_max)
        i = i + 2
        rmin = p.drw_rate[j]
    end

    freq_min = post.f_min
    for j in 1:post.nosc
        x[i] = Parameterizations.bounded_param(p.osc_rms[j], post.rms_min, post.rms_max)
        x[i+1] = Parameterizations.bounded_param(p.osc_freq[j], freq_min, post.f_max)
        x[i+2] = Parameterizations.bounded_param(p.osc_Q[j], post.Q_min, post.Q_max)
        i = i + 3
        freq_min = p.osc_freq[j]
    end

    x
end

function log_prior(post::MultiEpochPosterior, x::Array{Float64, 1})
    log_prior(post, to_params(post, x))
end

function bounded_logjac_value(x, low, high)
    Parameterizations.bounded_logjac(x, Parameterizations.bounded_param(x, low, high), low, high)
end

function log_prior(post::MultiEpochPosterior, p::MultiEpochParams)
    logp = 0.0

    for i in 1:size(post.ts, 1)
        mu = mean(post.ys[i])
        sigma = std(post.ys[i])

        logp += logpdf(Normal(mu, 10*sigma), p.mu[i])

        logp -= log(p.nu[i]) # flat in log(nu)
        logp += bounded_logjac_value(p.nu[i], 0.5, 2)
    end

    logp -= log(p.K)
    logp += bounded_logjac_value(p.K, post.K_min, post.K_max)

    logp -= log(p.P)
    logp += bounded_logjac_value(p.P, post.P_min, post.P_max)

    # Uniform prior on unit disk for (e*cos(omega), e*sin(omega))
    z = [p.e*cos(p.omega), p.e*sin(p.omega)]
    logp += Parameterizations.unit_disk_logjac(z, Parameterizations.unit_disk_param(z))

    # Uniform prior in chi
    logp += bounded_logjac_value(p.chi, 0.0, 1.0)

    rmin = post.rate_min
    for i in 1:post.ndrw
        logp -= log(p.drw_rms[i])
        logp += bounded_logjac_value(p.drw_rms[i], post.rms_min, post.rms_max)

        logp -= log(p.drw_rate[i])
        logp += bounded_logjac_value(p.drw_rate[i], rmin, post.rate_max)

        rmin = p.drw_rate[i]
    end

    fmin = post.f_min
    for i in 1:post.nosc
        logp -= log(p.osc_rms[i])
        logp += bounded_logjac_value(p.osc_rms[i], post.rms_min, post.rms_max)

        logp -= log(p.osc_freq[i])
        logp += bounded_logjac_value(p.osc_freq[i], fmin, post.f_max)
        fmin = p.osc_freq[i]

        logp -= log(p.osc_Q[i])
        logp += bounded_logjac_value(p.osc_Q[i], post.Q_min, post.Q_max)
    end

    logp
end

function draw_prior(post::MultiEpochPosterior, n)
    hcat([draw_prior(post) for i in 1:n]...)
end

function rand_flatlog(low, high)
    exp(log(low) + rand()*log(high/low))
end

function draw_prior(post::MultiEpochPosterior)
    nts = size(post.ts, 1)

    mus = zeros(nts)
    nus = zeros(nts)

    for i in eachindex(mus)
        mu = mean(post.ys[i])
        sigma = std(post.ys[i])

        mus[i] = mu + 10*sigma*randn()

        nus[i] = rand_flatlog(0.5, 2.0)
    end

    K = rand_flatlog(post.K_min, post.K_max)
    P = rand_flatlog(post.P_min, post.P_max)

    x = 1.0
    y = 1.0
    while x*x + y*y > 1.0
        x = rand()
        y = rand()
    end

    e = sqrt(x*x + y*y)
    omega = atan(y,x)

    chi = rand()

    drw_rmss = zeros(post.ndrw)
    drw_rates = zeros(post.ndrw)

    for i in eachindex(drw_rmss)
        drw_rmss[i] = rand_flatlog(post.rms_min, post.rms_max)
        drw_rates[i] = rand_flatlog(post.rate_min, post.rate_max)
    end
    inds = sortperm(drw_rates)
    drw_rmss = drw_rmss[inds]
    drw_rates = drw_rates[inds]

    osc_rmss = zeros(post.nosc)
    osc_freqs = zeros(post.nosc)
    osc_Qs = zeros(post.nosc)

    for i in eachindex(osc_rmss)
        osc_rmss[i] = rand_flatlog(post.rms_min, post.rms_max)
        osc_freqs[i] = rand_flatlog(post.f_min, post.f_max)
        osc_Qs[i] = rand_flatlog(post.Q_min, post.Q_max)
    end
    inds = sortperm(osc_freqs)
    osc_rmss = osc_rmss[inds]
    osc_freqs = osc_freqs[inds]
    osc_Qs = osc_Qs[inds]

    p = MultiEpochParams(mus, nus, K, P, e, omega, chi, drw_rmss, drw_rates, osc_rmss, osc_freqs, osc_Qs)

    to_array(post, p)
end

function log_likelihood(post::MultiEpochPosterior, x::Array{Float64, 1})
    log_likelihood(post, to_params(post, x))
end

function produce_ys_dys(post::MultiEpochPosterior, p::MultiEpochParams)
    n = size(post.allts, 1)

    ys = zeros(n)
    dys = zeros(n)

    for i in eachindex(post.ts)
        ys[post.inds[i]] = post.ys[i] .- p.mu[i]
        dys[post.inds[i]] = post.dys[i].*p.nu[i]
    end

    for i in eachindex(ys)
        ys[i] = ys[i] - Kepler.rv(post.allts[i], p.K, p.P, p.e, p.omega, p.chi)
    end

    ys, dys
end

function log_likelihood(post::MultiEpochPosterior, p::MultiEpochParams)
    ys, dys = produce_ys_dys(post, p)

    filt = Celerite.CeleriteKalmanFilter(0.0, p.drw_rms, p.drw_rate, p.osc_rms, p.osc_freq, p.osc_Q)

    Celerite.log_likelihood(filt, post.allts, ys, dys)
end

function residuals(post::MultiEpochPosterior, x::Array{Float64, 1})
    residuals(post, to_params(post, x))
end

function residuals(post::MultiEpochPosterior, p::MultiEpochParams)
    ys, dys = produce_ys_dys(post, p)

    filt = Celerite.CeleriteKalmanFilter(0.0, p.drw_rms, p.drw_rate, p.osc_rms, p.osc_freq, p.osc_Q)

    Celerite.residuals(filt, post.allts, ys, dys)
end

function psd(post::MultiEpochPosterior, x::Array{Float64, 1}, fs::Array{Float64, 1})
    psd(post, to_params(post, x), fs)
end

function psd(post::MultiEpochPosterior, p::MultiEpochParams, fs::Array{Float64, 1})
    filt = Celerite.CeleriteKalmanFilter(0.0, p.drw_rms, p.drw_rate, p.osc_rms, p.osc_freq, p.osc_Q)

    Celerite.psd(filt, fs)
end

function predict(post::MultiEpochPosterior, x::Array{Float64, 1}, ts::Array{Float64, 1})
    predict(post, to_params(post, x), ts)
end

function predict(post::MultiEpochPosterior, p::MultiEpochParams, ts::Array{Float64, 1})
    rvs = [Kepler.rv(t, p.K, p.P, p.e, p.omega, p.chi) for t in ts]

    ys, dys = produce_ys_dys(post, p)

    filt = Celerite.CeleriteKalmanFilter(0.0, p.drw_rms, p.drw_rate, p.osc_rms, p.osc_freq, p.osc_Q)

    ysp, vysp = Celerite.predict(filt, post.allts, ys, dys, ts)

    ysp = ysp .+ rvs

    ysp, vysp
end

function posterior_predictive(post::MultiEpochPosterior, x::Array{Float64, 1})
    posterior_predictive(post, to_params(post, x))
end

function posterior_predictive(post::MultiEpochPosterior, p::MultiEpochParams)
    filt = Celerite.CeleriteKalmanFilter(0.0, p.drw_rms, p.drw_rate, p.osc_rms, p.osc_freq, p.osc_Q)

    allnoise = Celerite.generate(filt, post.allts, zeros(size(post.allts, 1)))

    allrvnoise = zeros(size(allnoise, 1))
    for i in eachindex(allnoise)
        allrvnoise[i] = allnoise[i] + Kepler.rv(post.allts[i], p.K, p.P, p.e, p.omega, p.chi)
    end

    obs_ys = Array{Float64, 1}[]
    for i in eachindex(post.ts)
        ys = allrvnoise[post.inds[i]]
        push!(obs_ys, ys + p.nu[i]*randn(size(ys, 1)) + p.mu[i])
    end

    obs_ys
end

function write(f::Union{HDF5File, HDF5Group}, post::MultiEpochPosterior, ns::EnsembleNest.NestState)
    ns_group = g_create(f, "nest_state")
    write(ns_group, ns)

    xs, lnlikes = EnsembleNest.postsample(ns)

    ps = [to_params(post, xs[:,j]) for j in 1:size(xs,2)]

    f["lnlike", "compress", 3, "shuffle", ()] = lnlikes
    f["mu", "compress", 3, "shuffle", ()] = hcat([p.mu for p in ps]...)
    f["nu", "compress", 3, "shuffle", ()] = hcat([p.nu for p in ps]...)
    f["K", "compress", 3, "shuffle", ()] = [p.K for p in ps]
    f["P", "compress", 3, "shuffle", ()] = [p.P for p in ps]
    f["e", "compress", 3, "shuffle", ()] = [p.e for p in ps]
    f["omega", "compress", 3, "shuffle", ()] = [p.omega for p in ps]
    f["chi", "compress", 3, "shuffle", ()] = [p.chi for p in ps]

    f["drw_rms", "compress", 3, "shuffle", ()] = hcat([p.drw_rms for p in ps]...)
    f["drw_rate", "compress", 3, "shuffle", ()] = hcat([p.drw_rate for p in ps]...)

    f["osc_rms", "compress", 3, "shuffle", ()] = hcat([p.osc_rms for p in ps]...)
    f["osc_freq", "compress", 3, "shuffle", ()] = hcat([p.osc_freq for p in ps]...)
    f["osc_Q", "compress", 3, "shuffle", ()] = hcat([p.osc_Q for p in ps]...)
end

function read(f::Union{HDF5File, HDF5Group}, post::MultiEpochPosterior)
    mus = read(f, "mu")
    nus = read(f, "nu")
    Ks = read(f, "K")
    Ps = read(f, "P")
    es = read(f, "e")
    omegas = read(f, "omega")
    chis = read(f, "chi")
    drw_rmss = read(f, "drw_rms")
    drw_rates = read(f, "drw_rate")
    osc_rmss = read(f, "osc_rms")
    osc_freqs = read(f, "osc_freq")
    osc_Qs = read(f, "osc_Q")

    ps = [MultiEpochParams(mus[:,i], nus[:,i], Ks[i], Ps[i], es[i], omegas[i], chis[i], drw_rmss[:,i], drw_rates[:,i], osc_rmss[:,i], osc_freqs[:,i], osc_Qs[:,i]) for i in eachindex(Ks)]

    logl = x -> log_likelihood(post, x)
    logp = x -> log_prior(post, x)

    ns = EnsembleNest.NestState(f["nest_state"], logl=logl, logp=logp)

    (ps, ns)
end

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GNUMake_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GNUMake")
JLLWrappers.@generate_main_file("GNUMake", UUID("6c1a3432-6b93-5c89-98b5-f3caf0c092ba"))
end  # module GNUMake_jll


text:
module ErrorFreeTransforms

typealias IEEE754 Union(Float32, Float64)

include("add.jl")
include("mul.jl")
include("fma.jl")

export
  # add.jl
  err_add,
  err_fast_add,

  # mul.jl
  fsplit,
  err_mul,
  err_fast_mul,

  # fma.jl
  err_fma,
  err_fma_compr,
  err_fma_nearest,
  err_fma_approx

end


text:
# Benchmarking cudnn activation functions; ours are faster so not included in Knet.jl

# typedef enum
# {
#     CUDNN_ACTIVATION_SIGMOID      = 0,
#     CUDNN_ACTIVATION_RELU         = 1,
#     CUDNN_ACTIVATION_TANH         = 2,
#     CUDNN_ACTIVATION_CLIPPED_RELU = 3,
#     CUDNN_ACTIVATION_ELU          = 4,
#     CUDNN_ACTIVATION_IDENTITY     = 5
# } cudnnActivationMode_t;

cudnnsigm(x)=cudnnActivationForward(x,mode=0)
cudnnrelu(x)=cudnnActivationForward(x,mode=1)
cudnntanh(x)=cudnnActivationForward(x,mode=2)
cudnncrelu(x;coef=1)=cudnnActivationForward(x,mode=3,coef=coef)
cudnnelu(x;coef=1)=cudnnActivationForward(x,mode=4,coef=coef)
cudnnidentity(x)=cudnnActivationForward(x,mode=5)

# typedef enum{
#     CUDNN_NOT_PROPAGATE_NAN  = 0,
#     CUDNN_PROPAGATE_NAN      = 1,
# } cudnnNanPropagation_t;

mutable struct AD; ptr; end
Base.unsafe_convert(::Type{Cptr}, ad::AD)=ad.ptr
function AD(; mode=0, coef=0, reluNanOpt=false)
    d = Cptr[0]
    @cudnn(cudnnCreateActivationDescriptor,(Ptr{Cptr},),d)
    @cudnn(cudnnSetActivationDescriptor,(Cptr,Cint,Cint,Cdouble),d[1],mode,(reluNanOpt ? 1 : 0), coef)
    ad = AD(d[1])
    finalizer(x->@cudnn(cudnnDestroyActivationDescriptor,(Cptr,),x.ptr), ad)
    return ad
end

t4d(x::KnetArray{T,1}) where T = TD(T,(1,1,length(x),1))
t4d(x::KnetArray{T,2}) where T = TD(T,(1,1,size(x,1),size(x,2)))
t4d(x::KnetArray{T,3}) where T = TD(T,(1,size(x)...))
t4d(x::KnetArray{T,4}) where T = TD(T,size(x))
t4d(x::KnetArray{T,5}) where T = TD(T,size(x))
t4d(x::KnetArray{T,N}) where {T,N} = TD(T,(1,1,length(x),1)) # TODO: do something smarter here

function cudnnActivationForward(x::KnetArray{T}; handle=cudnnhandle(), alpha=1, o...) where T # TODO: do not rely on KnetArray
    beta=0
    ad = AD(;o...)
    td = t4d(x)
    y = similar(x)
    @cudnn(cudnnActivationForward,(Cptr,Cptr,Ptr{T},Cptr,Ptr{T},Ptr{T},Cptr,Ptr{T}), handle, ad, Ref(T(alpha)), td, x, Ref(T(beta)), td, y) # TODO: understand Ref
    return y
end

function cudnnActivationBackward(dy::KnetArray{T},y::KnetArray{T},x::KnetArray{T}; handle=cudnnhandle(), alpha=1, o...) where T
    beta = 0
    ad = AD(;o...)
    td = t4d(x)
    dx = similar(x)
    @cudnn(cudnnActivationBackward,
          (Cptr,Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Ptr{T},Cptr,Ptr{T}),
          handle,ad, Ref(T(alpha)),td,y,td,dy,td,x,Ref(T(beta)),td,dx)
    return dx
end

@primitive cudnnActivationForward(x;o...),dy,y  cudnnActivationBackward(value(dy),value(y),value(x);o...)


text:
using Mimi
using MimiFUND
using Plots
using Statistics
using Distributions

# get the model
m = MimiFUND.get_model()

update_param!(m, :climatedynamics, :climatesensitivity, 4.5)

discounting = [i for i in 0:0.005:0.05]

results_sc = []

for i in eachindex(discounting)
    push!(results_sc, MimiFUND.compute_scc(m, year = 2020, eta = 1., prtp = discounting[i]) * 1.68)
end

plot(discounting, results_sc, label = nothing, xlabel = "PRTP (in %)", ylabel = "Social Cost of Carbon (2020 USD)")

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
@testset "Active Set Methods" begin

    # Test active
    N = 10
    p = 3
    model = UnicycleGame(p=p)
    probsize = ProblemSize(N,model)
    game_con = GameConstraintValues(probsize)

    T = Float64
    radius = 1.0
    add_collision_avoidance!(game_con, radius)
    s0 = stampify(:v, :col, 1, 2, 12)
    @test active(game_con, s0)[1] == 0

    game_con.state_conval[1][1].λ[1][1] =  10.0
    game_con.state_conval[1][1].λ[2][1] = -20.0
    game_con.state_conval[1][1].λ[3][1] = -20.0

    game_con.state_conval[1][1].vals[3][1] =  1.0
    game_con.state_conval[1][1].vals[4][1] =  2.0
    game_con.state_conval[1][1].vals[5][1] = -2.0

    Altro.update_active_set!(game_con.state_conval[1][1], Altro.Val(0.0))
    s1 = stampify(:v, :col, 1, 2, 2)
    @test active(game_con, s1)[1] == 1
    s1 = stampify(:v, :col, 1, 2, 3)
    @test active(game_con, s1)[1] == 1
    s1 = stampify(:v, :col, 1, 2, 4)
    @test active(game_con, s1)[1] == 1
    s1 = stampify(:v, :col, 1, 2, 5)
    @test active(game_con, s1)[1] == 1
    s1 = stampify(:v, :col, 1, 2, 6)
    @test active(game_con, s1)[1] == 0

    # Test active_vertical_mask
    N = 10
    dt = 0.1
    p = 3
    model = UnicycleGame(p=p)
    probsize = ProblemSize(N,model)
    ascore = ActiveSetCore(probsize)
    game_con = GameConstraintValues(probsize)
    radius = 1e-8
    add_collision_avoidance!(game_con, radius)
    n = probsize.n
    m = probsize.m
    S = n*p*(N-1) + m*(N-1) + n*(N-1)
    pdtraj = PrimalDualTraj(probsize, dt)
    # All activated
    active_vertical_mask!(ascore, game_con)
    ascore.vmask == Vector{Int}(1:S+(N-1)*p*(p-1)/2)

    # All deactivated
    Random.seed!(100)
    init_traj!(pdtraj, f=rand, amplitude=1e3)
    update_active_set!(game_con, pdtraj.pr)
    active_vertical_mask!(ascore, game_con)
    ascore.vmask == Vector{Int}(1:S)

    # Test active_horizontal_mask
    N = 10
    dt = 0.1
    p = 3
    model = UnicycleGame(p=p)
    probsize = ProblemSize(N,model)
    ascore = ActiveSetCore(probsize)
    game_con = GameConstraintValues(probsize)
    radius = 1e-8
    add_collision_avoidance!(game_con, radius)
    n = probsize.n
    m = probsize.m
    S = n*p*(N-1) + m*(N-1) + n*(N-1)
    pdtraj = PrimalDualTraj(probsize, dt)
    # All activated
    active_horizontal_mask!(ascore, game_con)
    ascore.hmask == Vector{Int}(1:S+(N-1)*p*(p-1))

    # All deactivated
    Random.seed!(100)
    init_traj!(pdtraj, f=rand, amplitude=1e3)
    update_active_set!(game_con, pdtraj.pr)
    active_horizontal_mask!(ascore, game_con)
    ascore.hmask == Vector{Int}(1:S)

    # Test Nullspace
    T = Float64
    N = 10
    dt = 0.1
    p = 3
    model = UnicycleGame(p=p)
    probsize = ProblemSize(N,model)
    x0 = rand(SVector{model.n,T})
    opts = Options()

    Q = [Diagonal(rand(SVector{model.ni[i],T})) for i=1:p]
    R = [Diagonal(rand(SVector{model.mi[i],T})) for i=1:p]
    xf = [i*ones(SVector{model.ni[i],T}) for i=1:p]
    uf = [2i*ones(SVector{model.mi[i],T}) for i=1:p]
    game_obj = GameObjective(Q,R,xf,uf,N,model)

    game_con = GameConstraintValues(probsize)
    radius = 1.0
    add_collision_avoidance!(game_con, radius)
    prob = GameProblem(N, dt, x0, model, opts, game_obj, game_con)

    ascore = ActiveSetCore(probsize)
    residual!(ascore, prob, prob.pdtraj)
    residual_jacobian!(ascore, prob, prob.pdtraj)

    update_nullspace!(ascore, prob, prob.pdtraj)
    @test size(ascore.null.mat) == (probsize.S+(N-1)*p*(p-1), (N-1)*p)
    @test length(ascore.null.vec) == (N-1)*p
    @test length(ascore.null.vec[1]) == probsize.S+(N-1)*p*(p-1)

end


text:
module Land

using DocStringExtensions
using LinearAlgebra, StaticArrays

using CLIMAParameters
using CLIMAParameters.Planet:
    ρ_cloud_liq, ρ_cloud_ice, cp_l, cp_i, T_0, LH_f0, T_freeze, grav

using ..VariableTemplates
using ..MPIStateArrays
using ..BalanceLaws
import ..BalanceLaws:
    BalanceLaw,
    vars_state,
    flux_first_order!,
    flux_second_order!,
    source!,
    boundary_conditions,
    boundary_state!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    nodal_init_state_auxiliary!,
    init_state_prognostic!,
    nodal_update_auxiliary_state!
using ..DGMethods: LocalGeometry, DGModel
export LandModel

"""
    LandModel{PS, S, SRC, IS} <: BalanceLaw

A BalanceLaw for land modeling.
Users may over-ride prescribed default values for each field.

# Usage

    LandModel(
        param_set,
        soil,
        source
        init_state_prognostic
    )

# Fields
$(DocStringExtensions.FIELDS)
"""
struct LandModel{PS, S, SRC, IS} <: BalanceLaw
    "Parameter set"
    param_set::PS
    "Soil model"
    soil::S
    "Source Terms (Problem specific source terms)"
    source::SRC
    "Initial Condition (Function to assign initial values of state variables)"
    init_state_prognostic::IS
end

"""
    LandModel(
        param_set::AbstractParameterSet,
        soil::BalanceLaw;
        source::SRC = (),
        init_state_prognostic::IS = nothing
    ) where {SRC, IS}

Constructor for the LandModel structure.
"""
function LandModel(
    param_set::AbstractParameterSet,
    soil::BalanceLaw;
    source::SRC = (),
    init_state_prognostic::IS = nothing,
) where {SRC, IS}
    @assert init_state_prognostic ≠ nothing
    land = (param_set, soil, source, init_state_prognostic)
    return LandModel{typeof.(land)...}(land...)
end

function vars_state(land::LandModel, st::Prognostic, FT)
    @vars begin
        soil::vars_state(land.soil, st, FT)
    end
end

function vars_state(land::LandModel, st::Auxiliary, FT)
    @vars begin
        z::FT
        soil::vars_state(land.soil, st, FT)
    end
end

function vars_state(land::LandModel, st::Gradient, FT)
    @vars begin
        soil::vars_state(land.soil, st, FT)
    end
end

function vars_state(land::LandModel, st::GradientFlux, FT)
    @vars begin
        soil::vars_state(land.soil, st, FT)
    end
end

function nodal_init_state_auxiliary!(
    land::LandModel,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
)
    aux.z = geom.coord[3]
    land_init_aux!(land, land.soil, aux, geom)
end

function flux_first_order!(
    land::LandModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    directions,
) end

function compute_gradient_argument!(
    land::LandModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)

    compute_gradient_argument!(land, land.soil, transform, state, aux, t)
end

function compute_gradient_flux!(
    land::LandModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    compute_gradient_flux!(
        land,
        land.soil,
        diffusive,
        ∇transform,
        state,
        aux,
        t,
    )

end

function flux_second_order!(
    land::LandModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    flux_second_order!(
        land,
        land.soil,
        flux,
        state,
        diffusive,
        hyperdiffusive,
        aux,
        t,
    )

end

function nodal_update_auxiliary_state!(
    land::LandModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    land_nodal_update_auxiliary_state!(land, land.soil, state, aux, t)
end

function source!(
    land::LandModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    land_source!(land.source, land, source, state, diffusive, aux, t, direction)
end

function init_state_prognostic!(
    land::LandModel,
    state::Vars,
    aux::Vars,
    coords,
    t,
    args...,
)
    land.init_state_prognostic(land, state, aux, coords, t, args...)
end

include("land_bc.jl")
include("SoilWaterParameterizations.jl")
using .SoilWaterParameterizations
include("SoilHeatParameterizations.jl")
using .SoilHeatParameterizations
include("Runoff.jl")
using .Runoff
include("soil_model.jl")
include("soil_boundary_types.jl")
include("soil_water.jl")
include("soil_heat.jl")
include("soil_bc.jl")
include("source.jl")
end # Module

max_forks_repo_name: jm-c/CLIMA


text: println("Hello,World!")

max_stars_repo_name: amanoese/splash-julia

max_issues_repo_name: amanoese/splash-julia

max_forks_repo_name: amanoese/splash-julia


text:
using Spotify
using Spotify.Library
using Spotify.Tracks

function listallmine()
    offset = 0
    st, repeat_after_sec = library_get_saved_tracks(2, offset, "NO");
    st
    @info "Listing $(st.total) tracks. Pauses may occur due the the current API rate limit. It may be increased by calling the API with a wider scope."
    offset = 0
    while true
        while true
            repeat_after_sec !=0 && sleep(repeat_after_sec + 0.5)
            st, repeat_after_sec = library_get_saved_tracks(50, offset, "NO")
            repeat_after_sec > 5 && @info "Waiting for $repeat_after_sec s for library API. Scope may affect limits."
            repeat_after_sec == 0 && break
        end
        if length(st) != 0
            for (i, t) in enumerate(get(st, :items, []))
                tr = t.track
                s1 = lpad(i + offset, 5)
                s2 = lpad(tr.name, 100)
                id = tr.id
                te = 0
                af = Spotify.JSON3.Object()
                if Spotify.isid(id)
                    while true
                        repeat_after_sec !=0 && sleep(repeat_after_sec + 0.5)
                        af, repeat_after_sec =  tracks_get_audio_features(id)
                        if !(typeof(af)  <: Spotify.JSON3.Object)
                            @warn typeof(af)
                        end
                        repeat_after_sec > 5 && @info "\n Waiting for $repeat_after_sec s for tracks API. Scope may affect limits."
                        repeat_after_sec == 0 && break
                    end
                    te = round(Int64, get(af, :tempo, 0))
                end
                s3 = lpad(te, 4)
                println(s1, s2, s3)
            end
        end
        offset += length(st.items)
        offset >= st.total && break
    end
end

if !Spotify.has_ig_access_token()
    Spotify.get_implicit_grant()
    Spotify.wait_for_ig_access()
end
listallmine()


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule QCDNUM_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("QCDNUM")
JLLWrappers.@generate_main_file("QCDNUM", UUID("4f99f0f9-6eee-532a-a65d-722344444835"))
end  # module QCDNUM_jll


text:
import Soundcard
using Test
using Statistics
using SharedArrays
using Libaudio

let (ndev, digest) = Soundcard.device()
    @test ndev > 0
    @info "Soundcard device populate" digest
end

let x = 10^(-20/20) * randn(192000,1), ms = zeros(1,2)
    ms[1,2] = 1
    y = Soundcard.play(x, ms, 48000)
    @test y == nothing
end

let mm = zeros(8,1)
    mm[2,1] = 1
    y = Soundcard.record(48000, mm, 48000)
    @test typeof(y) == Matrix{Float32}
end

let x = 10^(-20/20) * randn(192000,1), ms = zeros(1,2), mm = zeros(8,1)
    ms[1,1] = 1
    mm[2,1] = 1
    y = Soundcard.playrecord(x, ms, mm, 48000)
    @info "loop gain" std(y)/std(x)
    @test std(y)/std(x) < 1e-4

    ms[1,2] = 1
    mm[2,1] = 1
    y = Soundcard.playrecord(x, ms, mm, 48000)
    @info "loop gain" std(y)/std(x)
    @test std(y)/std(x) > 0.1
end

let out = 10^(-20/20) * randn(192000,1), ms = zeros(1,2)
    ms[1,2] = 1
    x = Soundcard.mixer(convert(Matrix{Float32},out), ms)
    y = SharedArray{Float32,1}(Soundcard.interleave(x))
    r = Soundcard.play(size(x), y, 48000)
        # async tasks...
    # fetch(done)
    @info r
    @test r == 0
end

let out = 10^(-20/20) * randn(192000,1), mm = zeros(8,1)
    mm[2,1] = 1
    samples = 96000
    y = SharedArray{Float32,1}(zeros(Float32, size(mm,1)*samples))
    Soundcard.record(y, size(mm), samples, 48000)
        # async tasks...
    # fetch(done)
    recording = Soundcard.mixer(transpose(reshape(y,size(mm,1),samples)), mm)
    @test typeof(recording) == Matrix{Float32}
end

let out = 10^(-20/20) * randn(192000,1), ms = zeros(1,2), mm = zeros(8,1)
    ms[1,1] = 1
    mm[2,1] = 1
    samples = size(out,1)
    x = Soundcard.mixer(convert(Matrix{Float32},out), ms)
    y = SharedArray{Float32,1}(Soundcard.interleave(x))
    z = SharedArray{Float32,1}(zeros(Float32,size(mm,1)*samples))
    r = Soundcard.playrecord(size(x), y, z, size(mm), 48000)
        # async tasks...
    # fetch(done)
    recording = Soundcard.mixer(transpose(reshape(z,size(mm,1),samples)), mm)
    @info r
    @test r == 0
    @test typeof(recording) == Matrix{Float32}
end

let s = 10^(-10/20) * sinpi.(2*1000*(0:48000-1)/48000), ms = zeros(1,2), mm = zeros(8,1)
    x = s[:,:]
    ms[1,2] = 1
    mm[2,1] = 1
    y = Soundcard.playrecord(x, ms, mm, 48000)
    Libaudio.wavwrite("loop.wav", y, 48000, 32)
    # m = Libaudio.zerocrossingrate(s)
    # n = Libaudio.zerocrossingrate(y[:,1])
    # @info length(m[m.==1.0])
    # @info length(n[n.==1.0])
end


text:
module GTrends

using RCall: rimport, rcopy

export gtrends

gtrendsR = Ref{Module}()

function __init__()
    gtrendsR[] = rimport("gtrendsR")
end

#-----------------------------------------------------------------------------# Query
const UVS = Union{String, Vector{String}}

"""
    gtrends(q; kw...)

# Keyword Args

- `geo = "US"`: Geographic region, e.g. `"US"`, `["US", "CA"]`.
- `time="all"`: Date range - `now 1-H`, `now 4-H`, `now 1-d`, `today 1-m`, `today 3-m`, `today 12-m`, `today+5-y`, `all`, `Y-m-d Y-m-d`
- `gprop = "web"`: `web`, `news`, `images`, `froogle`, `youtube`.
- `category = 0`: I'm unsure of what categories are here.
- `hl = "en-US"`: ISO language code.
- `low_search_volume=false`: Should include regions with low search volumes .
- `tz = ""`:  See gtrendsR documentation.
- `onlyInterest=false`: If you only want the main interest-over-time data.
"""
function gtrends(q::UVS; geo::UVS = "US", time::String = "all", gprop::UVS = "web", category::Int=0,
                 hl::String = "en-US", low_search_volume::Bool=false, tz="", onlyInterest=false)
    rcopy(gtrendsR[].gtrends(q; geo, time, gprop, category, hl, low_search_volume, onlyInterest))
end

end