text:
export Conv1d

doc"""
    Conv1d(T, ksize, inchannel, outchannel, [padding=0, stride=1, dilation=1])

1-dimensional convolution function.

```julia
T = Float32
x = Var(rand(T,10,5))
f = Conv1d(T, 5, 10, 3, padding=2)
y = f(x)
```
"""
mutable struct Conv1d
    ksize::Int
    padding::Int
    stride::Int
    dilation::Int
    W::Var
    b::Var
end

function Conv1d(::Type{T}, ksize::Int, inchannel::Int, outchannel::Int;
    padding=0, stride=1, dilation=1, init_W=Xavier(), init_b=Fill(0)) where T

    W = init_W(T, ksize*inchannel, outchannel)
    b = init_b(T, outchannel)
    Conv1d(ksize, padding, stride, dilation, parameter(W), parameter(b))
end

function (f::Conv1d)(x, dims)
    h = window1d(x, dims, f.ksize, f.padding, f.stride, f.dilation)
    linear(h, f.W, f.b)
end


text:
using Flux3D
using Documenter
using AbstractPlotting

makedocs(;
    modules = [Flux3D],
    doctest = false,
    authors = "Nirmal P. Suthar <nirmalps@iitk.ac.in>",
    repo = "https://github.com/FluxML/Flux3D.jl/blob/{commit}{path}#L{line}",
    sitename = "Flux3D.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://fluxml.ai/Flux3D.jl",
        assets = String["assets/favicon.ico"],
        analytics = "UA-154580699-2",
    ),
    pages = [
        "Home" => "index.md",
        "Tutorials" => [
            "PointNet classication" => "tutorials/pointnet.md",
            "DGCNN classication" => "tutorials/dgcnn.md",
            "Supervised 3D reconstruction" => "tutorials/fit_mesh.md",
        ],
        "3D Structure" => [
            "PointCloud" => "rep/pointcloud.md",
            "TriMesh" => "rep/trimesh.md",
            "VoxelGrid" => "rep/voxels.md",
        ],
        "Datasets" =>
            ["ModelNet" => "datasets/modelnet.md", "Custom Dataset" => "datasets/utils.md"],
        "Transforms" => "api/transforms.md",
        "Metrics" => "api/metrics.md",
        "API Documentation" => [
            "Conversions" => "api/conversions.md",
            "Helper function" => "api/utils.md",
            "Visualization" => "api/visualize.md",
            "3D Models" => "api/models.md",
        ],
    ],
)

deploydocs(; repo = "github.com/FluxML/Flux3D.jl.git")

ext: jl


text:
"""
運動学
"""
module Kinematics

export calc_local
export calc_global

const c1 = 837019575
const c2 = 4133430
const c3 = 32805
const c4 = 486
const c5 = 18
const c6 = 55801305
const c7 = 688905
const c8 = 3645
const c9 = 81

const c10 = 279006525
const c11 = 1377810
const c12 = 10935
const c13 = 162

const c14 = 243
const c15 = 2066715

const r = 0.0125
const L0 = 0.15

const sq3 = sqrt(3)

"""Aを計算"""
function As(q::T) where T

    A1 = q[1]^2 + q[2]^2 + q[3]^2 - q[1]*q[2] - q[1]*q[3] - q[2]*q[3]
    A2 = 2*q[1] - q[2] - q[3]
    A3 = q[2] - q[3]
    A4 = 3*L0 + q[1] + q[2] + q[3]

    A1, A2, A3, A4
end

function P1(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    -(A2 * A1^4 * A4 * ξ^10) / ((c1 * r^9)) +
    (A2 * A1^3 * A4 * ξ^8) / (c2 * r^7) -
    (A2 * A1^2 * A4 * ξ^6) / (c3 * r^5) +
    (A2 * A1 * A4 * ξ^4) / (c4 * r^3)
    (A2 * A4 * ξ^2) / (c5 * r)
end

function P2(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    -(sq3 * A4 * A3 * A1^4 * ξ^10) / (c1 * r^9) +
    (sq3 * A4 * A3 * A1^3 * ξ^8) / (c2 * r^7) -
    (sq3 * A4 * A3 * A1^2 * ξ^6) / (c3 * r^5) +
    (sq3 * A4 * A1 * A2 * ξ^4) / (c4 * r^3) -
    (sq3 * A4 * A3 * ξ^2) / (c5 * r)
end

function P3(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    (2 * A1^4 * A4 * ξ^9) / (c6 * r^8) -
    (4 * A1^3 * A4 * ξ^7) / (c7 * r^6) +
    (2 * A1^2 * A4 * ξ^5) / (c8 * r^4) -
    (2 * A1 *A4 * ξ^3) / (c9 * r^2) +
    (A4 * ξ) / 3
end

function P(q::T, ξ::U) where {T, U}
    A1, A2, A3, A4 = As(q)

    [
        P1(A1, A2, A3, A4, ξ)
        P2(A1, A2, A3, A4, ξ)
        P3(A1, A2, A3, A4, ξ)
    ]
end

function R11(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    1 - (A2^2 * A1^4 * ξ^10) / (c1 * r^10) +
    (A2^2 * A1^3 * ξ^8) / (c1 * r^8) -
    (A2^2 * A1^2 * ξ^6) / (c3 * r^6) +
    (A1 * A2^2 * ξ^4) / (c4 * r^4) -
    (A2^2 * ξ^2) / (c5 * r^2)
end

function R12(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    (sq3 * A2 * A3 * A1^4 * ξ^10) / (c1 * r^10) +
    (sq3 * A2 * A3 * A1^3 * ξ^8) / (c2 * r^8) -
    (sq3 * A2 * A3 * A1^2 * ξ^6) / (c3 * r^6) +
    (sq3 * A2 * A3 * A1 * ξ^4) / (c4 * r^4) -
    (sq3 * A2 * A3 * ξ^2) / (c5 * r^2)
end

function R13(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    -(2 * A2 * A1^4 * ξ^9) / (c6 * r^9) +
    (4 * A2 * A1^3 * ξ^7) / (c7 * r^7) -
    (2 * A2 * A1^2 * ξ^5) / (c8 * r^5) +
    (2 * A2 * A1 * ξ^3) / (c9 * r^3) -
    (A2 * ξ) / (3 * r)
end

function R21(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    R12(A1, A2, A3, A4, ξ)
end

function R22(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    1 - (A3^2 * A1^4 * ξ^10) / (c10 * r^10) +
    (A3^2 * A1^3 * ξ^8) / (c11 * r^8) -
    (A3^2 * A1^2 * ξ^6) / (c12 * r^6) +
    (A3^2 * A1 * ξ^4) / (c13 * r^4) -
    (A3^2 * ξ^2) / (6 * r^2)
end

function R23(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    -(2*sq3 * A3 * A1^4 * ξ^9) / (c6 * r^9) +
    (4*sq3 * A3 * A1^3 * ξ^7) / (c7 * r^7) -
    (2*sq3 * A3 * A1^2 * ξ^5) / (c8 * r^5) +
    (2*sq3 * A3 * A1 * ξ^3) / (c9 * r^3) -
    (sq3 * A3 * ξ) / (3 * r)
end

function R31(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    -R13(A1, A2, A3, A4, ξ)
end

function R32(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    -R23(A1, A2, A3, A4, ξ)
end

function R33(A1::T, A2::T, A3::T, A4::T, ξ::T) where T
    1 - (2 * ξ^2 * A1) / (9 * r^2) +
    (2 * ξ^4 * A1^2) / (c14 * r^4) -
    (4 * ξ^6 * A1^3) / (c3 * r^6) +
    (2 * ξ^8 * A1^4) / (c15 * r^8) -
    (4 * ξ^10 * A1^5) / (c1 * r^10)
end

function R(q::T, ξ::U) where {T, U}
    A1, A2, A3, A4 = As(q)

    [
        R11(A1, A2, A3, A4, ξ) R12(A1, A2, A3, A4, ξ) R13(A1, A2, A3, A4, ξ)
        R21(A1, A2, A3, A4, ξ) R22(A1, A2, A3, A4, ξ) R23(A1, A2, A3, A4, ξ)
        R31(A1, A2, A3, A4, ξ) R32(A1, A2, A3, A4, ξ) R33(A1, A2, A3, A4, ξ)
    ]
end

function calc_local(N::Int64, Q::T, Ξ::T) where T

    Ps = Vector{Vector{typeof(T)}}(undef, N)
    Rs = Vector{Matrix{typeof(T)}}(undef, N)
    for i in 1:N
        Ps[i] = P(Q[3(i-1)+1:3i], Ξ[i])
        println(R(Q[3(i-1)+1:3i], Ξ[i]))
        Rs[i] = R(Q[3(i-1)+1:3i], Ξ[i])
    end

    Ps, Rs
end

function calc_global(N::Int64, Ps::Vector{T}, Rs::Vector{U}) where {T, U}

    Φs = Vector{T}(undef, N)
    Θs = Vector{U}(undef, N)

    Φs[1] = Ps[1]
    Θs[1] = Rs[1]

    for i in 2:N
        Θs[i] = Θs[i-1] * Rs[i]
        Φs[i] = Φs[i-1] + Θs[i-1]*Ps[i]
    end

    Φs, Θs
end

end


text:
using Revise
using MetaGen
using Gen
using StatProfilerHTML

#possible_objects = ["person","bicycle","car","motorcycle","airplane"]
possible_objects = [1, 2, 3, 4, 5]
#call it
#@profilehtml gt_trace,_ = Gen.generate(metacog, (possible_objects,))
#@profilehtml gt_trace,_ = Gen.generate(metacog, (possible_objects,))
num_frames = 3
num_videos = 2
gt_trace,_ = Gen.generate(metacog, (possible_objects, num_videos, num_frames))
#println(gt_trace)
gt_choices = get_choices(gt_trace)

#################################################################################
num_particles = 100
tr = unfold_particle_filter(num_particles);

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
module InitOptions

export defaultopts, fromcurrentvm, forjavahome, setfromcurrentvm!, unsetfromcurrentvm!,
       setjavahome!, unsetjavahome!, pushclasspath!, pushoptions!, JavaCallInitOptions

#=
Struct to configure the JavaCall package when initializing

fromcurrentvm specifies if the jvm that is currently running should be used
    when executing java calls. If set to true then the other parameters will
    be ignored as the jvm is already running.

javahome specifies a specific javahome to use, if unset than a javahome will be searched
    by examining the JAVA_HOME variable and other common directories for a Java installation

jvmclasspath defines a vector of locations to add to the classpath when creating the jvm

jvmopts defines a vector of options to send the jvm
=#
mutable struct JavaCallInitOptions
    fromcurrentvm::Bool
    javahome::Union{String, Nothing}
    jvmclasspath::Vector{String}
    jvmopts::Vector{String}
end

defaultopts()::JavaCallInitOptions = JavaCallInitOptions(false, nothing, [], [])

fromcurrentvm()::JavaCallInitOptions = JavaCallInitOptions(true, nothing, [], [])

forjavahome(javahome::String) = JavaCallInitOptions(false, javahome, [], [])

function setfromcurrentvm!(initopts::JavaCallInitOptions)
    initopts.fromcurrentvm = true
    return initopts
end

function unsetfromcurrentvm!(initopts::JavaCallInitOptions)
    initopts.fromcurrentvm = false
    return initopts
end

function setjavahome!(initopts::JavaCallInitOptions, javahome::String)
    initopts.javahome = javahome
    return initopts
end

function unsetjavahome!(initopts::JavaCallInitOptions)
    initopts.javahome = nothing
    return initopts
end

function pushclasspath!(initopts::JavaCallInitOptions, entries::String...)
    push!(initopts.jvmclasspath, entries...)
    return initopts
end

function pushoptions!(initopts::JavaCallInitOptions, opts::String...)
    push!(initopts.jvmopts, opts...)
    return initopts
end

function Base.show(io::IO, opts::JavaCallInitOptions)
    print(io, "JavaCallInitOptions("*
        "fromcurrentvm=$(opts.fromcurrentvm), "*
        "javahome=$(opts.javahome), "*
        "jvmclasspath=$(opts.jvmclasspath), "*
        "jvmopts=$(opts.jvmopts))")
end

function Base.copy(opts::JavaCallInitOptions)
    JavaCallInitOptions(
        opts.fromcurrentvm,
        opts.javahome,
        copy(opts.jvmclasspath),
        copy(opts.jvmopts))
end

end


text:
"""
 Tests for https://arxiv.org/abs/1303.5103 with test problems as in the paper.
 DRI1, RI1, RI3, RI5, RI6, RDI1WM, RDI2WM, RDI3WM, RDI4WM
"""

import Statistics # for mean values of trajectories
import LinearAlgebra # for the normn
using StochasticDiffEq
using Test
using Random
using DiffEqDevTools
#using DiffEqGPU

function prob_func(prob, i, repeat)
    remake(prob,seed=seeds[i])
end

"""
 Test Scalar SDEs (oop)
"""

@info "Scalar oop noise"

numtraj = Int(2e5) # in the paper they use 1e9
u₀ = 0.0
f(u,p,t) = 1//2*u+sqrt(u^2+1)
g(u,p,t) = sqrt(u^2+1)
dts = 1 .//2 .^(4:-1:1)
tspan = (0.0,2.0) # 2.0 in paper

h1(z) = z^3-6*z^2+8*z
#analytical_sol(t) = E(f(X(t))) = E(h1(arsinh(X(t))) = t^3-3*t^2+2*t
#analytical_sol(2) = 0 and analytical_sol(1)=0

seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

prob = SDEProblem(f,g,u₀,tspan)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h1(asinh(sol[end])),false),
        prob_func = prob_func
        )

sim = test_convergence(dts,ensemble_prob,DRI1(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("DRI1:", sim.𝒪est[:weak_final])

numtraj = Int(5e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI1(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RI1:", sim.𝒪est[:weak_final])

numtraj = Int(5e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI3(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RI3:", sim.𝒪est[:weak_final])

numtraj = Int(5e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI5(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RI5:", sim.𝒪est[:weak_final])

numtraj = Int(7e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI6(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RI6:", sim.𝒪est[:weak_final])

numtraj = Int(1e3)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
println("RDI1WM:", sim.𝒪est[:weak_final])

numtraj = Int(7e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI2WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RDI2WM:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RDI3WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RDI3WM:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RDI4WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.568 with 3e5 trajectories
println("RDI4WM:", sim.𝒪est[:weak_final])

"""
 Test Scalar SDEs (iip)
"""

@info "Scalar iip noise"

u₀ = [0.0]
f1!(du,u,p,t) = @.(du = 1//2*u+sqrt(u^2 +1))
g1!(du,u,p,t) = @.(du = sqrt(u^2 +1))
dts = 1 .//2 .^(4:-1:1)
tspan = (0.0,2.0)

h1(z) = z^3-6*z^2+8*z

prob = SDEProblem(f1!,g1!,u₀,tspan)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h1(asinh(sol[end][1])),false),
        prob_func = prob_func
        )

numtraj = Int(2e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,DRI1(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("DRI1:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,DRI1NM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("DRI1NM:", sim.𝒪est[:weak_final])

numtraj = Int(5e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI1(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RI1:", sim.𝒪est[:weak_final])

numtraj = Int(5e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI3(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RI3:", sim.𝒪est[:weak_final])

numtraj = Int(5e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI5(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RI5:", sim.𝒪est[:weak_final])

numtraj = Int(7e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI6(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RI6:", sim.𝒪est[:weak_final])

numtraj = Int(1e3)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
println("RDI1WM:", sim.𝒪est[:weak_final])

numtraj = Int(7e5)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI2WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RDI2WM:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RDI3WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RDI3WM:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RDI4WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=0.0
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3
println("RDI4WM:", sim.𝒪est[:weak_final])

"""
 Test Diagonal noise SDEs (iip), SIAM Journal on Numerical Analysis, 47 (2009), pp. 1713–1738
"""

@info "Diagonal noise"

u₀ = [0.1,0.1]
function f3!(du,u,p,t)
  du[1] = 3//2*u[1]
  du[2] = 3//2*u[2]
end
function g3!(du,u,p,t)
  du[1] = 1//10*u[1]
  du[2] = 1//10*u[2]
end
dts = 1 .//2 .^(3:-1:0)
tspan = (0.0,1.0)

h3(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and  == 1//100**exp(301//100*t) if h3(z) = z^2 )

prob = SDEProblem(f3!,g3!,u₀,tspan)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h3(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(5e4)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,DRI1(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.91
println("DRI1:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,DRI1NM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.91
println("DRI1NM:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RI1(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 3.05
println("RI1:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RI3(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.77
println("RI3:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RI5(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RI5:", sim.𝒪est[:weak_final])

numtraj = Int(1e5)
seed = 70
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RI6(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.45 # order is 1.55
println("R6:", sim.𝒪est[:weak_final])

numtraj = Int(1e3)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-1) < 0.45 # order is 1.44
println("RDI1WM:", sim.𝒪est[:weak_final])

numtraj = Int(1e5)
seed = 70
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI2WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.45 # order is 1.55
println("RDI2WM:", sim.𝒪est[:weak_final])

numtraj = Int(5e4)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RDI3WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.84
println("RDI3WM:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RDI4WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.91
println("RDI4WM:", sim.𝒪est[:weak_final])


text:
function sample(f, nsamples::Integer)
    T = typeof(f())
    samples = Vector{T}(undef, nsamples)  # avoid allocation during measurements
    for i in eachindex(samples)
        samples[i] = f()
    end
    return samples
end

function metadata()
    return Dict(
        :nthreads => Threads.nthreads(),
        :julia_version => string(VERSION),
        # ... what else?
    )
end


text:
include("../../../get_home_dir.jl")

using DataFrames
using CSV
#using Plots
using Random
using StatsBase

include(string(PROJECT_HOME,"/assim_code/mironov.jl"));
include(string(PROJECT_HOME,"/assim_code/tau_omega_funs.jl"));
include(string(PROJECT_HOME,"/assim_code/time_averaging.jl"))

df_raw = CSV.read(string(PROJECT_HOME,"/data/moflux_land_data_skipyear_hourly2.csv"), DataFrame);
df_raw[!,"RAIN"] *= 2;

N = 24*365
istart = 24*365*2 + 1;
soil0 = 0.30;

include(string(PROJECT_HOME,"/simulation_code/sim_vary_new_stomata_med.jl"));

deltaT = FT(60*60);
alpha = FT(50);
nsoil = FT(1.2);

function run_sim_2(vcmax_par::FT, k_frac::FT, k_plant::FT,  z_soil::FT, weibB::FT, vol_factor::FT, g1::FT)
	return run_sim_vary(vcmax_par, k_frac, weibB, FT(2), k_plant, FT(1e-4), z_soil, istart, N, soil0, vol_factor, 1e-5, 3, df_raw, g1, deltaT, alpha, nsoil);
end

sim_res1 = run_sim_2(FT(31),FT(0.5), FT(16), FT(2000),FT(5),FT(0.33),FT(506));

#vcmax_par, k_frac, k_plant, k_soil, z_soil, istart, N, smc0, slope_index,use_flux, weibB, weibC

noise_std = 1.3;
noise_var = noise_std^2

tsoil = sim_res1[1].T_SOIL .+ 273.15;
tair = sim_res1[1].T_AIR .+ 273.15;
tcan = (sim_res1[1].LW_OUT / 0.97 / K_STEFAN ) .^ 0.25;
laiM = sim_res1[1].LAI_modis;
canpot_true = mean(sim_res1[3],dims=2);

using Distributions
using LabelledArrays
using LinearAlgebra

vodA = 0.067;
vodB = 0.82;
vodC = 0.051;

trueTB = get_TB_2(sim_res1[2][:,1], tsoil, tcan, canpot_true, laiM, vodA, vodB, vodC);
obsH = trueTB[1] + noise_std*randn(Float64, length(trueTB[1]));
obsV = trueTB[2] + noise_std*randn(Float64, length(trueTB[1]));

obsDF = DataFrame(hpol=obsH, vpol=obsV);
CSV.write("obsTB_witherr_1.csv",obsDF);

obsH2 = trueTB[1] + noise_std*randn(Float64, length(trueTB[1]));
obsV2 = trueTB[2] + noise_std*randn(Float64, length(trueTB[1]));

obsDF = DataFrame(hpol=obsH2, vpol=obsV2);
CSV.write("obsTB_witherr_2.csv",obsDF);

obsH3 = trueTB[1] + noise_std*randn(Float64, length(trueTB[1]));
obsV3 = trueTB[2] + noise_std*randn(Float64, length(trueTB[1]));

obsDF = DataFrame(hpol=obsH3, vpol=obsV3);
CSV.write("obsTB_witherr_3.csv",obsDF);


text:
abstract AbstractSampleMethod
immutable SampleUniform        <: AbstractSampleMethod end
immutable SampleBinCenter      <: AbstractSampleMethod end
immutable SampleUniformZeroBin <: AbstractSampleMethod end

const SAMPLE_UNIFORM          = SampleUniform()
const SAMPLE_BIN_CENTER       = SampleBinCenter()
const SAMPLE_UNIFORM_ZERO_BIN = SampleUniformZeroBin()


text:
⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b)
⊗(a::AbstractArray, b::AbstractArray) = a * b'

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using LaTeXStrings
using ControlSystems
using Plots; gr()

tspan = 0:.01:10

s = tf("s")

plot()
for τ in [.25, .5, .8, 1, 1.35]
    sys = 1/(1+τ*s)
    # PLOT THE STEP RESPONSE
    plot!(step(sys, tspan), xlabel="t (sec)", label=L"\tau"*"=$τ", legend=:bottomright)
end
current()

max_stars_repo_name: a-mhamdi/Eg-Softs

max_issues_repo_name: a-mhamdi/Eg-Softs

max_forks_repo_name: a-mhamdi/Eg-Softs


text:
_isfinite(x) = isfinite(x)
_isfinite(x::VecTypes) = all(isfinite, x)
isfinite_rect(x::Rect) = all(isfinite.(minimum(x))) &&  all(isfinite.(maximum(x)))
scalarmax(x::Union{Tuple, AbstractArray}, y::Union{Tuple, AbstractArray}) = max.(x, y)
scalarmax(x, y) = max(x, y)
scalarmin(x::Union{Tuple, AbstractArray}, y::Union{Tuple, AbstractArray}) = min.(x, y)
scalarmin(x, y) = min(x, y)

extrema_nan(itr::Pair) = (itr[1], itr[2])
extrema_nan(itr::ClosedInterval) = (minimum(itr), maximum(itr))

function extrema_nan(itr)
    vs = iterate(itr)
    vs === nothing && return (NaN, NaN)
    v, s = vs
    vmin = vmax = v
    # find first finite value
    while vs !== nothing && !_isfinite(v)
        v, s = vs
        vmin = vmax = v
        vs = iterate(itr, s)
    end
    while vs !== nothing
        x, s = vs
        vs = iterate(itr, s)
        _isfinite(x) || continue
        vmax = scalarmax(x, vmax)
        vmin = scalarmin(x, vmin)
    end
    return (vmin, vmax)
end

function point_iterator(plot::Union{Scatter, MeshScatter, Lines, LineSegments})
    return plot.positions[]
end

# TODO?
function point_iterator(text::Text{<: Tuple{<: Union{GlyphCollection, AbstractVector{GlyphCollection}}}})
    if is_data_space(text.markerspace[])
        return decompose(Point, boundingbox(text))
    else
        if text.position[] isa VecTypes
            return [to_ndim(Point3f, text.position[], 0.0)]
        else
            return convert_arguments(PointBased(), text.position[])[1]
        end
    end
end

function point_iterator(text::Text)
    return point_iterator(text.plots[1])
end

point_iterator(mesh::GeometryBasics.Mesh) = decompose(Point, mesh)

function point_iterator(list::AbstractVector)
    Iterators.flatten((point_iterator(elem) for elem in list))
end

point_iterator(plot::Combined) = point_iterator(plot.plots)

point_iterator(plot::Mesh) = point_iterator(plot.mesh[])

function br_getindex(vector::AbstractVector, idx::CartesianIndex, dim::Int)
    return vector[Tuple(idx)[dim]]
end

function br_getindex(matrix::AbstractMatrix, idx::CartesianIndex, dim::Int)
    return matrix[idx]
end

function get_point_xyz(linear_indx::Int, indices, X, Y, Z)
    idx = indices[linear_indx]
    x = br_getindex(X, idx, 1)
    y = br_getindex(Y, idx, 2)
    z = Z[linear_indx]
    return Point(x, y, z)
end

function get_point_xyz(linear_indx::Int, indices, X, Y)
    idx = indices[linear_indx]
    x = br_getindex(X, idx, 1)
    y = br_getindex(Y, idx, 2)
    return Point(x, y, 0.0)
end

function point_iterator(plot::Surface)
    X = plot.x[]
    Y = plot.y[]
    Z = plot.z[]
    indices = CartesianIndices(Z)
    return (get_point_xyz(idx, indices, X, Y, Z) for idx in 1:length(Z))
end

function point_iterator(plot::Heatmap)
    X = plot.x[]
    Y = plot.y[]
    Z = plot[3][]
    zsize = size(Z) .+ 1
    indices = CartesianIndices(zsize)
    return (get_point_xyz(idx, indices, X, Y) for idx in 1:prod(zsize))
end

function point_iterator(plot::Image)
    X = plot.x[]
    Y = plot.y[]
    Z = plot[3][]
    zsize = size(Z)
    indices = CartesianIndices(zsize)
    return (get_point_xyz(idx, indices, X, Y) for idx in 1:prod(zsize))
end

function point_iterator(x::Volume)
    axes = (x[1], x[2], x[3])
    extremata = map(extrema∘to_value, axes)
    minpoint = Point3f(first.(extremata)...)
    widths = last.(extremata) .- first.(extremata)
    rect = Rect3f(minpoint, Vec3f(widths))
    return unique(decompose(Point, rect))
end

function foreach_plot(f, s::Scene)
    foreach_plot(f, s.plots)
    foreach(sub-> foreach_plot(f, sub), s.children)
end

foreach_plot(f, s::Figure) = foreach_plot(f, s.scene)
foreach_plot(f, s::FigureAxisPlot) = foreach_plot(f, s.figure)
foreach_plot(f, plot::Combined) = foreach_plot(f, plot.plots)
foreach_plot(f, list::AbstractVector) = foreach(f, list)

function foreach_transformed(f, point_iterator, model, trans_func)
    for point in point_iterator
        point_t = apply_transform(trans_func, point)
        point_m = project(model, point_t)
        f(point_m)
    end
    return
end

function foreach_transformed(f, plot)
    points = point_iterator(plot)
    t = transformation(plot)
    model = model_transform(t)
    trans_func = t.transform_func[]
    # use function barrier since trans_func is Any
    foreach_transformed(f, points, model, identity)
end

function update_boundingbox!(bb_ref, point)
    if all(isfinite, point)
        vec = to_ndim(Vec3f, point, 0.0)
        bb_ref[] = update(bb_ref[], vec)
    end
end

function update_boundingbox!(bb_ref, bb::Rect)
    # ref is uninitialized, so just set it to the first bb
    if !isfinite_rect(bb_ref[])
        bb_ref[] = bb
        return
    end
    # don't update if not finite
    !isfinite_rect(bb) && return
    # ok, update!
    bb_ref[] = union(bb_ref[], bb)
    return
end

function data_limits(plot::AbstractPlot)
    # Because of closure inference problems
    # we need to use a ref here which gets updated inplace
    bb_ref = Base.RefValue(Rect3f())
    foreach_transformed(plot) do point
        update_boundingbox!(bb_ref, point)
    end
    return bb_ref[]
end

function data_limits(scenelike, exclude=(p)-> false)
    bb_ref = Base.RefValue(Rect3f())
    foreach_plot(scenelike) do plot
        if !exclude(plot)
            update_boundingbox!(bb_ref, data_limits(plot))
        end
    end
    return bb_ref[]
end

# A few overloads for performance
function data_limits(plot::Surface)
    mini_maxi = extrema_nan.((plot.x[], plot.y[], plot.z[]))
    mini = first.(mini_maxi)
    maxi = last.(mini_maxi)
    return Rect3f(mini, maxi .- mini)
end

function data_limits(plot::Heatmap)
    mini_maxi = extrema_nan.((plot.x[], plot.y[]))
    mini = Vec3f(first.(mini_maxi)..., 0)
    maxi = Vec3f(last.(mini_maxi)..., 0)
    return Rect3f(mini, maxi .- mini)
end

function data_limits(plot::Image)
    mini_maxi = extrema_nan.((plot.x[], plot.y[]))
    mini = Vec3f(first.(mini_maxi)..., 0)
    maxi = Vec3f(last.(mini_maxi)..., 0)
    return Rect3f(mini, maxi .- mini)
end


text:
module Neovim

using Compat
using MsgPack
import MsgPack: pack, unpack
import Sockets: connect

export NvimClient, nvim_connect, nvim_env, nvim_spawn, nvim_child, start_host
export Buffer, Tabpage, Window
export reply_result, reply_error

# types that have api methods (nvim itself + api defined types)
abstract type NvimObject end

include("client.jl")
include("api_gen.jl")
include("interface.jl")
include("plugin_host.jl")

# too inconvenient api to supply handler here?
function nvim_connect(path::String, args...)
    s = connect(path)
    NvimClient(s, s, args...)
end

nvim_env(args...) = nvim_connect(ENV["NVIM_LISTEN_ADDRESS"], args...)

function nvim_spawn(args...; cmd=`nvim --embed`)
    # TODO(smolck): Maybe just use the Process handle instead of p.in and p.out
    p = open(cmd, "r+")
    (NvimClient(p.in, p.out, args...), p)
end

function nvim_child(args...)
    # make stdio private. Reversed since from nvim's perspective
    input, output = stdout, stdin
    if haskey(ENV, "NEOVIM_JL_DEBUG")
        debug = open(ENV["NEOVIM_JL_DEBUG"], "w")
        redirect_stdout(debug)
        redirect_stderr(debug)
    else
        redirect_stdout()
        redirect_stderr()
    end
    redirect_stdin()

    NvimClient(input, output, args...)
end

end


text:
function uniform_leftorth(A::PeriodicArray{T,1}; tol::Float64 = Defaults.tolgauge, maxiter::Int = Defaults.maxiter,
    cguess = PeriodicArray([isomorphism(Matrix{eltype(A[i])},space(A[i],1),space(A[i],1)) for i in 1:length(A)])) where T <: GenericMPSTensor{S,N1} where {S,N1}
    iteration=1;delta = 2*tol; len = length(A)

    (_,cnew) = leftorth(cguess[1], alg=TensorKit.QRpos())
    normalize!(cnew)

    Al = similar(A)
    Cs = copy(cguess);
    Cs[1] = cnew;

    while iteration<maxiter && delta>tol
        if iteration>10 #when qr starts to fail, start using eigs
            alg=Arnoldi(krylovdim = 30, tol = max(delta*delta,tol/10),maxiter=maxiter)
            #Projection of the current guess onto its largest self consistent eigenvector + isolation of the unitary part

            outp = eigsolve(Cs[1], 1, :LM,alg) do x
                    transfer_left(x,A,Al)
                end

            Cs[1] = leftorth!(outp[2][1],alg=TensorKit.QRpos())[2]
        end

        cold = Cs[1]

        for loc in 1:len
            Al[loc] = permute(Cs[loc]*permute(A[loc],(1,),ntuple(x->x+1,Val{N1}())),ntuple(x->x,Val{N1}()),(N1+1,))
            Al[loc], Cs[loc+1] = leftorth!(Al[loc], alg = QRpos())
            Cs[loc+1]/=norm(Cs[loc+1])
        end

        #update delta
        if domain(cold) == domain(Cs[1]) && codomain(cold) == codomain(Cs[1])
            delta = norm(cold-Cs[1])
        end

        iteration += 1
    end

    delta>tol && @info "leftorth failed to converge $(delta)"

    return Al, Cs, delta
end

function uniform_rightorth(A::PeriodicArray{T,1}; tol::Float64 = Defaults.tolgauge, maxiter::Int = Defaults.maxiter,
    cguess = PeriodicArray([isomorphism(Matrix{eltype(A[i])},space(A[i],1),space(A[i],1)) for i in 1:length(A)])) where T <: GenericMPSTensor{S,N1} where {S,N1}
    iteration=1; delta = 2*tol; len = length(A)

    (cnew,_) = rightorth(cguess[1], alg=TensorKit.LQpos())
    normalize!(cnew)

    Ar = similar(A)
    Cs = copy(cguess);
    Cs[1] = cnew

    while iteration<maxiter && delta>tol
        if iteration>10#when qr starts to fail, start using eigs
            alg=Arnoldi(krylovdim = 30, tol = max(delta*delta,tol/10),maxiter=maxiter)
            #Projection of the current guess onto its largest self consistent eigenvector + isolation of the unitary part

            outp = eigsolve(Cs[1], 1, :LM,alg) do x
                transfer_right(x,A,Ar)
            end
            Cs[1] = rightorth!(outp[2][1],alg=TensorKit.LQpos())[1]
        end

        cold = Cs[1]
        for loc in len:-1:1
            #@tensor Ar[loc][-1;-2 -3]=A[loc][-1,-2, 1]*Cs[loc][1,-3]
            Ar[loc] = A[loc]*Cs[loc+1]

            temp = permute(Ar[loc],(1,),ntuple(x->x+1,Val{N1}()));

            Cs[loc], temp = rightorth!(temp, alg=LQpos())
            Ar[loc] = permute(temp,ntuple(x->x,Val{N1}()),(N1+1,))
            Cs[loc]/=norm(Cs[loc])
        end

        #update counters and delta
        delta = norm(cold-Cs[1])

        #are we done ?
        iteration += 1
  end

  delta>tol && @info "rightorth failed to converge $(delta)"

  return Ar, Cs, delta
end

#=
#https://arxiv.org/abs/1909.06341
function uniform_leftorth(A::PeriodicArray{T,1}; tol::Float64 = Defaults.tolgauge, maxiter::Int = Defaults.maxiter,
    cguess = PeriodicArray([isomorphism(Matrix{eltype(A[i])},space(A[i],1),space(A[i],1)) for i in 1:length(A)])) where T <: GenericMPSTensor{S,N1} where {S,N1}

    AL = similar(A);
    C = similar(cguess);
    temp_C = similar(cguess);
    for (i,(cprev,a,cnext)) in enumerate(zip(cguess,A,circshift(cguess,-1)))
        (AL[i],temp_C[i+1]) = leftorth!(_permute_front(cprev*_permute_tail(a*inv(cnext))),alg=QRpos())
        C[i+1] = temp_C[i+1]*cnext
    end

    delta = 2*tol; iteration=1;
    while iteration<maxiter && delta>tol
        delta = 0;

        for (i,(a,c)) in enumerate(zip(copy(AL),copy(temp_C)))
            (AL[i],temp_C[i+1]) = leftorth!(_permute_front(c*_permute_tail(a)),alg = QRpos())
            C[i+1] = temp_C[i+1]*C[i+1]
            normalize!(C[i+1])

            delta = max(delta,norm(temp_C[i+1]-one(temp_C[i+1])*norm(temp_C[i+1],Inf)))
        end

        iteration+=1;

    end
    delta>tol && @info "leftorth failed to converge $(delta)"

    return AL, C, delta
end

function uniform_rightorth(A::PeriodicArray{T,1}; tol::Float64 = Defaults.tolgauge, maxiter::Int = Defaults.maxiter,
    cguess = PeriodicArray([isomorphism(Matrix{eltype(A[i])},space(A[i],1),space(A[i],1)) for i in 1:length(A)])) where T <: GenericMPSTensor{S,N1} where {S,N1}

    AR = PeriodicArray(_permute_tail.(A));
    C = similar(cguess);
    temp_C = similar(cguess);
    for (i,(cprev,a,cnext)) in reverse(collect(enumerate(zip(cguess,A,circshift(cguess,-1)))))
        (temp_C[i],AR[i]) = rightorth!(inv(cprev)*_permute_tail(a*cnext),alg=LQpos())
        C[i] = cprev*temp_C[i];
    end

    delta = 2*tol; iteration=1;
    while iteration<maxiter && delta>tol
        delta = 0;

        for (i,(a,c)) in enumerate(zip(copy(AR),circshift(temp_C,-1)))
            (temp_C[i],AR[i]) = rightorth!(_permute_tail(_permute_front(a)*c),alg = LQpos())
            C[i] = C[i]*temp_C[i]
            normalize!(C[i])

            delta = max(delta,norm(temp_C[i]-one(temp_C[i])*norm(temp_C[i],Inf))
        end

        iteration+=1;
    end
    delta>tol && @info "rightorth failed to converge $(delta)"

    return PeriodicArray(_permute_front.(AR)), C, delta
end
=#


text:
using Distributed
addprocs()

@everywhere module SFEM

    using Distributed
    using LinearAlgebra
    using DelimitedFiles
    using Pardiso
    using SparseArrays
    using SharedArrays

    include("assembleout.jl")
    include("assemblein.jl")

    export IOsfem
    export Ksfem
    export MKL

    function IOsfem(A, B, C, D, E)

        Cen = readdlm(A)
        Cen = Array{Float64,2}(Cen)

        Ele = readdlm(B)
        Ele = Array{Int64,2}(Ele)

        Face_in = readdlm(C)
        Face_in = Array{Int64,2}(Face_in)

        Face_out = readdlm(D)
        Face_out = Array{Int64,2}(Face_out)

        Node = readdlm(E)
        Node = Array{Float64,2}(Node)

        return Node, Ele, Cen, Face_out, Face_in
    end

    function Ksfem(Node, Ele, Cen, Face_out, Face_in, EE, NU, FF)

        Q = zeros(Float64,size(Node,1)*3,1)
        max = maximum(Node,dims = 1)
        force_point =  Vector{Int64}()

        @inbounds for i = 1:size(Node,1)
            if Node[i,1] == max[1]
                force_point = append!(force_point,i)
            end
        end

        @inbounds for i = force_point
            Q[i*3] = FF/size(force_point,1)
        end

        Node = SharedArray(Node); Ele = SharedArray(Ele); Cen = SharedArray(Cen)
        Face_in = SharedArray(Face_in); Face_out = SharedArray(Face_out)
        IK = SharedArray{Int64,1}(Int(size(Face_in, 1)*225)+Int(size(Face_out, 1)*144))
        JK = SharedArray{Int64,1}(Int(size(Face_in, 1)*225)+Int(size(Face_out, 1)*144))
        VK = SharedArray{Float64,1}(Int(size(Face_in, 1)*225)+Int(size(Face_out, 1)*144))
        IK, JK, VK = Assembleout(IK,JK,VK,Node, Ele, Cen, Face_out, EE, NU, FF)

        con = size(Face_out, 1)*144
        IK, JK, VK = Assemblein(IK,JK,VK,con,Node, Ele, Cen, Face_in, EE, NU, FF)
        KE = sparse(IK, JK, VK)
        min = minimum(Node, dims = 1)
        node_dex = Vector{Int64}()
        @inbounds for i = 1:size(Node, 1)
            if Node[i,1] == min[1]
                append!(node_dex, i)
            end
        end
        @inbounds for i = 1:size(node_dex, 1)
            KE[node_dex[i]*3, node_dex[i]*3] = KE[node_dex[i]*3, node_dex[i]*3]*1E20
            KE[node_dex[i]*3-1, node_dex[i]*3-1] = KE[node_dex[i]*3-1, node_dex[i]*3-1]*1E20
            KE[node_dex[i]*3-2, node_dex[i]*3-2] = KE[node_dex[i]*3-2, node_dex[i]*3-2]*1E20
        end

        return KE, Q

    end

    function MKL(KE,Q)

        ps = MKLPardisoSolver()
        set_nprocs!(ps,48)
        result = similar(Q)
        result = solve(ps, KE, Q)
#=
        set_matrixtype!(ps, 11)
        KE_pardiso = get_matrix(ps, KE, :T)
        set_phase!(ps, 13)
        result = similar(Q)
        printstyled("[STATUS] --> 开始解方程\n", color=:green)
        pardiso(ps, result, KE, Q)
        set_phase!(ps, -1)
=#
        Disp = zeros(Float64,Int(size(Q,1)/3),3)
        @inbounds for i = 1:Int(size(Q,1)/3)
            Disp[i,1] = result[i*3-2,1]
            Disp[i,2] = result[i*3-1,1]
            Disp[i,3] = result[i*3,1]
        end

        return Disp
    end

end

max_stars_repo_path: src/multi-core juSFEM/SFEM.jl

max_stars_repo_name: Elliothuo/juSFEM

max_issues_repo_path: src/multi-core juSFEM/SFEM.jl

max_issues_repo_name: Elliothuo/juSFEM

max_forks_repo_path: src/multi-core juSFEM/SFEM.jl

max_forks_repo_name: Elliothuo/juSFEM


text:
"""
IterManager: convergence based only on number of iterations

@author : Spencer Lyon <spencer.lyon@stern.nyu.edu>
@date : 2015-04-13 11:29:12

"""

struct IterManager <: IterationManager
    maxiter::Int
    verbose::Bool
    print_skip::Int
end

IterManager(maxiter::Int) = IterManager(maxiter, true, div(maxiter, 5))

function Base.show(io::IO, ::MIME"text/plain", im::IterManager)
    verb = im.verbose ? "verbose" : "not verbose"
    println(io, "IterManager: $(im.maxiter) max iterations ($(verb))")
end

finished(mgr::IterManager, istate::IterationState) =
    num_iter(istate) > mgr.maxiter


text:
using Pkg # hideall
Pkg.activate("_literate/D0-categorical/Project.toml")
Pkg.update()

# This tutorial follows loosely [the docs](https://juliadata.github.io/CategoricalArrays.jl/latest/using.html).
#
# ## Defining a categorical vector

using CategoricalArrays

v = categorical(["AA", "BB", "CC", "AA", "BB", "CC"])

# This declares a categorical vector, i.e. a Vector whose entries are expected to represent a group or category.
# You can retrieve the group labels using `levels`:

levels(v)

# which, by  default, returns the labels in lexicographic order.
#
# ## Working with categoricals
#
# ### Ordered categoricals
#
# You can specify that categories are *ordered* by specifying `ordered=true`, the order then follows that of the levels. If you wish to change that order, you  need to  use the `levels!` function.
# Let's see two examples.

v = categorical([1, 2, 3, 1, 2, 3, 1, 2, 3], ordered=true)

levels(v)

# Here the lexicographic order matches what we want so no  need to change it, since we've specified  that the categories are ordered we can do:

v[1] < v[2]

# Let's now consider another example

v = categorical(["high", "med", "low", "high", "med", "low"], ordered=true)

levels(v)

# The levels follow the lexicographic order which  is not what  we want:

v[1] < v[2]

# In order to re-specify the order we need to  use `levels!`:

levels!(v, ["low", "med", "high"])

# now things are properly ordered:

v[1] < v[2]

# ### Missing values
#
# You can also have a categorical vector with missing values:

v = categorical(["AA", "BB", missing, "AA", "BB", "CC"]);

# that doesn't change the levels:

levels(v)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# Add PATH (adapt!)
push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/Leaf/");
push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/Utils/");
#push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/tools/");

# Use Plots:
using Plots
#using PhotoStructs
#plotly()

using CanopyRTMod
using BenchmarkTools
using Statistics
using Parameters
using Revise

@unpack wl,wle,wlf, soil = CanopyRTMod;

leaf = leafbio{FT, length(wl), length(wle), length(wlf)}();
soilAlbedo = 0.07;
incomingLW = 400.0 #W/m2
arrayOfLeaves = Array{leafbio{FT}, 1}(undef, 1)
arrayOfLeaves[1]=leaf

leaf.fqe=0.01
CanopyRTMod.fluspect!(leaf, CanopyRTMod.optis)
CanopyRTMod.computeCanopyGeomProps!(canopy, angles,canOpt)
CanopyRTMod.computeCanopyMatrices!(arrayOfLeaves,canOpt);

using PhysCon
consts = PhysCon.phys{Float32}()

# Hack Thermal optical properties per layer here:
nl = 20
nwl = 1
iLAI    = 0.2;#canopy.LAI/nl;
sigf = canOpt.ddf*leaf.ρ_LW + canOpt.ddb*leaf.τ_LW
sigb = canOpt.ddb*leaf.ρ_LW + canOpt.ddf*leaf.τ_LW
τ_dd = (1 - (1-sigf)*iLAI)*ones(nwl,nl)
ρ_dd = (sigb*iLAI)*ones(nwl,nl)

1.0-τ_dd[1,1]-ρ_dd[1,1]

S⁺ = zeros(nwl,nl).+152*iLAI*0.98
S⁻ = zeros(nwl,nl).+152*iLAI*0.98
S⁺[1,9:11].=200*iLAI
S⁻[1,9:11].=200*iLAI
size(S⁺)

Emin,Eplu,netLW =  CanopyRTMod.RTM_diffuseS(τ_dd, ρ_dd,S⁻, S⁺, [0.0], [300.0], [0.06]);

#iLAI = canopy.LAI/canopy.nlayers;
sumLAI = [0:iLAI:nl*iLAI;];

plot(Emin',-sumLAI, label="Downwelling Thermal")
plot!(Eplu',-sumLAI, label="Upwelling Thermal")
plot!(1 .+Emin2',-sumLAI2, label="Downwelling Thermal fine")
plot!(1 .+Eplu2',-sumLAI2, label="Upwelling Thermal fine")
xlabel!("W/m2")
ylabel!("-Cum LAI")
#plot!(Eplu',1:1:21, label="Upwnwelling Thermal")

(0.98*consts.σ*290^4)

plot(netLW'-S⁺'-S⁻',-sumLAI[1:nl], label="net Thermal")
#plot!(Eplu',-sumLAI, label="Upwelling Thermal")
xlabel!("W/m2")
ylabel!("-Cum LAI")

# Hack Thermal optical properties per layer here:
nl = 200
nwl = 1
iLAI    = 0.02;#canopy.LAI/nl;
sigf = canOpt.ddf*leaf.ρ_LW + canOpt.ddb*leaf.τ_LW
sigb = canOpt.ddb*leaf.ρ_LW + canOpt.ddf*leaf.τ_LW
τ_dd = (1 .- (1-sigf)*iLAI)*ones(nwl,nl)
ρ_dd = (sigb*iLAI)*ones(nwl,nl)

1.0-τ_dd[1,1]-ρ_dd[1,1]
@show 1-exp(-sigf*iLAI)
@show (1-sigf)*iLAI

S⁺ = zeros(nwl,nl).+150*iLAI*0.98
S⁻ = zeros(nwl,nl).+150*iLAI*0.98
#S⁺[1,80:109].=200*iLAI
#S⁻[1,80:109].=200*iLAI
size(S⁺)
Emin2,Eplu2,netLW2 =  CanopyRTMod.RTM_diffuseS(τ_dd, ρ_dd,S⁻, S⁺, [0.0], [300.0], [0.06]);

sumLAI2 = [0:iLAI:nl*iLAI;];
plot(Emin2',-sumLAI2, label="Downwelling Thermal")
plot!(Eplu2',-sumLAI2, label="Upwelling Thermal")
xlabel!("W/m2")
ylabel!("-Cum LAI")

# Hack Thermal optical properties per layer here:
nl = 20
nwl = 1
iLAI    = 4/nl;0.002;#canopy.LAI/nl;
sigf = canOpt.ddf*leaf.ρ_LW + canOpt.ddb*leaf.τ_LW
sigb = canOpt.ddb*leaf.ρ_LW + canOpt.ddf*leaf.τ_LW
τ_dd = (1 - (exp(-sigf*iLAI))*ones(nwl,nl)
ρ_dd = (sigb*iLAI)*ones(nwl,nl)

#
#
#
#
# a parsing error here
# 1.0-τ_dd[1,1]-ρ_dd[1,1]
#
#
#
#
S⁺ = zeros(nwl,nl).+150*iLAI*0.98
S⁻ = zeros(nwl,nl).+150*iLAI*0.98
#S⁺[1,80:109].=200*iLAI
#S⁻[1,80:109].=200*iLAI
size(S⁺)
Emin3,Eplu3,netLW3 =  CanopyRTMod.RTM_diffuseS(τ_dd, ρ_dd,S⁻, S⁺, [0.0], [300.0], [0.06]);
sumLAI3 = [0:iLAI:nl*iLAI;];
plot(Emin2',-sumLAI2, label="Downwelling Thermal")
plot!(Eplu2',-sumLAI2, label="Upwelling Thermal")
xlabel!("W/m2")
ylabel!("-Cum LAI")

plot(Emin3',-sumLAI3, label="Downwelling Thermal")
plot!(Eplu3',-sumLAI3, label="Upwelling Thermal")
plot!(1 .+Emin2',-sumLAI2, label="Downwelling Thermal fine")
plot!(1 .+Eplu2',-sumLAI2, label="Upwelling Thermal fine")
xlabel!("W/m2")
ylabel!("-Cum LAI")

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module RocketLowercaseOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: lowercase()" begin

    println("Testing: operator lowercase()")

    run_proxyshowcheck("Lowercase", lowercase())

    run_testset([
        (
            source = from("Hello, world") |> lowercase(),
            values = @ts(['h', 'e', 'l', 'l', 'o', ',', ' ', 'w', 'o', 'r', 'l', 'd', c])
        ),
        (
            source = completed() |> lowercase(),
            values = @ts(c)
        ),
        (
            source      = faulted(String, "e") |> lowercase(),
            values      = @ts(e("e")),
            source_type = String
        ),
        (
            source = never() |> lowercase(),
            values = @ts()
        )
    ])

end

end


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools
using Optim

# Wavelength range [nm]
λ = 400:1000
# Angle of incidence [degrees]
θ = [5.]
# Polarisation (1.0 = p, 0.0 = s, between 0.0 and 1.0 = average)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Refractive indices of incident (0) and substrate (2)
incident = RIdb.air(beam.λ)
emergent = RIdb.silicon(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(incident),
    ModelFit(:looyenga; N=(ninc=incident, nhost=emergent)),
    ModelFit(:looyenga; N=(ninc=incident, nhost=emergent)),
    LayerTMMO(emergent),
]

# Set the order of the layers (build the system) to fit
order = [
    1 # incident medium
    repeat([2, 3], 9) # top DBR of the MC
    4 # substrate
]

# Absolute reflectance spectrum to fit stored in Utils
Rexp = SpectraDB.bragg_spectrum(beam.λ)
# plot(beam.λ, Rexp)
# gui()

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, show_trace=true, store_trace=true,
);

# Seeds for each ModelFit layer defined above (without alpha)
seed = [
    [69, 0.54], # layers[2]
    [119.0, 0.77], # layers[3]
]

solOptim = fit_tmm_optics(
    Reflectance(Rexp), seed, beam, layers;
    order=order, options=options, alg=LBFGS(),
)

plot(FitSpectrum(), solOptim.beam.λ, solOptim.spectrumExp, solOptim.spectrumFit)
gui()

# Seeds for each ModelFit layer defined above plus the alpha (with alpha)
seed2 = [
    [69, 0.54], # layers[2]
    [119.0, 0.77], # layers[3]
    [0.995], # alpha
]

solOptim2 = fit_tmm_optics(
    Reflectance(Rexp), seed2, beam, layers;
    order=order, options=options, alg=LBFGS(), alpha=true,
)

plot(FitSpectrum(), solOptim2.beam.λ, solOptim2.spectrumExp, solOptim2.spectrumFit)
gui()


text:
module download

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
"Values for units for strings found in option line."
const units = Dict{ String, Float64 }(
  "HZ"    =>      1e0,
  "KHZ"   =>      1e3,
  "MHZ"   =>      1e6,
  "GHZ"   =>      1e9
)
"The default unit, if not found in option line."
const default_units = "GHZ"

"Symbols for different parameter formats for strings found in option line"
const parameters = Dict{ String, Symbol }(
  "S" => :ScatteringParameters,
  "Y" => :AdmittanceParameters,
  "Z" => :ImpedanceParameters,
  "H" => :HybridHParameters,
  "G" => :HybridGParameters,
)
"The default parameter format, if not found in option line."
const default_parameter = "S"

"Symbols for data format for strings found in option line."
const formats = Dict{ String, Symbol }(
  "DB" => :DecibelAngle,
  "MA" => :MagnitudeAngle,
  "RI" => :RealImaginary
)
"The default data format, if not found in option line."
const default_format = "MA"
"The default characteristic impedance, if not found in option line."
const default_resistance = 50

"""
    is_commentline( line )

Checks, if a line is a comment line.
"""
is_comment_line( line ) = length( line ) > 0 && line[ 1 ] == '!'
"""
    parse_comment_line( line )

Returns the comment of a comment line.
"""
parse_comment_line( line ) = line[ 2:end ]

"""
    is_option_line( line )

Checks, if a line is a option line.
"""
is_option_line( line ) = length( line ) > 0 && line[ 1 ] == '#'

"""
    is_frequency_unit_option( option )

Checks, if a string is a valid unit in an option line.
"""
is_frequency_unit_option( option ) = haskey( units,  uppercase( option ) )

"""
    parse_frequency_unit_option( option )

Returns the multiplier for a unit string in an option line.
"""
parse_frequency_unit_option( option ) = units[ uppercase( option ) ]

"""
    is_parameter_option( option )

Checks, if a string is a valid parameter format in an option line.
"""
is_parameter_option( option ) = haskey( parameters, uppercase( option ) )

"""
    parse_parameter_option( option )

Returns the symbol for a valid parameter format string in an option line.
"""
parse_parameter_option( option ) = parameters[ uppercase( option ) ]

"""
    is_format_option( option )

Checks, if a string is a valid data format in an option line.
"""
is_format_option( option ) = haskey( formats, uppercase( option ) )

"""
    parse_format_option( option )

Returns the symbol for a valid data format string in an option line.
"""
parse_format_option( option ) = formats[ uppercase( option ) ]

"""
    is_resistance_option( option )

Checks, if a string is a valid resistance in an option line.
"""
is_resistance_option( option ) = uppercase( option ) == "R"

"""
    parse_resistance_option( option )

Returns the value of the characteristic impedance for a valid resistance string in an option line.
"""
parse_resistance_option( option ) = parse( Float64, option )

"""
    is_empy_line( line )

Checks, if a line is empty.
"""
is_empy_line( line ) = length( line ) == 0

"""
    parse_option_line( line )

Returns the option structure for a valid option line.
"""
function parse_option_line( line )
  line = uppercase( line )

  default = Options()
  unit        = default.unit
  parameter   = default.parameter
  format      = default.format
  resistance  = default.resistance

  entries = split( line[ 2:end ] )

  next_is_resistance = false
  for entry in entries
    if next_is_resistance
      resistance = parse_resistance_option( entry )
      next_is_resistance = false
    else
      if is_frequency_unit_option( entry )
        unit = parse_frequency_unit_option( entry )
      elseif is_parameter_option( entry )
        parameter = parse_parameter_option( entry )
      elseif is_format_option( entry )
        format = parse_format_option( entry )
      else
        next_is_resistance =  is_resistance_option( entry )
      end
    end
  end
  Options( unit, parameter, format, resistance )
end

"""
    ma2comp( m, a )

Converts a pair of values in magnitude / angle format to a complex number.
"""
ma2comp( m, a ) =  m * cis( deg2rad( a ) )

dB2Mag( dB ) = 10 ^ ( dB / 20 )

"""
    da2comp( m, a )

Converts a pair of values in dB / angle format to a complex number.
"""
da2comp( d, a ) = ma2comp( dB2Mag( d ), a )

"Holds conversion functions for parameter format symbols when parsing."
const ParseConversions = Dict{ Symbol, Function }(
  :RealImaginary  => complex,
  :MagnitudeAngle => ma2comp,
  :DecibelAngle   => da2comp,
)

"""
    parse_data( line, N, [ options ] )

Returns a data point structure for a valid line of a one-port Touchstone file.

For 3 or more ports, the line should be the concatenation of 3 or more lines from the Touchstone file.
"""
function parse_data( line::String, N, options = Options(), version = 1 )
  vals = map( s -> parse( Float64, s ), split( line ) )
  parse_data( vals, N, options, version )
end

"""
    parse_data( vals, N, [ options ] )

Returns a data point structure for a vector of numbers from a one-port Touchstone file.

For 3 or more ports, the array should contain exactly 2N² + 1 values.
"""
function parse_data( vals::Array{Float64}, N, options = Options(), version = 1, twoPortDataFlipped = true )
  if length( vals ) != 2N * N + 1
    error( "Wrong number of values!" )
  end
  freq = vals[ 1 ]  * options.unit
  pairs = zip( vals[ 2:2:end - 1 ], vals[ 3:2:end ] )
  converted = map( pair -> ParseConversions[ options.format ]( pair... ), pairs )
  # renormalize Z, Y, G and H parameters for version 1.0
  if version == 1
    if options.parameter == :ImpedanceParameters
      converted = map( x -> x * options.resistance, converted )
    end
    if options.parameter == :AdmittanceParameters
      converted = map( x -> x / options.resistance, converted )
    end
    if N == 2
      if options.parameter == :HybridGParameters
        converted[ 1 ] /= options.resistance
        converted[ 4 ] *= options.resistance
      end
      if options.parameter == :HybridHParameters
        converted[ 1 ] *= options.resistance
        converted[ 4 ] /= options.resistance
      end
    end
  end
  res = reshape( converted, N, N )
  if N > 2 || ( N == 2 && !twoPortDataFlipped )
    res = permutedims( res, [ 2, 1 ] )
  end
  return DataPoint( freq, res )
end

function parseNoiseData( vals::Array{Float64}, options = Options(), version = 1 )
  if length( vals ) != 5
    error( "Wrong number of noise data values!" )
  end
  frequency = vals[ 1 ]  * options.unit
  minNoiseFigure = dB2Mag( vals[ 2 ] )
  reflCoeff = ma2comp( vals[ 3 ], vals[ 4 ] )
  effNoiseRes = vals[ 5 ]
  if version == 1
    effNoiseRes *= options.resistance
  end
  return NoiseDataPoint( frequency, minNoiseFigure, reflCoeff, effNoiseRes )
end

function stripCommentFromLine( line )
  comment = ""
  strings = split( line, "!" )
  if length( strings ) > 1
    comment = join( strings[ 2:end ], "!" )
    line = strings[ 1 ]
  end
  return( line, comment)
end

function getValsAndPushComments!( comments, line )
  line, comment = stripCommentFromLine( line )
  if comment != ""
    push!( comments, comment )
  end
  vals = map( s -> parse( Float64, s ), split( line ) )
  return vals
end

"""
    parse_touchstone_stream( stream, [ ports ] )

Returns a TouchstoneData structure for a valid Touchstone data stream.
"""
function parse_touchstone_stream( stream::IO, ports::Integer = 1 )
  version = 1
  options = Options()
  option_line_found = false
  comments = Vector{ String }()
  data = Vector{ DataPoint }()
  noiseData = Vector{ NoiseDataPoint }()
  neededValues = 1 + 2ports * ports
  vals = Vector{ Float64 }()
  keywordparams = Dict{ Symbol, Any }()
  networkdata = false
  noiseDataExpected = false
  twoPortDataFlipped = true
  referenceNextLine = false
  endfound = false
  lastParamFrequency = -1
  for line in eachline( stream )
    if referenceNextLine
      keywordparams[ :Reference ] = parse_reference( split( line ) )
      referenceNextLine = false
    elseif is_empy_line( line )
      ;# nothing
    elseif is_keyword_line( line )
      version = 2
      keyword, params = parse_keyword_line( line )
      keywordparams[ keyword ] = params
      if option_line_found && keyword == :Version
        error( "V2.0: [Version] 2.0 before option line expected." )
      end
      if keyword in unOrderedKeywordSymbols
        keywordString = keywordSymbolStrings[ keyword ]
        if !haskey( keywordparams, :Version ) || keywordparams[ :Version ] != 2
          error( "V2.0: [Version] keyword before [$( keywordString )] keyword expected." )
        end
        if !option_line_found
          error( "V2.0: Option line before [$( keywordString )] keyword expected." )
        end
        if !haskey( keywordparams, :NumberOfPorts )
          error( "V2.0: [Number of Ports] keyword before [$( keywordString )] keyword expected." )
        end
        if haskey( keywordparams, :NetworkData )
          error( "V2.0: [Network Data] keyword after [$( keywordString )] keyword expected." )
        end
      end
      if keyword == :Reference
        if length( keywordparams[ :Reference ] ) == 0
          referenceNextLine = true
        end
      elseif keyword == :NetworkData
        networkdata = true
        if haskey( keywordparams, :NoiseData )
          error( "V2.0: [Noise Data] keyword after [Network Data] keyword expected." )
        end
        if !haskey( keywordparams, :NumberOfFrequencies )
          error( "V2.0: [Number of Frequencies] keyword before [Network Data] keyword expected." )
        end
      elseif keyword == :NumberOfPorts
        ports = keywordparams[ :NumberOfPorts ]
        neededValues = 1 + 2ports * ports
      elseif keyword == :TwoPortDataOrder
        dataOrder = keywordparams[ :TwoPortDataOrder ]
        if dataOrder == "12_21"
          twoPortDataFlipped = true
        elseif dataOrder == "21_12"
          twoPortDataFlipped = false
        end
      elseif keyword == :NoiseData
        noiseDataExpected   = true
        networkdata = false
      elseif keyword == :End
        networkdata = false
        noiseDataExpected   = false
        endfound = true
      end
    elseif is_comment_line( line )
      push!( comments, parse_comment_line( line ) )
    elseif !option_line_found && is_option_line( line )
      if haskey( keywordparams, :Version ) && keywordparams[ :Version ] == 2 && haskey( keywordparams, :NumberOfPorts )
        error( "V2.0: Option line before [Number of Ports] keyword expected." )
      end
      options = parse_option_line( line )
      option_line_found = true
    elseif noiseDataExpected
      vals = getValsAndPushComments!( comments, line )
      push!( noiseData, parseNoiseData( vals, options, version ) )
    elseif !haskey( keywordparams, :Version ) || ( keywordparams[ :Version ] == 2 && networkdata )
      newVals = getValsAndPushComments!( comments, line )
      if version == 1
        # V1.0 Noise Data starts, when frequency decreases.
        if length( vals ) == 0
          if newVals[ 1 ] <= lastParamFrequency
            push!( noiseData, parseNoiseData( newVals, options, 1 ) )
            vals = Vector{Float64}()
            noiseDataExpected = true
          end
          lastParamFrequency = newVals[ 1 ]
        end
        # V1.0: not more then 4 value pairs per line. (#7)
        # The first  line contains the frequency.
        maxVals = length( vals ) == 0 ? 9 : 8
        if length( newVals ) > maxVals
          error( "V1.0: more than four complex values in one parameter line." )
        end
      end
      vals = [ vals; newVals ]
      if length( vals ) >= neededValues
        push!( data, parse_data( vals, ports, options, version, twoPortDataFlipped  ) )
        vals = Vector{Float64}()
      end
    elseif endfound
      error( "V2.0: Non empty or comment line found after [End] keyword." )
    end
  end
  ret = TouchstoneData( data, noiseData, options, comments, keywordparams )
  return ret
end

"""
    parse_touchstone_string( in, [ ports ] )

Returns a TouchstoneData structure for a valid Touchstone data string.
"""
parse_touchstone_string( in::String, ports::Integer = 1 ) = parse_touchstone_stream( IOBuffer( in ), ports )

"""
    parse_touchstone_file( filename, [ ports ] )

Returns a TouchstoneData structure for a valid Touchstone data file.

When called without ports, tries to interpret the file extension, such as ".s2p", or ".S4P".
"""
function parse_touchstone_file( filename::String, ports::Integer = 0 )
  if ports == 0
    extorig =  splitext( filename )[ 2 ][ 2: end ]
    ext = uppercase( extorig )
    if ext[ 1 ] == 'S' && ext[ end ] == 'P'
      ports = parse( Int64, ext[ 2:end - 1 ] )
    else
      error( "File extension ( $( extorig ) ) not recognized, unknown number of ports." )
    end
  end
  open( filename ) do io
    parse_touchstone_stream( open( filename ), ports )
  end
end

# V2.0
"Symbols for different keywords"
const keywordStringSymbols = Dict{ String, Symbol }(
  "VERSION"         => :Version,
  "NUMBER OF PORTS" => :NumberOfPorts,
  "TWO-PORT DATA ORDER"         => :TwoPortDataOrder,
  "NUMBER OF FREQUENCIES"       => :NumberOfFrequencies,
  "NUMBER OF NOISE FREQUENCIES" => :NumberOfNoiseFrequencies,
  "REFERENCE"         => :Reference,
  "MATRIX FORMAT"     => :MatrixFormat,
  "MIXED-MODE ORDER"  => :MixedModeOrder,
  "BEGIN INFORMATION" => :BeginInformation,
  "END INFORMATION"   => :EndInformation,
  "NETWORK DATA"      => :NetworkData,
  "NOISE DATA"        => :NoiseData,
  "END"   => :End,
)

const keywordSymbolStrings = Dict{ Symbol, String }(
  :Version        => "Version",
  :NumberOfPorts  => "Number of Ports",
  :TwoPortDataOrder         => "Two-Port Data Order",
  :NumberOfFrequencies      => "Number of Frequencies",
  :NumberOfNoiseFrequencies => "Number of Noise Frequencies",
  :Reference        => "Reference",
  :MatrixFormat     => "Matrix Format",
  :MixedModeOrder   => "Mixed-Mode Order",
  :BeginInformation => "Begin Information",
  :EndInformation   => "End Information",
  :NetworkData      => "Network Data",
  :NoiseData  => "NoiseData",
  :End        => "End",
)

const unOrderedKeywordSymbols = [
  :TwoPortDataOrder,
  :NumberOfFrequencies,
  :NumberOfNoiseFrequencies,
  :Reference,
  :MatrixFormat,
  :MixedModeOrder,
  :BeginInformation,
  :EndInformation,
]

"""
    is_keyword_line( line )

Checks, if a line is a valid keyword line.
"""
is_keyword_line( line::String ) = length( line ) > 0 && line[ 1 ] == '['

function parse_N_params( params, N::Integer, tp::DataType = String )
  if length( params ) != N
    error( "Exactly $( N ) parameter expected." )
  else
    if tp != String
      params = map( x -> parse( tp, x ), params )
    end
    if N == 1
      params = params[ 1 ]
    end
    return params
  end
end

function parse_N_params( params, tp::DataType )
  N = length( params )
  if tp != String
    params = map( x -> parse( tp, x ), params )
  end
end

function parse_version_params( params )
  versionstring = parse_N_params( params, 1 )
  if versionstring == "2.0"
    return 2
  else
    error( "Unknown version $( versionstring )" )
  end
end

parse_number_of_ports( params ) = parse_N_params( params, 1, Int64 )
parse_two_port_data_order( params )  = parse_N_params( params, 1 )
parse_number_of_frequencies( params )       = parse_N_params( params, 1, Int64 )
parse_number_of_noise_frequencies( params ) = parse_N_params( params, 1, Int64 )
parse_reference( params )         = parse_N_params( params, Float64 )
parse_matrix_format( params )     = []
parse_mixed_mode_order( params )  = []
parse_begin_information( params ) = []
parse_end_information( params )   = []
parse_network_data( params )  = []
parse_noise_data( params )  = []
parse_end( params )  = []

const version_param_parsers = Dict{ Symbol, Function }(
  :Version        => parse_version_params,
  :NumberOfPorts  => parse_number_of_ports,
  :TwoPortDataOrder          => parse_two_port_data_order,
  :NumberOfFrequencies       => parse_number_of_frequencies,
  :NumberOfNoiseFrequencies  => parse_number_of_noise_frequencies,
  :Reference      => parse_reference,
  :MatrixFormat   => parse_matrix_format,
  :MixedModeOrder => parse_mixed_mode_order,
  :BeginInformation => parse_begin_information,
  :EndInformation   => parse_end_information,
  :NetworkData      => parse_network_data,
  :NoiseData  => parse_noise_data,
  :End        => parse_end,
)

"""
    parse_keyword_line( line )

Returns the keyword parameters for a keyword line.
"""
function parse_keyword_line( line )
  keyword, paramstring = split( line[ 2:end ], "]" )
  keyword = uppercase( keyword )
  params = split( paramstring )
  if haskey( keywordStringSymbols, keyword )
    kwSymbol = keywordStringSymbols[ keyword ]
    params = version_param_parsers[ kwSymbol ]( params )
  else
    error( "Not a valid keyword." )
  end
  return kwSymbol, params
end


text:
using BeliefPropagation
using Test
using Random, Statistics
using Erdos
using RRRMC
import OnlineStats

include("test_utils.jl")

@testset "FactorGraph" begin
    include("factor_graph.jl")
end

@testset "Matching" begin
    include("matching.jl")
end

@testset "B-Matching" begin
    include("bmatching.jl")
end

@testset "Ising" begin
    include("ising.jl")
end


text:
# plot/z-list.jl

include("jim.jl")


text:
module Cryptography

export  factors, fast_exponentiation, odd,
        next_prime, relative_prime,
        integer_to_base64_string, base64_string_to_integer,
        integer_to_uint8_array, uint8_array_to_integer

include("factors.jl")
include("fast_exponentiation.jl")
include("next_prime.jl")
include("rsa.jl")

end


text:
using Test, BioCCP

include("BioCCP.jl")

max_issues_repo_name: kirstvh/BioCCP

max_forks_repo_name: kirstvh/BioCCP


text:
module ODS

using ZipFile, LightXML, DataFrames

"""
    ODSFile

A handle for Open Document Format Spreadsheet
"""
mutable struct ODSFile

end

"""
    openODS(filename::AbstractString)

opens an ODS file and parse it as xml
"""
function openODS(filename::AbstractString)
    zipsheet = ZipFile.Reader(filename)
        for file in zipsheet.files
            if file.name == "content.xml"
                xroot = read(file,String) |> parse_string |> root
                xhseets = child nodes(xroot["body"][1])
                return xsheets
            end
        end

end

end


text:
function convert_FloatData(T :: DataType, fd_T0 :: QM_FloatData{T0}) where {T0<:Real}
    return QM_FloatData(SparseMatrixCSC{T, Int}(fd_T0.Q.m, fd_T0.Q.n,
                                                fd_T0.Q.colptr, fd_T0.Q.rowval, Array{T}(fd_T0.Q.nzval)),
                        SparseMatrixCSC{T, Int}(fd_T0.AT.m, fd_T0.AT.n,
                                                fd_T0.AT.colptr, fd_T0.AT.rowval, Array{T}(fd_T0.AT.nzval)),
                        Array{T}(fd_T0.b),
                        Array{T}(fd_T0.c),
                        T(fd_T0.c0),
                        Array{T}(fd_T0.lvar),
                        Array{T}(fd_T0.uvar))
end

function convert_types(T :: DataType, pt :: Point{T_old}, itd :: IterData{T_old}, res :: Residuals{T_old},
                       dda :: DescentDirectionAllocs{T_old}, pad :: PreallocatedData{T_old}, T0 :: DataType) where {T_old<:Real}

   pt = convert(Point{T}, pt)
   res = convert(Residuals{T}, res)
   itd = convert(IterData{T}, itd)
   pad = convertpad(PreallocatedData{T}, pad, T0)
   dda = convert(DescentDirectionAllocs{T}, dda)
   return pt, itd, res, dda, pad
end

function iter_and_update_T!(pt :: Point{T}, itd :: IterData{T}, fd_T :: Abstract_QM_FloatData{T}, id :: QM_IntData, res :: Residuals{T},
                            sc :: StopCrit{Tsc}, dda :: DescentDirectionAllocs{T}, pad :: PreallocatedData{T}, ϵ_T :: Tolerances{T},
                            ϵ :: Tolerances{T0}, cnts :: Counters, max_iter_T :: Int, T_next :: DataType,
                            display :: Bool) where {T<:Real, T0<:Real, Tsc<:Real}
    # iters T
    sc.max_iter = max_iter_T
    iter!(pt, itd, fd_T, id, res, sc, dda, pad, ϵ_T, cnts, T0, display)

    # convert to T_next
    pt, itd, res, dda, pad = convert_types(T_next, pt, itd, res, dda, pad, T0)
    sc.optimal = itd.pdd < ϵ.pdd && res.rbNorm < ϵ.tol_rb && res.rcNorm < ϵ.tol_rc
    sc.small_Δx, sc.small_μ = res.n_Δx < ϵ.Δx, itd.μ < ϵ.μ
    return pt, itd, res, dda, pad
end


text:
"""
Asynchronous client for Kafka brokers.
Workflow:

    # create KafkaClient using one of constructors
    kc = KafkaClient(...)
    # send requests obtaining one-time channels with a response
    channel = metadata(kc, ...)
    # take the result
    resp = take!(channel)

"""
mutable struct KafkaClient
    brokers::Dict{Int, TCPSocket}
    meta::Dict{Symbol, Any}
    last_cor_id::Int64
    inprogress::Dict{Int64, Type}
    results::Dict{Int64, Channel{Any}}
end

function KafkaClient(host::AbstractString, port::Int; resp_loop=true)
    sock = connect(host, port)
    meta = init_metadata(sock)
    brokers = Dict(b.node_id => connect(b.host, b.port) for b in meta[:brokers])
    inprogress = Dict{Int64, Type}()
    results = Dict{Int64, Channel{Any}}()
    kc = KafkaClient(brokers, meta, 0, inprogress, results)
    if resp_loop
        sockets = collect(values(brokers))
        for sock in sockets
            @async begin
                while isopen(sock)
                    Base.start_reading(sock)
                    while bytesavailable(sock) == 0 sleep(0.1) end
                    handle_response(kc, sock)
                end
            end
        end
    end
    return kc
end

function Base.show(io::IO, kc::KafkaClient)
    n_brokers = length(kc.brokers)
    broker_str = n_brokers == 1 ? "1 broker" : "$n_brokers brokers"
    print(io, "KafkaClient($broker_str)")
end

function register_request(kc::KafkaClient, ::Type{T}) where T
    kc.last_cor_id += 1
    id = kc.last_cor_id
    kc.inprogress[id] = T
    kc.results[id] = Channel(1)
    return id, kc.results[id]
end

function find_leader_id(kc::KafkaClient, topic::AbstractString, partition::Integer)
    node_id = kc.meta[:topics][topic][partition][:leader]
    return node_id
end

function find_leader(kc::KafkaClient, topic::AbstractString, partition::Integer)
    node_id = kc.meta[:topics][topic][partition][:leader]
    return kc.brokers[node_id]
end

function random_broker(kc::KafkaClient)
    sockets = collect(values(kc.brokers))
    return sockets[rand(1:length(sockets))]
end

function handle_response(kc::KafkaClient, sock::TCPSocket)
    header = recv_header(sock)
    cor_id = header.correlation_id
    T = kc.inprogress[cor_id]
    resp = readobj(sock, T)
    put!(kc.results[cor_id], resp)
    delete!(kc.inprogress, cor_id)
    delete!(kc.results, cor_id) # if nobody listens to a channel,
                                #  GC should delete it as well
end

"""
Ensure that KafkaClient has open connection to the leader for these
topic and partition. Reconnect if previous socket has been closed
"""
function ensure_leader(kc::KafkaClient, topic::AbstractString, partition::Integer)
    node_id = find_leader_id(kc, topic, partition)
    if !isopen(kc.brokers[node_id])
        idx = find(b -> b.node_id == node_id, kc.meta[:brokers])[1]
        host = kc.meta[:brokers][idx].host
        port = kc.meta[:brokers][idx].port
        kc.brokers[node_id] = connect(host, port)
    end
end


text:
using ImageInTerminal, ImageCore
using TestImages, ImageTransformations, CoordinateTransformations
using Rotations, OffsetArrays
using SparseArrays
using ImageQualityIndexes
using Test, ReferenceTests

function ensurecolor(f, args...)
    old_color = Base.have_color
    try
        Core.eval(Base, :(have_color = true))
        return f(args...)
    finally
        Core.eval(Base, :(have_color = $old_color))
    end
end

# ====================================================================

# define some test images
gray_square = colorview(Gray, N0f8[0. 0.3; 0.7 1])
gray_square_alpha = colorview(GrayA, N0f8[0. 0.3; 0.7 1], N0f8[1 0.7; 0.3 0])
gray_line = colorview(Gray, N0f8[0., 0.3, 0.7, 1])
gray_line_alpha = colorview(GrayA, N0f8[0., 0.3, 0.7, 1], N0f8[1, 0.7, 0.3, 0])
rgb_line = colorview(RGB, range(0, stop=1, length=20), zeroarray, range(1, stop=0, length=20))
rgb_line_4d = repeat(repeat(rgb_line', 1, 1, 1, 1), 1, 1, 2, 2)

camera_man = testimage("camera")
lighthouse = testimage("lighthouse")
toucan = testimage("toucan")
lena = testimage("lena_color_256")

# ====================================================================

tests = [
    "tst_colorant2ansi.jl",
    "tst_encodeimg.jl",
    "tst_imshow.jl",
    "tst_baseshow.jl",
]

if VERSION >= v"1.6"
    # Sixel is an additional test target and requires Julia at least v1.6
    # It's not installed with `Pkg.instantiate` because it's not listed in
    # `Project.toml`, but instead be manually added in
    # `.github/workflows/UnitTest.yml` in CI setup time.
    try
        @info "Install additional dependencies"
        using Pkg
        Pkg.add(PackageSpec(name="Sixel"))
        using Sixel
        push!(tests, "sixel.jl")
        @info "Sixel test: enabled."
    catch e
        @warn "Sixel test: disabled."
        @warn e
    end
end

ImageInTerminal.encoder_backend[1] = :ImageInTerminal # manually disable Sixel
for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
module ITensorsGPU

using CUDA
using CUDA.CUTENSOR
using CUDA.CUBLAS
using CUDA.CUSOLVER
using LinearAlgebra
using Random, Strided
using TimerOutputs
using StaticArrays
using ITensors
using ITensors.NDTensors
using Strided
import CUDA: CuArray, CuMatrix, CuVector, cu
import CUDA.CUTENSOR: cutensorContractionPlan_t, cutensorAlgo_t

import CUDA.Mem: pin
#=
const devs = Ref{Vector{CUDAdrv.CuDevice}}()
const dev_rows = Ref{Int}(0)
const dev_cols = Ref{Int}(0)
function __init__()
  voltas    = filter(dev->occursin("V100", CUDAdrv.name(dev)), collect(CUDAdrv.devices()))
  pascals    = filter(dev->occursin("P100", CUDAdrv.name(dev)), collect(CUDAdrv.devices()))
  devs[] = voltas[1:1]
  #devs[] = pascals[1:2]
  CUBLASMG.cublasMgDeviceSelect(CUBLASMG.mg_handle(), length(devs[]), devs[])
  dev_rows[] = 1
  dev_cols[] = 1
end
=#
import ITensors: randn!, compute_contraction_labels,
                 eigen, tensor,
                 scale!, unioninds, array, matrix, vector,
                 polar, tensors, truncate!, leftlim, rightlim,
                 permute, BroadcastStyle, Indices
import ITensors.NDTensors: similartype, ContractionProperties, contract!!, _contract!!, _contract!, contract!, contract,
                         contraction_output, UniformDiagTensor, CombinerTensor, contraction_output_type,
                         UniformDiag, Diag, DiagTensor, NonuniformDiag, NonuniformDiagTensor, zero_contraction_output,
                         outer!, outer!!, is_trivial_permutation, ind, permutedims!!, Dense, DenseTensor, Combiner,
                         Tensor, data, getperm, compute_contraction_properties!, Atrans, Btrans, Ctrans, _contract_scalar!, _contract_scalar_noperm!
import Base.*, Base.permutedims!
import Base: similar
include("tensor/cudense.jl")
include("tensor/dense.jl")
include("tensor/culinearalgebra.jl")
include("tensor/cutruncate.jl")
include("tensor/cucombiner.jl")
include("tensor/cudiag.jl")
include("cuitensor.jl")
include("mps/cumps.jl")
include("mps/cumpo.jl")

#const ContractionPlans = Dict{String, Tuple{cutensorAlgo_t, cutensorContractionPlan_t}}()
const ContractionPlans = Dict{String, cutensorAlgo_t}()

export cu,
       cpu,
       cuITensor,
       randomCuITensor,
       cuMPS,
       randomCuMPS,
       productCuMPS,
       randomCuMPO,
       cuMPO

end #module


text:
using Documenter, MDTable

makedocs(
    sitename = "MDTable.jl"
)

deploydocs(
    repo = "github.com/TheRoniOne/MDTable.jl.git",
)


text:
function mergeSort!(A)
    return mergeSortStep!(A,1,length(A))
end

function mergeSortStep!(A,p,r)
    if p < r
        q = Int(floor((p+r)/2));
        mergeSortStep!(A, p, q);
        mergeSortStep!(A, q+1, r);
        mergeRez!(A, p, q, r);
    end
    return A
end

function mergeRez!(A, p, q, r)
    L = copy(A[p:q])
    R = copy(A[q+1:r])
    push!(L, typemin(Int64))
    push!(R, typemin(Int64))
    i = 1;
    j = 1;
    for k = p:r
        if R[i] <= L[j]
            A[k] = L[j];
            j = j + 1;
        else
            A[k] = R[i];
            i = i + 1;
        end
    end
end


text:
# test for hasEdge of a QuartetNetwork
# Claudia January 2015
# Also, tests for net.ht, net.numht, qnet.indexht

#println("----- Case G ------")
include("../examples/case_g_example.jl");
# include(joinpath(dirname(pathof(PhyloNetworks)),  "..","examples","case_g_example.jl"))

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);
#parameters!(net)

extractQuartet!(net,d)
error1 = false
try
    q1.qnet.hasEdge == [true,true,true,true] ? nothing : error("q1 wrong hasEdge")
    q2.qnet.hasEdge == [true,false,true,true] ? nothing : error("q2 wrong hasEdge")
    q3.qnet.hasEdge == [true,false,true,true] ? nothing : error("q3 wrong hasEdge")
    q4.qnet.hasEdge == [false,true,true,false] ? nothing : error("q4 wrong hasEdge")
    q5.qnet.hasEdge == [true, true,true,false] ? nothing : error("q5 wrong hasEdge")

    net.ht == [0.1,0.2,0.1,1.0] ? nothing : error("net.ht not correct")
    net.numht == [7,3,6,9] ? nothing : error("net.numth not correct")

    q1.qnet.indexht == [1,2,3,4] ? nothing : error("q1.qnet.indexht not correct")
    q2.qnet.indexht == [1,3,4] ? nothing : error("q2.qnet.indexht not correct")
    q3.qnet.indexht == [1,3,4] ? nothing : error("q3.qnet.indexht not correct")
    q4.qnet.indexht == [2,3] ? nothing : error("q4.qnet.indexht not correct")
    q5.qnet.indexht == [1,2,3] ? nothing : error("q5.qnet.indexht not correct")

    q1.qnet.index == [7,3,6,9] ? nothing : error("q1.qnet.index not correct")
    q2.qnet.index == [5,4,7] ? nothing : error("q2.qnet.index not correct")
    q3.qnet.index == [5,4,7] ? nothing : error("q3.qnet.index not correct")
    q4.qnet.index == [3,4] ? nothing : error("q4.qnet.index not correct")
    q5.qnet.index == [7,3,6] ? nothing : error("q5.qnet.index not correct")
catch
    println("---- error in case G -----")
    global error1 = true
end

#println("----- Case F: bad diamond ------")
include("../examples/case_f_example.jl");
# include(joinpath(dirname(pathof(PhyloNetworks)),  "..","examples","case_f_example.jl"))

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);
parameters!(net)

extractQuartet!(net,d)
error1 = false
try
    q1.qnet.hasEdge == [false,true,true] ? nothing : error("q1 wrong hasEdge")
    q2.qnet.hasEdge == [true,false,false] ? nothing : error("q2 wrong hasEdge")
    q3.qnet.hasEdge == [true,false,true] ? nothing : error("q3 wrong hasEdge")
    q4.qnet.hasEdge == [true,true,false] ? nothing : error("q4 wrong hasEdge")
    q5.qnet.hasEdge == [false, true,true] ? nothing : error("q5 wrong hasEdge")

    all(map(approxEq,net.ht,[0.1,0.7*(1-exp(-0.2)),0.3*(1-exp(-0.1))])) ? nothing : error("net.ht not correct")
    net.numht == [9,21,22] ? nothing : error("net.numth not correct")

    q1.qnet.indexht == [2,3] ? nothing : error("q1.qnet.indexht not correct")
    q2.qnet.indexht == [1] ? nothing : error("q2.qnet.indexht not correct")
    q3.qnet.indexht == [1,3] ? nothing : error("q3.qnet.indexht not correct")
    q4.qnet.indexht == [1,2] ? nothing : error("q4.qnet.indexht not correct")
    q5.qnet.indexht == [2,3] ? nothing : error("q5.qnet.indexht not correct")

    q1.qnet.index == [1,3] ? nothing : error("q1.qnet.index not correct")
    q2.qnet.index == [4] ? nothing : error("q2.qnet.index not correct")
    q3.qnet.index == [5,3] ? nothing : error("q3.qnet.index not correct")
    q4.qnet.index == [5,3] ? nothing : error("q4.qnet.index not correct")
    q5.qnet.index == [1,3] ? nothing : error("q5.qnet.index not correct")
catch
    println("---- error in case F -----")
    global error1 = true
end

#println("----- Case I: bad diamondII ------")
include("../examples/case_i_example.jl");

q1 = Quartet(1,["6","7","4","8"],[0.5,0.4,0.1]);
q2 = Quartet(2,["6","7","10","8"],[0.5,0.4,0.1]);
q3 = Quartet(3,["10","7","4","8"],[0.5,0.4,0.1]);
q4 = Quartet(4,["6","10","4","8"],[0.5,0.4,0.1]);
q5 = Quartet(5,["6","7","4","10"],[0.5,0.4,0.1]);

d = DataCF([q1,q2,q3,q4,q5]);

extractQuartet!(net,d)
error1 = false
try
    q1.qnet.hasEdge == [true,false,true,false,true] ? nothing : error("q1 wrong hasEdge")
    q2.qnet.hasEdge == [true,true,true,true,true] ? nothing : error("q2 wrong hasEdge")
    q3.qnet.hasEdge == [true,true,true,true,true] ? nothing : error("q3 wrong hasEdge")
    q4.qnet.hasEdge == [true,true,true,true,true] ? nothing : error("q4 wrong hasEdge")
    q5.qnet.hasEdge == [true,false,true,false,true] ? nothing : error("q5 wrong hasEdge")

    net.ht == [0.1,2.,1.,1.,1.] ? nothing : error("net.ht not correct")
    net.numht == [9,4,6,9,10] ? nothing : error("net.numth not correct")

    q1.qnet.indexht == [1,3,5] ? nothing : error("q1.qnet.indexht not correct")
    q2.qnet.indexht == [1,2,3,4,5] ? nothing : error("q2.qnet.indexht not correct")
    q3.qnet.indexht == [1,2,3,4,5] ? nothing : error("q3.qnet.indexht not correct")
    q4.qnet.indexht == [1,2,3,4,5] ? nothing : error("q4.qnet.indexht not correct")
    q5.qnet.indexht == [1,3,5] ? nothing : error("q5.qnet.indexht not correct")

    q1.qnet.index == [7,4,8] ? nothing : error("q1.qnet.index not correct")
    q2.qnet.index == [8,4,6,8,9] ? nothing : error("q2.qnet.index not correct")
    q3.qnet.index == [8,4,6,8,9] ? nothing : error("q3.qnet.index not correct")
    q4.qnet.index == [8,4,5,8,9] ? nothing : error("q4.qnet.index not correct")
    q5.qnet.index == [7,4,8] ? nothing : error("q5.qnet.index not correct")

catch
    println("---- error in case I -----")
    global error1 = true
end

if error1
    throw("errors in has edge")
end


text:
# wrappers for LLVM-specific functionality

@inline trap() = ccall("llvm.trap", llvmcall, Cvoid, ())


text:
struct TestComponent <: InfrastructureSystemsType
    name::AbstractString
    val::Int
    _forecasts::Forecasts
    internal::InfrastructureSystemsInternal
end

function TestComponent(name, val)
    return TestComponent(name, val, Forecasts(), InfrastructureSystemsInternal())
end

get_val(component::TestComponent) = component.val

function get__forecasts(component::TestComponent)
    return component._forecasts
end

function JSON2.read(io::IO, ::Type{TestComponent})
    data = JSON2.read(io)
    return TestComponent(data.name, data.val, convert_type(Forecasts, data._forecasts),
                         JSON2.read(JSON2.write(data.internal),
                                    InfrastructureSystemsInternal))
end

function runtests(args...)
    test_prefix = "test_"
    for arg in args
        if !startswith(arg, test_prefix)
            arg = test_prefix * arg
        end
        push!(ARGS, arg)
    end

    try
        include("test/runtests.jl")
    finally
        empty!(ARGS)
    end
end


text:
struct UnknownIndexerError <: Exception end
struct BadSubgraphIpfsHashError <: Exception end


text:
using CauseMap
using PyCall

function getvectors(vec1::AbstractVector, vec2::AbstractVector, libsizemin::Int64, libsizemax::Int64, E::Int64,
    tau_s::Int64, tau_p::Int64, npred::Int64, pred_start::Int64, lib_start::Int64=1)

    shadowmat_dict_vec1, distmat_dict_vec1 = precalc_manif_dists(3:3, 1:1, vec1)
    shadowmat_dict_vec2, distmat_dict_vec2 = precalc_manif_dists(3:3, 1:1, vec2)

    librange12, yval_12 = calcCCM(vec1, vec2,
                                                        shadowmat_dict_vec1, distmat_dict_vec1,
                                                        E, tau_s, tau_p;
                                                        lib_start=lib_start, libsizemin=libsizemin, libsizemax=libsizemax,
                                                        npred=npred, pred_start=pred_start, pred_rand=false)
    librange21, yval_21 = calcCCM(vec2, vec1,
                                                        shadowmat_dict_vec2, distmat_dict_vec2,
                                                        E, tau_s, tau_p;
                                                        lib_start=lib_start, libsizemin=libsizemin, libsizemax=libsizemax,
                                                        npred=npred, pred_start=pred_start, pred_rand=false)
    return yval_12, yval_21
end

function validatevecs(yval_12::AbstractVector, yval_21::AbstractVector; margin=.01,
                                            fname="Para_xmap_Didi_values.txt", col12=2, col21=3)
    col12 = 2
    col21 = 3

    ds_validation = readdlm(fname, '\t', Float64, header=true)[1]
    sum1 = (yval_12 - ds_validation[:, col12])
    sum1 = dot(sum1, sum1)
    sum2 = (yval_21 - ds_validation[:, col21])
    sum2 = dot(sum2, sum2)
    fail1 =  (sum1 > (length(yval_12)*margin))
    fail2 = (sum2 > (length(yval_21)*margin))
    if fail1 | fail2
        if fail1
            println("Para didi xmap off by $sum1")
        end
        if fail2
            println("Didi para xmap off by $sum2")
        end
        error("Error greater than a margin of $margin")
    end
end

E = 3
tau_s = 1
tau_p = 0

libsizemin = 8
libsizemax = 58
pred_start = 13
npred = 60

ds   = readdlm("vr.raw_fixed.txt", '\t', Float64, header=false)
para = ds[:, 2]
didi = ds[:, 3]

para_didi, didi_para = getvectors(para, didi, libsizemin, libsizemax, E,
                                                    tau_s, tau_p, npred, pred_start)
println("Starting Para Didi test")
validatevecs(para_didi, didi_para)
println("Test passed!")


text:
########## Fold ##########

@enum Fold_type f_x=Int('x') f_y=Int('y')

struct Fold{Fold_type}
    line::Int
end
line(fold::Fold) = fold.line

Fold(entry::String) =
    Fold{Fold_type(Int(entry[12]))}(parse(Int, entry[14:end]))

########## Point ##########

mutable struct Point
    x::Int
    y::Int
end
x(point::Point) = point.x
y(point::Point) = point.y

Point(entry::String) =
    Point(parse.(Int, String.(split(entry, ",")))...)

Base.:(==)(point1::Point, point2::Point)::Bool =
    x(point1) == x(point2) && y(point1) == y(point2)

########## functions ##########

mirror(coord::Int, line::Int)::Int =
    coord - line > 0 ? 2 * line - coord : coord

fold!(point::Point, fold::Fold{f_x}) =
    point.x = mirror(x(point), line(fold))
fold!(point::Point, fold::Fold{f_y}) =
    point.y = mirror(y(point), line(fold))

function Base.unique(points::Vector{Point})
    result = Point[]
    for point in points
        !(point in result) && append!(result, [point])
    end
    return result
end

function fold!(points::Vector{Point}, fold::Fold)
    for point in points
        fold!(point, fold)
    end
    points = unique(points)
end

function code(points::Vector{Point})::Vector{String}
    width = maximum(x.(points))
    height = maximum(y.(points))
    return [join(row) for row in eachrow(
        [exist ? '#' : '.' for exist in
            [Point(x, y) in points for x = 0:width, y = 0:height]'
        ]
    )]
end

########## Script ##########

function main(input_name::String)

    # open file and store content into an input string
    input = open(join(["inputs//", input_name]), "r") do io
        read(io, String)
    end

    # split input into a the points part and the folding instructions part
    data = String.(split(input, "\r\n\r\n"))
    # parse data as a vector of points and a vector of fold instructions
    points = Point.(String.(split(data[1], "\r\n")))
    folds = Fold.(String.(split(data[2], "\r\n")))

    # part 1
    fold!(points, folds[1])
    println("Part 1: ", length(points))

    # part 2
    for fold in folds[2:end]
        fold!(points, fold)
    end
    println("Part 2: ")
    for row in eachrow(code(points))
        println(row)
    end

end

main("day13.txt");


text:
function Resistor(; name, R::Float64)
    @named component = VolumelessComponent()
    @unpack Q, ΔP = component

    ps = @parameters (R = R)

    eqs = [
        Q ~ ΔP / R
    ]

    extend(ODESystem(eqs, t, [Q, ΔP], ps; name), component)
end

function Valve(; name, R::Float64)
    @named component = VolumelessComponent()
    @unpack Q, ΔP = component

    ps = @parameters (R = R)

    eqs = [
        Q ~ (ΔP / R) * (ΔP > 0.0)
    ]

    extend(ODESystem(eqs, t, [Q, ΔP], ps; name), component)
end

"""
Vessel model (e.g. artery or vein). This Vessel has an elasticity but no inertia or dampening.

# Arguments
   - `Ees`: Elasticity.
   - `Vd`: Unstressed volume (volume with no pressure).
"""
function ElasticVessel(; name,
    # Parameters
    Ees::Float64,
    Vd::Float64,
    ext_pressure=0.)

    @named compartment = CompliantCompartment(;ext_pressure)
    @unpack V, P_tm = compartment

    ps = @parameters (Ees = Ees, Vd = Vd)

    eqs = [
        P_tm ~ max(0, Ees*(V-Vd))
    ]

    extend(ODESystem(eqs, t, [], ps; name), compartment)
end

"""
Vessel model (e.g. artery or vein). This Vessel has an elasticity and dampening.

# Arguments
   - `Ees`: Elasticity.
   - `Vd`: Unstressed volume (volume with no pressure).
   - `R_damp`: Damping resistance (vessel wall resistance to movement).
"""
function DampedElasticVessel(; name,
    # Parameters
    Ees::Float64,
    Vd::Float64,
    R_damp::Float64 = 0.,
    ext_pressure=0.)

    @named compartment = CompliantCompartment(;ext_pressure)
    @unpack V, P_tm, dV = compartment

    sts = @variables (P_elastic(t), P_damp(t) = 0.)

    ps = @parameters (Ees = Ees, Vd = Vd, R_damp = R_damp)

    eqs = [
        P_elastic ~ max(0, Ees*(V-Vd)),
        dV ~ P_damp / R_damp,
        P_tm ~ P_elastic + P_damp
    ]

    extend(ODESystem(eqs, t, sts, ps; name), compartment)
end

"""
Ventricle model

# Arguments
   - `Ees`: Maximum elasticity (systole).
   - `Vd`: (Unstressed volume) Volume where not contraction can occour.
   - `V0`: Volume with no passive elasticitiy.
   - `λ`: ?. Something diastolic
   - `P0`: ?. Something diastolic Unstressed
"""
function Ventricle(; name,
    # Parameters
    Ees::Float64,
    Vd::Float64,
    V0::Float64,
    λ::Float64,
    P0::Float64 = 0.,
    # Arguments that are passed along to PressurisedCompartment
    ext_pressure = 0.,
    )

    @named compartment = CompliantCompartment(;ext_pressure)
    @unpack V, P_tm = compartment

    ps = @parameters (
        Ees = Ees,
        Vd = Vd,
        V0 = V0,
        λ = λ,
        P0 = P0
    )

    sts = @variables (
        drv(t),
        Pes(t),
        Ped(t)
    )

    eqs = [
        Pes ~ Ees*(V-Vd),
        Ped ~ P0*(exp(λ*(V-V0))-1.0),
        P_tm ~ drv*Pes+Ped
    ]

    extend(ODESystem(eqs, t, [drv], ps; name), compartment)

end

function ConstantPressure(; name,
    # Parameters
    P::Float64)

    @named out = Con()

    ps = @parameters (P = P)

    eqs = [
        out.P ~ P
    ]

    compose(ODESystem(eqs, t, [], ps; name), out)

end

# Register cardiac driver
e_card(t; A = 1.0, B = 80.0, C = 0.27) = A*exp(-B*(t-C)^2)
e_card_rep(t; A = 1.0, B = 80.0, C = 0.27, cycle_len = 1.0) = e_card(t % cycle_len; A = A, B = B, C = C)
@register_symbolic e_card_rep(t, A, B, C, cycle_len)

function automaticDriver(; name,
    A::Float64 = 1.,
    B::Float64 = 80.,
    C::Float64 = 0.27,
    cycle_len = 1.0)

    ps = @parameters (A=A, B=B, C=C)

    sts = @variables contraction(t)

    eqs = [
        contraction ~ e_card_rep(t, A = A, B = B, C = C, cycle_len = cycle_len)
    ]

    ODESystem(eqs, t, [contraction], ps; name)
end

function e_card_qrs(t; A = 1.0, B = 80.0, C = 0.27, qrs_times::Vector{Float64})
    # Sum of the cardiac driver effect over all time_since_qrs.
    # This works since the driver is essetially zero shortly after each QRS.
    t_since_qrs = max.(0, t .- qrs_times)
    return sum(e_card.(t_since_qrs, A = A, B = B, C = C))
end

@register_symbolic e_card_qrs(t, A, B, C, qrs_times)

"""
QRS driver

# Arguments
    - `A`: Amplitude.
    - `B`: Acceleration.
    - `C`: Time (seconds) from activation (QRS) to peak.

"""
function qrsDriver(; name,
    A::Float64 = 1.,
    B::Float64 = 80.,
    C::Float64 = 0.27,
    qrs_times::Vector{Float64})

    ps = @parameters (A=A, B=B, C=C)

    sts = @variables contraction(t)

    eqs = [
        contraction ~ e_card_qrs(t,
                                 A = A, B = B, C = C,
                                 qrs_times = qrs_times)
    ]

    ODESystem(eqs, t, [contraction], ps; name)
end


text:
# @test_throw, with additional testing for the exception message
macro test_throws_message(f, typ, ex...)
    quote
        msg = ""
        @test_throws $(esc(typ)) try
            $(esc(ex...))
        catch err
            msg = sprint(showerror, err)
            rethrow()
        end

        if !$(esc(f))(msg)
            # @test should return its result, but doesn't
            @error "Failed to validate error message\n$msg"
        end
        @test $(esc(f))(msg)
    end
end

# @test_throw, peeking into the load error for testing macro errors
macro test_throws_macro(ty, ex)
    return quote
        Test.@test_throws $(esc(ty)) try
            $(esc(ex))
        catch err
            @test err isa LoadError
            @test err.file === $(string(__source__.file))
            @test err.line === $(__source__.line + 1)
            rethrow(err.error)
        end
    end
end

# NOTE: based on test/pkg.jl::capture_stdout, but doesn't discard exceptions
macro grab_output(ex)
    quote
        mktemp() do fname, fout
            ret = nothing
            open(fname, "w") do fout
                redirect_stdout(fout) do
                    ret = $(esc(ex))
                end
            end
            ret, read(fname, String)
        end
    end
end

# Run some code on-device, returning captured standard output
macro on_device(ex)
    @gensym kernel
    esc(quote
        let
            function $kernel()
                $ex
                return
            end

            @cuda $kernel()
            synchronize()
        end
    end)
end

# helper function for sinking a value to prevent the callee from getting optimized away
@inline sink(i::Int32) =
    Base.llvmcall("""%slot = alloca i32
                     store volatile i32 %0, i32* %slot
                     %value = load volatile i32, i32* %slot
                     ret i32 %value""", Int32, Tuple{Int32}, i)
@inline sink(i::Int64) =
    Base.llvmcall("""%slot = alloca i64
                     store volatile i64 %0, i64* %slot
                     %value = load volatile i64, i64* %slot
                     ret i64 %value""", Int64, Tuple{Int64}, i)

function julia_script(code, args=``)
    # FIXME: this doesn't work when the compute mode is set to exclusive
    script = """using CUDAnative, CUDAdrv
                const CuArray = CUDAnative.CuHostArray
                device!($(device()))

                $code"""
    cmd = Base.julia_cmd()
    if Base.JLOptions().project != C_NULL
        cmd = `$cmd --project=$(unsafe_string(Base.JLOptions().project))`
    end
    cmd = `$cmd --eval $script $args`

    out = Pipe()
    err = Pipe()
    proc = run(pipeline(cmd, stdout=out, stderr=err), wait=false)
    close(out.in)
    close(err.in)
    wait(proc)
    proc.exitcode, read(out, String), read(err, String)
end

# tests that are conditionall broken
macro test_broken_if(cond, ex...)
    quote
        if $(esc(cond))
            @test_broken $(map(esc, ex)...)
        else
            @test $(map(esc, ex)...)
        end
    end
end


text:
#=
src\DataStructures\heap.jl:

=#

using BenchmarkTools
i=j=0
@btime begin
    m = n = 10
    for i = 1:m, j = 1:n
        arr[j, i] = tuple[j, i][1] #no method matching getindex (typeof(tuple), Int64, Int64)
    end
end

typeof(tuple[j, i])

#@btime begin
function convertTuple(tuple::Tuple{Int64,Int64})
    m = n = 10
    @inbounds for i = 1:m, j = 1:n
        #convert tuple to array
        arr[j, i] = deepcopy(size(collect(minimum(tuple[j, i]))))
        arr2[j, i] = deepcopy(size(collect(minimum(tuple[i, j]))))

    end
    return arr, arr2
    end

#end

typeof(tuple[j, i])

@btime begin
    m = n = 10
    @inbounds for i = 1:m, j = 1:n
        #convert tuple to array
        arr[j, i] = tuple[j, i][1]

    end
end

typeof(tuple[j, i])

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner