text:
@testset "Tests: Minimize depth for controlled-NOT gate" begin

    params = Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 5,
    "elementary_gates" => ["H_1", "H_2", "CNot_1_2", "Identity"],
    "initial_gate" => "Identity",
    "target_gate" => QCO.CNotRevGate(),
    "set_cnot_lower_bound" => 1,
    "set_cnot_upper_bound" => 1,
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    model_options = Dict{Symbol, Any}(:model_type => "compact_formulation_1", # Testing incorrect model_type
                                      :all_valid_constraints => 2)            # Testing incorrect all_valid_constraints

    result_qc = QCO.run_QCModel(params, CBC; options = model_options)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 5, atol=tol_0)
    @test isapprox(result_qc["solution"]["z_onoff_var"][1,1], 1, atol=tol_0) || isapprox(result_qc["solution"]["z_onoff_var"][2,1], 1, atol=tol_0)
    @test isapprox(result_qc["solution"]["z_onoff_var"][1,2], 1, atol=tol_0) || isapprox(result_qc["solution"]["z_onoff_var"][2,2], 1, atol=tol_0)
    @test isapprox(result_qc["solution"]["z_onoff_var"][3,3], 1, atol=tol_0)
    @test isapprox(result_qc["solution"]["z_onoff_var"][1,4], 1, atol=tol_0) || isapprox(result_qc["solution"]["z_onoff_var"][2,4], 1, atol=tol_0)
    @test isapprox(result_qc["solution"]["z_onoff_var"][1,5], 1, atol=tol_0) || isapprox(result_qc["solution"]["z_onoff_var"][2,5], 1, atol=tol_0)

end

@testset "Tests: Minimum CNOT swap gate decomposition" begin

    params = Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 4,
        "elementary_gates" => ["CNot_1_2", "CNot_2_1", "Identity"],
        "target_gate" => QCO.SwapGate(),
        "set_cnot_lower_bound" => 2,
        "set_cnot_upper_bound" => 3,
        "objective" => "minimize_cnot",
        "decomposition_type" => "exact"
        )

    model_options = Dict{Symbol, Any}(:model_type => "balas_formulation",
                                      :all_valid_constraints => 1,
                                      :unit_magnitude_constraints => true)

    result_qc = QCO.run_QCModel(params, CBC; options = model_options)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 3, atol = tol_0)
    @test isapprox(sum(result_qc["solution"]["z_onoff_var"][1:2,:]), 3, atol=tol_0)

end

@testset "Tests: Minimum depth U3(0,0,π/4) gate" begin

    params = Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 2,
    "elementary_gates" => ["U3_1", "U3_2", "Identity"],
    "target_gate" => QCO.kron_single_qubit_gate(2, QCO.U3Gate(0,0,π/4), "q1"),
    "U3_θ_discretization" => [0, π/2],
    "U3_ϕ_discretization" => [0],
    "U3_λ_discretization" => [0, π/4],
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    result_qc = QCO.run_QCModel(params, CBC)

    data = QCO.get_data(params)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 1, atol = tol_0)
    @test "Identity" in data["gates_dict"]["3"]["type"]
    @test isapprox(sum(result_qc["solution"]["z_onoff_var"][3,:]), 1, atol = tol_0)
    @test isapprox(sum(result_qc["solution"]["z_onoff_var"][4,:]), 1, atol = tol_0)
    @test data["gates_dict"]["4"]["qubit_loc"] == "qubit_1"
end

@testset "Tests: Minimum depth CU3(0,0,π/4) gate" begin

    params = Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 2,
    "elementary_gates" => ["CU3_1_2", "CU3_2_1", "Identity"],
    "target_gate" => QCO.CU3Gate(0, 0, π/4),
    "CU3_θ_discretization" => [0, π/2],
    "CU3_ϕ_discretization" => [0],
    "CU3_λ_discretization" => [0, π/4],
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    result_qc = QCO.run_QCModel(params, CBC)

    data = QCO.get_data(params)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 1, atol = tol_0)
    @test "Identity" in data["gates_dict"]["3"]["type"]
    @test isapprox(sum(result_qc["solution"]["z_onoff_var"][3,:]), 1, atol = tol_0)
    @test isapprox(sum(result_qc["solution"]["z_onoff_var"][4,:]), 1, atol = tol_0)
    @test data["gates_dict"]["4"]["qubit_loc"] == "qubit_1_2"
end

@testset "Tests: Minimum depth Rev CU3(0,0,π/4) gate" begin

    params = Dict{String, Any}(
    "num_qubits" => 3,
    "maximum_depth" => 2,
    "elementary_gates" => ["CU3_3_1", "CU3_1_3", "Identity"],
    "target_gate" => QCO.get_full_sized_gate("CU3_3_1", 3, angle = [0, 0, pi/4]),
    "CU3_θ_discretization" => [0, π/2],
    "CU3_ϕ_discretization" => [0],
    "CU3_λ_discretization" => [0, π/4],
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    result_qc = QCO.run_QCModel(params, CBC)

    data = QCO.get_data(params)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 1, atol = tol_0)
    @test "Identity" in data["gates_dict"]["3"]["type"]
    @test isapprox(sum(result_qc["solution"]["z_onoff_var"][3,:]), 1, atol = tol_0)
    @test isapprox(sum(result_qc["solution"]["z_onoff_var"][4,:]), 1, atol = tol_0)
    @test data["gates_dict"]["4"]["qubit_loc"] == "qubit_3_1"
end

@testset "Tests: Minimum depth RX, RY, RZ gate decomposition" begin

    params = Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 3,
    "elementary_gates" => ["RX_1", "RY_2", "RZ_1", "Identity"],
    "target_gate" => QCO.kron_single_qubit_gate(2, QCO.RXGate(π/4), "q1") * QCO.kron_single_qubit_gate(2, QCO.RYGate(π/4), "q2") * QCO.kron_single_qubit_gate(2, QCO.RZGate(π/4), "q1"),
    "RX_discretization" => [0, π/4],
    "RY_discretization" => [π/4],
    "RZ_discretization" => [π/2, π/4],
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    model_options = Dict{Symbol, Any}(:model_type => "balas_formulation",
                                      :commute_gate_constraints => true)

    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 3, atol = tol_0)

end

@testset "Tests: Minimum depth CRX, CRY, CRZ gate decomposition" begin

    params = Dict{String, Any}(
    "num_qubits" => 3,
    "maximum_depth" => 3,
    "elementary_gates" => ["CRX_1_2", "CRY_2_3", "CRZ_3_1", "Identity"],
    "target_gate" => QCO.get_full_sized_gate("CRX_1_2", 3, angle = pi/4) * QCO.get_full_sized_gate("CRY_2_3", 3, angle=pi/4) * QCO.get_full_sized_gate("CRZ_3_1", 3, angle=pi/4),
    "CRX_discretization" => [0, π/4],
    "CRY_discretization" => [π/4],
    "CRZ_discretization" => [π/2, π/4],
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    model_options = Dict{Symbol, Any}(:model_type => "balas_formulation",
                                      :commute_gate_constraints => true)

    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 3, atol = tol_0)

end

@testset "Tests: Minimum depth Rev CRX, CRY, CRZ gate decomposition" begin

    params = Dict{String, Any}(
    "num_qubits" => 3,
    "maximum_depth" => 3,
    "elementary_gates" => ["CRX_3_1", "CRY_3_1", "CRZ_1_3", "Identity"],
    "target_gate" => QCO.get_full_sized_gate("CRX_3_1", 3, angle=pi/4) * QCO.get_full_sized_gate("CRY_3_1", 3, angle=pi/4) * QCO.get_full_sized_gate("CRZ_1_3", 3, angle = pi/2),
    "CRX_discretization" => [0, π/4],
    "CRY_discretization" => [π/4],
    "CRZ_discretization" => [π/2, π/4],
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    model_options = Dict{Symbol, Any}(:model_type => "balas_formulation",
                                      :commute_gate_constraints => true)

    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 3, atol = tol_0)

end

@testset "Tests: 3-qubit RX gate decomposition" begin

    params = Dict{String, Any}(
    "num_qubits" => 3,
    "maximum_depth" => 2,
    "elementary_gates" => ["U3_1", "U3_2", "U3_3", "Identity"],
    "target_gate" => QCO.kron_single_qubit_gate(3, QCO.RXGate(π/4), "q3"),
    "U3_θ_discretization" => [0, π/4],
    "U3_ϕ_discretization" => [0, -π/2],
    "U3_λ_discretization" => [0, π/2],
    "objective" => "minimize_cnot",
    "decomposition_type" => "exact"
    )

    model_options = Dict{Symbol, Any}(:model_type => "balas_formulation")

    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)

    data = QCO.get_data(params)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 1, atol = tol_0)
    if isapprox(sum(result_qc["solution"]["z_onoff_var"][14,:]), 1, atol=tol_0)
        @test data["gates_dict"]["14"]["qubit_loc"] == "qubit_3"
        @test isapprox(rad2deg(data["gates_dict"]["14"]["angle"]["θ"]),  45, atol=tol_0)
        @test isapprox(rad2deg(data["gates_dict"]["14"]["angle"]["ϕ"]), -90, atol=tol_0)
        @test isapprox(rad2deg(data["gates_dict"]["14"]["angle"]["λ"]),  90, atol=tol_0)
    end

end

@testset "Tests: 3-qubit CRX gate decomposition" begin

    params = Dict{String, Any}(
        "num_qubits" => 3,
        "maximum_depth" => 2,
        "elementary_gates" => ["CU3_1_2", "CU3_2_3", "CU3_1_3", "Identity"],
        "target_gate" => QCO.get_full_sized_gate("CRX_1_3", 3, angle = pi/4),
        "CU3_θ_discretization" => [0, π/4],
        "CU3_ϕ_discretization" => [0, -π/2],
        "CU3_λ_discretization" => [0, π/2],
        "objective" => "minimize_cnot",
        "decomposition_type" => "exact"
    )

    model_options = Dict{Symbol, Any}(:model_type => "balas_formulation")

    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)

    data = QCO.get_data(params)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 1, atol = tol_0)
    if isapprox(sum(result_qc["solution"]["z_onoff_var"][18,:]), 1, atol=tol_0)
        @test data["gates_dict"]["18"]["qubit_loc"] == "qubit_1_3"
        @test isapprox(rad2deg(data["gates_dict"]["18"]["angle"]["θ"]),  45, atol=tol_0)
        @test isapprox(rad2deg(data["gates_dict"]["18"]["angle"]["ϕ"]), -90, atol=tol_0)
        @test isapprox(rad2deg(data["gates_dict"]["18"]["angle"]["λ"]),  90, atol=tol_0)
    end

end

@testset "Tests: feasibility problem" begin
    params = Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 3,
        "elementary_gates" => ["H_1", "H_2"],
        "target_gate" => QCO.kron_single_qubit_gate(2, QCO.HGate(), "q1"),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact"
        )

    model_options = Dict{Symbol, Any}(:all_valid_constraints => -1)

    result_qc = QCO.run_QCModel(params, CBC; options =  model_options)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    data = QCO.get_data(params)

    for i in keys(data["gates_dict"])
        if data["gates_dict"][i]["type"] == "H_1"
            @test isapprox(sum(result_qc["solution"]["z_onoff_var"][parse(Int64, i),:]), 3 , atol = tol_0)
        end
    end

end

@testset "Tests: JuMP set_start_value for z_onoff_var variables" begin
    function input_circuit()
        # [(depth, gate)]
        return [(1, "CNot_2_1"),
                (2, "S_1"),
                (3, "H_2"),
                (4, "S_2")
                ]
    end

    params = Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 4,
    "elementary_gates" => ["S_1", "S_2", "H_1", "H_2", "CNot_1_2", "CNot_2_1", "Identity"],
    "target_gate" => QCO.MGate(),
    "input_circuit" => input_circuit(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    result_qc = QCO.run_QCModel(params, CBC)
    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(sum(result_qc["solution"]["z_onoff_var"][6,:]), 1, atol=tol_0)
    @test isapprox(sum(result_qc["solution"]["z_onoff_var"][5,:]), 0, atol=tol_0)

end

@testset "Tests: Involutory gate constraints" begin

    params = Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 4,
    "elementary_gates" => ["S_1", "S_2", "H_1", "H_2", "CNot_1_2", "CNot_2_1", "Identity"],
    "target_gate" => QCO.MGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    gates_dict = QCO.get_data(params)["gates_dict"]
    involutory_gates = QCO.get_involutory_gates(gates_dict)
    @test length(involutory_gates) == 4 #excluding Identity gate

    for i in keys(gates_dict)
        if ("S_1" in gates_dict[i]["type"]) || ("S_2" in gates_dict[i]["type"])
            @test !(parse(Int, i) in involutory_gates)
        end
        if ("H_1" in gates_dict[i]["type"]) || ("H_2" in gates_dict[i]["type"]) || ("CNot_1_2" in gates_dict[i]["type"]) || ("CNot_2_1" in gates_dict[i]["type"])
            @test (parse(Int, i) in involutory_gates)
        end
    end

    result_qc = QCO.run_QCModel(params, CBC)
    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 4, atol = tol_0)

end

@testset "Tests: TIME_LIMIT for building results dict and log" begin

    params = Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 5,
    "elementary_gates" => ["S_1", "S_2", "H_1", "H_2", "CNot_1_2", "CNot_2_1", "Identity"],
    "target_gate" => QCO.MGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact"
    )

    model_options = Dict{Symbol, Any}(:time_limit => 1)

    result_qc = QCO.run_QCModel(params, CBC; options = model_options)
    @test result_qc["termination_status"] == MOI.TIME_LIMIT
    @test result_qc["primal_status"] == MOI.NO_SOLUTION

end

@testset "Tests: constraint_redundant_gate_product_pairs" begin

    function target_gate()
        T1 = QCO.get_full_sized_gate("U3_2", 2, angle = [0,π/2,π])
        T2 = QCO.get_full_sized_gate("U3_1", 2, angle = [π/2,π/2,-π/2])
        return T1*T2
    end

    params = Dict{String, Any}(
               "num_qubits" => 2,
               "maximum_depth" => 2,
               "elementary_gates" => ["U3_1", "U3_2", "Identity"],
               "target_gate" => target_gate(),
               "U3_θ_discretization" => [0, π/2],
               "U3_ϕ_discretization" => [π/2],
               "U3_λ_discretization" => [-π/2, π],
               "objective" => "minimize_depth")

    data = QCO.get_data(params)
    redundant_pairs = QCO.get_redundant_gate_product_pairs(data["gates_dict"])
    @test length(redundant_pairs) == 2
    @test redundant_pairs[1] == (1,4)
    @test redundant_pairs[2] == (5,7)
    @test isapprox(data["gates_dict"]["1"]["matrix"] * data["gates_dict"]["4"]["matrix"], data["gates_dict"]["2"]["matrix"], atol = tol_0)
    @test isapprox(data["gates_dict"]["5"]["matrix"] * data["gates_dict"]["7"]["matrix"], data["gates_dict"]["6"]["matrix"], atol = tol_0)

    result_qc = QCO.run_QCModel(params, CBC)
    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 2.0, atol = tol_0)

end

@testset "Tests: constraint_idempotent_gates" begin
    params = Dict{String, Any}(
    "num_qubits" => 2,
    "maximum_depth" => 3,
    "elementary_gates" => ["U3_1", "U3_2", "CNot_1_2", "Identity"],
    "target_gate" => QCO.CZGate(),
    "U3_θ_discretization" => [-π/2, 0, π/2],
    "U3_ϕ_discretization" => [0, π/2],
    "U3_λ_discretization" => [0, π/2])

    data = QCO.get_data(params)
    idempotent_pairs = QCO.get_idempotent_gates(data["gates_dict"])
    @test length(idempotent_pairs) == 2

    model_options = Dict{Symbol, Any}(:idempotent_gate_constraints => true)

    result_qc = QCO.run_QCModel(params, CBC; options = model_options)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 3.0, atol = tol_0)

    for i in keys(data["gates_dict"])
        if "Identity" in data["gates_dict"][i]["type"]
            @test isapprox(sum(result_qc["solution"]["z_onoff_var"][parse(Int64, i),:]), 0, atol = tol_0)
        end
    end

end

@testset "Tests: constraint_convex_hull_complex_gates" begin

    function target_gate()
        num_qubits = 4
        CV_1_4 = QCO.get_full_sized_gate("CV_1_4", num_qubits);
        CV_2_4 = QCO.get_full_sized_gate("CV_2_4", num_qubits);
        CV_3_4 = QCO.get_full_sized_gate("CV_3_4", num_qubits);
        CVdagger_3_4 = QCO.get_full_sized_gate("CVdagger_3_4", num_qubits);
        CNot_1_2 = QCO.get_full_sized_gate("CNot_1_2", num_qubits);
        CNot_2_3 = QCO.get_full_sized_gate("CNot_2_3", num_qubits);

        return CV_2_4 * CNot_1_2 * CV_3_4 * CV_1_4 * CNot_2_3 * CVdagger_3_4
    end

    params = Dict{String, Any}(
    "num_qubits" => 4,
    "maximum_depth" => 2,
    "elementary_gates" => ["CV_1_4", "CV_2_4", "CV_3_4", "CVdagger_3_4", "CNot_1_2", "CNot_2_3", "Identity"],
    "target_gate" => target_gate()
    )

    model_options = Dict{Symbol, Any}(:relax_integrality => true,
                                      :convex_hull_gate_constraints => true,
                                      :optimizer_log => false)

    result_qc = QCO.run_QCModel(params, CBC; options = model_options)

    @test result_qc["termination_status"] == MOI.OPTIMAL
    @test result_qc["primal_status"] == MOI.FEASIBLE_POINT
    @test isapprox(result_qc["objective"], 0.8275862068965518, atol = tol_0)
end


text:
import Base.transpose, Base.ctranspose
export transpose, ctranspose

transpose(x::Constant) = Constant(transpose(x.value))
ctranspose(x::Constant) = transpose(x)

function transpose(x::AffineExpr)
  vec_sz = x.size[1] * x.size[2]
  num_rows = x.size[1]
  num_cols = x.size[2]

  I = Array{Int64}(vec_sz)
  J = Array{Int64}(vec_sz)

  k = 1
  for r = 1:num_rows
    for c = 1:num_cols
      J[k] = (c - 1) * num_rows + r
      I[k] = (r - 1) * num_cols + c
      k += 1
    end
  end

  perm_constant = Constant(sparse(I, J, 1.0))

  vars_to_coeffs_map = Dict{UInt64, Constant}()
  for (v, c) in x.vars_to_coeffs_map
    vars_to_coeffs_map[v] = perm_constant * c
  end
  constant = perm_constant * x.constant

  this = AffineExpr(:transpose, (x,), vars_to_coeffs_map, constant, (num_cols, num_rows))
  this.evaluate = ()->x.evaluate()'
  return this
end

ctranspose(x::AffineExpr) = transpose(x::AffineExpr)

ext: jl


text:
using Base.Test

@testset "EarthMoversDistance test suite" begin
    include("test_emd.jl")
end


text:
using ProfileView
using JLD, ProfileView
@load "profdata.jld"
ProfileView.view(li, lidict=lidict)
readline(stdin)

max_stars_repo_name: probcomp/DepthRenderer

max_issues_repo_name: probcomp/DepthRenderer

max_forks_repo_name: probcomp/DepthRenderer


text:
immutable Record{Tf<:Tuple, To}
    fields::Tf
end

function Record{T<:Tuple}(t::T)
    To = Tuple{map(fieldtype, t)...}
    #Tov = Tuple{map(s->Vector{s},map(fieldtype, t))...}
    Record{T, To}(t)
end

# for dispatch on N
if VERSION >= v"0.6.0-dev.2123" # new type system, julia PR #18457
    include_string("const RecN{N,U} = Record{T,U} where T<:NTuple{N, Any}")
else
    include_string("typealias RecN{N,U} Record{NTuple{N}, U}")
end

@generated function tryparsenext{N, To}(r::RecN{N, To}, str, i, len, opts=default_opts)
    quote
        R = Nullable{To}
        i > len && @goto error

        Base.@nexprs $N j->begin
            @chk2 (val_j, i) = tryparsenext(r.fields[j], str, i, len, opts)
        end

        @label done
        return R(Base.@ntuple $N val), i

        @label error
        R(), i
    end
end

const PARSE_SUCCESS = 0x00
const PARSE_ERROR   = 0x01
const POOL_CROWDED  = 0x02
const POOL_OVERFLOW = 0x03

@generated function tryparsesetindex{N,To}(r::RecN{N,To}, str::AbstractString, i::Int, len::Int, columns::Tuple, row::Int, opts)
    quote
        R = Result{Int, Tuple{Int,Int,Int,UInt8}}
        err_field = 1
        ii = i
        err_code = PARSE_ERROR
        i > len && @goto error

        Base.@nexprs $N j->begin
            err_field = j
            @chk2 val_j, ii = tryparsenext(r.fields[j], str, i, len, opts)
            err = setcell!(columns[j], row, val_j, str)
            if err != PARSE_SUCCESS
                err_code = err
                @goto error
            end
            i = ii
        end

        @label done
        return R(true, i)

        @label error
        R(false, (ii, i, err_field, err_code)) # error char, start of error field, error field
    end
end

@inline function setcell!(col, i, val, str)
    col[i] = val
    PARSE_SUCCESS
end

@inline function setcell!(col::DataValueArray{Union{}}, i, val, str)
    PARSE_SUCCESS
end

@inline function setcell!(col::DataValueArray, i, val::Nullable, str)
    col[i] = DataValue(val)
    PARSE_SUCCESS
end

const MAX_POOL_FRACTION = 0.05
const ROWS_BEFORE_CROWDING = 510
@inline function setcell!{R}(col::PooledArray{String,R}, i, val::StrRange, str)
    if i > ROWS_BEFORE_CROWDING && length(col.pool) > i * MAX_POOL_FRACTION
        return POOL_CROWDED
    elseif length(col.pool) >= typemax(R)
        return POOL_OVERFLOW
    else
        nonallocating_setindex!(col, i, val, str)
        return PARSE_SUCCESS
    end
end

@inline function setcell!(col::Array{String,1}, i, val::StrRange, str)
    col[i] = alloc_string(str, val)
    PARSE_SUCCESS
end

# Weird hybrid of records and fields

immutable UseOne{T,R<:Record,use} <: AbstractToken{T}
    record::R
end

fieldtype{T}(::UseOne{T}) = T

function UseOne(fields::Tuple, use)
    r = Record(fields)
    UseOne{fieldtype(fields[use]), typeof(r), use}(r)
end
getthing{n}(x, ::Type{Val{n}}) = x[n]
function tryparsenext{T,S,use}(f::UseOne{T,S,use}, str, i, len, opts=default_opts)
    R = Nullable{T}
    @chk2 xs, i = tryparsenext(f.record, str, i, len, opts)

    @label done
    return R(getthing(xs, Val{use})), i

    @label error
    return R(), i
end

immutable Repeated{F, T, N}
    field::F
end

Repeated{F}(f::F, n) = Repeated{F, fieldtype(f), n}(f)

fieldtype{F,T,N}(::Repeated{F,T,N}) = NTuple{N,T}

@generated function tryparsenext{F,T,N}(f::Repeated{F,T,N}, str, i, len, opts=default_opts)
    quote
        R = Nullable{NTuple{N,T}}
        i > len && @goto error

        # pefect candidate for #11902
        Base.@nexprs $N j->begin
            @chk2 (val_j, i) = tryparsenext(f.field, str, i, len, opts)
        end

        @label done
        return R(Base.@ntuple $N val), i

        @label error
        R(), i
    end
end


text:
using Dates
using CodeTracking
using JLD2
using FileIO
using Logging

# Config files
using Pkg.TOML
using JSON

"""
    Experiment(dir, file, module_name, func_name, args_iter, config=nothing; modify_save_path=true)

This is a struct which encapsulates the idea of an experiment.
# Arguments
- `dir`: The directory associated w/ this experiment
- `file`: The path to the file containing the experiment module
- `module_name`: The module encapsulating your experiment (either a string or symbol)
- `func_name`: The function inside the module which runs your experiment
- `args_iter`: The settings your want to iterate over. This can be any iterator which returns `(job_id, params)`
- `config`: [Optional] The path of the config file your experiment was built from. See [`Experiment(config_path, save_path="")`](@ref)
- `modify_save_path`: If true and you are using ArgIterator or ArgLooper the arg_iter added to the experiment will have a new static_arg "save_dir" added.

"""
struct Experiment{I}
    dir::String
    file::String
    module_name::Symbol
    func_name::Symbol
    args_iter::I
    hash::UInt64
    config::Union{String, Nothing}
end

function Experiment(dir, file, module_name, func_name, arg_iter, config=nothing; modify_save_path=true)
    arg_iter = if modify_save_path
        cp_arg_iter = deepcopy(arg_iter)
        set_save_dir!(cp_arg_iter, joinpath(dir, "data"))
        cp_arg_iter
    else
        arg_iter
    end
    Experiment(dir, file, Symbol(module_name), Symbol(func_name), arg_iter, hash(string(arg_iter)), config)
end

"""
    Experiment(config_path, save_path="")

Build and experiment from either a toml file or a json.
"""
function Experiment(config_path, save_path="")

    ext = splitext(config_path)[end]
    dict = if ext == ".toml"
        TOML.parsefile(config_path)
    elseif ext == ".json"
        JSON.Parser.parsefile(config_path)
    else
        throw(ErrorException("Experiment currently doesn't support $(ext) files."))
    end

    cdict = dict["config"]
    save_dir = joinpath(save_path, cdict["save_dir"])
    exp_file = cdict["exp_file"]
    exp_module_name = cdict["exp_module_name"]
    exp_func_name = cdict["exp_func_name"]

    iter_type = cdict["arg_iter_type"]

    args = if iter_type == "iter"
        arg_order = get(cdict, "arg_list_order", nothing)

        @assert arg_order isa Nothing || all(sort(arg_order) .== sort(collect(keys(dict["sweep_args"]))))

        static_args_dict = get(dict, "static_args", Dict{String, Any}())
        static_args_dict["save_dir"] = joinpath(save_dir, "data")
        sweep_args_dict = dict["sweep_args"]

        for key ∈ keys(sweep_args_dict)
            if sweep_args_dict[key] isa String
                sweep_args_dict[key] = eval(Meta.parse(sweep_args_dict[key]))
            end
        end

        ArgIterator(sweep_args_dict,
                    static_args_dict,
                    arg_order=arg_order)

    elseif iter_type == "looper"
        static_args_dict = get(dict, "static_args", Dict{String, Any}())
        static_args_dict["save_dir"] = joinpath(save_dir, "data")
        args_dict_list = if "loop_args" ∈ keys(dict)
            [dict["loop_args"][k] for k ∈ keys(dict["loop_args"])]
        elseif "arg_file" ∈ keys(cdict)
            d = FileIO.load(cdict["arg_file"])
            d["args"]
        end

        run_param = cdict["run_param"]
        run_list = cdict["run_list"]
        if run_list isa String
            run_list = eval(Meta.parse(cdict["run_list"]))
        end
        ArgLooper(args_dict_list, static_args_dict, run_list, run_param)
    else
        throw("$(iter_type) not supported.")
    end

    set_save_dir!(args, joinpath(save_dir, "data"))

    experiment = Experiment(save_dir,
                            exp_file,
                            exp_module_name,
                            exp_func_name,
                            args,
                            config_path)

end

function _safe_mkdir(exp_dir)
    if !isdir(exp_dir)
        try
            mkdir(exp_dir)
        catch ex
            @info "Somebody else created directory... Waiting"
            if isa(ex, SystemError) && ex.errnum == 17
                sleep(0.1) # Other Process Made folder. Waiting...
            else
                throw(ex)
            end
        end
    end
end

function _safe_mkpath(exp_dir)
    if !isdir(exp_dir)
        try
            mkpath(exp_dir)
        catch ex
            @info "Somebody else created directory... Waiting"
            if isa(ex, SystemError) && ex.errnum == 17
                sleep(0.1) # Other Process Made folder. Waiting...
            else
                throw(ex)
            end
        end
    end
end

"""
    pre_experiment(exp::Experiment; kwargs...)

Experiment setup phase. This helps deal with all the setup that needs to occur to setup an experiment folder.
"""
function pre_experiment(exp::Experiment; kwargs...)
    create_experiment_dir(exp; kwargs...)
    add_experiment(exp)
end

"""
    create_experiment_dir(exp_dir; org_file=false, replace=false, tldr="")
    create_experiment_dir(exp::Experiment; kwargs...)

This creates an experiment directory and sets up a notes file (if org_file is true). It is recommended to use [`pre_experiment(exp::Experiment; kwargs...)`](@ref).

# Arguments
- `org_file=false`: Indicator on whether to create an org file w/ notes on the job (experimental).
- `replace=false`: If true the directory will be deleted. This is not recommended when calling create_experiment_dir from multiple jobs.
- `tldr=""`: A TLDR which will be put into the org file.

"""
function create_experiment_dir(exp_dir;
                               org_file=false,
                               replace=false,
                               tldr="")
    mkdir_fn = occursin("/", exp_dir) ? _safe_mkpath : _safe_mkdir
    if isdir(exp_dir)
        if !replace
            @info "directory already created - told to not replace..."
            return
        else
            @info "directory already created - told to replace..."
            rm(exp_dir; force=true, recursive=true)
            mkdir_fn(exp_dir)
        end
    else
        @info "creating experiment directory"
        mkdir_fn(exp_dir)
    end

    if isdir(joinpath(exp_dir, "data"))
        if !replace
            @info "data directory already created - told to not replace..."
            return
        else
            @info "directory already created - told to replace..."
            # rm(exp_dir; force=true, recursive=true)
            _safe_mkdir(joinpath(exp_dir, "data"))
        end
    else
        @info "creating experiment directory"
        _safe_mkdir(joinpath(exp_dir, "data"))
    end

    if org_file
        open(joinpath(exp_dir, "notes.org"), "w") do f
            write(f, "#+title: Experimental Notes for $(exp_dir)\n\n")
            write(f, "TL;DR: $(tldr)\n\n")
        end
    end

    return
end

function create_experiment_dir(exp::Experiment;
                               kwargs...)
    create_experiment_dir(exp.dir; kwargs...)
    return
end

"""
    append_experiment_notes_file(exp_dir, experiment_file, exp_module_name, exp_func_name, settings_file, config_file, args_iter, config)

Write into an org file which contains the notes for the run experiments

"""
function append_experiment_notes_file(
    exp_dir,
    experiment_file,
    exp_module_name,
    exp_func_name,
    settings_file,
    config_file,
    args_iter,
    config)

    if isfile(joinpath(exp_dir, "notes.org"))

        date_str = Dates.format(now(), dateformat"<yyyy-mm-dd e HH:MM:SS>")
        tab = "    "
        make_args_str = "nothing"

        open(joinpath(exp_dir, "notes.org"), "a") do f
            exp_str = "* " * date_str * "\n\n" *
                tab*"Git-head: $(git_head())\n" *
                tab*"Git-branch: $(git_branch())\n" *
                tab*"experiment file: $(experiment_file)\n" *
                tab*"experiment module: $(string(exp_module_name))\n" *
                tab*"experiment function: $(string(exp_func_name))\n\n" *
                tab*"settings file: $(basename(settings_file))\n" *
                (!(config isa Nothing) ? tab*"config file: $(basename(config_file))\n\n" : "\n") *
                tab*"#+BEGIN_src julia\n" *
                (args_iter isa ArgIterator ? tab*"dict = $(args_iter.dict)\n" : tab*"runs_iter=$(args_iter.runs_iter)\n") *
                (args_iter isa ArgIterator ? tab*"arg_order = $(args_iter.arg_order)\n" : tab*"arg_list = $(args_iter.dict_list)\n") *
                tab*"static_arg = $(args_iter.static_args)\n\n" *
                tab*"# Make Arguments\n" *
                tab*make_args_str*"\n" *
                tab*"#+END_src\n\n"
            write(f, exp_str)
        end
    end

end

"""
    add_experiment(exp_dir, experiment_file, exp_module_name, exp_func_name, args_iter, hash, config)
    add_experiment(exp::Experiment)

Set up details of experiment in `exp_dir`. This includes saving settings and config files, etc...
It is recommended to use [`pre_experiment(exp::Experiment; kwargs...)`](@ref).
"""
function add_experiment(exp_dir,
                        experiment_file,
                        exp_module_name,
                        exp_func_name,
                        args_iter,
                        hash::UInt64,
                        config=nothing)

    if "SLURM_ARRAY_TASK_ID" in keys(ENV)
        if parse(Int64, ENV["SLURM_ARRAY_TASK_ID"]) != 1
            job_id = parse(Int64, ENV["SLURM_ARRAY_TASK_ID"])
            @info "Told to not add all experiments... job_id : $(job_id) $(job_id == 1)"
            return
        end
    end

    @info "Adding Experiment to $(exp_dir)"

    settings_dir = "settings"
    settings_dir = joinpath(exp_dir, settings_dir)
    _safe_mkdir(settings_dir)

    settings_file = joinpath(settings_dir, "settings_0x"*string(hash, base=16)*".jld2")
    config_file = if config isa Nothing
        nothing
    else
        joinpath(settings_dir, "config_0x"*string(hash, base=16)*splitext(config)[end])
    end

    append_experiment_notes_file(
        exp_dir,
        experiment_file,
        exp_module_name,
        exp_func_name,
        settings_file,
        config_file,
        args_iter,
        config)

    make_args_str = "nothing"

    jldopen(settings_file, "w") do file
        file["args_iter"]=args_iter
        file["make_args_str"]=make_args_str
    end

    if !(config isa Nothing)
        cp(config, config_file; force=true)
    end

end

function add_experiment(exp::Experiment)
    add_experiment(exp.dir,
                   exp.file,
                   String(exp.module_name),
                   String(exp.func_name),
                   exp.args_iter,
                   exp.hash,
                   exp.config)
end

function post_experiment(exp_dir::AbstractString, canceled_jobs::Array{Int64, 1})

    if "SLURM_ARRAY_TASK_ID" in keys(ENV)
        @info "Post_experiment doesn't work with slurm job arrays."
        return
    end

    if isfile(joinpath(exp_dir, "notes.org"))
        tab = "\t"
        date_str = Dates.format(now(), dateformat"<yyyy-mm-dd e HH:MM:SS>")

        open(joinpath(exp_dir, "notes.org"), "a") do f

            post_exp_str = tab*"Post Experiment: \n" *
                tab*"Canceled Jobs: $(canceled_jobs)\n" *
                tab*"Ended: $(date_str)\n"
            write(f, post_exp_str)
        end
    end
end

function post_experiment(exp_dir::AbstractString, finished_job::Bool)
    @info "Post_experiment not supported with task jobs."
end

function post_experiment(exp::Experiment, job_ret)
    post_experiment(exp.dir, job_ret)
end

function exception_file(exc_file::AbstractString, job_id, exception, trace)

    if isfile(exc_file)
        @warn "$(exc_file) already exists. Overwriting..."
    end

    open(exc_file, "w") do f
        exception_string =
            "Exception for job_id: $(job_id)\n\n" * string(exception) * "\n\n"

        write(f, exception_string)
        Base.show_backtrace(f, trace)
    end

    return
end


text:
""" Algebras of operads of wiring diagrams.
"""
module WiringDiagramAlgebras
export oapply, query

using Compat: isnothing
import TypedTables

using ...Theories, ...CategoricalAlgebra, ...CategoricalAlgebra.FinSets
using ...CSetDataStructures: make_table
using ..UndirectedWiringDiagrams

""" Compose morphisms according to UWD.

The morphisms corresponding to the boxes, and optionally also the objects
corresponding to the junctions, are given by dictionaries indexed by
box/junction attributes. The default attributes are those compatible with the
`@relation` macro.
"""
function oapply(composite::UndirectedWiringDiagram, hom_map::AbstractDict,
                ob_map::Union{AbstractDict,Nothing}=nothing;
                hom_attr::Symbol=:name, ob_attr::Symbol=:variable)
  homs = [ hom_map[name] for name in subpart(composite, hom_attr) ]
  obs = isnothing(ob_map) ? nothing :
    [ ob_map[name] for name in subpart(composite, ob_attr) ]
  oapply(composite, homs, obs)
end

# UWD algebras of multi(co)spans
################################

function oapply(composite::UndirectedWiringDiagram,
                spans::AbstractVector{<:Multispan}; Ob=nothing, Hom=nothing)
  @assert nboxes(composite) == length(spans)
  # FIXME: This manual type inference is hacky and bad. The right solution is to
  # extend `Multi(co)span` with type parameters that allow abstract types.
  if isnothing(Ob); Ob = typejoin(mapreduce(typeof∘apex, typejoin, spans),
                                  mapreduce(eltype∘feet, typejoin, spans)) end
  if isnothing(Hom); Hom = mapreduce(eltype∘legs, typejoin, spans) end
  junction_feet = Vector{Ob}(undef, njunctions(composite))

  # Create bipartite free diagram whose vertices of types 1 and 2 are the UWD's
  # boxes and junctions, respectively.
  diagram = BipartiteFreeDiagram{Ob,Hom}()
  add_vertices₁!(diagram, nboxes(composite), ob₁=map(apex, spans))
  add_vertices₂!(diagram, njunctions(composite))
  for (b, span) in zip(boxes(composite), spans)
    for (p, leg) in zip(ports(composite, b), legs(span))
      j = junction(composite, p)
      add_edge!(diagram, b, j, hom=leg)
      foot = codom(leg)
      if isassigned(junction_feet, j)
        foot′ = junction_feet[j]
        foot == foot′ || error("Feet of spans are not equal: $foot != $foot′")
      else
        junction_feet[j] = foot
      end
    end
  end
  all(isassigned(junction_feet, j) for j in junctions(composite)) ||
    error("Limits with isolated junctions are not supported")
  diagram[:ob₂] = junction_feet

  # The composite multispan is given by the limit of this diagram.
  lim = limit(diagram)
  outer_legs = map(junction(composite, outer=true)) do j
    e = first(incident(diagram, j, :tgt))
    legs(lim)[src(diagram, e)] ⋅ hom(diagram, e)
  end
  Multispan(ob(lim), outer_legs)
end

function oapply(composite::UndirectedWiringDiagram,
                cospans::AbstractVector{<:StructuredMulticospan{L}},
                junction_feet::Union{AbstractVector,Nothing}=nothing) where L
  @assert nboxes(composite) == length(cospans)
  if isnothing(junction_feet)
    junction_feet = Vector{first(dom(L))}(undef, njunctions(composite))
  else
    @assert njunctions(composite) == length(junction_feet)
  end

  # Create bipartite free diagram whose vertices of types 1 and 2 are the UWD's
  # junctions and boxes, respectively.
  diagram = BipartiteFreeDiagram{codom(L)...}()
  add_vertices₁!(diagram, njunctions(composite))
  add_vertices₂!(diagram, nboxes(composite), ob₂=map(apex, cospans))
  for (b, cospan) in zip(boxes(composite), cospans)
    for (p, leg, foot) in zip(ports(composite, b), legs(cospan), feet(cospan))
      j = junction(composite, p)
      add_edge!(diagram, j, b, hom=leg)
      if isassigned(junction_feet, j)
        foot′ = junction_feet[j]
        foot == foot′ || error("Feet of cospans are not equal: $foot != $foot′")
      else
        junction_feet[j] = foot
      end
    end
  end
  diagram[:ob₁] = map(L, junction_feet)

  # Find, or if necessary create, an outgoing edge for each junction. The
  # existence of such edges is an assumption for colimits of bipartite diagrams.
  # The edges are also needed to construct inclusions for the outer junctions.
  junction_edges = map(junctions(composite)) do j
    out_edges = incident(diagram, j, :src)
    if isempty(out_edges)
      x = ob₁(diagram, j)
      v = add_vertex₂!(diagram, ob₂=x)
      add_edge!(diagram, j, v, hom=id(x))
    else
      first(out_edges)
    end
  end

  # The composite multicospan is given by the colimit of this diagram.
  colim = colimit(diagram)
  outer_js = junction(composite, outer=true)
  outer_legs = map(junction_edges[outer_js]) do e
    hom(diagram, e) ⋅ legs(colim)[tgt(diagram, e)]
  end
  outer_feet = junction_feet[outer_js]
  StructuredMulticospan{L}(Multicospan(ob(colim), outer_legs), outer_feet)
end

# Queries via UWD algebras
##########################

""" Evaluate a conjunctive query on an attributed C-set.

The query is a undirected wiring diagram whose boxes and ports are assumed to be
named through attributes `:name` and `:port_name`/`:outer_port_name`. To define
such a diagram, use the named form of the [`@relation`](@ref) macro.

This function straightforwardly wraps the [`oapply`](@ref) method for
multispans, which implements the UWD algebra of multispans.
"""
function query(X::AbstractACSet, diagram::UndirectedWiringDiagram,
               params=NamedTuple(); table_type::Type=TypedTables.Table)
  # For each box in the diagram, extract span from ACSet.
  spans = map(boxes(diagram), subpart(diagram, :name)) do b, name
    apex = FinSet(nparts(X, name))
    legs = map(subpart(diagram, ports(diagram, b), :port_name)) do port_name
      FinDomFunction(X, port_name == :_id ? name : port_name)
    end
    Multispan(apex, legs)
  end

  # Add an extra box and corresponding span for each parameter.
  if !isempty(params)
    diagram = copy(diagram)
    spans = vcat(spans, map(collect(pairs(params))) do (var, value)
      box = add_part!(diagram, :Box, name=:_const)
      add_part!(diagram, :Port, port_name=:_value,
                box=box, junction=incident(diagram, var, :variable))
      SMultispan{1}(ConstantFunction(value, FinSet(1)))
    end)
  end

  # Call `oapply` and make a table out of the resulting span.
  outer_names = subpart(diagram, :outer_port_name)
  outer_span = oapply(diagram, spans, Ob=SetOb, Hom=FinDomFunction{Int})
  table = NamedTuple{Tuple(outer_names)}(Tuple(map(collect, outer_span)))
  make_table(table_type, table)
end

end


text:
module RocketWithLatestOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: with_latest()" begin

    println("Testing: operator with_latest()")

    run_testset([
        (
            source = of(0) |> with_latest(of(0.0)),
            values = @ts([ (0, 0.0), c ]),
            source_type = Tuple{Int, Float64}
        ),
        (
            source = of(0) |> with_latest(of(0.0), of('0')),
            values = @ts([ (0, 0.0, '0'), c ]),
            source_type = Tuple{Int, Float64, Char}
        ),
        (
            source = from(1:5) |> with_latest(of(1)),
            values = @ts([ (1, 1), (2, 1), (3, 1), (4, 1), (5, 1), c ])
        ),
        (
            source = from(1:2) |> async(0) |> with_latest(from(1:2) |> async(0)),
            values = @ts([ (1, 1) ] ~ [ (2, 2) ] ~ c)
        ),
        (
            source = from(1:2) |> with_latest(completed()),
            values = @ts(c)
        ),
        (
            source = completed(Int) |> with_latest(of(1)),
            values = @ts(c)
        ),
        (
            source = completed(Int) |> with_latest(faulted("e")),
            values = @ts(e("e"))
        ),
        (
            source = completed() |> with_latest(completed()),
            values = @ts(c)
        ),
        (
            source = completed() |> with_latest(never()),
            values = @ts(c)
        ),
        (
            source = completed(Int) |> with_latest(of(1) |> async(0)),
            values = @ts(c)
        ),
        (
            source = never(Int) |> with_latest(of(1)),
            values = @ts()
        ),
        (
            source = never() |> with_latest(completed()),
            values = @ts(c)
        ),
        (
            source = never() |> with_latest(never()),
            values = @ts()
        ),
        (
            source = never() |> with_latest(faulted("e")),
            values = @ts(e("e"))
        ),
        (
            source = faulted("e") |> with_latest(of(1), of(2)),
            values = @ts(e("e"))
        ),
        (
            source = faulted("e") |> with_latest(faulted("e2")),
            values = @ts(e("e2"))
        ),
        (
            source = faulted("e") |> with_latest(completed()),
            values = @ts(c)
        ),
        (
            source = faulted("e") |> with_latest(never()),
            values = @ts(e("e"))
        )
    ])

end

end


text:
using Formatting
using ArgParse, JSON
using NCDatasets

println("""
This program output the forcing file with QFLX_TEMP and QFLX_SALT from WKRSTT and WKRSTS
 variables in history files. It also needs the data folder generated by program
main_concat_and_convert_units.jl that has monthly TEMP and SALT profile.
""")

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

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

function makeTimeFile(output_file, years)

    dom = [31.0, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    sum_dom = sum(dom)
    sum_dom == 365 || throw(ErrorException("Sum of dom is $(sum_dom) rather than 365."))

    _t    = zeros(Float64, length(dom))
    _bnds = zeros(Float64, 2, length(dom))

    for m=1:length(dom)
        #bnds[m, 1] = beg of month  m
        #bnds[m, 2] = end of month  m
        if m==1
            _bnds[1, m] = 0.0
        else
            _bnds[1, m] = _bnds[2, m-1]
        end

        _bnds[2, m] = _bnds[1, m] + dom[m]

        _t[m] = (_bnds[1, m] + _bnds[2, m]) / 2.0
    end

    t    = zeros(Float64, 12 * years)
    bnds = zeros(Float64, 2, 12 * years)

    for y = 1:years
        i_offset = (y-1)*12
        t_offset = (y-1)*sum_dom
        t[i_offset+1:i_offset+12]       .+= _t    .+ t_offset
        bnds[:, i_offset+1:i_offset+12] .+= _bnds .+ t_offset
    end

    Dataset(output_file, "c") do ds

        defDim(ds, "time", Inf)
        defDim(ds, "d2", 2)

        defVar(ds, "time", t, ("time", ), ; attrib = Dict(
            "long_name" => "time",
            "bounds"    => "time_bound",
            "calendar"  =>  "noleap",
            "units"     => "days since 0001-01-01 00:00:00",
        ))

        defVar(ds, "time_bound", bnds, ("d2", "time"), ; attrib = Dict(
            "long_name" => "boundaries for time-averaging interval",
            "units"     => "days since 0001-01-01 00:00:00",
        ))
    end
end

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--hist-dir"
            help = "Casename"
            arg_type = String
            required = true

        "--data-dir"
            help = "Casename"
            arg_type = String
            required = true

        "--year-rng"
            help = "Casename"
            arg_type = Int64
            nargs    = 2
            required = true

        "--output-file"
            help = "Casename"
            arg_type = String
            default = "forcing_cyclic.nc"

    end

    return parse_args(s)
end

parsed = parse_commandline()

JSON.print(parsed,4)

beg_yr, end_yr = parsed["year-rng"]
println("Beg year: $(beg_yr)")
println("End year: $(end_yr)")

yr_rng_str  = format( "{:04d}-{:04d}", beg_yr, end_yr )
yr_rng_eval = format( "{:04d}..{:04d}", beg_yr, end_yr )

output_file = parsed["output-file"]

println("Make time file: tmp_time.nc")
makeTimeFile("tmp_time.nc", 1)

println("Making mean profile")

mkpath("tmp")
for m = 1:12
    m_str = format( "{:02d}", m)
    pleaseRun(`bash -c "ncra -v WKRSTT,WKRSTS,dz_cT,lat_sT,mask_sT,area_sT -O $(parsed["hist-dir"])/*.h0.*.{$(yr_rng_eval)}-$(m_str).nc tmp/monthly_mean_$(m_str).nc"`)
end

println("Output file : $(output_file)")
pleaseRun(`bash -c "ncrcat -O tmp/monthly_mean_{01..12}.nc $output_file"`)
pleaseRun(`ncks -O -3 $output_file $output_file`)
pleaseRun(`ncrename -d Nx,nlon -d Ny,nlat -d Nz,z_t -v WKRSTT,QFLXT -v WKRSTS,QFLXS $output_file`)
#pleaseRun(`ncap2 -O -s 'QFLX_TEMP=QFLX_TEMP*3996*1026;' $output_file $output_file`)
pleaseRun(`ncks -A -v SALT,z_w_top,z_w_bot  $(parsed["data-dir"])/monthly/SALT_monthly.nc $output_file`)
pleaseRun(`ncks -A -v TEMP                  $(parsed["data-dir"])/monthly/TEMP_monthly.nc $output_file`)
pleaseRun(`ncks -A -v HMXL                  $(parsed["data-dir"])/monthly/HMXL_monthly.nc $output_file`)
pleaseRun(`ncks -A -v time                  tmp_time.nc $output_file`)

rm("tmp_time.nc", force=true)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
#= Embedded code

Provides syntax cases of macro-embedded code in julia.
=#

using Cxx
using PyCall
using RCall
using Markdown

# Expected: the contents of the strings should look like R code.
R"x <- 3"

R"""
 y <- fit(x, q)
"""

# Expected: the contents should look like markdown. Embedded code blocks
# should look like blocks in their respective languages.
md"""
# This is a header

Here is some markdown text

 - This is a list item
 - This is another item

This should look like a block of julia code:
```julia
x = x + y
```

And this should look like a block of C++:
```c++
auto x = 42;
```
"""

# Expected: the string delimiters should look like strings, and the string
# contents should look like C++
cxx"""
#include <cxx>
"""

# Expected: this should look like an ordinary julia expression
q = 42

# Expected: the macro captures should be highlighted specially
icxx"f(($x));"

icxx"""
  auto x = $y;
  auto z = q($x);
  return x;
"""

str = "the sky is $(g)"

# Expected: should look like javascript inside
js"""
alert("Hello")
"""

# Expected: should look like python inside
py"""
import numpy as np
print 'Hello, world!'
"""

# Expected: should look like an ordinary julia expression
x = 2 + 2

max_stars_repo_name: jarrison/julia-syntax-test-cases

max_issues_repo_name: jarrison/julia-syntax-test-cases

max_forks_repo_name: twadleigh/julia-syntax-test-cases


text:
abstract AbstractDistanceField
abstract AbstractUnsignedDistanceField <: AbstractDistanceField
abstract AbstractSignedDistanceField <: AbstractDistanceField
"""
Abstract to categorize geometry primitives of dimensionality `N` and
the numeric element type `T`.
"""
abstract AbstractGeometry{N, T}
abstract AbstractMesh{VertT, FaceT} <: AbstractGeometry
abstract GeometryPrimitive{N, T} <: AbstractGeometry{N, T}

"""
Abstract to classify Simplices. The convention for N starts at 1, which means
a Simplex has 1 point. A 2-simplex has 2 points, and so forth. This convention
is not the same as most mathematical texts.
"""
abstract AbstractSimplex{N,T} <: FixedVector{N,T}

"""
A `Simplex` is a generalization of an N-dimensional tetrahedra and can be thought
of as a minimal convex set containing the specified points.

* A 0-simplex is a point.
* A 1-simplex is a line segment.
* A 2-simplex is a triangle.
* A 3-simplex is a tetrahedron.

Note that this datatype is offset by one compared to the traditional
mathematical terminology. So a one-simplex is represented as `Simplex{2,T}`.
This is for a simpler implementation.

It applies to infinite dimensions. The structure of this type is designed
to allow embedding in higher-order spaces by parameterizing on `T`.
"""
immutable Simplex{N,T} <: AbstractSimplex{N,T}
    _::NTuple{N,T}
end

immutable Normal{N, T} <: FixedVector{N, T}
    _::NTuple{N, T}
end
immutable TextureCoordinate{N, T} <: FixedVector{N, T}
    _::NTuple{N, T}
end

"""
A Face is typically used when constructing subtypes of `AbstractMesh` where
the `Face` should not reproduce the vertices, but rather index into them.
Face is parameterized by:

* `N` - The number of vertices in the face.
* `T` - The type of the indices in the face, a subtype of Integer.
* `O` - The offset relative to Julia arrays. This helps reduce copying when
communicating with 0-indexed systems such ad OpenGL.
"""
immutable Face{N, T, IndexOffset} <: FixedVector{N, T}
    _::NTuple{N, T}
end

"""
A `HyperRectangle` is a generalization of a rectangle into N-dimensions.
Formally it is the cartesian product of intervals, which is represented by the
`origin` and `width` fields, whose indices correspond to each of the `N` axes.
"""
immutable HyperRectangle{N, T} <: GeometryPrimitive{N, T}
    origin::Vec{N, T}
    widths::Vec{N, T}
end

immutable HyperCube{N, T} <: GeometryPrimitive{N, T}
    origin::Vec{N, T}
    width::T
end

"""
A `HyperSphere` is a generalization of a sphere into N-dimensions.
A `center` and radius, `r`, must be specified.
"""
immutable HyperSphere{N, T} <: GeometryPrimitive{N, T}
    center::Point{N, T}
    r::T
end

immutable SimpleRectangle{T} <: GeometryPrimitive{2, T}
    x::T
    y::T
    w::T
    h::T
end

# TODO remove before 0.2.0 tag
const Rectangle = SimpleRectangle

"""
A rectangle in 3D space.
"""
immutable Quad{T} <: GeometryPrimitive{3, T}
    downleft::Vec{3, T}
    width   ::Vec{3, T}
    height  ::Vec{3, T}
end

immutable Pyramid{T} <: GeometryPrimitive{3, T}
    middle::Point{3, T}
    length::T
    width ::T
end

immutable Particle{N, T} <: GeometryPrimitive{N, T}
    position::Point{N, T}
    velocity::Vec{N, T}
end

"""
A `SignedDistanceField` is a uniform sampling of an implicit function.
The `bounds` field corresponds to the sampling space intervals on each axis.
The `data` field represents the value at each point whose exact location
can be rationalized from `bounds`.
The type is parameterized by:

* `N` - The dimensionality of the sampling space.
* `SpaceT` - the type of the space where we will uniformly sample.
* `FieldT` - the type resulting from evaluation of the implicit function.

Note that decoupling the space and field types is useful since geometry can
be formulated with integers and distances can be measured with floating points.
"""
type SignedDistanceField{N,SpaceT,FieldT} <: AbstractSignedDistanceField
    bounds::HyperRectangle{N,SpaceT}
    data::Array{FieldT,N}
end

"""
The `HomogenousMesh` type describes a polygonal mesh that is useful for
computation on the CPU or on the GPU.
All vectors must have the same length or must be empty, besides the face vector
Type can be void or a value, this way we can create many combinations from this
one mesh type.
This is not perfect, but helps to reduce a type explosion (imagine defining
every attribute combination as a new type).
"""
immutable HomogenousMesh{VertT, FaceT, NormalT, TexCoordT, ColorT, AttribT, AttribIDT} <: AbstractMesh{VertT, FaceT}
    vertices            ::Vector{VertT}
    faces               ::Vector{FaceT}
    normals             ::Vector{NormalT}
    texturecoordinates  ::Vector{TexCoordT}
    color               ::ColorT
    attributes          ::AttribT
    attribute_id        ::Vector{AttribIDT}
end

"""
AbstractFlexibleGeometry{T}

AbstractFlexibleGeometry refers to shapes, which are somewhat mutable.
"""
abstract AbstractFlexibleGeometry{T}
typealias AFG AbstractFlexibleGeometry

"""
FlexibleConvexHull{T}

Represents the convex hull of finitely many points. The number of points is not fixed.
"""
immutable FlexibleConvexHull{T} <: AFG{T}
    _::Vector{T}
end

"""
FlexibleSimplex{T}

Represents a Simplex whos number of vertices is not fixed.
"""
immutable FlexibleSimplex{T} <: AFG{T}
    _::Vector{T}
end

"""
AbstractConvexHull

Groups all geometry types, that can be described as the convex hull of finitely
many points.
"""
typealias AbstractConvexHull Union{Simplex, FlexibleConvexHull, FlexibleSimplex,
HyperCube, HyperRectangle} # should we parametrize ACH by the type of points T?


text:
using ChainRulesCore

@non_differentiable generate_horizontal_rules()
@non_differentiable generate_vertical_rules()
@non_differentiable hamiltonian(model::HamiltonianModel)
@non_differentiable swapgate(n1,n2)
@non_differentiable Hubbard_hand(model::HamiltonianModel)
@non_differentiable hop_pair_hand(t,γ)

function ChainRulesCore.rrule(::typeof(parity_conserving),T::Union{Array,CuArray})
	result = parity_conserving(T)
	function pullback_parity_conserving(ΔT)
		return (NoTangent(),parity_conserving(ΔT))
	end
	return result,pullback_parity_conserving
end


text:
using DataStructures: BinaryMinHeap

mutable struct Heap{T}
        h::BinaryMinHeap{T}
        maxsz::Int
end

Heap{T}(maxsz::Int) where {T} = Heap{T}(BinaryMinHeap{T}(), maxsz)

import Base.length
length(h::Heap{T}) where {T} = length(h.h)

import Base.isempty
isempty(h::Heap{T}) where {T} = isempty(h.h)

import Base.push!
function push!(h::Heap{T}, v::T) where {T}
        if length(h) == h.maxsz && top(h.h) < v; pop!(h.h); end
        if length(h) < h.maxsz; push!(h.h, v); end
end

import Base.empty!
empty!(h::Heap{T}) where {T} = empty!(h.h.valtree)

import Base.iterate
iterate(h::Heap{T}, i=1) where {T} = i > length(h) ? nothing : (h.h.valtree[i], i + 1)


text:
## BlinkDisplay
##
## Use Blink for a standalone display
## The following window management is
## shamelessly copied from Tim Holy's GadflyDisplay in Immerse.jl

## closefig with mouse needs a callback...
## closefig should decrement list.

export figure, scf, gcf, closefig, closeall
export figure_size, figure_position, figure_raise, figure_floating, figure_title, figure_title!

type BlinkDisplay <: Display
    figs::Dict{Int,Blink.AtomShell.Window}
    fig_order::Vector{Int}
    current_fig::Int
    next_fig::Int

    BlinkDisplay() = new(Dict{Int,Blink.AtomShell.Window}(), Int[], 0, 1)
end

const _display = BlinkDisplay()

"""
`figure(;name="Figure \$n", width=400, height=400)` creates a new
figure window for displaying plots.

`figure(n)` raises the `n`th figure window and makes it the current
default plotting window, and returns the
"""
function figure(;name::AbstractString="Figure $(nextfig(_display))",
                 width::Integer=660,    # TODO: make configurable
                height::Integer=500)
    i = nextfig(_display)
    w = Blink.Window(Blink.@d(:show => false))
    empty_page(w)
    size(w, width, height)
    title(w, name)
    addfig(_display, i, w)
    front(w)
    w
end

function figure(i::Integer; displayfig::Bool = true)
    switchfig(_display, i)
    w = curfig(_display)
    front(w)
    w
end

function addfig(d::BlinkDisplay, i::Int, fig)
    @assert !haskey(d.figs,i)
    d.figs[i] = fig
    push!(d.fig_order, i)
    while haskey(d.figs,d.next_fig)
        d.next_fig += 1
    end
    d.current_fig = i
end

hasfig(d::BlinkDisplay, i::Int) = haskey(d.figs,i)

function switchfig(d::BlinkDisplay, i::Int)
    haskey(d.figs,i) && (d.current_fig = i)
end

function getfig(d::BlinkDisplay, i::Int)
    haskey(d.figs,i) ? d.figs[i] : error("no figure with index $i")
end

function curfig(d::BlinkDisplay)
    cur = d.current_fig
    if !active(d.figs[cur])
        cur = cur + 1
        cur > length(d.figs) && return figure()
    end
    d.figs[cur]
end

## keep trying, but not if not alive
function nextfig(d::BlinkDisplay)
    i = d.next_fig
    !haskey(d.figs, i) && return i
    while(!active(d.figs[i]))
        i = i + 1
        i > length(d.figs) && break
    end
    i
end

function dropfig(d::BlinkDisplay, i::Int)
    haskey(d.figs,i) || return
    delete!(d.figs, i)
    splice!(d.fig_order, findfirst(d.fig_order,i))
    d.next_fig = min(d.next_fig, i)
    d.current_fig = isempty(d.fig_order) ? 0 : d.fig_order[end]
end

gcf() = _display.current_fig

"`scf()` (\"show current figure\") raises (makes visible) the current figure"
scf() = gcf() > 0 && front(figure(gcf()))

"""

`closefig(n)` closes the `n`th figure window.

`closefig()` closes the current figure window.
"""
closefig() = closefig(_display.current_fig)

function closefig(i::Integer)
    ## XXX Should closefig remove from _display? Clear up Resources XXX
    fig = getfig(_display,i)
    close(fig)
    _display.current_fig = _display.next_fig
    _display.next_fig = _display.current_fig  + 1
end

"`closeall()` closes all existing figure windows."
closeall() = (map(closefig, keys(_display.figs)); nothing)

## Some Blink specific things we can do to windows
"Adjust size of figure, defaulting to current one"
figure_size(w::Integer, h::Integer, figure_number=gcf()) = size(figure(figure_number), w,h)

"set position of figure"
figure_position(x::Integer, y::Integer, figure_number=gcf()) = position(figure(figure_number), x, y)

"Raise figure"
figure_raise(figure_number=gcf()) = front(figure(figure_number))

"Set if figure is always on top."
figure_floating(flag::Bool=true, figure_number=gcf()) = floating(figure(figure_number), flag)

"Set title of figure"
figure_title(figure_number=gcf()) = title(figure(figure_number))
figure_title!(str::AbstractString, figure_number=gcf()) = title(figure(figure_number), str)

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

function Base.display(d::BlinkDisplay, plt)
    isempty(d.figs) && figure()
    ## need to sleep or make asychronous...
    w = curfig(d)
    scf()

    code = "Plotly.newPlot(PLOT, $(Plots.get_series_json(plt)), $(Plots.get_plot_json(plt)));"
    Blink.js(w, DThree.JSString(code))
end

## Make a web page
## from Plots.jl
function open_browser_window(filename::AbstractString)
    @osx_only   return run(`open $(filename)`)
    @linux_only return run(`xdg-open $(filename)`)
    @windows_only return run(`$(ENV["COMSPEC"]) /c start $(filename)`)
    warn("Unknown OS... cannot open browser window.")
end


text:
export cg

function cg(A::Union{SparseMatrixCSC{TF,Int},Matrix{TF}},b::Vector{TF}; kwargs...) where {TF<:Real}
	x = zeros(TF,size(A,2)) # pre-allocate
	return cg(v -> mul!(x, transpose(A), v, TF(1.0), TF(0.0)),b; kwargs...) # multiply with transpose of A for efficiency
end

cg(A,b::Vector;kwargs...)= cg(x -> A*x,b;kwargs...)

"""
This is a slightly modified version from the Julia Package KrylovMethods

x,flag,err,iter,resvec = cg(A,b,tol=1e-2,maxIter=100,M=1,x=[],out=0)

(Preconditioned) Conjugate Gradient applied to the linear system A*x = b, where A is assumed
to be symmetric positive semi definite.

Input:

  A       - matrix or function computing A*x
  b       - right hand side vector
  tol     - error tolerance
  maxIter - maximum number of iterations
  M       - preconditioner, a function that computes M\\x
  x       - starting guess
  out     - flag for output (-1: no output, 0: only errors, 1: final status, 2: residual norm at each iteration)

Output:

  x       - approximate solution
  flag    - exit flag ( 0 : desired tolerance achieved,
  					 -1 : maxIter reached without converging
  					 -2 : Matrix A is not positive definite
  					 -3 : cg stalled, i.e. two consecutive residuals have same norm
  					 -9 : right hand side was zero)
  err     - norm of relative residual, i.e., norm(A*x-b)/norm(b)
  iter    - number of iterations
  resvec  - norm of relative residual at each iteration
"""
# function cg(A::Function,b::Vector; tol::Real=1e-2,maxIter::Integer=100,M::Function=identity,x::Vector=[],out::Integer=0,
# 	storeInterm::Bool=false)
function cg(A::Function,b::Vector{TF}; tol::Real=1e-2,maxIter::Integer=100,M::Function=identity,x::Vector=[],out::Integer=0) where {TF <: Real}
	n = length(b)

	if norm(b)==0; return zeros(eltype(b),n), -9, 0f0, 0, [0f0]; end
	if isempty(x)
		x = zeros(TF,n)#x = zeros(eltype(b),n)
		r = copy(b)
	else
		r = b .- A(x)::Vector{TF}
	end
	z = M(r)::Vector{TF}
	p = copy(z)::Vector{TF}

    # if storeInterm
    #     X = zeros(n,maxIter)	# allocate space for intermediates
    # end

	nr0  = norm(b)

	if out==2
		println("=== cg ===")
		println(@sprintf("%4s\t%7s","iter","relres"))
	end

	resvec = zeros(TF,maxIter)
	iter   = 1 # makes iter available outside the loop
	flag   = -1

	if norm(r)/nr0<=tol
		flag = 0;
		return x,flag,resvec[iter],iter,resvec[1:iter]
	end

	#TF    = typeof(x[1])
 	alpha = TF(0)
	beta  = TF(0)
	lastIter = 0
	for iter=1:maxIter
		lastIter = iter

		Ap 	  = A(p)::Vector{TF}
		gamma = dot(r,z)
		#gamma = BLAS.dot(n, r,1,z,1)#
		alpha = gamma / dot(p,Ap)
		#alpha = gamma / BLAS.dot(n, p,1,Ap,1)#dot(p,Ap)

		if alpha==Inf || alpha<0
			flag = -2; break
		end

		BLAS.axpy!(n,alpha,p,1,x,1) # x = alpha*p+x
		#if storeInterm; X[:,iter] = x; end
		BLAS.axpy!(n,-alpha,Ap,1,r,1) # r -= alpha*Ap

		#resvec[iter] = BLAS.nrm2(n, r, 1) / nr0#
		resvec[iter]  = norm(r)/nr0
		if out==2
			println(iter,resvec[iter])
		end
		if resvec[iter] <= tol
			flag = 0; break
		end

		z    = M(r)::Vector{TF}
		#beta = BLAS.dot(n, z,1,r,1) / gamma
		beta  = dot(z,r) / gamma
		# the following two lines are equivalent to p = z + beta*p
		#p = BLAS.scal!(n,beta,p,1)::Vector{TF}
		#p = BLAS.axpy!(n,TF(1.0),z,1,p,1)::Vector{TF}
		p  = BLAS.axpby!(TF(1.0), z, beta, p)::Vector{TF} # Overwrite Y with X*a + Y*b, where a and b are scalars. Return Y
	end

	if out>=0
		if flag==-1
			println("cg iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",
																								maxIter,resvec[lastIter],tol)
		elseif flag==-2
			println("Matrix A in cg has to be positive definite.")
		elseif flag==0 && out>=1
			println("cg achieved desired tolerance at iteration %d. Residual norm is %1.2e.",lastIter,resvec[lastIter])
		end
	end
    return x,flag,resvec[lastIter],lastIter,resvec[1:lastIter]
end

#specifically for CDS/DIA storage matrices
# function cg(R::Matrix{TF},offset::Vector{TI},b::Vector{TF},tol::Real,maxIter::Integer,x::Vector{TF},out::Integer) where {TF <: Real, TI <: Integer}
# 	n  	   = length(b)
# 	ndiags = size(R,2)
# 	r  	   = zeros(TF,n)

# 	if norm(b)==0; return zeros(eltype(b),n),-9,0.0,0,[0.0]; end

# 	CDS_MVp_MT(n,ndiags,R,offset,x,r)
# 	r .= -r .+ b
# 		#r = b .- A(x)::Vector{TF}

# 	z = copy(r)

# 	p = copy(r)::Vector{TF}

# 	nr0  = norm(b)

# 	resvec = zeros(TF,maxIter)
# 	iter   = 1 # makes iter available outside the loop
# 	flag   = -1

# 	if norm(r)/nr0<=tol
# 		flag = 0;
# 		return x,flag,resvec[iter],iter,resvec[1:iter]
# 	end

#  	alpha  	 = TF(0)
# 	beta     = TF(0)
# 	lastIter = 0
# 	Ap  	 = zeros(TF,n)
# 	for iter=1:maxIter
# 		lastIter = iter

# 		#Ap 	  = A(p)::Vector{TF}
# 		fill!(Ap,TF(0))
# 		CDS_MVp_MT(n,ndiags,R,offset,p,Ap)
# 		gamma = dot(r,z)
# 		#gamma = BLAS.dot(n, r,1,z,1)#
# 		alpha = gamma / dot(p,Ap)
# 		#alpha = gamma / BLAS.dot(n, p,1,Ap,1)#dot(p,Ap)

# 		if alpha==Inf || alpha<0
# 			flag = -2; break
# 		end

# 		BLAS.axpy!(n,alpha,p,1,x,1) # x = alpha*p+x
# 		BLAS.axpy!(n,-alpha,Ap,1,r,1) # r -= alpha*Ap

# 		#resvec[iter] = BLAS.nrm2(n, r, 1) / nr0#
# 		resvec[iter]  = norm(r)/nr0
# 		if resvec[iter] <= tol
# 			flag = 0; break
# 		end

# 		#z    = M(r)::Vector{TF}
# 		z = copy(r)
# 		#beta = BLAS.dot(n, z,1,r,1) / gamma
# 		beta  = dot(z,r) / gamma
# 		# the following two lines are equivalent to p = z + beta*p
# 		#p = BLAS.scal!(n,beta,p,1)::Vector{TF}
# 		#p = BLAS.axpy!(n,TF(1.0),z,1,p,1)::Vector{TF}
# 		p  = BLAS.axpby!(TF(1.0), z, beta, p)::Vector{TF} # Overwrite Y with X*a + Y*b, where a and b are scalars. Return Y
# 	end
#
# 	if out>=0
# 		if flag==-1
# 			println("cg iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",
# 																								maxIter,resvec[lastIter],tol)
# 		elseif flag==-2
# 			println("Matrix A in cg has to be positive definite.")
# 		elseif flag==0 && out>=1
# 			println("cg achieved desired tolerance at iteration %d. Residual norm is %1.2e.",lastIter,resvec[lastIter])
# 		end
# 	end

#     return x,flag,resvec[lastIter],lastIter,resvec[1:lastIter]

# end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
@testset "Scalar1PBasis" begin

##
using ACE
using Printf, Test, LinearAlgebra, StaticArrays
using ACE: evaluate, evaluate_d, Rn1pBasis, Ylm1pBasis,
      EuclideanVectorState, Product1pBasis
using Random: shuffle
using ACEbase.Testing: dirfdtest, fdtest, print_tf

##

maxdeg = 10
r0 = 1.0
rcut = 3.0
trans = trans = PolyTransform(1, r0)
bscal = ACE.scal1pbasis(:x, :k, maxdeg, trans, rcut)

ACE.evaluate_d(bscal, 1.0)

##

end


text:
"""
Calculate value, gradient, and hessian of the variational ELBO using a PSF
image rather than a mixture of Gaussians.
"""

module DeterministicVIImagePSF

using StaticArrays, DiffBase

import ..DeterministicVI:
    ElboArgs, ElboIntermediateVariables, maximize_f,
    StarPosParams, GalaxyPosParams, CanonicalParams, VariationalParams,
    SourceBrightness, GalaxyComponent, SkyPatch,
    load_source_brightnesses, add_elbo_log_term!,
    accumulate_source_pixel_brightness!, subtract_kl_all_sources!, KL_HELPER_POOL

import ..Model:
    populate_gal_fsm!, getids, ParamSet, linear_world_to_pix, lidx,
    BvnComponent, GalaxyCacheComponent, GalaxySigmaDerivs,
    get_bvn_cov, galaxy_prototypes, linear_world_to_pix,
    Image, eval_psf, CatalogEntry, init_source

import ..SensitiveFloats:
    SensitiveFloat, zero_sensitive_float_array,
    multiply_sfs!, add_scaled_sfs!, clear!

import ..Infer: load_active_pixels!, get_sky_patches

import ..PSF: get_psf_at_point, trim_psf

import WCS

include("deterministic_vi_image_psf/sensitive_float_fft.jl")
include("deterministic_vi_image_psf/kernels.jl")
include("deterministic_vi_image_psf/fsm_matrices.jl")
include("deterministic_vi_image_psf/elbo_image_psf.jl")

export elbo_likelihood_with_fft!, FSMSensitiveFloatMatrices,
       initialize_fsm_sf_matrices!, initialize_fft_elbo_parameters

function infer_source_fft(images::Vector{Image},
                          neighbors::Vector{CatalogEntry},
                          entry::CatalogEntry)
   cat_local = vcat([entry], neighbors)
   vp = Vector{Float64}[init_source(ce) for ce in cat_local]
   patches = get_sky_patches(images, cat_local)
   load_active_pixels!(images, patches)

   ea_fft, fsm_mat = initialize_fft_elbo_parameters(
       images, vp, patches, [1], use_raw_psf=false)
   elbo_fft_opt = get_fft_elbo_function(ea_fft, fsm_mat)
   maximize_f(elbo_fft_opt, ea_fft)

   vp[1]
end

end


text:
import Random
export Objective, PSO, update_velocity!, update_position!, optimize!

"""
    Objective

Stores the objective function and the related `search_space` the dimensionality.
"""
struct Objective{T}
    fun::Function
    search_space::AbstractVector{<:AbstractVector{<:T}} # TODO: other  object for search_space, maybe allow ellipses and more, maybe constraints?
    function Objective(fun::Function, search_space::AbstractVector{<:AbstractVector})
        isempty(search_space) && error("`search_space` is empty")
        for e in search_space
            length(e) != 2 && error("$e has not exactly two elements")
            !(e[1] < e[2]) && error("$(e[1]) is not less than $(e[2])")
        end
        return new{eltype(eltype(search_space))}(fun, search_space)
    end
end
function Objective(fun::Function, search_space::AbstractVector, number_of_dimensions::Integer)
    length(search_space) != 2 && error("search must contain the lower and the upper limit")
    search_space = [search_space for i in 1:number_of_dimensions]
    Objective(fun, search_space)
end
function Objective(fun::Function, lower_limits::AbstractVector, upper_limits::AbstractVector)
    length(lower_limits) == length(upper_limits) || error("`lower_limits` and `upper_limits` must have the same length")
    search_space = [[lower_limits[i], upper_limits[i]] for i in 1:number_of_dimensions]
    Objective(fun, search_space)
end

const LOWER_BOUND_IDX = 1
const UPPER_BOUND_IDX = 2

"""
    PSO

Packs all data for the pso.
"""
mutable struct PSO{T, C, R, I<:Integer}
    objective::Objective{T}

    coefficients::AbstractVector{C}

    iterations::I

    position::AbstractVector{<:AbstractVector{<:T}}
    velocity::AbstractVector{<:AbstractVector{<:T}}

    pos_best::AbstractVector{<:AbstractVector{<:T}}
    results_best::AbstractVector{R}

    random_mat::AbstractArray{C, 3}
    rand1::AbstractVector{<:AbstractVector{C}}
    rand2::AbstractVector{<:AbstractVector{C}}

    update_velocity!::Function
    update_position!::Function
    confinement!::Function
    get_localbest::Function

    neighbours::Neighbourhood{I}

    compare::Function
end
function PSO(objective::Objective{T}, neighbours::Neighbourhood{I}; compare::Function=<,
             coefficients::AbstractVector{C}=[1.0/(2.0*log(2.0)), 0.5 + log(2), 0.5 + log(2)]) where T where C where I<:Integer

    number_of_particles = length(neighbours)
    number_of_dimensions = length(objective.search_space)

    # Initialization - SPSO 2011 default
    position = [Vector{T}(undef, number_of_dimensions) for i in 1:number_of_particles]
    velocity = [Vector{T}(undef, number_of_dimensions) for i in 1:number_of_particles]

    for p in 1:number_of_particles, d in 1:number_of_dimensions
        serach_type = eltype(objective.search_space[d])
        # TODO: use Distributions.jl instead of project
        position[p][d] = project(rand(serach_type), zero(serach_type), oneunit(serach_type), objective.search_space[d]...)
        velocity[p][d] = project(rand(serach_type), zero(serach_type), oneunit(serach_type),
                                          objective.search_space[d][LOWER_BOUND_IDX].-position[p][d],
                                          objective.search_space[d][UPPER_BOUND_IDX].-position[p][d])
    end

    results_best = objective.fun.(position)

    pos_best = deepcopy(position)

    rearrange!(neighbours, results_best, compare)

    random_mat = zeros(C, number_of_dimensions, length(position), 2)
    rand1 = [view(random_mat, :, i, 1) for i in 1:length(position)]
    rand2 = [view(random_mat, :, i, 2) for i in 1:length(position)]

    return PSO(objective, coefficients, zero(I), position, velocity,
        pos_best, results_best, random_mat, rand1, rand2, update_velocity!,
        update_position!, confinement!, get_localbest, neighbours, compare)
end

# update functions
"""
    update_velocity!(position::AbstractVector{T}, velocity::AbstractVector{T}, personal_best::AbstractVector{T},
                     local_best::AbstractVector{T}, rand_a::AbstractVector{T}, rand_b::AbstractVector{T}, w::Number,
                     c1::Number, c2::Number) where T<:Number

Default equation for calculating the velocity for the next step.
"""
@inline function update_velocity!(position::AbstractVector{<:T}, velocity::AbstractVector{<:T},
                                  personal_best::AbstractVector{<:T}, local_best::AbstractVector{<:T},
                                  rand_a::AbstractVector{C}, rand_b::AbstractVector{C},
                                  coefficients::AbstractVector{C}) where T where C
    w, c1, c2 = coefficients
    @inbounds for pidx in eachindex(velocity)
        velocity[pidx] = w*velocity[pidx] + c1*rand_a[pidx]*(personal_best[pidx] - position[pidx]) + c2*rand_b[pidx]*(local_best[pidx] - position[pidx])
    end
    return
end

# @inline function update_velocity_geometric!(position::AbstractVector{<:T}, velocity::AbstractVector{<:T}, personal_best::AbstractVector{<:T}, local_best::AbstractVector{<:T}, rand_a::AbstractVector{<:T}, rand_b::AbstractVector{<:T}, w::Number, c1::Number, c2::Number) where T<:Number
#     velocity .= w.*velocity .+ c1.*rand_a.*(personal_best .- position) .+ c2.*rand_b.*(local_best .- position)
# end

"""
    function update_position!(position::AbstractVector{<:T}, velocity::AbstractVector{<:T},
                              personal_best::AbstractVector{<:T}, local_best::AbstractVector{<:T}) where T<:Number

Default update of the next position.
"""
@inline function update_position!(position::AbstractVector{<:T}, velocity::AbstractVector{<:T},
                                  personal_best::AbstractVector{<:T}, local_best::AbstractVector{<:T}) where T
    position .+= velocity
    return
end

"""
    confinement!(results_best::AbstractVector{<:Number}, neighbours::AbstractVector{<:Integer}, compare::Function)

Checks, if a particle is over the edge of search space. If yes for a dimension, it will be set on the edge and the velocity for this dimension will be inverted and decreased by the half.
"""
@inline function confinement!(position::AbstractVector{<:T}, velocity::AbstractVector{<:T}, limits::AbstractVector{<:AbstractVector{<:T}}) where T
    @inbounds for d in eachindex(position)
        if position[d] < limits[d][LOWER_BOUND_IDX]
            velocity[d] *= -0.5
            position[d] = limits[d][LOWER_BOUND_IDX]
        end
        if position[d] > limits[d][UPPER_BOUND_IDX]
            velocity[d] *= -0.5
            position[d] = limits[d][UPPER_BOUND_IDX]
        end
    end
    return
end

"""
    get_localbest(results_best::AbstractVector{<:Number}, neighbours::AbstractVector{<:Integer}, compare::Function)

Returns the neighbour with the best evaluation of an particle.
"""
@inline function get_localbest(results_best::AbstractVector{<:Number}, neighbours::AbstractVector{<:Integer},
                               compare::Function)
    min_idx = one(eltype(neighbours))
    @inbounds min_val = results_best[neighbours[min_idx]]
    @inbounds for neig_idx in (min_idx+min_idx):lastindex(neighbours)
        if compare(results_best[neighbours[neig_idx]], min_val)
            min_idx = neig_idx
            min_val = results_best[neighbours[neig_idx]]
        end
    end
    return neighbours[min_idx]
end

"""
    evaluate!(position::AbstractVector{T}, results_best::AbstractVector{<:Number}, pos_best::AbstractVector{T}, idx::Int, cmp::Function, obj::Function) where T<:Number

Evaluates the a particle `position` and copies it to the particles memory if the compare function `cmp` says true. To the compare funciton the result from  objective function is applied as first argument, the current best evaluation value as second argument.
"""
@inline function evaluate!(position::AbstractVector{<:T}, results_best::AbstractVector, pos_best::AbstractVector{<:T}, idx::Integer, cmp::Function, obj::Function) where T
    result = obj(position)
    @inbounds if cmp(result, results_best[idx])
        copyto!(pos_best, position)
        results_best[idx] = result
    end
    return
end

"""
    optimize!(pso::PSO, number_of_iterations::Int; additional_arguments::Dict{Symbol, <:Any}=Dict{Symbol, Any}())

Iterates an PSO object for `number_of_iterations`. `additional_arguments` would be applied to the objective function
as named arguments.

# Examples

```jldoctest
julia> objective = Objective(TestFunctions.sphere, 2, (-10.0, 10.0))
julia> neighbours = GlobalNeighbourhood(20)
julia> pso = PSO(objective, neighbours)
julia> optimize!(pso, 10_000)

```
"""
function optimize!(pso::PSO{T, C, R, I}, number_of_iterations::I) where {T, C, R, I<:Integer}
    num_particles::Int = pso.neighbours.particle_number

    for iter in one(I):number_of_iterations
        Random.rand!(pso.random_mat)

        @inbounds for particle_idx in 1:num_particles
            l_best = pso.get_localbest(pso.results_best,  pso.neighbours[particle_idx], pso.compare)
            pso.update_velocity!(pso.position[particle_idx], pso.velocity[particle_idx], pso.pos_best[particle_idx],
                             pso.pos_best[l_best], pso.rand1[particle_idx], pso.rand2[particle_idx],
                             pso.coefficients)
            pso.update_position!(pso.position[particle_idx], pso.velocity[particle_idx], pso.pos_best[particle_idx],
                             pso.pos_best[l_best])
            pso.confinement!(pso.position[particle_idx], pso.velocity[particle_idx], pso.objective.search_space)
            evaluate!(pso.position[particle_idx], pso.results_best, pso.pos_best[particle_idx], particle_idx, pso.compare, pso.objective.fun)
        end

        rearrange!(pso.neighbours, pso.results_best, pso.compare)
        pso.iterations += one(I)
    end
    return
end

"""
    getoptimum(pso::PSO)

Returns a tuple with the n-dimensional position of the current found optimum and the function value on this position.

# Examples

```jldoctest
julia> getoptimum()

```
"""
@inline function getoptimum(pso::PSO{T, C, R, I})  where {T, C, R, I<:Integer}
    all_best = get_localbest(pso.results_best, one(I):pso.neighbours.particle_number, pso.compare)
    return pso.pos_best[all_best], pso.results_best[all_best]
end


text:
# Section 5.3.1 in Ericson, C. (2004). Real-time collision detection
function intersect(l::LineSegment{Dim, T}, plane::Hyperplane{Dim, T}) where {Dim, T}
    r = (plane.d - (plane.𝗻 ⋅ l.𝘅₁.coord))/(plane.𝗻 ⋅l.𝘂)
    (r ≥ 0 && r ≤ 1) && return true, l(r)
    return false, nan(Point{Dim,T})
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using ChaoticEncryption
using TestImages
using Test

@testset "ChaoticEncryption.jl" begin

    @testset "Logistic Map" begin
        @test typeof(logistic_key(0.01, 3.97, 20)) == Array{Int64, 1}
        @test length(logistic_key(0.01, 3.97, 20)) == 20
    end

    @testset "Lorenz System of Differential Equations" begin
        @test typeof(lorenz_key(0.01, 0.02, 0.03, 20)) == Tuple{Array{Int64, 1}, Array{Int64, 1}, Array{Int64, 1}}
        @test length(lorenz_key(0.01, 0.02, 0.03, 20)) == 3
        @test length(lorenz_key(0.01, 0.02, 0.03, 20)[1]) == 20
        @test length(lorenz_key(0.01, 0.02, 0.03, 20)[2]) == 20
        @test length(lorenz_key(0.01, 0.02, 0.03, 20)[3]) == 20
    end

    @testset "Substitution Encryption" begin
        img = testimage("mandrill")
        height, width = size(img)
        keys = logistic_key(0.01, 3.97, height * width)
        substitution_encryption(img, keys; path_for_result="../test_images/encrypted.png")
        @test isfile("../test_images/encrypted.png")

        keys = logistic_key(0.01, 3.97, 20)
        @test_throws ArgumentError("Number of keys must be equal to height * width of image.") substitution_encryption(img, keys; path_for_result="../test_images/encrypted.png")
    end

    @testset "Substitution Decryption" begin
        img = testimage("mandrill")
        height, width = size(img)
        keys = logistic_key(0.01, 3.97, height * width)

        substitution_decryption("../test_images/encrypted.png", keys; path_for_result="../test_images/decrypted.png")
        @test isfile("../test_images/decrypted.png")

        substitution_decryption(img, keys; path_for_result="../test_images/decrypted.png")
        @test isfile("../test_images/decrypted.png")

        @test_throws ArgumentError("image must be of the type ::String or ::Array{RGB{N0f8},2}") substitution_decryption(5, keys; path_for_result="../test_images/decrypted.png")

        keys = logistic_key(0.01, 3.97, 20)
        @test_throws ArgumentError("Number of keys must be equal to height * width of image.") substitution_decryption(img, keys; path_for_result="../test_images/decrypted.png")

        rm("../test_images", recursive=true)
    end
end


text:
# Creates a operator that permutes rows, in blocks of size
# length(perm)
struct PermutationOperator{T,DS,RS} <: Operator{T}
    perm::Vector{Int}
    domainspace::DS
    rangespace::RS
end
PermutationOperator{T}(prm, ds::DS, rs::RS) where {T, DS<:Space,RS<:Space} =
    PermutationOperator{T,DS,RS}(prm, ds, rs)
PermutationOperator(prm, ds, rs) = PermutationOperator{Int}(prm, ds, rs)
PermutationOperator(prm) = PermutationOperator(prm, ℓ⁰, ℓ⁰)

domainspace(P::PermutationOperator) = P.domainspace
rangespace(P::PermutationOperator) = P.rangespace

convert(::Type{Operator{T}},P::PermutationOperator) where {T} =
    PermutationOperator{T}(P.perm, P.domainspace, P.rangespace)

function bandwidths(P::PermutationOperator)
    dfs=P.perm-(1:length(P.perm))
    -minimum(dfs),maximum(dfs)
end

function getindex(P::PermutationOperator{T},k::Integer,j::Integer) where T
    n=length(P.perm)
    if (k-1)÷n == (j-1)÷n  # diagonal blocks
        k=mod(k-1,n)+1
        j=mod(j-1,n)+1
        convert(T,P.perm[k]==j)
    else
        zero(T)
    end
end

# the permutation that rearranges a to be b
multiplyperm(b,a) = Int[a[bk] for bk in b]
perm(a::Vector,b::Vector) = multiplyperm(invperm(sortperm(b)),sortperm(a))
perm(a::Tuple,b::Tuple) = perm(collect(a),collect(b))

struct NegateEven{T,DS,RS} <: Operator{T}
    domainspace::DS
    rangespace::RS
end

NegateEven{T}(ds::DS, rs::RS) where {DS<:Space,RS<:Space,T} =
    NegateEven{T,DS,RS}(ds, rs)
NegateEven(ds, rs) = NegateEven{Int}(ds, rs)
NegateEven() = NegateEven(ℓ⁰,ℓ⁰)

domainspace(P::NegateEven) = P.domainspace
rangespace(P::NegateEven) = P.rangespace

convert(::Type{Operator{T}},P::NegateEven) where {T} =
    NegateEven{T}(P.domainspace, P.rangespace)

bandwidths(P::NegateEven) = (0,0)

getindex(P::NegateEven{T},k::Integer,j::Integer) where {T} =
    k==j ? (iseven(k) ? -one(T) : one(T)) : zero(T)


text:
########## Compute the GBP marginal vector ##########
function marginal(gbp::ContinuousModel)
    @inbounds for i in gbp.graph.iterateMarginal
        Mcol = gbp.graph.meanDirect[i]; Wcol = gbp.graph.weightDirect[i]

        for j in gbp.graph.colptrMarginal[i]
            Mcol += gbp.inference.meanFactorVariable[j] / gbp.inference.varianceFactorVariable[j]
            Wcol += 1 / gbp.inference.varianceFactorVariable[j]
        end
        gbp.inference.variance[i] = 1 / Wcol
        gbp.inference.mean[i] = Mcol * gbp.inference.variance[i]
    end
end

########## Compute the GBP marginal vector for the tree factor graph ##########
function marginal(gbp::ContinuousTreeModel)
    @inbounds for i = 1:gbp.graph.Nvariable
        Mcol = gbp.graph.meanDirect[i]; Wcol = gbp.graph.weightDirect[i]

        for j in gbp.graph.incomingToVariable[i]
            Mcol += gbp.inference.meanFactorVariable[j] / gbp.inference.varianceFactorVariable[j]
            Wcol += 1 / gbp.inference.varianceFactorVariable[j]
        end
        gbp.inference.variance[i] = 1 / Wcol
        gbp.inference.mean[i] = Mcol * gbp.inference.variance[i]
    end
end

########## Compute the BP marginal vector for the tree factor graph ##########
function marginal(bp::DiscreteTreeModel)
    @inbounds for i = 1:bp.graph.Nvariable
        bp.inference.marginal[i] = copy(bp.graph.messageDirect[i])
        for k in bp.graph.incomingToVariable[i]
            bp.inference.marginal[i] .*= bp.inference.messageFactorVariable[k]
        end
        normalizeConstant = sum(bp.inference.marginal[i])
        bp.inference.marginal[i] = bp.inference.marginal[i] ./ normalizeConstant
    end
end

########## Compute the BP unnormalized marginal vector for the tree factor graph ##########
function marginalUnnormalized(bp::DiscreteTreeModel)
    @inbounds for i = 1:bp.graph.Nvariable
        bp.inference.marginal[i] = copy(bp.graph.messageDirect[i])
        for k in bp.graph.incomingToVariable[i]
            bp.inference.marginal[i] .*= bp.inference.messageFactorVariable[k]
        end
    end
end

########## Dynamic the GBP update ##########
@inline function dynamicFactor!(gbp::ContinuousModel; factor = 0::Int64, mean = 0, variance = 0)
    if (gbp.system.jacobianTranspose.colptr[factor + 1] - gbp.system.jacobianTranspose.colptr[factor]) == 1
        idx = gbp.system.jacobianTranspose.colptr[factor]
        variable = trunc(Int, gbp.system.jacobianTranspose.rowval[idx])

        gbp.graph.meanDirect[variable] -= gbp.system.observation[factor] * gbp.system.jacobianTranspose[variable, factor] / gbp.system.variance[factor]
        gbp.graph.weightDirect[variable] -= gbp.system.jacobianTranspose[variable, factor]^2 / gbp.system.variance[factor]

        gbp.graph.meanDirect[variable] += mean * gbp.system.jacobianTranspose[variable, factor] / variance
        gbp.graph.weightDirect[variable] += gbp.system.jacobianTranspose[variable, factor]^2 / variance
    else
        k = gbp.graph.dynamic[factor]
        @inbounds for i in gbp.graph.rowptr[k]
            gbp.graph.meanIndirect[i] = mean
            gbp.graph.varianceIndirect[i] = variance
        end
    end
    gbp.system.observation[factor] = mean
    gbp.system.variance[factor] = variance
end

######### Ageing the GBP update ##########
@inline function ageingVariance!(gbp::ContinuousModel; factor = 0::Int64, initial = 0, limit = 0, model = 0::Int64, a = 0, b = 0, tau = 0)
    if gbp.system.variance[factor] < limit
        Nvariable = gbp.system.jacobianTranspose.colptr[factor + 1] - gbp.system.jacobianTranspose.colptr[factor]
        if Nvariable == 1
            idx = gbp.system.jacobianTranspose.colptr[factor]
            variable = trunc(Int, gbp.system.jacobianTranspose.rowval[idx])

            gbp.graph.meanDirect[variable] -= gbp.system.observation[factor] * gbp.system.jacobianTranspose[variable, factor] / gbp.system.variance[factor]
            gbp.graph.weightDirect[variable] -= gbp.system.jacobianTranspose[variable, factor]^2 / gbp.system.variance[factor]
        end

        if model == 1
            gbp.system.variance[factor] = a * tau + initial
        elseif model == 2
            gbp.system.variance[factor] = a * log10((tau + 1 + b) / (1 + b)) + initial
        elseif model == 3
            gbp.system.variance[factor] = initial * (1 + b)^(a * tau)
        end

        if gbp.system.variance[factor] > limit
            gbp.system.variance[factor] =  limit
        end

        if Nvariable == 1
            gbp.graph.meanDirect[variable] += gbp.system.observation[factor] * gbp.system.jacobianTranspose[variable, factor] / gbp.system.variance[factor]
            gbp.graph.weightDirect[variable] += gbp.system.jacobianTranspose[variable, factor]^2 / gbp.system.variance[factor]
        else
            k = gbp.graph.dynamic[factor]
            @inbounds for i in gbp.graph.rowptr[k]
                gbp.graph.varianceIndirect[i] = gbp.system.variance[factor]
            end
        end
    end
end

########## Set damping parameters ##########
function damping!(gbp::ContinuousModel; prob::Float64 = 0.6, alpha::Float64 = 0.4)
    gbp.graph.alphaNew = fill(1.0, gbp.graph.Nlink)
    gbp.graph.alphaOld = fill(0.0, gbp.graph.Nlink)
    bernoulliSample = randsubseq(collect(1:gbp.graph.Nlink), prob)
    @inbounds for i in bernoulliSample
        gbp.graph.alphaNew[i] = 1.0 - alpha
        gbp.graph.alphaOld[i] = alpha
    end
end


text:
#use imshow(image_name) in order to see the images if you're running from this file

using Images, ImageView, ImageFiltering

misty = load("misty.jpg")
#= You can find misty.jpg in the resources folder,
make sure to put it in the same directory you run the code =#

cropped_misty = misty[500:800, 500:800] # we can crop the image by specifying a range of indices

reduced_misty = misty[1:5:end, 1:5:end] # this reduces the image size by only taking every 5th pixel.

red_square = zeros(RGB, 50, 50) #this creates a black 50x50 image
red_square[1:50, 1:50] .= RGB(1.0, 0.0, 0.0) #this makes every pixel value red since 1, 0, 0 is red in RGB

#similarly, we can modify the pixels of misty
squaerd_misty = misty[50:100, 50:100] .= RGB(0.0, 0.5, 0.8)
imshow(misty) #we need this because for some reason jupyter only displays the modified pixels (disregard the output below)

invert(color::AbstractRGB) = RGB(1-color.r, 1-color.g, 1-color.b)
# this function inverts r, g, b values of any image it is passed
inverted_misty = invert.(misty)

gray_misty = Gray.(misty)
#we can also convert misty to grayscale like this

#this was i hope a basic introduction to image modification in julia.
#please check other notebooks for more cool examples!

max_stars_repo_name: gunalpms/imageprojl

max_issues_repo_name: gunalpms/imageprojl

max_forks_repo_name: gunalpms/imageprojl