text:
# This file is part of Kpax3. License is MIT.

#=
#R = [1; 1; 1; 1; 1; 1]
#R = [1; 1; 1; 1; 2; 2]
#R = [1; 1; 2; 2; 3; 3]
#R = [1; 3; 2; 2; 1; 4]
k = length(unique(R))

priorR = Kpax3.EwensPitman(settings.α, settings.θ)
priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)

C = zeros(UInt8, n, m)

cs = ((UInt8[1], UInt8[2], UInt8[3], UInt8[4]),
      (UInt8[1; 1], UInt8[2; 2], UInt8[3; 3], UInt8[3; 4], UInt8[4; 3],
       UInt8[4; 4]),
      (UInt8[1; 1; 1], UInt8[2; 2; 2], UInt8[3; 3; 3], UInt8[3; 3; 4],
       UInt8[3; 4; 3], UInt8[4; 3; 3], UInt8[3; 4; 4], UInt8[4; 3; 4],
       UInt8[4; 4; 3], UInt8[4; 4; 4]),
      (UInt8[1; 1; 1; 1], UInt8[2; 2; 2; 2], UInt8[3; 3; 3; 3],
       UInt8[3; 3; 3; 4], UInt8[3; 3; 4; 3], UInt8[3; 4; 3; 3],
       UInt8[4; 3; 3; 3], UInt8[3; 3; 4; 4], UInt8[3; 4; 3; 4],
       UInt8[3; 4; 4; 3], UInt8[4; 3; 3; 4], UInt8[4; 3; 4; 3],
       UInt8[4; 4; 3; 3], UInt8[3; 4; 4; 4], UInt8[4; 3; 4; 4],
       UInt8[4; 4; 3; 4], UInt8[4; 4; 4; 3], UInt8[4; 4; 4; 4]),
      (UInt8[1; 1; 1; 1; 1], UInt8[2; 2; 2; 2; 2], UInt8[3; 3; 3; 3; 3],
       UInt8[3; 3; 3; 3; 4], UInt8[3; 3; 3; 4; 3], UInt8[3; 3; 4; 3; 3],
       UInt8[3; 4; 3; 3; 3], UInt8[4; 3; 3; 3; 3], UInt8[3; 3; 3; 4; 4],
       UInt8[3; 3; 4; 3; 4], UInt8[3; 3; 4; 4; 3], UInt8[3; 4; 3; 3; 4],
       UInt8[3; 4; 3; 4; 3], UInt8[3; 4; 4; 3; 3], UInt8[4; 3; 3; 3; 4],
       UInt8[4; 3; 3; 4; 3], UInt8[4; 3; 4; 3; 3], UInt8[4; 4; 3; 3; 3],
       UInt8[3; 3; 4; 4; 4], UInt8[3; 4; 3; 4; 4], UInt8[3; 4; 4; 3; 4],
       UInt8[3; 4; 4; 4; 3], UInt8[4; 3; 3; 4; 4], UInt8[4; 3; 4; 3; 4],
       UInt8[4; 4; 3; 3; 4], UInt8[4; 3; 4; 4; 3], UInt8[4; 4; 3; 4; 3],
       UInt8[4; 4; 4; 3; 3], UInt8[3; 4; 4; 4; 4], UInt8[4; 3; 4; 4; 4],
       UInt8[4; 4; 3; 4; 4], UInt8[4; 4; 4; 3; 4], UInt8[4; 4; 4; 4; 3],
       UInt8[4; 4; 4; 4; 4]),
      (UInt8[1; 1; 1; 1; 1; 1], UInt8[2; 2; 2; 2; 2; 2],
       UInt8[3; 3; 3; 3; 3; 3], UInt8[3; 3; 3; 3; 3; 4],
       UInt8[3; 3; 3; 3; 4; 3], UInt8[3; 3; 3; 4; 3; 3],
       UInt8[3; 3; 4; 3; 3; 3], UInt8[3; 4; 3; 3; 3; 3],
       UInt8[4; 3; 3; 3; 3; 3], UInt8[3; 3; 3; 3; 4; 4],
       UInt8[3; 3; 3; 4; 3; 4], UInt8[3; 3; 3; 4; 4; 3],
       UInt8[3; 3; 4; 3; 3; 4], UInt8[3; 3; 4; 3; 4; 3],
       UInt8[3; 3; 4; 4; 3; 3], UInt8[3; 4; 3; 3; 3; 4],
       UInt8[3; 4; 3; 3; 4; 3], UInt8[3; 4; 3; 4; 3; 3],
       UInt8[3; 4; 4; 3; 3; 3], UInt8[4; 3; 3; 3; 3; 4],
       UInt8[4; 3; 3; 3; 4; 3], UInt8[4; 3; 3; 4; 3; 3],
       UInt8[4; 3; 4; 3; 3; 3], UInt8[4; 4; 3; 3; 3; 3],
       UInt8[3; 3; 3; 4; 4; 4], UInt8[3; 3; 4; 3; 4; 4],
       UInt8[3; 3; 4; 4; 3; 4], UInt8[3; 3; 4; 4; 4; 3],
       UInt8[3; 4; 3; 3; 4; 4], UInt8[3; 4; 3; 4; 3; 4],
       UInt8[3; 4; 4; 3; 3; 4], UInt8[3; 4; 3; 4; 4; 3],
       UInt8[3; 4; 4; 3; 4; 3], UInt8[3; 4; 4; 4; 3; 3],
       UInt8[4; 3; 3; 3; 4; 4], UInt8[4; 3; 3; 4; 3; 4],
       UInt8[4; 3; 4; 3; 3; 4], UInt8[4; 4; 3; 3; 3; 4],
       UInt8[4; 3; 3; 4; 4; 3], UInt8[4; 3; 4; 3; 4; 3],
       UInt8[4; 4; 3; 3; 4; 3], UInt8[4; 3; 4; 4; 3; 3],
       UInt8[4; 4; 3; 4; 3; 3], UInt8[4; 4; 4; 3; 3; 3],
       UInt8[3; 3; 4; 4; 4; 4], UInt8[3; 4; 3; 4; 4; 4],
       UInt8[3; 4; 4; 3; 4; 4], UInt8[3; 4; 4; 4; 3; 4],
       UInt8[3; 4; 4; 4; 4; 3], UInt8[4; 3; 3; 4; 4; 4],
       UInt8[4; 3; 4; 3; 4; 4], UInt8[4; 3; 4; 4; 3; 4],
       UInt8[4; 3; 4; 4; 4; 3], UInt8[4; 4; 3; 3; 4; 4],
       UInt8[4; 4; 3; 4; 3; 4], UInt8[4; 4; 3; 4; 4; 3],
       UInt8[4; 4; 4; 3; 3; 4], UInt8[4; 4; 4; 3; 4; 3],
       UInt8[4; 4; 4; 4; 3; 3], UInt8[3; 4; 4; 4; 4; 4],
       UInt8[4; 3; 4; 4; 4; 4], UInt8[4; 4; 3; 4; 4; 4],
       UInt8[4; 4; 4; 3; 4; 4], UInt8[4; 4; 4; 4; 3; 4],
       UInt8[4; 4; 4; 4; 4; 3], UInt8[4; 4; 4; 4; 4; 4]))

logpnew = -Inf
logpold = -Inf
l = 0
for c1 in cs[state.k], c2 in cs[state.k], c3 in cs[state.k], c4 in cs[state.k]
  l += 1

  tmp = hcat(c1, c2, c3, c4)

  for l in 1:state.k
    state.C[state.cl[l], :] = tmp[l, :]
  end

  logpnew = Kpax3.logcondpostC(state.C, state.cl, state.k, state.v, state.n1s, priorC)

  if logpnew > logpold
    logpold = logpnew
    C[1:state.k, :] = copy(state.C[state.cl[1:state.k], :])
  end
end

C
=#

function test_local_mode()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, gamma=[0.3; 0.35; 0.35])

  data = UInt8[1 1 1 1 0 0;
               0 0 1 1 0 0;
               1 0 0 0 1 1;
               1 0 1 0 1 0]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # k = 1
  R = [1; 1; 1; 1; 1; 1]
  k = 1
  C = zeros(UInt8, 6, 4)
  C[1, :] .= 0x01
  cl = [1; 0; 0; 0; 0; 0]
  v = [6; 0; 0; 0; 0; 0]
  n1s = zeros(Float64, 6, 4)
  n1s[1, :] .= vec(sum(float(data), dims=2))

  s = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(s.v, s.n1s, s.C, s.cl, s.k, s.logpC, priorC)

  @test s.C == C
  @test isapprox(s.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(s.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  t = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(t, priorC)

  @test t.C == C
  @test isapprox(t.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(t.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  # k = 2
  R = [1; 1; 1; 1; 2; 2]
  k = 2
  C = zeros(UInt8, 6, 4)
  C[1:2, :] = UInt8[4 1 2 1; 3 1 2 1]
  cl = [1; 2; 0; 0; 0; 0]
  v = [4; 2; 0; 0; 0; 0]
  n1s = zeros(Float64, 6, 4)
  n1s[1:2, :] .= permutedims(hcat(sum(float(data[:, R .== 1]), dims=2), sum(float(data[:, R .== 2]), dims=2)))

  s = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(s.v, s.n1s, s.C, s.cl, s.k, s.logpC, priorC)

  @test s.C == C
  @test isapprox(s.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(s.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  t = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(t, priorC)

  @test t.C == C
  @test isapprox(t.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(t.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  # k = 3
  R = [1; 1; 2; 2; 3; 3]
  k = 3
  C = zeros(UInt8, 6, 4)
  C[1:3, :] = UInt8[4 3 2 1; 4 4 2 1; 3 3 2 1]
  cl = [1; 2; 3; 0; 0; 0]
  v = [2; 2; 2; 0; 0; 0]
  n1s = zeros(Float64, 6, 4)
  n1s[1:3, :] .= permutedims(hcat(sum(float(data[:, R .== 1]), dims=2), sum(float(data[:, R .== 2]), dims=2), sum(float(data[:, R .== 3]), dims=2)))

  s = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(s.v, s.n1s, s.C, s.cl, s.k, s.logpC, priorC)

  @test s.C == C
  @test isapprox(s.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(s.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  t = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(t, priorC)

  @test t.C == C
  @test isapprox(t.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(t.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  # k = 4
  R = [1; 3; 2; 2; 1; 4]
  k = 4
  C = zeros(UInt8, 6, 4)
  C[1:4, :] = UInt8[1 3 4 2; 1 4 3 2; 1 3 3 2; 1 3 4 2]
  cl = [1; 2; 3; 4; 0; 0]
  v = [2; 2; 1; 1; 0; 0]
  n1s = zeros(Float64, 6, 4)
  n1s[1:4, :] .= permutedims(hcat(sum(float(data[:, R .== 1]), dims=2), sum(float(data[:, R .== 2]), dims=2), sum(float(data[:, R .== 3]), dims=2), sum(float(data[:, R .== 4]), dims=2)))

  s = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(s.v, s.n1s, s.C, s.cl, s.k, s.logpC, priorC)

  @test s.C == C
  @test isapprox(s.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(s.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  t = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(t, priorC)

  @test t.C == C
  @test isapprox(t.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(t.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  # k = 5
  R = [1; 2; 3; 3; 4; 5]
  k = 5
  C = zeros(UInt8, 6, 4)
  C[1:5, :] = UInt8[2 3 2 1; 2 3 2 1; 2 4 2 1; 2 3 2 1; 2 3 2 1]
  cl = [1; 2; 3; 4; 5; 0]
  v = [1; 1; 2; 1; 1; 0]
  n1s = zeros(Float64, 6, 4)
  n1s[1:5, :] .= permutedims(hcat(sum(float(data[:, R .== 1]), dims=2), sum(float(data[:, R .== 2]), dims=2), sum(float(data[:, R .== 3]), dims=2), sum(float(data[:, R .== 4]), dims=2), sum(float(data[:, R .== 5]), dims=2)))

  s = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(s.v, s.n1s, s.C, s.cl, s.k, s.logpC, priorC)

  @test s.C == C
  @test isapprox(s.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(s.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  t = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(t, priorC)

  @test t.C == C
  @test isapprox(t.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(t.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  # k = 6
  R = [1; 2; 3; 4; 5; 6]
  k = 6
  C = UInt8[1 1 2 2; 1 1 2 2 ; 1 1 2 2; 1 1 2 2; 1 1 2 2; 1 1 2 2]
  cl = [1; 2; 3; 4; 5; 6]
  v = [1; 1; 1; 1; 1; 1]
  n1s = permutedims(float(data))

  s = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(s.v, s.n1s, s.C, s.cl, s.k, s.logpC, priorC)

  @test s.C == C
  @test isapprox(s.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(s.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  t = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  Kpax3.computelocalmode!(t, priorC)

  @test t.C == C
  @test isapprox(t.logpC[1], Kpax3.logpriorC(C, cl, k, priorC), atol=ε)
  @test isapprox(t.logpC[2], Kpax3.logcondpostC(C, cl, k, v, n1s, priorC), atol=ε)

  nothing
end

test_local_mode()


text:
using Test
using MianalyzerBackend
const KMeans = MianalyzerBackend.KMeans

@testset "KMeans" begin
  @testset ".call" begin
    actual = KMeans.call([[1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4], [5, 5, 5]], 1)
    @test actual[:clusters] == [[[1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4], [5, 5, 5]]]
    @test actual[:centers] == [[3.0, 3.0, 3.0]]

    actual = KMeans.call([[1, 1], [2, 2]], 2)
    @test actual[:clusters] == [[[1, 1]], [[2, 2]]]
    @test actual[:centers] == [[1.0, 1.0], [2.0, 2.0]]
  end
end

ext: jl


text:
function synthetic_sample(no_obs::Integer, obsl, bhmm_vec::AbstractVector{<:BHMM}, bhmm_dist::Categorical, spikes::AbstractVector{<:AbstractMatrix{<:AbstractFloat}}, spike_instruct::AbstractVector{<:Tuple{<:Bool, <:Tuple}})
    !(typeof(obsl) <:UnitRange || typeof(obsl) <:Integer) && throw(ArgumentError("obsl must be Integer or UnitRange"))
    length(spike_instruct)!=length(spikes) && throw(ArgumentError("spike_instruct must be as long as spikes"))

    obs, hmm_truth=obs_array_from_bhmms(no_obs,obsl,bhmm_vec,bhmm_dist)

    spike_truth = spike_obs!(obs, spikes, spike_instruct)

    bg_scores = score_synthetic(obs, bhmm_vec, hmm_truth)

    return obs, bg_scores, hmm_truth, spike_truth
end

function obs_array_from_bhmms(no_obs, obsl, bhmm_vec, bhmm_dist)
    obs=zeros(UInt8,max(obsl...)+1, no_obs)
    truthvec=zeros(UInt8,no_obs)
    for o in 1:no_obs
        hmm_idx=rand(bhmm_dist)
        truthvec[o]=hmm_idx
        bhmm=bhmm_vec[hmm_idx]
        typeof(obsl)<:UnitRange ? (l=rand(obsl)) : (l=obsl)
        obs[1:l,o]=rand(bhmm, l)
    end
    return obs, truthvec
end

function spike_obs!(obs, spikes, spike_instruct)
    truthmat=falses(size(obs,2),length(spikes))
    for (s,spike) in enumerate(spikes)
        structural,ins=spike_instruct[s]
        structural ? (truthmat[:,s] = spike_struc!(obs, spike, ins...)) : (truthmat[:,s] = spike_irreg!(obs, spike, ins...))
    end
    return truthmat
end

function spike_struc!(obs, spike, frac_obs, periodicity)
    truth=falses(size(obs,2))
    for o in 1:size(obs,2)
        if rand() < frac_obs
            truth[o]=true
            rand()<.5 ? (source=revcomp_pwm(spike)) : (source=spike)
            pos=rand(1:periodicity)
            oidx=findfirst(iszero,obs[:,o])
            while pos<oidx
                pos_ctr=pos
                pwm_ctr=1
                while pos_ctr<oidx&&pwm_ctr<=size(source,1)
                    obs[pos_ctr,o]=rand(Categorical(source[pwm_ctr,:]))
                    pos_ctr+=1
                    pwm_ctr+=1
                end
                pos+=periodicity+pwm_ctr
            end
        end
    end
    return truth
end

function spike_irreg!(obs, spike, frac_obs, recur)
    truth=falses(size(obs,2))
    for o in 1:size(obs,2)
        oidx=findfirst(iszero,obs[:,o])
        if rand()<frac_obs
            truth[o]=true
            for r in 1:rand(recur)
                rand()<.5 ? (source=revcomp_pwm(spike)) : source=spike
                pos=rand(1:oidx-1)
                pwm_ctr=1
                while pos<oidx && pwm_ctr<=size(source,1)
                    obs[pos,o]=rand(Categorical(source[pwm_ctr,:]))
                    pos+=1
                    pwm_ctr+=1
                end
            end
        end
    end
    return truth
end

function score_synthetic(obs, bhmm_vec, hmm_truth)
    lh_mat=zeros(size(obs,1)-1, size(obs,2))
    for o in 1:size(obs,2)
        oidx=findfirst(iszero,obs[:,o])
        lh_mat[1:oidx-1,o]=BioBackgroundModels.get_BGHMM_symbol_lh(transpose(obs[1:oidx,o:o]),bhmm_vec[hmm_truth[o]])
    end
    return lh_mat
end


text:
using NLopt, Test

count = 0 # keep track of # function evaluations

function myfunc(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end

    global count
    count::Int += 1
    println("f_$count($x)")

    sqrt(x[2])
end

function myconstraint(x::Vector, grad::Vector, a, b)
    if length(grad) > 0
        grad[1] = 3a * (a*x[1] + b)^2
        grad[2] = -1
    end
    (a*x[1] + b)^3 - x[2]
end

opt = Opt(:LD_MMA, 2)
opt.lower_bounds = [-Inf, 0.]
opt.xtol_rel = 1e-4
opt.min_objective = myfunc
opt.inequality_constraint = (x,g) -> myconstraint(x,g,2,0)
opt.inequality_constraint = (x,g) -> myconstraint(x,g,-1,1)

(minf,minx,ret) = optimize(opt, [1.234, 5.678])
println("got $minf at $minx after $count iterations (returned $ret)")

@test minx[1] ≈ 1/3 rtol=1e-5
@test minx[2] ≈ 8/27 rtol=1e-5
@test minf ≈ sqrt(8/27) rtol=1e-5
@test ret == :XTOL_REACHED
@test opt.numevals == count


text:
# Process QIN Breast Data using NRRM and CNRRM

# Run this script by launching julia then entering: include("path/to/e03_4_qinProcessNLSQ.jl")

# Estimated runtime: ~500 seconds, on 4 cores

# Choice of parameters for the reference tissue. Should match e03_05
# It's best to leave this as it is
refName = "refY"

# First, figure out the current file path
curFilePath = Base.source_path() # Path to current file
workDir = dirname(curFilePath) # Our working directory (directory of current file)
## If the above doesn't work, then user should enter the location of the CLRRM directory
## by manually modifying this script

# Load pre-requisite packages
using DCEMRI
using MAT

# Get the directories and files of interest
matDir = workDir * "/data/QINBreast/ForJulia"
warmupFile = workDir * "/data/QIBA-ToftsV6/QIBAv6-Mini4Jl-warmup.mat"
auxCode = workDir * "/jlfiles/auxCode.jl"
refRegCode = workDir * "/jlfiles/refRegionFunction.jl"

outDir = workDir * "/data/QINBreast/FromJulia"

if !(isdir(outDir))
  mkdir(outDir)
end

############### Warmup run
results = fitdata(datafile=warmupFile, models=[2])
rm("output.mat") # Delete output since it isn't needed
##############

# This will load the functions for fitting the models
include(auxCode)
include(refRegCode)

# Define the reference tissue parameters
refKTrans = 0.1;
refVe = 0.1;
if (refName == "refWS")
  refKTrans = 0.07;
  refVe = 0.14;
elseif (refName == "refP")
  refKTrans = 0.137;
  refVe = 0.115;
end
refKep = refKTrans/refVe

############## Main run
matFiles = readdir(matDir)
numFiles = length(matFiles)

# Loop through the .mat files
for q=1:numFiles
  curName = matFiles[q]
  curFile = matDir * "/$curName"

  # Load data
  matData = matread(curFile)
  Ct = matData["Ct"]
  (nT, nV) = size(Ct)
  Ct = reshape(Ct, nT, nV, 1)
  Cp = vec(matData["Cp"])
  t = vec(matData["t"])
  tResFactor = round(Int64,matData["tResFactor"])
  (nT, nV) = size(Ct)
  mask = ones(nV,1) .> 0

  # Simulate reference tissue curve
  Crr = DCEMRI.toftskety(t, [refKTrans, refKep], Cp)
  Ct = reshape(Ct,nT,nV)
  # Downsample so that Crr, Cp, and t have same temporal resolution as Ct
  Crr = downsample(Crr,tResFactor)
  Cp = downsample(Cp,tResFactor)
  t = downsample(t,tResFactor)

  # Non-Linear Tofts Model
  pkParamsT, residT = DCEMRI.fitdce(Ct, mask, t, Cp, models=[2])

  # Re-arrange the output from Tofts Model fit
  kt = pkParamsT[1,:,1]
  ve = pkParamsT[2,:,1]
  kep = kt./ve

  # Parameter mask using threshold on fitted Tofts Model parameters
  pMask = (kt.>0.001) & (kt.<1) & (ve.>0.001) & (ve.<1) & (kep.<35)

  kt = kt[pMask]
  ve = ve[pMask]
  kep = kep[pMask]

  # NRRM and CNRRM
  (pkParamsCN, residCN, estKepRR, pkParamsN, residN, runtimeN, runtimeCN)=fitCNRRM(Ct, Crr, t, estKepRR=0.0, doTime=true)

  # Output results
  results = Dict()
  results["pkParamsN"] = pkParamsN
  results["residN"] = residN
  results["pkParamsCN"] = pkParamsCN
  results["residCN"] = residCN
  results["runtimeN"] = runtimeN
  results["runtimeC"] = runtimeCN
  results["estKepRRN"] = estKepRR
  results["pkParamsT"] = pkParamsT
  results["residT"] = residT
  results["kt"] = kt
  results["ve"] = ve
  results["kep"] = kep
  results["Crr"] = Crr
  results["Cp"] = Cp
  results["t"] = t
  results["pMask"] = convert(Array{Bool,1}, pMask)
  matwrite(outDir * "/$curName", results)
end

max_stars_repo_name: MPUmri/CRRM

max_issues_repo_name: MPUmri/CRRM

max_forks_repo_name: MPUmri/CRRM


text:
####### Mixing factor functions #########

export mixFactorToFreq, mixFactorToFreqIdx, mixingFactors

"""
This function returns the multiplicator `k` corresponding to the frequency
given by the mixing factors `mx`, `my`, and `mz`. I.e. the integer value, for
which k*F = mx*fx + my*fy* + mz*fz, where F is the measurement cycle frequency
and fx, fy, and fz are the excitation frequencies for the x,y, and z channel.
"""
function mixFactorToFreq(b::MPIFile,mx,my,mz=0)
  mxyz, mask, freqNumber = calcPrefactors(b)
  k = (mx*mxyz[1]+my*mxyz[2]+mz*mxyz[3])
  return k
end

"""
This function returns the index `freqidx` of the frequency given by the mixing
factors `mx`, `my`, and `mz` with respect to the frequency list
`freq = frequencies(bSF)`.
"""
function mixFactorToFreqIdx(b::MPIFile,mx,my,mz=0)
  freqidx = mixFactorToFreq(b,mx,my,mz)
  return freqidx + 1
end

function calcPrefactors(b::MPIFile)
  freqNumber = rxNumFrequencies(b)
  mask = collect((dfStrength(b)[1,:,1] .>= 0.0000001))
  divider = vec(dfDivider(b))
  #mxyz = round.(Int64,divider.*mask./gcd(divider.*mask))
  mxyz_ = 2 .*rxBandwidth(b)*dfCycle(b)./divider
  mxyz = max.(1,round.(Int64,mxyz_.*mask))

  return mxyz, mask, freqNumber
end

"""
This function returns a lookup table with columns
`freqidx` | `mx` | `my` | `mz` | `mixingOrder`
for all frequencies in `freq = frequencies(bSF)`, where only the lowest order
mixing coefficients `mx`, `my`, and `mz` are listed.
"""
function mixingFactors(b::MPIFile)
  mxyz, mask, freqNumber = calcPrefactors(b)
  MoList = zeros(Int64,freqNumber,4)
  MoList[:,4] .= -1 # set all mixing orders to -1 initially to change them later
  if length(mxyz) == 3
    Nx,Ny,Nz = round.(Int64,freqNumber./mxyz.*mask)

    return _mixingFactors(MoList, mxyz, Nx,Ny,Nz, freqNumber)
  else
    MoList[:,1] = MoList[:,4] = 0:(freqNumber-1)
    return MoList
  end
end

function _mixingFactors(MoList, mxyz, Nx,Ny,Nz, freqNumber)
  for mx = -Nx:Nx
    #for my = abs(mx)-n0:n0-abs(mx)
    for my = -Ny:Ny
      for mz = -Nz:Nz
        k = (mx*mxyz[1]+my*mxyz[2]+mz*mxyz[3])+1

          if k>=1 &&
             k<=freqNumber &&
             (MoList[k,4]<0 || MoList[k,4]>=abs(mx)+abs(my)+abs(mz))
           MoList[k,1] = mx
           MoList[k,2] = my
           MoList[k,3] = mz
           MoList[k,4] = abs(mx)+abs(my)+abs(mz)
          end
       end
    end
  end
  return MoList
end


text:
#=
callback.jl:

Author: Peter Wolf (opus111@gmail.com)

Port of the FastAI V2 Callback API to Julia

This code is inspired by FastAI, but differs from it in important ways

The original source is here

https://github.com/fastai/fastai2/blob/master/fastai2/callback/core.py

The documentation is copied from here

https://github.com/fastai/fastai2/blob/master/docs/callback.core.html
=#

"called before train and validate, ideal for initial setup."
function before_fit(cb::AbstractCallback,lrn::AbstractLearner, epoch_count, batch_size) end
"called at the after train and validate, for final clean-up."
function after_fit(cb::AbstractCallback,lrn::AbstractLearner) end
"called after cancelling train and validate."
function after_cancel_fit(cb::AbstractCallback,lrn::AbstractLearner) end

"called before each epoch"
function before_epoch(cb::AbstractCallback,lrn::AbstractLearner, epoch) end
"called after each epoch"
function after_epoch(cb::AbstractCallback,lrn::AbstractLearner, epoch) end
"called after cancelling epoch"
function after_cancel_epoch(cb::AbstractCallback,lrn::AbstractLearner, epoch) end

"called at the beginning of the training part of an epoch."
function before_epoch_train(cb::AbstractCallback,lrn::AbstractLearner,epoch) end
"called at the end of the training phase of an epoch."
function after_epoch_train(cb::AbstractCallback,lrn::AbstractLearner,epoch) end
"called after cancelling the training phase of an epoch."
function after_cancel_epoch_train(cb::AbstractCallback,lrn::AbstractLearner,epoch) end

"called at the beginning of the validate part of each epoch."
function before_epoch_validate(cb::AbstractCallback,lrn::AbstractLearner,epoch) end
"called at the end of the validation phase of each epoch."
function after_epoch_validate(cb::AbstractCallback,lrn::AbstractLearner,epoch) end
"called after cancelling validation of an epoch."
function after_cancel_epoch_validate(cb::AbstractCallback,lrn::AbstractLearner,epoch) end

"called at the beginning of each training each batch, just after drawing said batch. It can be used to do any setup necessary for the batch (like hyper-parameter scheduling) or to change the input/target before it goes in the model (change of the input with techniques like mixup for instance)."
function before_batch_train(cb::AbstractCallback,lrn::AbstractLearner,epoch, batch) end
"called at the end of training a batch, for any clean-up before the next one."
function after_batch_train(cb::AbstractCallback,lrn::AbstractLearner,epoch, batch) end
"report the loss of training a batch"
function batch_train_loss(cb::AbstractCallback,lrn::AbstractLearner,loss,epoch, batch) end
"called after cancelling training a batch."
function after_cancel_batch_train(cb::AbstractCallback,lrn::AbstractLearner,epoch, batch) end

"called at the beginning of each batch validate."
function before_batch_validate(cb::AbstractCallback,lrn::AbstractLearner,epoch, batch) end
"called at the end of a batch validate."
function after_batch_validate(cb::AbstractCallback,lrn::AbstractLearner,epoch, batch) end
"called to report the total loss of a batch validate"
function batch_validate_loss(cb::AbstractCallback,lrn::AbstractLearner,loss,epoch, batch) end
"called after cancelling a batch validate."
function after_cancel_batch_validate(cb::AbstractCallback,lrn::AbstractLearner,epoch, batch) end

struct DummyCallback <: AbstractCallback end

function before_fit(cb::DummyCallback,lrn::AbstractLearner, epoch_count, batch_size) println("Before Fit") end
function after_fit(cb::DummyCallback,lrn::AbstractLearner) println("After Fit") end
function after_cancel_fit(cb::DummyCallback,lrn::AbstractLearner) println("After Cancel Fit") end
function before_epoch(cb::DummyCallback,lrn::AbstractLearner, epoch) println("Before Epoch") end
function after_epoch(cb::DummyCallback,lrn::AbstractLearner, epoch) println("After Epoch") end
function after_cancel_epoch(cb::DummyCallback,lrn::AbstractLearner, epoch) println("After Cancel Epoch") end
function before_epoch_train(cb::DummyCallback,lrn::AbstractLearner,epoch) println("\nBefore Epoch Train") end
function after_epoch_train(cb::DummyCallback,lrn::AbstractLearner,epoch) println("\nAfter Epoch Train") end
function after_cancel_epoch_train(cb::DummyCallback,lrn::AbstractLearner,epoch) println("\nAfter Cancel Epoch Train") end
function before_epoch_validate(cb::DummyCallback,lrn::AbstractLearner,epoch) println("\nBefore Epoch Validate") end
function after_epoch_validate(cb::DummyCallback,lrn::AbstractLearner,epoch) println("\nAfter Epoch Validate") end
function after_cancel_epoch_validate(cb::DummyCallback,lrn::AbstractLearner,epoch) println("\nAfter Cancel Epoch Validate") end
function before_batch_train(cb::DummyCallback,lrn::AbstractLearner,epoch, batch) print("Before Batch Train,") end
function after_batch_train(cb::DummyCallback,lrn::AbstractLearner,epoch, batch) print("After Batch Train,")  end
function batch_train_loss(cb::DummyCallback,lrn::AbstractLearner,loss, epoch, batch) print("Batch Train Loss = $(loss)")  end
function after_cancel_batch_train(cb::DummyCallback,lrn::AbstractLearner,epoch, batch) print("After Cancel Batch Train")  end
function before_batch_validate(cb::DummyCallback,lrn::AbstractLearner,epoch, batch) print("Before Batch Validate,")  end
function after_batch_validate(cb::DummyCallback,lrn::AbstractLearner,epoch, batch) print("After Batch Validate") end
function batch_validate_loss(cb::DummyCallback,lrn::AbstractLearner,loss,epoch, batch) println("Batch Validate Loss = $(loss)") end
function after_cancel_batch_validate(cb::DummyCallback,lrn::AbstractLearner,epoch, batch) println("After Cancel Batch Validate") end

"`Callback` that tracks the epoch and batch and calculates progress (fraction done)"
mutable struct ProgressCallback <: AbstractCallback
    epoch::Int
    batch::Int
    n_batch::Int
    n_epoch::Int
end

ProgressCallback() = ProgressCallback(0,0,0,0)

function before_fit(tecb::ProgressCallback,learn::AbstractLearner, n_epoch)
    tecb.n_epoch = n_epoch
    tecb.n_batch = 0
    tecb.epoch = 0
    tecb.batch = 0
end

function after_epoch_train(tecb::ProgressCallback,lrn::AbstractLearner, epoch)
    tecb.n_batch = batch
    tecb.epoch = epoch
end

function after_batch_train(tecb::ProgressCallback,lrn::AbstractLearner, epoch, batch)
    tecb.batch = batch
end

function after_batch_validate_(tecb::ProgressCallback,lrn::AbstractLearner, loss, epoch, batch)
    println("Loss=$(loss) Amount trained = $(progress(tecb))")
end

function progress(tecb::ProgressCallback)
    num = ((tecb.n_epoch-1)*tecb.n_batch + tecb.batch)
    deom = float(tecb.n_epoch*tecb.n_batch)
    return denom > 0.0 ? num / denom : 0.0
end


text:
const branin = let res_init=zeros(2), jac_init=zeros(2,2), x_init=[2.0, 0.0]

    function res(x, r)
        r[1] = 4.0 * (x[1] + x[2])
        r[2] = r[1] + (x[1] - x[2]) * ((x[1] - 2.0) ^ 2 + x[2] ^ 2 - 1.0)
        return r
    end

    function jac(x, jac)
        jac[1, 1] = 4.0
        jac[1, 2] = 4.0
        jac[2, 1] = 3.0 + (x[1] - 2.0) * (3.0 * x[1] - 2.0 * x[2] - 2.0) + x[2] * x[2]
        jac[2, 2] = 1.0 + 2.0 * (2.0*x[1] - x[2] * x[2]) - (x[1] - x[2])^2
        return jac
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("branin", res, jac,
                                                            res_init, jac_init, x_init;
                                                            scale=scale, verbose=verbose,
                                                            print_steps=print_steps)
end


text:
module Core
    using DataStructures

    include("../io/internal.jl")
    include("../io/output.jl")

    mutable struct MOEUM
        """
        An Tabular Struct for Massive Data Processing

        **Common methods**
        """

        ###parameters
        #arguments
        dat_sori::Array
        hol_sori::Array
        name::String
        #hidden parameters
        _sok_dict::DataStructures.OrderedDict
        _sok_arr::Array

        ###functions
        #modify
        insert::Function
        append::Function

        #simple selection and query
        select::Function
        where::Function
        equation::Function

        #descriptive statistics
        describe::Function

        #output
        to_string::Function
        to_dataframe::Function
        to_dict::Function
        to_csv::Function
        to_db::Function

        #init
        function MOEUM(;dat_sori = [], hol_sori = [], name = "moeum.MOEUM")
            #init instance
            instance = new(dat_sori, hol_sori, name)

            instance._sok_dict, instance._sok_arr = internal.set_sok(instance, instance.dat_sori, instance.hol_sori)

            #modify
            instance.insert = function(;dat_name::Any = "", hol_value::Array = [], dat_pos::Int = length(instance.dat_sori) + 1)
                return internal.insert(instance, dat_pos, dat_name, hol_value)
            end

            instance.append = function(;dat_name::Any = "", hol_value::Array = [])
                return internal.append(instance, dat_name, hol_value)
            end

            #simple selection and query
            instance.select = function(key)
                return internal.select(instance, key)
            end

            instance.where = function(where_string; drop_null::Bool = true)
                return internal.where(instance, where_string, drop_null)
            end

            instance.equation = function(equation_string; inplace::Bool = true)
                return internal.equation(instance, equation_string, inplace)
            end

            #descriptive statistics
            instance.describe = function()
                return internal.describe(instance)
            end

            #output functions
            instance.to_string = function()
                return output.to_string(instance)
            end

            instance.to_dataframe = function()
                return output.to_dataframe(instance)
            end

            instance.to_dict = function(;orientation::String = "Dict")
                return output.to_dict(instance, orientation)
            end

            instance.to_csv = function(;csv_path::String = joinpath(pwd(), string(instance.name, ".csv")))
                return output.to_csv(instance, csv_path)
            end

            return instance
        end
    end
    function Base.show(io::Base.IO, instance::MOEUM)
        println(io, instance.to_string())
    end
end


text:
module ArrayInterface

using ArrayInterfaceCore
import ArrayInterfaceCore: allowed_getindex, allowed_setindex!, aos_to_soa, buffer,
    parent_type, fast_matrix_colors, findstructralnz, has_sparsestruct,
    issingular, isstructured, matrix_colors, restructure, lu_instance,
    safevec, zeromatrix, ColoringAlgorithm, fast_scalar_indexing, parameterless_type,
    ndims_index, ndims_shape, is_splat_index, is_forwarding_wrapper, IndicesInfo

# ArrayIndex subtypes and methods
import ArrayInterfaceCore: ArrayIndex, MatrixIndex, VectorIndex, BidiagonalIndex, TridiagonalIndex
# managing immutables
import ArrayInterfaceCore: ismutable, can_change_size, can_setindex
# constants
import ArrayInterfaceCore: MatAdjTrans, VecAdjTrans, UpTri, LoTri
# device pieces
import ArrayInterfaceCore: AbstractDevice, AbstractCPU, CPUPointer, CPUTuple, CheckParent,
    CPUIndex, GPU, can_avx

import ArrayInterfaceCore: known_first, known_step, known_last

using Static
using Static: Zero, One, nstatic, eq, ne, gt, ge, lt, le, eachop, eachop_tuple,
    permute, invariant_permutation, field_type, reduce_tup, find_first_eq

using IfElse

using Base.Cartesian
using Base: @propagate_inbounds, tail, OneTo, LogicalIndex, Slice, ReinterpretArray,
    ReshapedArray, AbstractCartesianIndex

using Base.Iterators: Pairs
using LinearAlgebra

import Compat

n_of_x(::StaticInt{N}, x::X) where {N,X} = ntuple(Compat.Returns(x), Val{N}())

_add1(@nospecialize x) = x + oneunit(x)
_sub1(@nospecialize x) = x - oneunit(x)

@generated function merge_tuple_type(::Type{X}, ::Type{Y}) where {X<:Tuple,Y<:Tuple}
    Tuple{X.parameters...,Y.parameters...}
end

const CanonicalInt = Union{Int,StaticInt}
canonicalize(x::Integer) = Int(x)
canonicalize(@nospecialize(x::StaticInt)) = x

abstract type AbstractArray2{T,N} <: AbstractArray{T,N} end

Base.size(A::AbstractArray2) = map(Int, ArrayInterface.size(A))
Base.size(A::AbstractArray2, dim) = Int(ArrayInterface.size(A, dim))

function Base.axes(A::AbstractArray2)
    is_forwarding_wrapper(A) && return ArrayInterface.axes(parent(A))
    throw(ArgumentError("Subtypes of `AbstractArray2` must define an axes method"))
end
function Base.axes(A::AbstractArray2, dim::Union{Symbol,StaticSymbol})
    axes(A, to_dims(A, dim))
end

function Base.strides(A::AbstractArray2)
    defines_strides(A) && return map(Int, ArrayInterface.strides(A))
    throw(MethodError(Base.strides, (A,)))
end
Base.strides(A::AbstractArray2, dim) = Int(ArrayInterface.strides(A, dim))

function Base.IndexStyle(::Type{T}) where {T<:AbstractArray2}
    is_forwarding_wrapper(T) ? IndexStyle(parent_type(T)) : IndexCartesian()
end

function Base.length(A::AbstractArray2)
    len = known_length(A)
    if len === nothing
        return Int(prod(size(A)))
    else
        return Int(len)
    end
end

@propagate_inbounds Base.getindex(A::AbstractArray2, args...) = getindex(A, args...)
@propagate_inbounds Base.getindex(A::AbstractArray2; kwargs...) = getindex(A; kwargs...)

@propagate_inbounds function Base.setindex!(A::AbstractArray2, val, args...)
    return setindex!(A, val, args...)
end
@propagate_inbounds function Base.setindex!(A::AbstractArray2, val; kwargs...)
    return setindex!(A, val; kwargs...)
end

@inline static_first(x) = Static.maybe_static(known_first, first, x)
@inline static_last(x) = Static.maybe_static(known_last, last, x)
@inline static_step(x) = Static.maybe_static(known_step, step, x)

@inline function _to_cartesian(a, i::CanonicalInt)
    @inbounds(CartesianIndices(ntuple(dim -> indices(a, dim), Val(ndims(a))))[i])
end
@inline function _to_linear(a, i::Tuple{CanonicalInt,Vararg{CanonicalInt}})
    _strides2int(offsets(a), size_to_strides(size(a), static(1)), i) + static(1)
end

"""
    has_parent(::Type{T}) -> StaticBool

Returns `static(true)` if `parent_type(T)` a type unique to `T`.
"""
has_parent(x) = has_parent(typeof(x))
has_parent(::Type{T}) where {T} = _has_parent(parent_type(T), T)
_has_parent(::Type{T}, ::Type{T}) where {T} = False()
_has_parent(::Type{T1}, ::Type{T2}) where {T1,T2} = True()

"""
    device(::Type{T}) -> AbstractDevice

Indicates the most efficient way to access elements from the collection in low-level code.
For `GPUArrays`, will return `ArrayInterface.GPU()`.
For `AbstractArray` supporting a `pointer` method, returns `ArrayInterface.CPUPointer()`.
For other `AbstractArray`s and `Tuple`s, returns `ArrayInterface.CPUIndex()`.
Otherwise, returns `nothing`.
"""
device(A) = device(typeof(A))
device(::Type) = nothing
device(::Type{<:Tuple}) = CPUTuple()
device(::Type{T}) where {T<:Array} = CPUPointer()
device(::Type{T}) where {T<:AbstractArray} = _device(has_parent(T), T)
function _device(::True, ::Type{T}) where {T}
    if defines_strides(T)
        return device(parent_type(T))
    else
        return _not_pointer(device(parent_type(T)))
    end
end
_not_pointer(::CPUPointer) = CPUIndex()
_not_pointer(x) = x
_device(::False, ::Type{T}) where {T<:DenseArray} = CPUPointer()
_device(::False, ::Type{T}) where {T} = CPUIndex()

"""
    is_lazy_conjugate(::AbstractArray) -> Bool

Determine if a given array will lazyily take complex conjugates, such as with `Adjoint`. This will work with
nested wrappers, so long as there is no type in the chain of wrappers such that `parent_type(T) == T`

Examples

    julia> a = transpose([1 + im, 1-im]')
    2×1 transpose(adjoint(::Vector{Complex{Int64}})) with eltype Complex{Int64}:
     1 - 1im
     1 + 1im

    julia> ArrayInterface.is_lazy_conjugate(a)
    True()

    julia> b = a'
    1×2 adjoint(transpose(adjoint(::Vector{Complex{Int64}}))) with eltype Complex{Int64}:
     1+1im  1-1im

    julia> ArrayInterface.is_lazy_conjugate(b)
    False()

"""
is_lazy_conjugate(::T) where {T<:AbstractArray} = _is_lazy_conjugate(T, False())
is_lazy_conjugate(::AbstractArray{<:Real}) = False()

function _is_lazy_conjugate(::Type{T}, isconj) where {T<:AbstractArray}
    Tp = parent_type(T)
    if T !== Tp
        _is_lazy_conjugate(Tp, isconj)
    else
        isconj
    end
end

function _is_lazy_conjugate(::Type{T}, isconj) where {T<:Adjoint}
    Tp = parent_type(T)
    if T !== Tp
        _is_lazy_conjugate(Tp, !isconj)
    else
        !isconj
    end
end

"""
    insert(collection, index, item)

Returns a new instance of `collection` with `item` inserted into at the given `index`.
"""
Base.@propagate_inbounds function insert(collection, index, item)
    @boundscheck checkbounds(collection, index)
    ret = similar(collection, length(collection) + 1)
    @inbounds for i = firstindex(ret):(index-1)
        ret[i] = collection[i]
    end
    @inbounds ret[index] = item
    @inbounds for i = (index+1):lastindex(ret)
        ret[i] = collection[i-1]
    end
    return ret
end

function insert(x::Tuple{Vararg{Any,N}}, index, item) where {N}
    @boundscheck if !checkindex(Bool, StaticInt{1}():StaticInt{N}(), index)
        throw(BoundsError(x, index))
    end
    return unsafe_insert(x, Int(index), item)
end

@inline function unsafe_insert(x::Tuple, i::Int, item)
    if i === 1
        return (item, x...)
    else
        return (first(x), unsafe_insert(tail(x), i - 1, item)...)
    end
end

"""
    deleteat(collection, index)

Returns a new instance of `collection` with the item at the given `index` removed.
"""
Base.@propagate_inbounds function deleteat(collection::AbstractVector, index)
    @boundscheck if !checkindex(Bool, eachindex(collection), index)
        throw(BoundsError(collection, index))
    end
    return unsafe_deleteat(collection, index)
end
Base.@propagate_inbounds function deleteat(collection::Tuple{Vararg{Any,N}}, index) where {N}
    @boundscheck if !checkindex(Bool, StaticInt{1}():StaticInt{N}(), index)
        throw(BoundsError(collection, index))
    end
    return unsafe_deleteat(collection, index)
end

function unsafe_deleteat(src::AbstractVector, index)
    dst = similar(src, length(src) - 1)
    @inbounds for i in indices(dst)
        if i < index
            dst[i] = src[i]
        else
            dst[i] = src[i+1]
        end
    end
    return dst
end

@inline function unsafe_deleteat(src::AbstractVector, inds::AbstractVector)
    dst = similar(src, length(src) - length(inds))
    dst_index = firstindex(dst)
    @inbounds for src_index in indices(src)
        if !in(src_index, inds)
            dst[dst_index] = src[src_index]
            dst_index += one(dst_index)
        end
    end
    return dst
end

@inline function unsafe_deleteat(src::Tuple, inds::AbstractVector)
    dst = Vector{eltype(src)}(undef, length(src) - length(inds))
    dst_index = firstindex(dst)
    @inbounds for src_index in static(1):length(src)
        if !in(src_index, inds)
            dst[dst_index] = src[src_index]
            dst_index += one(dst_index)
        end
    end
    return Tuple(dst)
end

@inline unsafe_deleteat(x::Tuple{T}, i) where {T} = ()
@inline unsafe_deleteat(x::Tuple{T1,T2}, i) where {T1,T2} =
    isone(i) ? (x[2],) : (x[1],)
@inline function unsafe_deleteat(x::Tuple, i)
    if i === one(i)
        return tail(x)
    elseif i == length(x)
        return Base.front(x)
    else
        return (first(x), unsafe_deleteat(tail(x), i - one(i))...)
    end
end

include("array_index.jl")
include("ranges.jl")
include("axes.jl")
include("size.jl")
include("dimensions.jl")
include("indexing.jl")
include("stridelayout.jl")
include("broadcast.jl")

end


text:
# First build the animals directory - mistake included to show how easy to fix
# julia> mkpath("data/animals/invertebrates/arthropods/crustaceans")
# "data/animals/invertebrates/arthropods/crustaceans"

# julia> mkpath("data/animals/invertebrates/arthropods/insects")
# "data/animals/invertebrates/arthropods/insects"

# julia> mkpath("data/animals/invertebrates/flatworms")
# "data/animals/invertebrates/flatworms"

# julia> mkpath("data/animals/invertebrates/molluscs")
# "data/animals/invertebrates/molluscs"

# julia> mkpath("data/animals/invertebrates/amphibians")
# "data/animals/invertebrates/amphibians"

# julia> rm("data/animals/invertebrates/amphibians")

# julia> mkpath("data/animals/vertebrates/amphibians")
# "data/animals/vertebrates/amphibians"

# julia> mkpath("data/animals/vertebrates/birds")
# "data/animals/vertebrates/birds"

# julia> mkpath("data/animals/vertebrates/fish")
# "data/animals/vertebrates/fish"

# julia> mkpath("data/animals/vertebrates/mammals")
# "data/animals/vertebrates/mammals"

# julia> for (root, dirs, files) in walkdir(an)
#            println("Directories in $root")
#            for dir in dirs
#                println(joinpath(root, dir)) # path to directories
#            end
#            println("Files in $root")
#            for file in files
#                println(joinpath(root, file)) # path to files
#            end
#        end
# Directories in data/animals
# data/animals/invertebrates
# data/animals/vertebrates
# Files in data/animals
# Directories in data/animals/invertebrates
# data/animals/invertebrates/arthropods
# data/animals/invertebrates/flatworms
# data/animals/invertebrates/molluscs
# Files in data/animals/invertebrates
# Directories in data/animals/invertebrates/arthropods
# data/animals/invertebrates/arthropods/crustaceans
# data/animals/invertebrates/arthropods/insects
# Files in data/animals/invertebrates/arthropods
# Directories in data/animals/invertebrates/arthropods/crustaceans
# Files in data/animals/invertebrates/arthropods/crustaceans
# Directories in data/animals/invertebrates/arthropods/insects
# Files in data/animals/invertebrates/arthropods/insects
# Directories in data/animals/invertebrates/flatworms
# Files in data/animals/invertebrates/flatworms
# Directories in data/animals/invertebrates/molluscs
# Files in data/animals/invertebrates/molluscs
# Directories in data/animals/vertebrates
# data/animals/vertebrates/amphibians
# data/animals/vertebrates/birds
# data/animals/vertebrates/fish
# data/animals/vertebrates/mammals
# Files in data/animals/vertebrates
# Directories in data/animals/vertebrates/amphibians
# Files in data/animals/vertebrates/amphibians
# Directories in data/animals/vertebrates/birds
# Files in data/animals/vertebrates/birds
# Directories in data/animals/vertebrates/fish
# Files in data/animals/vertebrates/fish
# Directories in data/animals/vertebrates/mammals
# Files in data/animals/vertebrates/mammals

# Now create the file hierarchy
# julia> ver = "data/animals/vertebrates"
# "data/animals/vertebrates"

# julia> cd(joinpath(ver,"amphibians")) do
#          foreach(touch, ["frog", "salamander"])
#        end

# julia> cd(joinpath(ver,"birds")) do
#          foreach(touch, ["crow", "mockingjay","seagull"])
#        end

# julia> cd(joinpath(ver,"mammals")) do
#          foreach(touch, ["cow", "human"])
#        end

# julia> for (root, dirs, files) in walkdir(an)
#            println("Directories in $root")
#            for dir in dirs
#                println(joinpath(root, dir)) # path to directories
#            end
#            println("Files in $root")
#            for file in files
#                println(joinpath(root, file)) # path to files
#            end
#        end
# Directories in data/animals
# data/animals/invertebrates
# data/animals/vertebrates
# Files in data/animals
# Directories in data/animals/invertebrates
# data/animals/invertebrates/arthropods
# data/animals/invertebrates/flatworms
# data/animals/invertebrates/molluscs
# Files in data/animals/invertebrates
# Directories in data/animals/invertebrates/arthropods
# data/animals/invertebrates/arthropods/crustaceans
# data/animals/invertebrates/arthropods/insects
# Files in data/animals/invertebrates/arthropods
# Directories in data/animals/invertebrates/arthropods/crustaceans
# Files in data/animals/invertebrates/arthropods/crustaceans
# Directories in data/animals/invertebrates/arthropods/insects
# Files in data/animals/invertebrates/arthropods/insects
# Directories in data/animals/invertebrates/flatworms
# Files in data/animals/invertebrates/flatworms
# Directories in data/animals/invertebrates/molluscs
# Files in data/animals/invertebrates/molluscs
# Directories in data/animals/vertebrates
# data/animals/vertebrates/amphibians
# data/animals/vertebrates/birds
# data/animals/vertebrates/fish
# data/animals/vertebrates/mammals
# Files in data/animals/vertebrates
# Directories in data/animals/vertebrates/amphibians
# Files in data/animals/vertebrates/amphibians
# data/animals/vertebrates/amphibians/frog
# data/animals/vertebrates/amphibians/salamander
# Directories in data/animals/vertebrates/birds
# Files in data/animals/vertebrates/birds
# data/animals/vertebrates/birds/crow
# data/animals/vertebrates/birds/mockingjay
# data/animals/vertebrates/birds/seagull
# Directories in data/animals/vertebrates/fish
# Files in data/animals/vertebrates/fish
# Directories in data/animals/vertebrates/mammals
# Files in data/animals/vertebrates/mammals
# data/animals/vertebrates/mammals/cow
# data/animals/vertebrates/mammals/human

# It seems Engheim implements visit files so as not to get into the deatails of walkdir
# however walkdir serves this exect purpose and can do the same ans so much more
# function visitfiles(fn, root::AbstractString)
#     if isfile(root)
#         fn(root)
#         return
#     end

#     cd(root) do
#         for file in readdir()
#             visitfiles(fn, file)
#         end
#     end
# end

# Example

# julia> visitfiles(an) do file
#     println(joinpath(pwd(), file))
# end
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/amphibians/frog
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/amphibians/salamander
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/birds/crow
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/birds/mockingjay
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/birds/seagull
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/mammals/cow
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/mammals/human

# julia> for (root, dirs, files) in walkdir(an)
#            for file in files
#                println(joinpath(pwd(), joinpath(root,file))) # path to files
#            end
#        end
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/amphibians/frog
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/amphibians/salamander
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/birds/crow
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/birds/mockingjay
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/birds/seagull
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/mammals/cow
# /Users/aronet/Code/FourM/Study/Julia/JuliaforBeginners/data/animals/vertebrates/mammals/human
function replace_animal(animal::AbstractString)
    rm(animal)
    mkdir(animal)
    cd(animal) do
        foreach(touch, ["description.txt", "looks.jpg"])
    end
end

# Example
# julia> for (root, dirs, files) in walkdir(an)
#            for file in files
#                replace_animal(joinpath(pwd(), joinpath(root,file)))
#            end
#        end

# julia> for (root, dirs, files) in walkdir(an)
#            println("Directories in $root")
#            for dir in dirs
#                println(joinpath(root, dir)) # path to directories
#            end
#            println("Files in $root")
#            for file in files
#                println(joinpath(root, file)) # path to files
#            end
#        end
# Directories in data/animals
# data/animals/invertebrates
# data/animals/vertebrates
# Files in data/animals
# Directories in data/animals/invertebrates
# data/animals/invertebrates/arthropods
# data/animals/invertebrates/flatworms
# data/animals/invertebrates/molluscs
# Files in data/animals/invertebrates
# Directories in data/animals/invertebrates/arthropods
# data/animals/invertebrates/arthropods/crustaceans
# data/animals/invertebrates/arthropods/insects
# Files in data/animals/invertebrates/arthropods
# Directories in data/animals/invertebrates/arthropods/crustaceans
# Files in data/animals/invertebrates/arthropods/crustaceans
# Directories in data/animals/invertebrates/arthropods/insects
# Files in data/animals/invertebrates/arthropods/insects
# Directories in data/animals/invertebrates/flatworms
# Files in data/animals/invertebrates/flatworms
# Directories in data/animals/invertebrates/molluscs
# Files in data/animals/invertebrates/molluscs
# Directories in data/animals/vertebrates
# data/animals/vertebrates/amphibians
# data/animals/vertebrates/birds
# data/animals/vertebrates/fish
# data/animals/vertebrates/mammals
# Files in data/animals/vertebrates
# Directories in data/animals/vertebrates/amphibians
# data/animals/vertebrates/amphibians/frog
# data/animals/vertebrates/amphibians/salamander
# Files in data/animals/vertebrates/amphibians
# Directories in data/animals/vertebrates/amphibians/frog
# Files in data/animals/vertebrates/amphibians/frog
# data/animals/vertebrates/amphibians/frog/description.txt
# data/animals/vertebrates/amphibians/frog/looks.jpg
# Directories in data/animals/vertebrates/amphibians/salamander
# Files in data/animals/vertebrates/amphibians/salamander
# data/animals/vertebrates/amphibians/salamander/description.txt
# data/animals/vertebrates/amphibians/salamander/looks.jpg
# Directories in data/animals/vertebrates/birds
# data/animals/vertebrates/birds/crow
# data/animals/vertebrates/birds/mockingjay
# data/animals/vertebrates/birds/seagull
# Files in data/animals/vertebrates/birds
# Directories in data/animals/vertebrates/birds/crow
# Files in data/animals/vertebrates/birds/crow
# data/animals/vertebrates/birds/crow/description.txt
# data/animals/vertebrates/birds/crow/looks.jpg
# Directories in data/animals/vertebrates/birds/mockingjay
# Files in data/animals/vertebrates/birds/mockingjay
# data/animals/vertebrates/birds/mockingjay/description.txt
# data/animals/vertebrates/birds/mockingjay/looks.jpg
# Directories in data/animals/vertebrates/birds/seagull
# Files in data/animals/vertebrates/birds/seagull
# data/animals/vertebrates/birds/seagull/description.txt
# data/animals/vertebrates/birds/seagull/looks.jpg
# Directories in data/animals/vertebrates/fish
# Files in data/animals/vertebrates/fish
# Directories in data/animals/vertebrates/mammals
# data/animals/vertebrates/mammals/cow
# data/animals/vertebrates/mammals/human
# Files in data/animals/vertebrates/mammals
# Directories in data/animals/vertebrates/mammals/cow
# Files in data/animals/vertebrates/mammals/cow
# data/animals/vertebrates/mammals/cow/description.txt
# data/animals/vertebrates/mammals/cow/looks.jpg
# Directories in data/animals/vertebrates/mammals/human
# Files in data/animals/vertebrates/mammals/human
# data/animals/vertebrates/mammals/human/description.txt
# data/animals/vertebrates/mammals/human/looks.jpg
function replace_github_math(filename)
    out = IOBuffer()
    open(filename) do io
        # need to change Engheim's version to handle non-markdown text at the end of the file
        # as in our example book.txt
        while true
            # read and then write out anything till first '$' you find and swallow '$'
            s = readuntil(io, '$')
            if eof(io)
                write(out, s)
                break
            end
            write(out, s)
            # this tests if the line you just read is Markua
            # considering that the exercise is converting Pandoc to Markua this is confusing.
            # why would there be markua and Pandoc in same file?
            # Perhaps Engheim wanted to illustrate the continue
            # So I am commenting this out and giving two alternate functions
            #if s[end] == '`'
            #    write(out, '$')
            #    continue
            #end

            mark(io)
            s = readuntil(io, '$')

            if isempty(s)
                reset(io)
                # need to use raw since $ is for string interpolation
                s = readuntil(io, raw"$$")
                write(out, '$')
                write(out, s)
                write(out, raw"$$")
                continue
            end

            write(out, '`')
            write(out, s)
            write(out, raw"`$")
        end
    end
    write_file(out, "markau", filename)
end

function replace_markau_math(filename)
    out = IOBuffer()
    open(filename) do io
        # need to change Engheim's version to handle non-markdown text at the end of the file
        # as in our example book.txt
        while true
            # read and then write out anything till first '$' you find and swallow '$'
            s = readuntil(io, '`')
            if eof(io)
                write(out, s)
                break
            end
            write(out, s)
            s = readuntil(io, '`')
            write(out, '$')
            write(out, s)
        end
    end
    write_file(out, "github", filename)
end

function write_file(out::IOBuffer, filetype, filename)
    # write out buffer
    seekstart(out)
    # write out to a different filename, don't overwrite original
    f = splitext(filename)
    rfilename = f[1]*filetype*f[2]
    s = read(out, String)
    open(rfilename, "w") do io
        write(io, s)
    end
end

# Example
# julia> include("bash.jl")
# write_file (generic function with 2 methods)
# Convert to Markau
# julia> replace_github_math("data/book.txt")
# 104
# Convert back to GITHUB
# julia> replace_markau_math("data/bookmarkau.txt")
# 103
# If it works properly the last version should b identical to the original
# (base) ➜  JuliaforBeginners git:(main) ✗ diff data/book.txt data/bookmarkaugithub.txt
# (base) ➜  JuliaforBeginners git:(main) ✗ $ Yay! It is!

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
export tarjan

mutable struct Node
    id::Int
    index::Int
    lowlink::Int
    onstack::Bool
    tails::Vector{Node}
end

doc"""
    tarjan(heads::Vector{Int})

Computes Tarjan's algorithm for finding strongly connected components (cycles) of a graph.

# Example
```julia
heads = [3, 4, 2, 1, 4, 5]
tarjan(heads)
```
"""
function tarjan(heads::Vector{Int})
    for h in heads
        1 <= h <= length(heads) || throw("Invalid input: $h")
    end

    nodes = [Node(i,0,0,false,Node[]) for i=1:length(heads)]
    for i = 1:length(heads)
        push!(nodes[heads[i]].tails, nodes[i])
    end
    stack = Node[]
    cycles = Vector{Int}[]
    index = 1
    function strongconnect!(node::Node)
        node.index = index
        node.lowlink = index
        index += 1
        push!(stack, node)
        node.onstack = true

        for t in node.tails
            if t.index == 0
                # Successor t has not yet been visited
                strongconnect!(t)
                node.lowlink = min(node.lowlink, t.lowlink)
            elseif t.onstack
                # Successor t is in stack and hence in the current SCC
                node.lowlink = min(node.lowlink, t.index)
            end
        end

        # If node is a root, pop the stack and generate an SCC
        if node.lowlink == node.index
            cycle = Int[]
            while true
                t = pop!(stack)
                t.onstack = false
                push!(cycle, t.id)
                t == node && break
            end
            length(cycle) > 1 && push!(cycles,sort!(cycle))
        end
    end

    for node in nodes
        node.index == 0 && strongconnect!(node)
    end
    cycles
end


text:
import Pkg
println(pwd())
if !isfile("Project.toml")
    println("Project.toml が見つからない。")
    exit(1)
end
Pkg.activate(".")
using HtmlTable
html = fromCSV("./test/data.csv", title="From CSV")
write("./test/result_csv.html", html)
println("test/result_csv.html が更新されました。")

max_stars_repo_name: makandat/julia

max_issues_repo_name: makandat/julia

max_forks_repo_name: makandat/julia


text:
using LinearAlgebra, SimpleDifferentialOperators, DiffEqOperators, Test, Parameters

#parameters
params = @with_kw (
    μ = -0.1, # constant negative drift
    σ = 0.1,
    ρ = 0.05,
    M = 3, # size of grid (interior points)
    x̄ = range(0.0, 1.0, length = (M+2)),
    x = interiornodes(x̄), # i.e., x̄[2:end-1]
    S = 3.0
)
p = params();

#----------------------------------
#Testing For Negative Drifts
#----------------------------------
# payoff function
π(x) = x^2
# SimpleDifferentialOperators setup
function SDO_negative_drift(π, params)
    bc = (Reflecting(), Reflecting())
    Lₓ = params.μ*L₁₋bc(params.x̄, bc) + params.σ^2 / 2 * L₂bc(params.x̄, bc)
    L_bc = I * params.ρ - Lₓ
    v = L_bc \ π.(params.x);
    return (v = v, L_bc = L_bc, Lₓ = Lₓ, L₁₋bc = L₁₋bc(params.x̄, bc), L₂bc = L₂bc(params.x̄, bc))
end
# DiffEqOperators Setup
function DEO_negative_drift(π, params)
    dx = params.x[2] - params.x[1]
    # discretize L = ρ - μ D_x - σ^2 / 2 D_xx
    # subject to reflecting barriers at 0 and 1
    L1 = UpwindDifference(1,1,dx,params.M,1.)
    L2 = CenteredDifference(2,2,dx,params.M)
    Q = Neumann0BC(dx, 1)
    L₁₋bc = -1. .* Array(UpwindDifference(1,1,dx,params.M,-1.) * Q)[1]
    # Here Array(A::GhostDerivativeOperator) will return a tuple of the linear part
    # and the affine part of the operator A, hence we index Array(µ*L1*Q).
    # The operators in this example are purely linear, so we don't worry about Array(µ*L1*Q)[2]
    Lₓ = Array(params.μ*L1*Q)[1] + Array(params.σ^2/2*L2*Q)[1] # when scalar coefficient is properly implemented
    # Lₓ = Array(L1*Q)[1] + Array(params.σ^2/2*L2*Q)[1] # when scalar coefficient is properly implemented
    L_bc = I * params.ρ - Lₓ

    # solve the value function
    v = L_bc \ π.(params.x);
    return (v = v, L_bc = L_bc, Lₓ = Lₓ, L₁₋bc = L₁₋bc, L₂bc = Array(L2*Q)[1] )
end
@testset "Constant Negative Drifts" begin
    @test DEO_negative_drift(π, p).v ≈ SDO_negative_drift(π, p).v
    @test DEO_negative_drift(π, p).L_bc ≈ SDO_negative_drift(π, p).L_bc
    @test DEO_negative_drift(π, p).Lₓ ≈ SDO_negative_drift(π, p).Lₓ
    @test DEO_negative_drift(π, p).L₁₋bc ≈ SDO_negative_drift(π, p).L₁₋bc
    @test DEO_negative_drift(π, p).L₂bc ≈ SDO_negative_drift(π, p).L₂bc
end
p = params(μ = 0.1);
#----------------------------------
#Testing For Positive Drifts
#----------------------------------
# SimpleDifferentialOperators setup
function SDO_positive_drift(π, params)
    bc = (Reflecting(), Reflecting())
    Lₓ = params.μ*L₁₊bc(params.x̄, bc) + params.σ^2 / 2 * L₂bc(params.x̄, bc)
    L_bc = I * params.ρ - Lₓ
    v = L_bc \ π.(params.x);
    return (v = v, L_bc = L_bc, Lₓ = Lₓ, L₁₊bc = L₁₊bc(params.x̄, bc), L₂bc = L₂bc(params.x̄, bc))
end
# DiffEqOperators Setup
function DEO_positive_drift(π, params)
    dx = params.x[2] - params.x[1]
    # discretize L = ρ - μ D_x - σ^2 / 2 D_xx
    # subject to reflecting barriers at 0 and 1
    L1 = UpwindDifference(1,1,dx,params.M,1.0)
    L2 = CenteredDifference(2,2,dx,params.M)
    Q = Neumann0BC(dx, 1)
    # Here Array(A::GhostDerivativeOperator) will return a tuple of the linear part
    # and the affine part of the operator A, hence we index Array(µ*L1*Q).
    # The operators in this example are purely linear, so we don't worry about Array(µ*L1*Q)[2]
    Lₓ = Array(params.µ*L1*Q)[1] + Array(params.σ^2/2*L2*Q)[1]
    L_bc = I * params.ρ - Lₓ

    # solve the value function
    v = L_bc \ π.(params.x);
    return (v = v, L_bc = L_bc, Lₓ = Lₓ, L₁₊bc = Array(L1*Q)[1], L₂bc = Array(L2*Q)[1] )
end

@testset "Constant Positive Drifts" begin
    @test DEO_positive_drift(π, p).v ≈ SDO_positive_drift(π, p).v
    @test DEO_positive_drift(π, p).L_bc ≈ SDO_positive_drift(π, p).L_bc
    @test DEO_positive_drift(π, p).Lₓ ≈ SDO_positive_drift(π, p).Lₓ
    @test DEO_positive_drift(π, p).L₁₊bc ≈ SDO_positive_drift(π, p).L₁₊bc
    @test DEO_positive_drift(π, p).L₂bc ≈ SDO_positive_drift(π, p).L₂bc
end

#----------------------------------
#Testing For State Dependent Drifts
#----------------------------------
μ(x) = -x
# SimpleDifferentialOperators setup
function SDO_state_dependent_drift(π, μ, params)
    bc = (Reflecting(), Reflecting())
    L₁ = Diagonal(min.(μ.(params.x), 0.0)) * L₁₋bc(params.x̄, bc) + Diagonal(max.(μ.(params.x), 0.0)) * L₁₊bc(params.x̄, bc)
    Lₓ = L₁ - params.σ^2 / 2 * L₂bc(params.x̄, bc)
    L_bc = I * params.ρ - Lₓ
    v = L_bc \ π.(params.x)
    return (v = v, L_bc = L_bc,  Lₓ = Lₓ, L₁ = L₁, L₂bc = L₂bc(params.x̄, bc))
end
# DiffEqOperators Setup
function DEO_state_dependent_drift(π, μ, params)
    dx = params.x[2] - params.x[1]
    # discretize L = ρ - μ D_x - σ^2 / 2 D_xx
    # subject to reflecting barriers at 0 and 1
    drift = μ.(params.x)
    L1 = UpwindDifference(1,1,dx,params.M,drift)
    L2 = CenteredDifference(2,2,dx,params.M)
    Q = Neumann0BC(dx, 1)
    # Here Array(A::GhostDerivativeOperator) will return a tuple of the linear part
    # and the affine part of the operator A, hence we index Array(µ*L1*Q).
    # The operators in this example are purely linear, so we don't worry about Array(µ*L1*Q)[2]
    L₁ = Array(L1*Q)[1]
    Lₓ = L₁ - Array(params.σ^2/2*L2*Q)[1]
    L_bc = I * params.ρ - Lₓ

    # solve the value function
    v = L_bc \ π.(params.x);
    return (v = v, L_bc = L_bc, Lₓ = Lₓ, L₁ = L₁, L₂bc = Array(L2*Q)[1] )
end

@testset "State Dependent Drifts" begin
    @test DEO_state_dependent_drift(π, μ, p).v ≈ SDO_state_dependent_drift(π, μ, p).v
    @test DEO_state_dependent_drift(π, μ, p).L_bc ≈ SDO_state_dependent_drift(π, μ, p).L_bc
    @test DEO_state_dependent_drift(π, μ, p).Lₓ ≈ SDO_state_dependent_drift(π, μ, p).Lₓ
    @test DEO_state_dependent_drift(π, μ, p).L₁ ≈ SDO_state_dependent_drift(π, μ, p).L₁
    @test DEO_state_dependent_drift(π, μ, p).L₂bc ≈ SDO_state_dependent_drift(π, μ, p).L₂bc
end

#----------------------------------
#Testing For Absorbing BC
#----------------------------------

# payoff function
π(x) = x^2
# SimpleDifferentialOperators setup
function SDO_absorbing_bc(π, params)
    bc = (NonhomogeneousAbsorbing(params.S), Reflecting())
    Lₓbc = params.μ*L₁₋bc(params.x̄, bc) + params.σ^2 / 2 * L₂bc(params.x̄ , bc)
    L_bc = I * params.ρ - Lₓbc

    # construct the RHS with affine boundary
    π_star = π.(params.x) + L₀affine(params.x̄, π.(params.x), bc)

    # solve the value function
    v = L_bc \ π_star

    return (v = v, L_bc = L_bc, Lₓbc = Lₓbc, L₁₋bc = L₁₋bc(params.x̄, bc), L₂bc = L₂bc(params.x̄, bc),
        π_star = π_star, π = π.(params.x))
end

# DiffEqOperators Setup
function DEO_absorbing_bc(π, params)
    dx = params.x[2] - params.x[1]

    L1 = UpwindDifference(1,1,dx,params.M,params.μ)
    L2 = CenteredDifference(2,2,dx,params.M)
    # RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::T, order = 1)
    # The variables in l are [αl, βl, γl], and correspond to a BC of the form αl*u(0) + βl*u'(0) = γl
    # imposed on the lower index boundary. The variables in r are [αr, βr, γr],
    # and correspond to an analagous boundary on the higher index end.
    l = (1.0, 0.0, p.S)
    r = (0.0, 1.0, 0.0)
    Q = RobinBC(l, r, dx, 1)

    Lₓbc = Array(L1*Q)[1] + Array(params.σ^2/2*L2*Q)[1]
    L_bc = I * params.ρ - Lₓbc

    # solve the value function
    v = L_bc \ π.(params.x)
    return (v = v, L_bc = L_bc, Lₓbc = Lₓbc, L₁₋bc = Array(L1*Q)[1] ./ params.μ, L₂bc = Array(L2*Q)[1],
        π_star = π.(params.x), π = π.(params.x))
end

p = params(x̄ = range(0.0, 1.0, length = (p.M+2)))

@testset "Absorbing BC" begin
    @test DEO_absorbing_bc(π, p).v ≈ SDO_absorbing_bc(π, p).v
    @test DEO_absorbing_bc(π, p).L_bc ≈ SDO_absorbing_bc(π, p).L_bc
    @test DEO_absorbing_bc(π, p).Lₓbc ≈ SDO_absorbing_bc(π, p).Lₓbc
    @test DEO_absorbing_bc(π, p).L₁₋bc ≈ SDO_absorbing_bc(π, p).L₁₋bc
    @test DEO_absorbing_bc(π, p).L₂bc ≈ SDO_absorbing_bc(π, p).L₂bc
    @test DEO_absorbing_bc(π, p).π_star ≈ SDO_absorbing_bc(π, p).π_star
    @test DEO_absorbing_bc(π, p).π ≈ SDO_absorbing_bc(π, p).π
end

#----------------------------------
#Testing For Solving KFE
#----------------------------------

function SDO_Solve_KFE(params)
    # ξ values for mixed boundary conditions
    ξ_lb = ξ_ub = -2 * params.μ / params.σ^2
    # define the corresponding mixed boundary conditions
    # note that the direction on the lower bound is backward (default is forward)
    # as the drift μ is negative.
    bc = (Mixed(ξ = ξ_lb, direction = :backward), Mixed(ξ = ξ_ub))

    # use SimpleDifferentialOperators.jl to construct the operator on the interior
    L_KFE_with_drift = Array(-params.μ*L₁₊bc(params.x̄, bc) + params.σ^2 / 2 * L₂bc(params.x̄, bc))
    L_KFE_without = params.σ^2 / 2 * L₂bc(params.x̄, bc)
    return (L_KFE_with_drift = L_KFE_with_drift, L_KFE_without = L_KFE_without)
end

function SDO_Solve_KFE_forward(params)
    # ξ values for mixed boundary conditions
    ξ_lb = ξ_ub = -2 * params.μ / params.σ^2
    # define the corresponding mixed boundary conditions
    # but with forward finite difference derivatives at boundaries
    bc = (Mixed(ξ = ξ_lb, direction = :forward), Mixed(ξ = ξ_ub))

    # use SimpleDifferentialOperators.jl to construct the operator on the interior
    L_KFE_with_drift = Array(-params.μ*L₁₊bc(params.x̄, bc) + params.σ^2 / 2 * L₂bc(params.x̄, bc))
    L_KFE_without = params.σ^2 / 2 * L₂bc(params.x̄, bc)
    return (L_KFE_with_drift = L_KFE_with_drift, L_KFE_without = L_KFE_without)
end

function DEO_Solve_KFE(params)
    dx = params.x[2] - params.x[1]

    L1 = UpwindDifference(1,1,dx,params.M,-params.μ)
    # L2l = UpwindDifference(2,2,dx,params.M,
    #                          vcat(-1.,zeros(params.M-1)))
    # L2r = UpwindDifference(2,2,dx,params.M,
    #                          vcat(zeros(params.M-1),1.))
    L2 = CenteredDifference(2,2,dx,params.M)

    ξ_lb = ξ_ub = -2 * params.μ / params.σ^2

    l = (ξ_lb, 1., 0.)
    r = (ξ_ub, 1., 0.)
    Q = RobinBC(l, r, dx, 1)

    # use SimpleDifferentialOperators.jl to construct the operator on the interior
    # Only difference is handling of L2, we may need to upwind for boundary
    L_KFE_with_drift = Array(L1*Q)[1] + Array(params.σ^2/2 *L2*Q)[1]
    L_KFE_without = (params.σ^2 / 2) .* Array(L2 * Q)[1]
    # L_KFE_without = (params.σ^2 / 2) .* ([0., 1., 0.] .* Array(L2 * Q)[1] +
    #     -(Array(L2l * Q)) + Array(L2r * Q))

    return (L_KFE_with_drift = L_KFE_with_drift, L_KFE_without = L_KFE_without)
end

p = params(x̄ = range(0.0, 1.0, length = (p.M+2)))
@testset "Solving KFE" begin
    @test_broken DEO_Solve_KFE(p).L_KFE_with_drift ≈ SDO_Solve_KFE(p).L_KFE_with_drift # concretization of Robin conditions appears broken
    @test_broken DEO_Solve_KFE(p).L_KFE_without ≈ SDO_Solve_KFE(p).L_KFE_without # Hard to check what's wrong since KFE uses forward/backward at boundaries rather than central differences when w/mixed bcs
    @test DEO_Solve_KFE(p).L_KFE_with_drift ≈ SDO_Solve_KFE_forward(p).L_KFE_with_drift
    @test DEO_Solve_KFE(p).L_KFE_without ≈ SDO_Solve_KFE_forward(p).L_KFE_without
end

nothing


text:
using Test

# write your own tests here
@test 1 == 1

using DataFrames
using ExoplanetsSysSim

function run_constructor_tests()
  ExoplanetsSysSim.SimulationParameters.test_sim_param_constructors()
  sim_param = ExoplanetsSysSim.setup_sim_param_demo()
  ExoplanetsSysSim.test_orbit_constructors()
  ExoplanetsSysSim.test_planet_constructors(sim_param)
  ExoplanetsSysSim.test_star_constructors(sim_param)
  ExoplanetsSysSim.test_planetary_system_constructors(sim_param)

  ExoplanetsSysSim.test_target(sim_param)
  ExoplanetsSysSim.test_transit_observations(sim_param)
  (cat_phys, cat_obs) = ExoplanetsSysSim.test_catalog_constructors(sim_param)
  ExoplanetsSysSim.test_summary_statistics(cat_obs, cat_phys, sim_param)
  ExoplanetsSysSim.test_abc_distance(cat_obs, cat_phys, sim_param)
  return 0
end
@test run_constructor_tests() == 0  # Just tests that the basic elements compile and run  # TODO: Write tests that will be useful in diagnosing any bugs

#Test CORBITS moved from ExoplanetsSysSim.test_corbits() to
using CORBITS
include(joinpath(dirname(pathof(CORBITS)),"..","test","runtests.jl"))


text:
# Accelerated row matrix product

C = zeros(m,p)

for j=1:m
    C[j,:] = BLAS.gemv(’T’,B,A[j,:])
end

max_stars_repo_path: Code/Numerical Analysis & Machine Learning/Numerical Linear Algebra/RowMatrixProduct.jl

max_stars_repo_name: BambooFlower/Math-Scripts

max_issues_repo_name: BambooFlower/Math-Scripts

max_forks_repo_name: BambooFlower/Math-Scripts


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

JLLWrappers.@generate_main_file_header("FastTransforms")
JLLWrappers.@generate_main_file("FastTransforms", UUID("34b6f7d7-08f9-5794-9e10-3819e4c7e49a"))
end  # module FastTransforms_jll


text:
using LinearAlgebra

function gen_combs(options)

    na = options["na"]
    nb = options["nb"]
    ne = options["ne"]
    nd = options["nd"]
    nk = nb+na+ne+1

    # Repeat powers
    comb = reshape(collect(0:nd), (1,nd+1))

    # Start combinations array
    combs = reshape(collect(0:nd), (1,nd+1))
    for ii = 2:nb+na+1

        # Current width
        width = size(combs,2)

        # Increment combinations array
        combs = [repeat(combs,1,nd+1); kron(comb,ones(1,width))]

        # remove combinations which have degree higher than nd
        ndComb = sum(combs,dims=1)
        combs = combs[:, vec(ndComb .<= nd)]
    end

    if options["noiseCrossTerms"]
        for ii = nb+na+2:nk

            # Current width
            width = size(combs,2)

            # Add noise cross terms
            combs = [repeat(combs,1,nd+1); kron(comb,ones(1,width))]

            # remove combinations which have degree higher than nd
            ndComb = sum(combs, dims=1)
            combs = combs[:, vec(ndComb .<= nd)]
        end
    else
        for ii = nb+na+2:nk
    #         noisecomb = [zeros(ii-1,1); 1]; % only linear terms
            noisecomb = [zeros(ii-1, nd); reshape(collect(1:nd), 1,nd)]
            combs = [[combs; zeros(1,size(combs,2))] noisecomb]
        end
    end

    if !options["crossTerms"]
        combs = combs[:, vec(sum(combs,dims=1) .> maximum(combs,dims=1))]
    end

    if !options["dc"]
        combs = combs[:,2:end]
    end

    return combs
end

function tmean(x::AbstractArray, dims::Int64; tr::Real=0.2)
    """`tmean(x; tr=0.2)`
    Trimmed mean of real-valued array `x`.
    Find the mean of `x`, omitting the lowest and highest `tr` fraction of the data.
    This requires `0 <= tr <= 0.5`. The amount of trimming defaults to `tr=0.2`.
    """

    nrows,ncols = size(x)

    if dims==1
        tmeans = zeros(1,ncols)
        for n = 1:nrows
            tmeans[1,n] = tmean!(copy(x[:,n]), tr=tr)
        end
    elseif dims==2
        tmeans = zeros(nrows,1)
        for n = 1:ncols
            tmeans[n,1] = tmean!(copy(x[n,:]), tr=tr)
        end
    end
    return tmeans
end

function tmean(x::AbstractArray; tr::Real=0.2)
    """`tmean(x; tr=0.2)`
    Trimmed mean of real-valued array `x`.
    Find the mean of `x`, omitting the lowest and highest `tr` fraction of the data.
    This requires `0 <= tr <= 0.5`. The amount of trimming defaults to `tr=0.2`.
    """
    if size(x,1) == 0
        return NaN
    else
        tmean!(copy(x), tr=tr)
    end
end

function tmean!(x::AbstractArray; tr::Real=0.2)
    """`tmean!(x; tr=0.2)`
    Trimmed mean of real-valued array `x`, which sorts the vector `x` in place.
    Find the mean of `x`, omitting the lowest and highest `tr` fraction of the data.
    This requires `0 <= tr <= 0.5`. The trimming fraction defaults to `tr=0.2`.
    """
    if tr < 0 || tr > 0.5
        error("tr cannot be smaller than 0 or larger than 0.5")
    elseif tr == 0
        return mean(x)
    elseif tr == .5
        return median!(x)
    else
        n   = length(x)
        lo  = floor(Int64, n*tr)+1
        hi  = n+1-lo
        return mean(sort!(x)[lo:hi])
    end
end

function trimse(x::AbstractArray; tr::Real=0.2)
    """`trimse(x; tr=0.2)`
    Estimated standard error of the mean for Winsorized real-valued array `x`.
    See `winval` for what Winsorizing (clipping) signifies.
    """
    if size(x,1) == 0
        return NaN
    else
        return sqrt(winvar(x,tr=tr))/((1-2tr)*sqrt(length(x)))
    end
end

function winval(x::AbstractArray; tr::Real=0.2)
    """`winval(x; tr=0.2)`
    Winsorize real-valued array `x`.
    Return a copy of `x` in which extreme values (that is, the lowest and highest
    fraction `tr` of the data) are replaced by the lowest or highest non-extreme
    value, as appropriate. The trimming fraction defaults to `tr=0.2`.
    """
    n = length(x)
    if n == 0
        return NaN
    else
        xcopy   = sort(x)
        ibot    = floor(Int64, tr*n)+1
        itop    = n-ibot+1
        xbot, xtop = xcopy[ibot], xcopy[itop]
        return  [x[i]<=xbot ? xbot : (x[i]>=xtop ? xtop : x[i]) for i=1:n]
    end
end

function winmean(x::AbstractArray; tr=0.2)
    """`winmean(x; tr=0.2)`
    Winsorized mean of real-valued array `x`.
    See `winval` for what Winsorizing (clipping) signifies.
    """
    return mean(winval(x, tr=tr))
end

function winvar(x::AbstractArray; tr=0.2)
    """`winvar(x; tr=0.2)`
    Winsorized variance of real-valued array `x`.
    See `winval` for what Winsorizing (clipping) signifies.
    """
    return var(winval(x, tr=tr))
end

function winstd(x::AbstractArray; tr=0.2)
    """`winstd(x; tr=0.2)`
    Winsorized standard deviation of real-valued array `x`.
    See `winval` for what Winsorizing (clipping) signifies.
    """
    return std(winval(x, tr=tr))
end

function wincov(x::AbstractArray, y::AbstractArray; tr::Real=0.2)
    """`wincov(x, y; tr=0.2)`
    Compute the Winsorized covariance between `x` and `y`.
    See `winval` for what Winsorizing (clipping) signifies.
    """
    xvec = winval(x, tr=tr)
    yvec = winval(y, tr=tr)
    return cov(xvec, yvec)
end


text:
using CharDisplay

@info "Initialization..."

d = DisplayP4(
    2, # RS
    3, # RW
    4, # E
    5, # D7
    6, # D6
    7, # D5
    8; # D4
    max_rate = 25000
)

set_cursor(d, 0, 0)
write(d, "Julia +")
set_cursor(d, 0, 1)
write(d, "Raspberry Pi")

display_shift(d)
sleep(1)
display_shift(d)
sleep(1)
display_shift(d)
sleep(1)

@info "STOP!"


text:
# ## Пример №1, подробный, статическая сетка
#
# Сопряженная задача определяется для сеточного решения некоторой прямой задачи,
# и некоторых априорных данных. Поэтому, сначала необходимо сгенерировать
# априорные данные и некоторую матрицу `u`.
# Сделаем это на увеличенном числе интервалов, т.к. решение сопряженной задачи
# менее гладкое.
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: phidetermination, Φ;
using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;

u_l(t) = -8
u_r(t) =  4
q(x) = 4*sin(3 * π * x);        # Коэффициент линейного усиления, который в обратной
                                # задаче необходимо определить, но при генерации априорной
                                # информации мы задаем некоторый коэффициент, который,
                                # собственно, после имея априорную информацию и будем определять.
ε = 0.2;                        # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0, 0.28;                # Область по T
N, M = 50, 80;                  # Увеличенное Кол-во разбиений по X, T
h = (b-a)/N;                    # шаг по X
τ = (T-t₀)/M;                   # шаг по T
Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
qₙ =      q.(Xₙ);               # Сеточные значения коэффициента линейного усиления
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
y₀ = u_init.(Xₙ);               # Начальные условия
#
u, XX, TP = solve(y₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
nothing #hide

## Генерация априорной информации
ϕl      = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                  # Левый вырожденный корень
ϕr      = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true); # Нужно подать инвертированную сетку
ϕr      = reverse(ϕr, dims=1);                                      # А после — инвертировать решение по X
ϕ       = Φ(ϕl, ϕr, N, M);                                          # Серединный корень
f1_data = f1(ϕ, u, Xₙ, N, M);                                       # Положение переходного слоя
f2_data = f2(f1_data, u, Xₙ, N, M);                                 # Значение функции на переходном слое
nothing #hide

## Решение сопряженной задачи
Uₙₘ = u;                        # Сохраним старую матрицу
y₀ = [0.0 for i in 1:N+1];      # Нулевые начальные условия
ψl = [0.0 for i in 1:M+1];      # Нулевые ГУ
ψr = [0.0 for i in 1:M+1];      # Нулевые ГУ

ψ = solve_adjoint(y₀, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, Uₙₘ, f1_data, f2_data, w = 0.0005)
nothing #hide


text:
"""
    TermSet <: AbstractSet{AbstractTerm}

Wrapped `Set{AbstractTerm}` that specifies a collection of terms.
Commonly used methods for `Set` work in the same way for `TermSet`.

Compared with `StatsModels.TermOrTerms`,
it does not maintain order of terms
but is more suitable for dynamically constructed terms.
"""
struct TermSet <: AbstractSet{AbstractTerm}
    terms::Set{AbstractTerm}
    TermSet(set::Set{AbstractTerm}) = new(set)
end

"""
    TermSet([itr])
    TermSet(ts::Union{Int, Symbol, AbstractTerm}...)

Construct a [`TermSet`](@ref) from a collection of terms.
Instead of passing an iterable collection,
one may pass the terms as arguments directly.
In the latter case, any `Int` or `Symbol` will be converted to a `Term`.
See also [`termset`](@ref), which is an alias for the constructor.
"""
TermSet() = TermSet(Set{AbstractTerm}())
TermSet(ts) = TermSet(Set{AbstractTerm}(ts))
TermSet(ts::Union{Int, Symbol, AbstractTerm}...) =
    TermSet(Set{AbstractTerm}(t isa AbstractTerm ? t : term(t) for t in ts))

"""
    termset([itr])
    termset(ts::Union{Int, Symbol, AbstractTerm}...)

Construct a [`TermSet`](@ref) from a collection of terms.
Instead of passing an iterable collection,
one may pass the terms as arguments directly.
In the latter case, any `Int` or `Symbol` will be converted to a `Term`.
"""
termset() = TermSet(Set{AbstractTerm}())
termset(ts) = TermSet(Set{AbstractTerm}(ts))
termset(ts::Union{Int, Symbol, AbstractTerm}...) =
    TermSet(Set{AbstractTerm}(t isa AbstractTerm ? t : term(t) for t in ts))

Base.isempty(ts::TermSet) = isempty(ts.terms)
Base.length(ts::TermSet)  = length(ts.terms)
Base.in(x, ts::TermSet) = in(x, ts.terms)
Base.push!(ts::TermSet, x) = push!(ts.terms, x)
Base.pop!(ts::TermSet, x) = pop!(ts.terms, x)
Base.pop!(ts::TermSet, x, default) = pop!(ts.terms, x, default)
Base.delete!(ts::TermSet, x) = delete!(ts.terms, x)
Base.empty!(ts::TermSet) = empty!(ts.terms)

Base.eltype(::Type{TermSet}) = AbstractTerm
Base.iterate(ts::TermSet, i...) = iterate(ts.terms, i...)
Base.emptymutable(::TermSet, ::Type{<:AbstractTerm}) = termset()

const Terms{N} = NTuple{N, AbstractTerm} where N

"""
    eachterm(t)

Return an iterable collection of terms in `t`.
"""
eachterm(@nospecialize(t::AbstractTerm)) = (t,)
eachterm(@nospecialize(t::Terms)) = t

==(@nospecialize(a::Terms), @nospecialize(b::Terms)) =
    length(a)==length(b) && all(t->t in b, a)
==(@nospecialize(a::InteractionTerm), @nospecialize(b::InteractionTerm)) =
    a.terms==b.terms
==(@nospecialize(a::FormulaTerm), @nospecialize(b::FormulaTerm)) =
    a.lhs==b.lhs && a.rhs==b.rhs

"""
    TreatmentTerm{T<:AbstractTreatment} <: AbstractTerm

A term that contains specifications on treatment and parallel trends assumption.
See also [`treat`](@ref).

# Fields
- `sym::Symbol`: the column name of data representing treatment status.
- `tr::T`: a treatment type that specifies the causal parameters of interest.
- `pr::P`: a parallel type that specifies the parallel trends assumption.
"""
struct TreatmentTerm{T<:AbstractTreatment,P<:AbstractParallel} <: AbstractTerm
    sym::Symbol
    tr::T
    pr::P
end

"""
    treat(s::Symbol, t::AbstractTreatment, p::AbstractParallel)

Construct a [`TreatmentTerm`](@ref) with fields set by the arguments.
"""
treat(s::Symbol, t::AbstractTreatment, p::AbstractParallel) = TreatmentTerm(s, t, p)

"""
    treat(s::Term, t::FunctionTerm{FT}, p::FunctionTerm{FP})

Construct a [`TreatmentTerm`](@ref) with fields set by `s.sym`,
`FT.instance(t.args_parsed...)` and `FP.instance(p.args_parsed...)`.
This method is called by [`parse_treat(formula)`](@ref)
when the formula is constructed by `@formula`.
"""
treat(s::Term, t::FunctionTerm{FT}, p::FunctionTerm{FP}) where {FT,FP} =
    treat(s.sym, FT.instance(t.args_parsed...), FP.instance(p.args_parsed...))

"""
    hastreat(t)

Determine whether the term `t` contains an instance of [`TreatmentTerm`](@ref).
"""
hastreat(::TreatmentTerm) = true
hastreat(::FunctionTerm{typeof(treat)}) = true
hastreat(t::InteractionTerm) = any(hastreat(x) for x in t.terms)
hastreat(::AbstractTerm) = false
hastreat(@nospecialize(t::FormulaTerm)) = any(hastreat(x) for x in eachterm(t.rhs))

"""
    parse_treat(formula::FormulaTerm)

Extract terms related to treatment specifications from `formula`.

# Returns
- `TreatmentTerm`: the unique `TreatmentTerm` contained in the `formula`.
- `TermSet`: a set of terms that are interacted with the `TreatmentTerm`.
- `TermSet`: a set of remaining terms in `formula.rhs`.

Error will be raised if either existence or uniqueness of the `TreatmentTerm` is violated.
"""
function parse_treat(@nospecialize(formula::FormulaTerm))
    mettreat = false
    treatterm = nothing
    ints = TermSet()
    xterms = TermSet()
    for term in eachterm(formula.rhs)
        if term isa TreatmentTerm
            if !mettreat
                mettreat = true
                treatterm = term
            else
                throw(ArgumentError("cannot accept more than one TreatmentTerm"))
            end
        elseif term isa FunctionTerm{typeof(treat)}
            if !mettreat
                mettreat = true
                treatterm = treat(term.args_parsed...)
            else
                throw(ArgumentError("cannot accept more than one TreatmentTerm"))
            end
        elseif term isa InteractionTerm
            if hastreat(term)
                !mettreat ||
                    throw(ArgumentError("cannot accept more than one TreatmentTerm"))
                for t in term.terms
                    if t isa TreatmentTerm
                        if !mettreat
                            mettreat = true
                            treatterm = t
                        else
                            throw(ArgumentError("cannot accept more than one TreatmentTerm"))
                        end
                    elseif t isa FunctionTerm{typeof(treat)}
                        if !mettreat
                            mettreat = true
                            treatterm = treat(t.args_parsed...)
                        else
                            throw(ArgumentError("cannot accept more than one TreatmentTerm"))
                        end
                    else
                        push!(ints, t)
                    end
                end
            else
                push!(xterms, term)
            end
        else
            push!(xterms, term)
        end
    end
    mettreat || throw(ArgumentError("no TreatmentTerm is found"))
    return treatterm::TreatmentTerm, ints, xterms
end

isintercept(t::AbstractTerm) = t in (InterceptTerm{true}(), ConstantTerm(1))
isomitsintercept(t::AbstractTerm) =
    t in (InterceptTerm{false}(), ConstantTerm(0), ConstantTerm(-1))

"""
    parse_intercept(ts::TermSet)

Remove any `ConstantTerm` or `InterceptTerm`
and return Boolean values indicating whether terms explictly requiring
including/excluding the intercept exist.

This function is useful for obtaining a unique way of specifying the intercept
before going through the `schema`--`apply_schema` pipeline defined in `StatsModels`.
"""
function parse_intercept!(ts::TermSet)
    hasintercept = false
    hasomitsintercept = false
    for t in ts
        if isintercept(t)
            delete!(ts, t)
            hasintercept = true
        end
        if isomitsintercept(t)
            delete!(ts, t)
            hasomitsintercept = true
        end
    end
    return hasintercept, hasomitsintercept
end

schema(ts::TermSet, data, hints=nothing) =
    Schema(t=>concrete_term(t, VecColumnTable(data), hints) for t in ts)

concrete_term(t::Term, dt::VecColumnTable, hint) =
    concrete_term(t, getproperty(dt, t.sym), hint)
concrete_term(t::Term, dt::VecColumnTable, hints::Dict{Symbol}) =
    concrete_term(t, getproperty(dt, t.sym), get(hints, t.sym, nothing))
concrete_term(::Term, ::VecColumnTable, hint::AbstractTerm) = hint

termvars(ts::TermSet) = mapreduce(termvars, union, ts, init=Symbol[])


text:
include("DataDictionary.jl")
include("FluxDriver.jl")
include("FVA.jl")

using GLPK
using DelimitedFiles

max_stars_repo_name: varnerlab/SEML

max_issues_repo_name: varnerlab/SEML

max_forks_repo_name: varnerlab/SEML


text:
@test_throws DimensionMismatch splitobs((X, rand(149)), at=0.7)

# ## These tests pass on julia 1.6 but fail on higher versions
# @testset "typestability" begin
#     @testset "Int" begin
#         @test typeof(@inferred(splitobs(10, at=0.))) <: NTuple{2}
#         @test eltype(@inferred(splitobs(10, at=0.))) <: UnitRange
#         @test typeof(@inferred(splitobs(10, at=1.))) <: NTuple{2}
#         @test eltype(@inferred(splitobs(10, at=1.))) <: UnitRange
#         @test typeof(@inferred(splitobs(10, at=0.7))) <: NTuple{2}
#         @test eltype(@inferred(splitobs(10, at=0.7))) <: UnitRange
#         @test typeof(@inferred(splitobs(10, at=0.5))) <: NTuple{2}
#         @test typeof(@inferred(splitobs(10, at=(0.5,0.2)))) <: NTuple{3}
#         @test eltype(@inferred(splitobs(10, at=0.5))) <: UnitRange
#         @test eltype(@inferred(splitobs(10, at=(0.5,0.2)))) <: UnitRange
#         @test eltype(@inferred(splitobs(10, at=0.5))) <: UnitRange
#         @test eltype(@inferred(splitobs(10, at=(0.,0.2)))) <: UnitRange
#     end
#     for var in vars
#         @test typeof(@inferred(splitobs(var, at=0.7))) <: NTuple{2}
#         @test eltype(@inferred(splitobs(var, at=0.7))) <: SubArray
#         @test typeof(@inferred(splitobs(var, at=0.5))) <: NTuple{2}
#         @test typeof(@inferred(splitobs(var, at=(0.5,0.2)))) <: NTuple{3}
#         @test eltype(@inferred(splitobs(var, at=0.5))) <: SubArray
#         @test eltype(@inferred(splitobs(var, at=(0.5,0.2)))) <: SubArray
#     end
#     for tup in tuples
#         @test typeof(@inferred(splitobs(tup, at=0.5))) <: NTuple{2}
#         @test typeof(@inferred(splitobs(tup, at=(0.5,0.2)))) <: NTuple{3}
#         @test eltype(@inferred(splitobs(tup, at=0.5))) <: Tuple
#         @test eltype(@inferred(splitobs(tup, at=(0.5,0.2)))) <: Tuple
#     end
# end

@testset "Int" begin
    @test splitobs(10, at=0.7) == (1:7,8:10)
    @test splitobs(10, at=0.5) == (1:5,6:10)
    @test splitobs(10, at=(0.5,0.3)) == (1:5,6:8,9:10)
    @test splitobs(150, at=0.7) == splitobs(150, at=0.7)
    @test splitobs(150, at=0.5) == splitobs(150, at=0.5)
    @test splitobs(150, at=(0.5,0.2)) == splitobs(150, at=(0.5,0.2))
    @test numobs.(splitobs(150, at=0.7)) == (105,45)
    @test numobs.(splitobs(150, at=(.2,.3))) == (30,45,75)
    @test numobs.(splitobs(150, at=(.1,.2,.3))) == (15,30,45,60)
    # tests if all obs are still present and none duplicated
    @test sum(sum.(getobs.(splitobs(150, at=0.7)))) == 11325
    @test sum(sum.(splitobs(150,at=.1))) == 11325
    @test sum(sum.(splitobs(150,at=(.2,.1)))) == 11325
    @test sum(sum.(splitobs(150,at=(.1,.4,.2)))) == 11325
    @test sum.(splitobs(150, at=0.7)) == (5565, 5760)
end

@testset "Array, SparseArray, and SubArray" begin
    for var in (Xs, ys, vars...)
        @test numobs.(splitobs(var, at=10)) == (10,5)
        @test numobs.(splitobs(var, at=0.7)) == (10, 5)
        @test numobs.(splitobs(var, at=(.2,.3))) == (3,4, 8)
        @test numobs.(splitobs(var, at=(.11,.2,.3))) == (2,3,4,6)
        @test numobs.(splitobs(var, at=(2,3,4))) == (2,3,4,6)
    end
    # tests if all obs are still present and none duplicated
    @test sum(sum.(splitobs(X1,at=.1))) == 1200
    @test sum(sum.(splitobs(X1,at=(.2,.1)))) == 1200
    @test sum(sum.(splitobs(X1,at=(.1,.4,.2)))) == 1200
end

@testset "Tuple of Array, SparseArray, and SubArray" begin
    for tup in ((Xs,ys), (X,ys), (Xs,y), (Xs,Xs), (XX,X,ys), (X,yv), (Xv,y), tuples...)
        @test_throws MethodError splitobs(tup..., at=0.5)
        @test all(map(x->(typeof(x)<:Tuple), splitobs(tup,at=0.5)))
        @test numobs.(splitobs(tup, at=0.6)) == (9, 6)
        @test numobs.(splitobs(tup, at=(.2,.3))) == (3,4,8)
    end
end

@testset "shuffle" begin
    s = splitobs(X, at=0.1)[1]
    @test s == X[:,1:2]
    s = splitobs(X, at=0.1, shuffle=true)[1]
    @test size(s) == size(X[:,1:2])
    @test s[:,1] !== s[:,2]
    @test s != X[:,1:2]
    @test any(s[:,1] == x for x in eachcol(X))
    @test any(s[:,2] == x for x in eachcol(X))
end

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase

max_forks_repo_name: JuliaML/MLBase


text:
function get_iodc(obj::GPSOrbitalElements)
    return jcall(obj, "getIODC", jint, ())
end

function get_iode(obj::GPSOrbitalElements)
    return jcall(obj, "getIODE", jint, ())
end

function get_tgd(obj::GPSOrbitalElements)
    return jcall(obj, "getTGD", jdouble, ())
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct RouteFilterRulesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a route in the specified route filter.
Param: resourceGroupName::String (required)
Param: routeFilterName::String (required)
Param: ruleName::String (required)
Param: routeFilterRuleParameters::RouteFilterRule (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: RouteFilterRule
"""
function _swaggerinternal_routeFilterRulesCreateOrUpdate(_api::RouteFilterRulesApi, resourceGroupName::String, routeFilterName::String, ruleName::String, routeFilterRuleParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", RouteFilterRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/routeFilters/{routeFilterName}/routeFilterRules/{ruleName}", ["azure_auth"], routeFilterRuleParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "routeFilterName", routeFilterName)  # type String
    Swagger.set_param(_ctx.path, "ruleName", ruleName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function routeFilterRulesCreateOrUpdate(_api::RouteFilterRulesApi, resourceGroupName::String, routeFilterName::String, ruleName::String, routeFilterRuleParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_routeFilterRulesCreateOrUpdate(_api, resourceGroupName, routeFilterName, ruleName, routeFilterRuleParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function routeFilterRulesCreateOrUpdate(_api::RouteFilterRulesApi, response_stream::Channel, resourceGroupName::String, routeFilterName::String, ruleName::String, routeFilterRuleParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_routeFilterRulesCreateOrUpdate(_api, resourceGroupName, routeFilterName, ruleName, routeFilterRuleParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified rule from a route filter.
Param: resourceGroupName::String (required)
Param: routeFilterName::String (required)
Param: ruleName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_routeFilterRulesDelete(_api::RouteFilterRulesApi, resourceGroupName::String, routeFilterName::String, ruleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/routeFilters/{routeFilterName}/routeFilterRules/{ruleName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "routeFilterName", routeFilterName)  # type String
    Swagger.set_param(_ctx.path, "ruleName", ruleName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function routeFilterRulesDelete(_api::RouteFilterRulesApi, resourceGroupName::String, routeFilterName::String, ruleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_routeFilterRulesDelete(_api, resourceGroupName, routeFilterName, ruleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function routeFilterRulesDelete(_api::RouteFilterRulesApi, response_stream::Channel, resourceGroupName::String, routeFilterName::String, ruleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_routeFilterRulesDelete(_api, resourceGroupName, routeFilterName, ruleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified rule from a route filter.
Param: resourceGroupName::String (required)
Param: routeFilterName::String (required)
Param: ruleName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: RouteFilterRule
"""
function _swaggerinternal_routeFilterRulesGet(_api::RouteFilterRulesApi, resourceGroupName::String, routeFilterName::String, ruleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", RouteFilterRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/routeFilters/{routeFilterName}/routeFilterRules/{ruleName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "routeFilterName", routeFilterName)  # type String
    Swagger.set_param(_ctx.path, "ruleName", ruleName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function routeFilterRulesGet(_api::RouteFilterRulesApi, resourceGroupName::String, routeFilterName::String, ruleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_routeFilterRulesGet(_api, resourceGroupName, routeFilterName, ruleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function routeFilterRulesGet(_api::RouteFilterRulesApi, response_stream::Channel, resourceGroupName::String, routeFilterName::String, ruleName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_routeFilterRulesGet(_api, resourceGroupName, routeFilterName, ruleName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all RouteFilterRules in a route filter.
Param: resourceGroupName::String (required)
Param: routeFilterName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: RouteFilterRuleListResult
"""
function _swaggerinternal_routeFilterRulesListByRouteFilter(_api::RouteFilterRulesApi, resourceGroupName::String, routeFilterName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", RouteFilterRuleListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/routeFilters/{routeFilterName}/routeFilterRules", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "routeFilterName", routeFilterName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function routeFilterRulesListByRouteFilter(_api::RouteFilterRulesApi, resourceGroupName::String, routeFilterName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_routeFilterRulesListByRouteFilter(_api, resourceGroupName, routeFilterName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function routeFilterRulesListByRouteFilter(_api::RouteFilterRulesApi, response_stream::Channel, resourceGroupName::String, routeFilterName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_routeFilterRulesListByRouteFilter(_api, resourceGroupName, routeFilterName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export routeFilterRulesCreateOrUpdate, routeFilterRulesDelete, routeFilterRulesGet, routeFilterRulesListByRouteFilter