text:
abstract type TierrianInstruction end

struct NOP_0 <: TierrianInstruction end
struct NOP_1 <: TierrianInstruction end

struct OR_1 <: TierrianInstruction end
struct SH_L <: TierrianInstruction end
struct ZERO <: TierrianInstruction end

struct IF_0 <: TierrianInstruction end

struct SUB_AB <: TierrianInstruction end
struct SUB_AC <: TierrianInstruction end

struct INC_A <: TierrianInstruction end
struct INC_B <: TierrianInstruction end
struct INC_C <: TierrianInstruction end
struct DEC_C <: TierrianInstruction end

struct PUSH_A <: TierrianInstruction end
struct PUSH_B <: TierrianInstruction end
struct PUSH_C <: TierrianInstruction end
struct PUSH_D <: TierrianInstruction end

struct POP_A <: TierrianInstruction end
struct POP_B <: TierrianInstruction end
struct POP_C <: TierrianInstruction end
struct POP_D <: TierrianInstruction end

struct JMP <: TierrianInstruction end
struct JMP_BACK <: TierrianInstruction end

struct CALL <: TierrianInstruction end
struct RET <: TierrianInstruction end

struct MOV_AB <: TierrianInstruction end
struct MOV_CD <: TierrianInstruction end
struct MOVI_BA <: TierrianInstruction end

struct ADR_A <: TierrianInstruction end
struct ADR_BACK_A <: TierrianInstruction end
struct ADR_FORW_A <: TierrianInstruction end

struct MALLOC <: TierrianInstruction end
struct DIVIDE <: TierrianInstruction end

const INT_TO_INS = Dict{UInt8, TierrianInstruction}(
    00 => NOP_0(),
    01 => NOP_1(),
    02 => OR_1(),
    03 => SH_L(),
    04 => ZERO(),
    05 => IF_0(),
    06 => SUB_AB(),
    07 => SUB_AC(),
    08 => INC_A(),
    09 => INC_B(),
    10 => INC_C(),
    11 => DEC_C(),
    12 => PUSH_A(),
    13 => PUSH_B(),
    14 => PUSH_C(),
    15 => PUSH_D(),
    16 => POP_A(),
    17 => POP_B(),
    18 => POP_C(),
    19 => POP_D(),
    20 => JMP(),
    21 => JMP_BACK(),
    22 => CALL(),
    23 => RET(),
    24 => MOV_AB(),
    25 => MOV_CD(),
    26 => MOVI_BA(),
    27 => ADR_A(),
    28 => ADR_BACK_A(),
    29 => ADR_FORW_A(),
    30 => MALLOC(),
    31 => DIVIDE()
)

const INS_TO_INT = Dict{TierrianInstruction, UInt8}(
    NOP_0() => 00,
    NOP_1() => 01,
    OR_1() => 02,
    SH_L() => 03,
    ZERO() => 04,
    IF_0() => 05,
    SUB_AB() => 06,
    SUB_AC() => 07,
    INC_A() => 08,
    INC_B() => 09,
    INC_C() => 10,
    DEC_C() => 11,
    PUSH_A() => 12,
    PUSH_B() => 13,
    PUSH_C() => 14,
    PUSH_D() => 15,
    POP_A() => 16,
    POP_B() => 17,
    POP_C() => 18,
    POP_D() => 19,
    JMP() => 20,
    JMP_BACK() => 21,
    CALL() => 22,
    RET() => 23,
    MOV_AB() => 24,
    MOV_CD() => 25,
    MOVI_BA() => 26,
    ADR_A() => 27,
    ADR_BACK_A() => 28,
    ADR_FORW_A() => 29,
    MALLOC() => 30,
    DIVIDE() => 31
)

convert_instruction(instruction::Integer) = INT_TO_INS[instruction % 32]
convert_instruction(instruction::TierrianInstruction) = INS_TO_INT[instruction]

convert_instructions(instructions::Vector{Integer}) = [INT_TO_INS[i] for i in instructions]
convert_instructions(instructions::Vector{TierrianInstruction}) = [INS_TO_INT[i] for i in instructions]

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module Quads
using Pkg
Pkg.add("Random")
Pkg.add("Combinatorics")
using Random
using Combinatorics

include("PokerCard.jl")
include("PokerDeck.jl")
include("Lookup.jl")
include("Evaluator.jl")

Pkg.resolve()

PokerCard.newCard("Ah")

export PokerCard
export PokerDeck
export Lookup
export Evaluatorus

end # module

ext: jl

max_stars_repo_name: HigginTown/Quads

max_issues_repo_name: HigginTown/Quads

max_forks_repo_name: HigginTown/Quads


text:
function trailing_edge_wing(settings::PyObject, ac::PyObject, freq::Array{Float64, 1}, M_0, c_0, rho_0, mu_0, theta, phi)

    ### ---------------- Wing trailing-edge noise ----------------
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 5
    delta_w_star = 0.37 * (ac.af_S_w / ac.af_b_w^2) * (rho_0 * M_0 * c_0 * ac.af_S_w / (mu_0 * ac.af_b_w))^(-0.2)

    # Determine configuration constant and the sound power
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 7
    if ac.af_clean_w == 0
        K_w = 4.464e-5
    elseif ac.af_clean_w == 1
        K_w = 7.075e-6
    end

    Pi_star_w = K_w * M_0^5 * delta_w_star

    # Determine directivity function
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 8
    D_w = 4. * cos(phi * π / 180.)^2 * cos(theta / 2 * π / 180.)^2

    # Determine spectral distribution function
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 10-11-12
    S_w = freq * delta_w_star * ac.af_b_w / (M_0 * c_0) * (1 - M_0 * cos(theta * π / 180.))
    if ac.af_delta_wing == 1
        F_w = 0.613 * (10 * S_w).^4 .* ((10 * S_w).^1.35 .+ 0.5).^(-4)
    elseif ac.af_delta_wing == 0
        F_w = 0.485 * (10 * S_w).^4 .* ((10 * S_w).^1.5 .+ 0.5).^(-4)
    end

    # Determine msap
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 1
    r_s_star_af = settings.r_0 / ac.af_b_w
    msap_w = 1 / (4 * π * r_s_star_af^2) / (1 - M_0 * cos(theta * π / 180.))^4 * (Pi_star_w * D_w * F_w)

    #Write output
    return msap_w
end
function trailing_edge_horizontal_tail(settings::PyObject, ac::PyObject, freq::Array{Float64, 1}, M_0, c_0, rho_0, mu_0, theta, phi)

    # ---------------- Horizontal tail trailing-edge noise ----------------
    # Trailing edge noise of the horizontal tail
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 5
    delta_h_star = 0.37 * (ac.af_S_h / ac.af_b_h^2) * (rho_0 * M_0 * c_0 * ac.af_S_h / (mu_0 * ac.af_b_h))^(-0.2)

    # Determine configuration constant and the sound power
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 7
    if ac.af_clean_h == 0
        K_h = 4.464e-5
    elseif ac.af_clean_h == 1
        K_h = 7.075e-6
    end
    Pi_star_h = K_h * M_0^5 * delta_h_star * (ac.af_b_h / ac.af_b_w)^2

    # Determine the directivity function
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 8
    D_h = 4 * cos(phi * π / 180.)^2 * cos(theta / 2 * π / 180.)^2

    # Determine spectral distribution function
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 10-11-12
    S_h = freq * delta_h_star * ac.af_b_h / (M_0 * c_0) * (1 - M_0 * cos(theta * π / 180.))
    F_h = 0.485 * (10 * S_h).^4 .* ((10 * S_h).^1.5 .+ 0.5).^(-4)

    # Determine msap
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 1
    r_s_star_af = settings.r_0 / ac.af_b_w
    msap_h = 1. / (4. * π * r_s_star_af^2) / (1 - M_0 * cos(theta * π / 180.))^4 * (Pi_star_h * D_h * F_h)

    return msap_h
end
function trailing_edge_vertical_tail(settings::PyObject, ac::PyObject, freq::Array{Float64, 1}, M_0, c_0, rho_0, mu_0, theta, phi)

    ### ---------------- Vertical tail trailing-edge noise ----------------
    delta_v_star = 0.37 * (ac.af_S_v / ac.af_b_v^2) * (rho_0 * M_0 * c_0 * ac.af_S_v / (mu_0 * ac.af_b_v))^(-0.2)

    # Trailing edge noise of the vertical tail
    if ac.af_clean_v == 0
        K_v = 4.464e-5
    elseif ac.af_clean_v == 1
        K_v = 7.075e-6
    end

    Pi_star_v = K_v * M_0^5 * delta_v_star * (ac.af_b_v / ac.af_b_w)^2

    # Determine directivity function
    D_v = 4 * sin(phi * π / 180.)^2 * cos(theta / 2 * π / 180.)^2

    # Determine spectral distribution function
    S_v = freq * delta_v_star * ac.af_b_v / (M_0 * c_0) * (1 - M_0 * cos(theta * π / 180.))

    if ac.af_delta_wing == 1
        F_v = 0.613 * (10 * S_v).^4 .* ((10 * S_v).^1.35 .+ 0.5).^(-4)
    elseif ac.af_delta_wing == 0
        F_v = 0.485 * (10 * S_v).^4 .* ((10 * S_v).^1.35 .+ 0.5).^(-4)
    end

    # Determine msap
    r_s_star_af = settings.r_0 / ac.af_b_w
    msap_v = 1. / (4 * π * r_s_star_af^2) / (1 - M_0 * cos(theta * π / 180.))^4 * (Pi_star_v * D_v * F_v)

    return msap_v
end
function leading_edge_slat(settings::PyObject, ac::PyObject, freq::Array{Float64, 1}, M_0, c_0, rho_0, mu_0, theta, phi)

    ### ---------------- Slat noise ----------------
    delta_w_star = 0.37 * (ac.af_S_w / ac.af_b_w^2) * (rho_0 * M_0 * c_0 * ac.af_S_w / (mu_0 * ac.af_b_w))^(-0.2)

    # Noise power
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 4
    Pi_star_les1 = 4.464e-5 * M_0^5 * delta_w_star  # Slat noise
    Pi_star_les2 = 4.464e-5 * M_0^5 * delta_w_star  # Added trailing edge noise

    # Determine the directivity function
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 8
    D_les = 4 * cos(phi * π / 180.)^2 * cos(theta / 2 * π / 180.)^2

    # Determine spectral distribution function
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 10-12-13
    S_les = freq * delta_w_star * ac.af_b_w / (M_0 * c_0) * (1 - M_0 * cos(theta * π / 180.))

    F_les1 = 0.613 * (10 * S_les).^4 .* ((10. * S_les).^1.5 .+ 0.5).^(-4)
    F_les2 = 0.613 * (2.19 * S_les).^4 .* ((2.19 * S_les).^1.5 .+ 0.5).^(-4)

    # Calculate msap
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 1
    r_s_star_af = settings.r_0 / ac.af_b_w
    msap_les = 1 / (4 * π * r_s_star_af^2) / (1 - M_0 * cos(theta * π / 180.))^4 * (Pi_star_les1 * D_les * F_les1 + Pi_star_les2 * D_les * F_les2)

    return msap_les
end
function trailing_edge_flap(settings::PyObject, ac::PyObject, freq::Array{Float64, 1}, M_0, c_0, theta, phi, theta_flaps)

    ### ---------------- Flap noise ----------------
    # Calculate noise power
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 14-15
    if ac.af_s < 3
        Pi_star_tef = 2.787e-4 * M_0^6 * ac.af_S_f / ac.af_b_w^2 * sin(theta_flaps * π / 180.)^2
    elseif ac.af_s == 3
        Pi_star_tef = 3.509e-4 * M_0^6 * ac.af_S_f / ac.af_b_w^2 * sin(theta_flaps * π / 180.)^2
    end

    # Calculation of the directivity function
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 16
    D_tef = 3 * (sin(theta_flaps * π / 180.) * cos(theta * π / 180.) + cos(theta_flaps * π / 180.) * sin(theta * π / 180.) * cos(phi * π / 180.))^2
    # Strouhal number
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 19
    S_tef = freq * ac.af_S_f / (M_0 * ac.af_b_f * c_0) * (1 - M_0 * cos(theta * π / 180.))

    # Calculation of the spectral function
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 17-18
    if ac.af_s < 3
        F_tef = zeros(eltype(S_tef), (settings.N_f, ))
        for i in range(1, settings.N_f, step=1)
            if S_tef[i] .< 2.
                F_tef[i] = 0.0480 * S_tef[i]
            elseif 2. .< S_tef[i] .< 20.
                F_tef[i] = 0.1406 * S_tef[i].^(-0.55)
            else
                F_tef[i] = 216.49 * S_tef[i].^(-3)
            end
        end
    elseif ac.af_s == 3
        F_tef = zeros(eltype(S_tef), (settings.N_f, ))
        for i in range(1, settings.N_f, step=1)
            if S_tef[i] .< 2.
                F_tef[i] = 0.0257 * S_tef[i]
            elseif 2. .< S_tef[i] .< 20.
                F_tef[i] = 0.0536 * S_tef[i].^(-0.0625)
            else
                F_tef[i] = 17078 * S_tef[i].^(-3)
            end
        end
    end

    # Calculate msap
    # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 1
    r_s_star_af = settings.r_0 / ac.af_b_w
    msap_tef = 1. / (4 * π * r_s_star_af^2) / (1 - M_0 * cos(theta * π / 180.))^4 * (Pi_star_tef * D_tef * F_tef)

    return msap_tef
end
function landing_gear(settings::PyObject, ac::PyObject, freq::Array{Float64, 1}, M_0, c_0, theta, phi, I_landing_gear)

    ### ---------------- Landing-gear noise ----------------
    if I_landing_gear == 1
        # Calculate nose-gear noise
        # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 29
        S_ng = freq * ac.af_d_ng / (M_0 * c_0) * (1 - M_0 * cos(theta * π / 180.))

        # Calculate noise power and spectral distribution function
        # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 20-21-22-25-26-27-28
        if ac.af_n_ng == 1 || ac.af_n_ng == 2
            Pi_star_ng_w = 4.349e-4 * M_0^6 * ac.af_n_ng * (ac.af_d_ng / ac.af_b_w)^2
            Pi_star_ng_s = 2.753e-4 * M_0^6 * (ac.af_d_ng / ac.af_b_w)^2 * (ac.af_l_ng / ac.af_d_ng)
            F_ng_w = 13.59 * S_ng.^2 .* (12.5 .+ S_ng.^2).^(-2.25)
            F_ng_s = 5.32 * S_ng.^2 .* (30 .+ S_ng.^8).^(-1)
        elseif ac.af_n_ng == 4
            Pi_star_ng_w = 3.414 - 4 * M_0^6 * ac.af_n_ng * (ac.af_d_ng / ac.af_b_w)^2
            Pi_star_ng_s = 2.753e-4 * M_0^6 * (ac.af_d_ng / ac.af_b_w)^2 * (ac.af_l_ng / ac.af_d_ng)
            F_ng_w = 0.0577 * S_ng.^2 .* (1 .+ 0.25 * S_ng.^2).^(-1.5)
            F_ng_s = 1.28 * S_ng.^3 .* (1.06 .+ S_ng.^2).^(-3)
        end

        # Calculate main-gear noise
        # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 29
        S_mg = freq * ac.af_d_mg / (M_0 * c_0) * (1 - M_0 * cos(theta * π / 180.))
        # Calculate noise power and spectral distribution function
        # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 20-21-22-25-26-27-28
        if ac.af_n_mg == 1 || ac.af_n_mg == 2
            Pi_star_mg_w = 4.349e-4 * M_0^6 * ac.af_n_mg * (ac.af_d_mg / ac.af_b_w)^2
            Pi_star_mg_s = 2.753e-4 * M_0^6 * (ac.af_d_mg / ac.af_b_w)^2 * (ac.af_l_ng / ac.af_d_mg)
            F_ng_w = 13.59 * S_mg.^2 .* (12.5 .+ S_mg.^2).^(-2.25)
            F_ng_s = 5.32 * S_mg.^2 .* (30 .+ S_mg.^8).^(-1)
        elseif ac.af_n_mg == 4
            Pi_star_mg_w = 3.414e-4 * M_0^6 * ac.af_n_mg * (ac.af_d_mg / ac.af_b_w)^2
            Pi_star_mg_s = 2.753e-4 * M_0^6 * (ac.af_d_mg / ac.af_b_w)^2 * (ac.af_l_ng / ac.af_d_mg)
            F_mg_w = 0.0577 * S_mg.^2 .* (1 .+ 0.25 * S_mg.^2).^(-1.5)
            F_mg_s = 1.28 * S_mg.^3 .* (1.06 .+ S_mg.^2).^(-3)
        end

        # Directivity function
        # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 23-24
        D_w = 1.5 * sin(theta * π / 180.)^2
        D_s = 3 * sin(theta * π / 180.)^2 * sin(phi * π / 180.)^2
        # Calculate msap
        # Source: Zorumski report 1982 part 2. Chapter 8.8 Equation 1
        # If landing gear is down
        r_s_star_af = settings.r_0 / ac.af_b_w
        msap_lg = 1 / (4 * π * r_s_star_af^2) / (1 - M_0 * cos(theta * π / 180.))^4 * (
                              ac.af_N_ng * (Pi_star_ng_w * F_ng_w * D_w + Pi_star_ng_s * F_ng_s * D_s) +
                              ac.af_N_mg * (Pi_star_mg_w * F_mg_w * D_w + Pi_star_mg_s * F_mg_s * D_s))

    # If landing gear is up
    else
        msap_lg = zeros(eltype(I_landing_gear), 1)
    end

    return msap_lg
end

function airframe(settings::PyObject, data::PyObject, ac::PyObject, n_t::Int64, M_0, mu_0, c_0, rho_0, theta, phi, theta_flaps, I_landing_gear)

    # Initialize solution
    T = eltype(theta_flaps)
    msap_af = zeros(T, (n_t, settings.N_f))

    for i in range(1, n_t, step=1)

        # Airframe noise only when aircraft is moving
        if M_0[i] != 0.

            # Apply HSR-era airframe suppression/calibration levels
            # Source: validation noise assessment data set of NASA STCA (Berton et al., 2019)
            if settings.hsr_calibration
                # Suppression data
                data_angles = data.supp_af_angles
                data_freq = data.supp_af_freq
                data_supp = data.supp_af
                f_supp = LinearInterpolation((data_freq, data_angles), data_supp)
                supp = f_supp.(data.f, theta[i]*ones(T, (settings.N_f, )))
            else
                supp = ones(T, (settings.N_f, ))
            end

            # Add airframe noise components
            msap_j = zeros(T, (settings.N_f, ))
            if "wing" in ac.comp_lst
                msap_w = trailing_edge_wing(settings, ac, data.f, M_0[i], c_0[i], rho_0[i], mu_0[i], theta[i], phi[i])
                msap_j = msap_j .+ msap_w .* supp
            end
            if "tail_v" in ac.comp_lst
                msap_v = trailing_edge_vertical_tail(settings, ac, data.f, M_0[i], c_0[i], rho_0[i], mu_0[i], theta[i], phi[i])
                msap_j = msap_j .+ msap_v .* supp
            end
            if "tail_h" in ac.comp_lst
                msap_h = trailing_edge_horizontal_tail(settings, ac, data.f, M_0[i], c_0[i], rho_0[i], mu_0[i], theta[i], phi[i])
                msap_j = msap_j .+ msap_h .* supp
            end
            if "les" in ac.comp_lst
                msap_les = leading_edge_slat(settings, ac, data.f, M_0[i], c_0[i], rho_0[i], mu_0[i], theta[i], phi[i])
                msap_j = msap_j .+ msap_les .* supp
            end
            if "tef" in ac.comp_lst
                msap_tef = trailing_edge_flap(settings, ac, data.f, M_0[i], c_0[i], theta[i], phi[i], theta_flaps[i])
                msap_j = msap_j .+ msap_tef .* supp
            end
            if "lg" in ac.comp_lst
                msap_lg = landing_gear(settings, ac, data.f, M_0[i], c_0[i], theta[i], phi[i], I_landing_gear[i])
                msap_j = msap_j .+ msap_lg .* supp
            end

            # Normalize msap by reference pressure
            msap_af[i,:] = msap_j/settings.p_ref^2

        else
            msap_af[i,:] =  1e-99 * ones(T, (settings.N_f, ))
        end
    end

    return msap_af
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
module FractalFlow

using StaticArrays   # SVector
using IterTools      # Imap
using LightXML       # Parse XML, VTU
using LightGraphs    #
using Combinatorics  # Combinations
using DataStructures # Deque
using Polynomials
using WriteVTK
using PyCall        # need networkx on python side

import Base: length, write, size, split, start, next, done, scale!
import LightGraphs: nv, ne, Graph, degree
# Topology
export Curve, VertexFunction, EdgeFunction, is_consistent, entity_dim
export nv, ne, Graph, degree_f, num_loops, branch_map, num_branches, num_colors, color_f
# Geometry
export FractalDim, fractal_dim, estimate_fd, min_branch, max_branch
export min_segment, max_segment
# Generation, fractals
export koch_t_curve, koch_q1_curve, koch_q2_curve
# Generation, random
export gmsh_code
# IO
export write

# -------------------------------------------------------------------

"""
Curve is a collection of segments. Each segment is defined by 2 indices
corresponding to two points in the points array.
"""
struct Curve{D, T}
    # NOTE: The inner preconditioner assumes sensible data that is points
    # in the vector are unique and the segments are unique as well. We do
    # no check if this is true here
    points::Vector{SVector{D, T}}
    segments::Vector{Tuple{Int, Int}}
end
geometrical_dim{D, T}(::Curve{D, T}) = D

"""|AB|"""
struct Segment{D, T}
    A::SVector{D, T}
    B::SVector{D, T}
end

# Constructors
include("curve.jl")

abstract type CurveFunction end

geometrical_dim(f::CurveFunction) = geometrical_dim(f.curve)

"""Values in nodes"""
struct VertexFunction{R <: Real} <: CurveFunction
    curve::Curve
    data::Vector{R}

    function VertexFunction{R}(c::Curve, d::Vector{R}) where R <: Real
        @assert length(c.points) == length(d)
        new(c, d)
    end
end
VertexFunction(c::Curve, d::Vector) = VertexFunction{eltype(d)}(c, d)

topological_dim(::VertexFunction) = 0

"""Values in segments"""
struct EdgeFunction{R <: Real} <: CurveFunction
    curve::Curve
    data::Vector{R}

    function EdgeFunction{R}(c::Curve, d::Vector{R}) where R <: Real
        @assert length(c.segments) == length(d)
        new(c, d)
    end
end
EdgeFunction(c::Curve, d::Vector) = EdgeFunction{eltype(d)}(c, d)

topological_dim(::EdgeFunction) = 1

# IO - let curves be 'readable' from files
include("vtk_io.jl")
include("xml_io.jl")
include("msh_io.jl")

"""Load from VTU or XML file"""
function Curve(path::AbstractString)
    name, ext = splitext(path)

    ext == ".vtu" && return Curve(read_curve_vtu(path)...)

    ext == ".xml" && return Curve(read_curve_xml(path)...)

    ext == ".msh" && return Curve(read_curve_msh(path)...)
end

# Various properties of curves
include("topology.jl")  # Connectivity

include("geometry.jl")  # Distance, fractal properties

# Generation
include("generation.jl")  # Fractals and random curves
include("gmsh_code.jl")  #

end # module


text:
function main()
	N = 600_851_475_143
	for i = isqrt(N):-1:1
		if 0 == N % i && isprime(i)
			println(i)
			break
		end
	end
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
#=

  Hamming distance in ConstraintSolver.jl

  I.e. the number of bits differing in two (binary) arrays.
  See http://en.wikipedia.org/wiki/Hamming_distance

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function hamming_distance(a,print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    n = length(a)
    println("a:$a n:$n")
    @variable(model, x[1:n], Bin)
    @variable(model, 0 <= diffs <= n, Int)

    # This don't work in the current version (v0.6.3)
    # @constraint(model, diffs == sum(x .!= a))
    # So we have to rewrite it...

    count_ctr2(model,x, :(!=), a, diffs)
    #=
    @variable(model, b[1:n], Bin)
    for i in 1:n
        @constraint(model, b[i] := {x[i] != a[i]})
    end
    @constraint(model, diffs == sum(b))
    =#

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                diffs_val = convert.(Integer,JuMP.value.(diffs; result=sol))
                println("x:$x_val diffs:$diffs_val")

            end
        end
    else
        println("status:$status")
    end
    println("num_sols:$num_sols\n")
    return status, num_sols
end

a = ones(Int,4)
@time hamming_distance(a)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
@osx_only begin
	using Homebrew
	Homebrew.add("terminal-notifier")
end


text:
using BitempWebApp, Test
# implement your tests here
push!(LOAD_PATH,"app/resources/insurancecontracts")
include("testsCreateContract.jl")

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
################################################################################
#
# (q/p)adic integers
#
# complete enough to support hnf
################################################################################
# CHECK precision!!!

struct QadicRing{T} <: Generic.Ring
  Q::T
end

function Base.show(io::IO, Q::QadicRing)
  println("Integers of ", Q.Q)
end

function MaximalOrder(Q::FlintQadicField)
  return QadicRing{FlintQadicField}(Q)
end
#integers(Q::FlintQadicField) = ring_of_integers(Q)

function MaximalOrder(Q::FlintPadicField)
  return QadicRing{FlintPadicField}(Q)
end
#integers(Q::FlintPadicField) = ring_of_integers(Q)

struct QadicRingElem{S} <: RingElem
  x::S
  P::QadicRing
  function QadicRingElem(a::qadic, P::QadicRing)
    r = new{qadic}(a, P)
  end
  function QadicRingElem(a::padic, P::QadicRing)
    r = new{padic}(a, P)
  end
end

function Base.show(io::IO, a::QadicRingElem)
  print(io, a.x)
end

*(a::QadicRingElem, b::QadicRingElem) = QadicRingElem(a.x*b.x, a.P)
+(a::QadicRingElem, b::QadicRingElem) = QadicRingElem(a.x+b.x, a.P)
-(a::QadicRingElem, b::QadicRingElem) = QadicRingElem(a.x-b.x, a.P)
-(a::QadicRingElem) = QadicRingElem(-a.x, a.P)
^(a::QadicRingElem, b::QadicRingElem) = QadicRingElem(a.x^b.x, a.P)
^(a::T, b::QadicRingElem{T}) where {T} = a^b.x

function inv(a::QadicRingElem)
  valuation(a.x) == 0 || error("non unit")
  return QadicRingElem(inv(a.x), a.P)
end

==(a::QadicRingElem, b::QadicRingElem) = a.x == b.x

function divexact(a::QadicRingElem, b::QadicRingElem)
  @assert !iszero(b.x)
  iszero(a) && return a
  valuation(a.x) >= valuation(b.x) || error("division not exact")
  return QadicRingElem(a.x//b.x, a.P)
end

function divrem(a::QadicRingElem, b::QadicRingElem)
  if valuation(a.x) < valuation(b.x)
    return setprecision(a.P(0), precision(a)), a
  end
  q = divexact(a, b)
  return q, a-q*b
end

function Base.div(a::QadicRingElem, b::QadicRingElem)
  if valuation(a.x) < valuation(b.x)
    return setprecision(a.P(0), precision(a))
  end
  q = divexact(a, b)
  return q
end

parent(a::QadicRingElem) = a.P
elem_type(::Type{QadicRing{FlintPadicField}}) = QadicRingElem{padic}
elem_type(::Type{QadicRing{FlintQadicField}}) = QadicRingElem{qadic}
parent_type(::Type{QadicRingElem{padic}}) = QadicRing{FlintPadicField}
parent_type(::Type{QadicRingElem{qadic}}) = QadicRing{FlintQadicField}

zero(Q::QadicRing) = QadicRingElem(Q.Q(0), Q)
one(Q::QadicRing) = QadicRingElem(Q.Q(1), Q)

(Q::QadicRing)(a::qadic) = QadicRingElem(a, Q)
(Q::QadicRing)(a::padic) = QadicRingElem(a, Q)
(Q::QadicRing)(a::QadicRingElem) = QadicRingElem(a.x, a.P)
(Q::QadicRing)(a::Int) = QadicRingElem(Q.Q(a), Q)
(Q::QadicRing)() = QadicRingElem(Q.Q(), Q)
(Q::FlintQadicField)(a::QadicRingElem{qadic}) = a.x
(Q::FlintPadicField)(a::QadicRingElem{padic}) = a.x
valuation(a::QadicRingElem) = valuation(a.x)
isunit(a::QadicRingElem) = valuation(a) == 0
(Q::FlintQadicField)(a::padic) =  _map(Q, a) #TODO: do properly

function _map(Q::FlintQadicField, a::padic)
  K = parent(a)
  v = valuation(a)
  if v >= 0
    q = Q(lift(a))
    return setprecision(q, a.N)
  else
    d = uniformizer(K)^-v
    n = a*d
    n1 = divexact(Q(lift(n)), Q(lift(d)))
    return n1
  end
end

function Base.deepcopy_internal(a::QadicRingElem, dict::IdDict)
  return QadicRingElem(a.x, a.P)
end

function canonical_unit(a::QadicRingElem)
  iszero(a.x) && return setprecision(a.P(1), precision(a))
  v = valuation(a.x)
  return QadicRingElem(inv(a.x//prime(a.P.Q)^v), a.P)
end

function gcdx(a::QadicRingElem, b::QadicRingElem)
  if iszero(a)
    c = canonical_unit(b)
    return b*c, a, c
  end
  if iszero(b)
    c = canonical_unit(a)
    return a*c, c, b
  end
  if valuation(a.x) < valuation(b.x)
    c = canonical_unit(a)
    return a*c, c, setprecision(a.P(0), precision(a))
  else
    c = canonical_unit(b)
    return b*c, setprecision(b.P(0), precision(b)), c
  end
end

function mul_red!(a::QadicRingElem, b::QadicRingElem, c::QadicRingElem, f::Bool = false)
  return b*c
end

function mul!(a::QadicRingElem, b::QadicRingElem, c::QadicRingElem)
  return b*c
end

function add!(a::QadicRingElem, b::QadicRingElem, c::QadicRingElem)
  return b+c
end

function addeq!(a::QadicRingElem, b::QadicRingElem)
  return a+b
end

Base.iszero(a::QadicRingElem) = iszero(a.x)
Base.isone(a::QadicRingElem) = isone(a.x)

Base.precision(Q::QadicRing) = precision(Q.Q)
Base.precision(a::QadicRingElem) = precision(a.x)

function setprecision!(Q::QadicRing, n::Int)
  setprecision!(Q.Q, n)
end

function Base.setprecision(a::QadicRingElem, n::Int)
  return a.P(setprecision(a.x, n))
end

function setprecision!(a::QadicRingElem, n::Int)
  setprecision!(a.x, n)
end

function Base.setprecision(a::Generic.MatSpaceElem{QadicRingElem{qadic}}, n::Int)
  return matrix(map(x -> setprecision(x, n), a.entries))
end

coefficient_ring(Q::QadicRing) = integers(coefficient_ring(Q.Q))


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

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

"""

Gets a list of resource navigation links for a subnet.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: subnetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ResourceNavigationLinksListResult
"""
function _swaggerinternal_resourceNavigationLinksList(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ResourceNavigationLinksListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/subnets/{subnetName}/ResourceNavigationLinks", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subnetName", subnetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function resourceNavigationLinksList(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_resourceNavigationLinksList(_api, resourceGroupName, virtualNetworkName, subnetName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function resourceNavigationLinksList(_api::DefaultApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, subnetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_resourceNavigationLinksList(_api, resourceGroupName, virtualNetworkName, subnetName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets a list of service association links for a subnet.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: subnetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ServiceAssociationLinksListResult
"""
function _swaggerinternal_serviceAssociationLinksList(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ServiceAssociationLinksListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/subnets/{subnetName}/ServiceAssociationLinks", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subnetName", subnetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function serviceAssociationLinksList(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_serviceAssociationLinksList(_api, resourceGroupName, virtualNetworkName, subnetName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function serviceAssociationLinksList(_api::DefaultApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, subnetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_serviceAssociationLinksList(_api, resourceGroupName, virtualNetworkName, subnetName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Prepares a subnet by applying network intent policies.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: subnetName::String (required)
Param: prepareNetworkPoliciesRequestParameters::PrepareNetworkPoliciesRequest (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_subnetsPrepareNetworkPolicies(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, prepareNetworkPoliciesRequestParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/subnets/{subnetName}/PrepareNetworkPolicies", ["azure_auth"], prepareNetworkPoliciesRequestParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subnetName", subnetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function subnetsPrepareNetworkPolicies(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, prepareNetworkPoliciesRequestParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_subnetsPrepareNetworkPolicies(_api, resourceGroupName, virtualNetworkName, subnetName, prepareNetworkPoliciesRequestParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function subnetsPrepareNetworkPolicies(_api::DefaultApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, subnetName::String, prepareNetworkPoliciesRequestParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_subnetsPrepareNetworkPolicies(_api, resourceGroupName, virtualNetworkName, subnetName, prepareNetworkPoliciesRequestParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Unprepares a subnet by removing network intent policies.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: subnetName::String (required)
Param: unprepareNetworkPoliciesRequestParameters::UnprepareNetworkPoliciesRequest (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_subnetsUnprepareNetworkPolicies(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, unprepareNetworkPoliciesRequestParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/subnets/{subnetName}/UnprepareNetworkPolicies", ["azure_auth"], unprepareNetworkPoliciesRequestParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subnetName", subnetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function subnetsUnprepareNetworkPolicies(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, subnetName::String, unprepareNetworkPoliciesRequestParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_subnetsUnprepareNetworkPolicies(_api, resourceGroupName, virtualNetworkName, subnetName, unprepareNetworkPoliciesRequestParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function subnetsUnprepareNetworkPolicies(_api::DefaultApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, subnetName::String, unprepareNetworkPoliciesRequestParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_subnetsUnprepareNetworkPolicies(_api, resourceGroupName, virtualNetworkName, subnetName, unprepareNetworkPoliciesRequestParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Checks whether a private IP address is available for use.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: ipAddress::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: IPAddressAvailabilityResult
"""
function _swaggerinternal_virtualNetworksCheckIPAddressAvailability(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, ipAddress::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IPAddressAvailabilityResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/CheckIPAddressAvailability", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "ipAddress", ipAddress)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworksCheckIPAddressAvailability(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, ipAddress::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksCheckIPAddressAvailability(_api, resourceGroupName, virtualNetworkName, ipAddress, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworksCheckIPAddressAvailability(_api::DefaultApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, ipAddress::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksCheckIPAddressAvailability(_api, resourceGroupName, virtualNetworkName, ipAddress, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists usage stats.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualNetworkListUsageResult
"""
function _swaggerinternal_virtualNetworksListUsage(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualNetworkListUsageResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/usages", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworksListUsage(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksListUsage(_api, resourceGroupName, virtualNetworkName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworksListUsage(_api::DefaultApi, response_stream::Channel, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworksListUsage(_api, resourceGroupName, virtualNetworkName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export resourceNavigationLinksList, serviceAssociationLinksList, subnetsPrepareNetworkPolicies, subnetsUnprepareNetworkPolicies, virtualNetworksCheckIPAddressAvailability, virtualNetworksListUsage


text:
module gan
    mutable struct GAN
        layers::Array{Any}
        initialize::Any
        activate::Any
        activate_Generator::Any
        activate_Discriminator::Any
        update::Any

        noise_generator::Any

        add_Glayer::Any
        add_Dlayer::Any
        num_Glayer::Int64
        num_Dlayer::Int64

        loss::Float64

        default_input_shape::Tuple

        G_range::UnitRange{Int64}
        D_range::UnitRange{Int64}

        function GAN(noise_generator)
            new(Any[Hidden_Input_Layer(Float32[])], init_GAN, activate_GAN, activate_Generator, activate_Discriminator, update_GAN, noise_generator, GAN_add_Glayer, GAN_add_Dlayer, 0, 0, 0.0, ())
        end
    end

    mutable struct Hidden_Input_Layer
        output::Array{Float32}
    end
    mutable struct Hidden_Output_Layer
        δ::Array{Float32}
    end

    function GAN_add_Glayer(model::GAN, layer::Any;args...)
        push!(model.layers, layer(;input_shape=model.default_input_shape, args...))
        model.num_Glayer += 1
        model.default_input_shape = model.layers[end].output_shape
    end
    function GAN_add_Dlayer(model::GAN, layer::Any;args...)
        push!(model.layers, layer(;input_shape=model.default_input_shape, args...))
        model.num_Dlayer += 1
        model.default_input_shape = model.layers[end].output_shape
    end

    function init_GAN(model::GAN, mini_batch::Int64)
        if length(model.layers) == model.num_Glayer+model.num_Dlayer+1
            push!(model.layers, Hidden_Output_Layer(zeros(Float32, model.layers[end].output_shape..., mini_batch)))
        end
        for i in 2:length(model.layers)-1
            model.layers[i].initialize(model.layers[i], mini_batch)
        end
        model.G_range = 2:model.num_Glayer+1
        model.D_range = model.num_Glayer+2:model.num_Glayer+model.num_Dlayer+1
    end

    function activate_Generator(model::GAN)
        model.layers[1].output = model.noise_generator()
        for i in model.G_range
            model.layers[i].activate(model.layers[i], model.layers[i-1].output)
        end
        return model.layers[model.num_Glayer+1].output
    end
    function activate_Discriminator(model::GAN, data::Array{Float32})
        model.layers[1+model.num_Glayer].output = data
        for i in model.D_range
            model.layers[i].activate(model.layers[i], model.layers[i-1].output)
        end
    end
    function activate_GAN(model::GAN)
        model.layers[1].output = model.noise_generator()
        for i in 2:length(model.layers)-1
            model.layers[i].activate(model.layers[i], model.layers[i-1].output)
        end
    end

    function update_GAN(model::GAN, current_input_data::Array{Float32}, current_output_data::Array{Float32}, loss_function::Any, monitor::Any, optimizer::String, α::Float64, parameters...)
        model.activate_Discriminator(model, current_input_data)
        loss_function.prop!(model.layers[end].δ, model.layers[end-1].output, reshape([0.0f0; 1.0f0], (2,1)))
        for i in Iterators.reverse(model.D_range)
            model.layers[i].update(model.layers[i], optimizer, model.layers[i-1].output, model.layers[i+1].δ, α, parameters)
        end
        model.loss += monitor.func(model.layers[end-1].output, current_output_data)

        model.activate(model)
        loss_function.prop!(model.layers[end].δ, model.layers[end-1].output, reshape([0.0f0; 1.0f0], (2,1)))
        for i in Iterators.reverse(model.D_range)
            model.layers[i].update(model.layers[i], optimizer, model.layers[i-1].output, model.layers[i+1].δ, α, parameters, -1)
        end
        for i in Iterators.reverse(model.G_range)
            model.layers[i].update(model.layers[i], optimizer, model.layers[i-1].output, model.layers[i+1].δ, α, parameters)
        end
    end
end


text:
function metropolisBase(curr::Real, log_prob::Function, stepSD::Real)
  cand = curr + randn() * stepSD
  logU = log(rand())
  logP = log_prob(cand) - log_prob(curr)
  accept = logP > logU
  draw = accept ? cand : curr
  return (draw, accept)
end

function metropolis(curr::Real, log_prob::Function, stepSD::Real)
  return metropolisBase(curr, log_prob, stepSD)[1]
end

function metropolis(curr::AbstractVector{<:Real}, log_prob::Function,
                    propCov::AbstractMatrix{<:Real})
  cand = rand(MvNormal(curr, propCov))
  p = log_prob(cand) - log_prob(curr)
  return p > log(rand()) ? cand : curr
end

# TODO: Check.
# See section 2 of http://probability.ca/jeff/ftpdir/adaptex.pdf
function metropolisAdaptive(curr::AbstractVector{<:Real}, log_prob::Function,
                            tuner::MvTuner)
  # Update summary stats.
  update!(curr, tuner)

  # Proposal covariance.
  propCov = if tuner.iter <= 2*tuner.d || tuner.beta > rand()
    0.01 * eye(tuner.d) / tuner.d
  else
    5.6644 * tuner.current_cov / tuner.d
  end

  # Metropolis step
  return metropolis(curr, log_prob, propCov)
end

"""
Adaptive metropolis (within Gibbs). See section 3 of the paper below:
  http://probability.ca/jeff/ftpdir/adaptex.pdf

Another useful website:
  https://m-clark.github.io/docs/ld_mcmc/index_onepage.html
"""
function metropolisAdaptive(curr::Real, log_prob::Function,
                            tuner::TuningParam;
                            update::Function=update_tuning_param_default)
  draw, accept = metropolisBase(curr, log_prob, tuner.value)
  update(tuner, accept)
  return draw
end

"""
Log absolute jacobian term to be added when taking the log of a postiviely
supported random variable (x), so as to produce a random variable with
support on the real line. When added to the target log density, the resulting
expression can be used in a metropolis step conveniently.
"""
function logabsjacobian_logx(real_x::Real)
  return real_x
end

"""
Log absolute jacobian term to be added when taking the logit of a
lower-and-upper-bounded random variable (x), so as to produce a random variable
with support on the real line. When added to the target log density, the
resulting expression can be used in a metropolis step conveniently.

# Arguments
- `real_x::Float64`: the unconstrained value of the positively-supported
                     parameter x (i.e. log x)
- `a::Float64`: lower bound for unconstrained parameter (default: 0)
- `b::Float64`: upper bound for unconstrained parameter (default: 1)
"""
function logabsjacobian_logitx(real_x::Real)
  return logpdf(Logistic(), real_x)
end

function metLogAdaptive(curr::Real, log_prob::Function, tuner::TuningParam;
                        update::Function=update_tuning_param_default)

  function log_prob_plus_logabsjacobian(log_x::Real)
    x = exp(log_x)
    return log_prob(x) + logabsjacobian_logx(log_x)
  end

  log_x = metropolisAdaptive(log(curr), log_prob_plus_logabsjacobian, tuner,
                             update=update)

  return exp(log_x)
end

function metLogitAdaptive(curr::Real, log_prob::Function,
                          tuner::TuningParam;
                          update::Function=update_tuning_param_default)

  function log_prob_plus_logabsjacobian(logit_x::Real)
    x = logistic(logit_x)
    return log_prob(x) + logabsjacobian_logitx(logit_x)
  end

  logit_x = metropolisAdaptive(logit(curr),
                               log_prob_plus_logabsjacobian,
                               tuner, update=update)

  return logistic(logit_x)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
if !@isdefined(CGS)
  include("../../../src/CGS.jl")
  using Main.CGS
end

using JLD2; using PyPlot
using DelimitedFiles
using Statistics

#rhostar = 51.14; sig_rho=0.67 # In units of solar density
#rhostar = 53.7; sig_rho=0.85 # In units of solar density
#rhostar = 53.308; sig_rho=0.0208 # In units of solar density
#rhostar = 51.14; sig_rho=0.001 # In units of solar density

# Artificially increase density of star to see what it
# would take to make an Earth-like density:
#rhostar *= 1.23

# Mass of star from Mann et al. (2019):  https://github.com/awmann/M_-M_K-
mstar = 0.0898; smstar = 0.0023  # Solar units
#mstar = 0.0898; smstar = 0.0001  # Solar units

# Latest masses of planets assuming a stellar mass of 0.09 M_sun.
# Take from plot_masses.jl in
# /Users/ericagol/Observing/Spitzer/Cycle13/TRAPPIST1/campaign03/April2019
# Includes delta-omega prior, includes mixture model for handling outliers.

#mp= [1.37024,1.31417,0.364033,0.652497,0.997529,1.28074,0.340946]
#sp= [0.0848279,0.0602858,0.022669,0.0230492,0.0202153,0.022538,0.0355592]
#mp = [1.3785, 1.3201, 0.3911, 0.7042, 1.0463, 1.3237, 0.3283]
#sp = [0.0530, 0.0397, 0.0076, 0.0116, 0.0128, 0.0146, 0.0190]

#mp = [1.37927,1.31014,0.388479,0.693422,1.04092,1.32341,0.325985]
#sp = [0.0596173,0.0465207,0.00748253,0.0127809,0.0154444,0.0168958,0.0185088]

# Load in the posterior from the transit-timing analysis:
#@load "/Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v09/Hyak/T1_hmc_total_02132020.jld2" state_total
@load "../../../data/T1_hmc_total_02212020.jld2" state_total

# Load in the posterior from the photodynamic analysis:
@load "../../../data/T1_pd_MCMC_run_000_001.jld2" flatchain
flatchain1 = flatchain[:,35001:200000]
#@load "/Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v06/Hyak/T1_pd_MCMC_run_000_002.jld2" flatchain
@load "../../../data/T1_pd_MCMC_run_000_002.jld2" flatchain
flatchain12 = hcat(flatchain1,flatchain[:,35001:200000])
#@load "/Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v06/Hyak/T1_pd_MCMC_run_000_003.jld2" flatchain
@load "../../../data/T1_pd_MCMC_run_000_003.jld2" flatchain
flatchain = hcat(flatchain12,flatchain[:,35001:200000])

# Make an array to hold posterior of mass-ratios:
nplanet = 7
mass_ratio_total = zeros(nplanet,size(state_total)[2])
for ip=1:nplanet
  mass_ratio_total[ip,:] .= state_total[(ip-1)*5+1,:]
end

periods = [1.51087637,2.42180746,4.049959,6.099043,9.205585,12.354473,18.767953]
t0 = [7322.51654,7282.80879,7670.14227,7660.37910,7671.39470,7665.35084,7662.55467]

logplot = false
include("mass_radius2.jl")

#merror = mass_radius2(mass_ratio_total,mstar,smstar,flatchain,false)

scatter([1.0,0.815],[1.0,0.95])

# Plot new Dorn et al. models:
dmr2 = readdlm("../../../data/MR_trappist_Solar.ddat")
plot(dmr2[:,1],dmr2[:,2],label="Earth-like (Dorn et al. 2018)",linestyle="--")

legend(loc = "lower right",fontsize=10)

tight_layout()
#savefig("T1_mass_radius_Earth_carbon.pdf",bbox_inches="tight")
#read(stdin,Char)

## Semi-amplitudes:
fac = (2*pi*GRAV/(24*3600)/MSUN^2)^(1//3)*MEARTH
Kp = fac.*mavg./(mstar.^2 .*periods).^(1//3)
# Plot RV amplitude:
duration = 300
nrv = duration * 24
RV = zeros(nplanet+1,nrv)
time = zeros(nrv)
t0 .-= 7600.0
for i=1:nrv
  time[i] = i/24.0
  for j=1:nplanet
    RV[j,i] = Kp[j].*sin.((time[i].-t0[j])./periods[j])
  end
  RV[nplanet+1,i] = sum(RV[1:nplanet,i])
#    RV[j,i] = Kp[j].*sin.((time[i].-t0[j])./periods[j]))
end

clf()
plot(time,RV[nplanet+1,:]./100,label="Sum",color="k")
#plot(time,RV[nplanet+1,:],label="Sum",color="k")
xlabel("Time [days]")
#ylabel("RV [cm/sec]")
ylabel("RV [m/sec]")
#axis([0,duration,-1200,1200])
axis([0,duration,-12.00,12.00])
for i=1:7
#  plot(time,RV[i,:],label=planet_name[i],color=cname[i])
  plot(time,RV[i,:]./100,label=planet_name[i],color=cname[i],alpha=0.5)
#  errorbar([120+10*i],[-1000],Kp[i]*merror[i]/mavg[i],color=cname[i])
  errorbar([100+10*i],[-11],Kp[i]*merror[i]/mavg[i]/100,color=cname[i])
end
#text(40,-1020,"Equivalent RV precision:",fontsize=10)
text(5,-11.20,"Equivalent RV precision:",fontsize=10)
text(180,-11.2,"Measured RV precision: ")
errorbar([280],[-9.0],2.5,label="Hirano et al. (2020)",fmt="o")
legend(ncol=5)
tight_layout()
savefig("../Equivalent_RV_precision.pdf",bbox_inches="tight")

println("Semi-amplitudes: ",Kp," cm/s")
println("Precision: ",merror./mavg.*Kp," cm/s")
Kpsun = fac/(365.25)^(1//3)
println("Earth-Sun amplitude: ",Kpsun)
println("Equivalent Solar precision: ",Kpsun.*merror)

#include("/Users/ericagol/JWST/Exomoon/plot_moons.jl")

#@save "T1_mass_radius_cmf_000.jld2" nsamp msamp rsamp cmf rho_samp mstar_samp rstar_samp nr nm mrgrid rgrid mgrid level_mr conf m1 m2 r1 r2 mavg ravg merror rerror


text:
using Statistics: mean

"""
    LOF(k = 5,
        metric = Euclidean(),
        algorithm = :kdtree,
        leafsize = 10,
        reorder = true,
        parallel = false)

Calculate an anomaly score based on the density of an instance in comparison to its neighbors. This algorithm introduced
the notion of local outliers and was developed by Breunig et al., see [1].

Parameters
----------
$_k_param

$_knn_shared

Examples
--------
$(_score_unsupervised("LOF"))

References
----------
[1] Breunig, Markus M.; Kriegel, Hans-Peter; Ng, Raymond T.; Sander, Jörg (2000): LOF: Identifying Density-Based Local
Outliers.
"""
MMI.@mlj_model mutable struct LOF <: UnsupervisedDetector
    k::Integer = 5::(_ > 0)
    metric::DI.Metric = DI.Euclidean()
    algorithm::Symbol = :kdtree::(_ in (:kdtree, :balltree))
    leafsize::Integer = 10::(_ ≥ 0)
    reorder::Bool = true
    parallel::Bool = false
end

struct LOFModel <: Model
    # For efficient prediction, we need to store the learned tree, the distances of each training sample to its
    # k-nearest neighbors, as well as the training lrds.
    tree::NN.NNTree
    ndists::AbstractMatrix
    lrds::AbstractVector
end

function fit(detector::LOF, X::Data)::Fit
    # create the specified tree
    tree = buildTree(X, detector.metric, detector.algorithm, detector.leafsize, detector.reorder)

    # use tree to calculate distances
    idxs, dists = knn_others(tree, X, detector.k)

    # transform dists (vec of vec) to matrix to allow faster indexing later
    ndists = reduce(hcat, dists)

    # pre calculate lrds for later prediction use
    lrds = _calculate_lrd(idxs, dists, ndists, detector.k)

    # reduce distances to outlier score
    scores = _lof_from_lrd(idxs, lrds)

    Fit(LOFModel(tree, ndists, lrds), scores)
end

@score function score(detector::LOF, model::Fit, X::Data)::Result
    if detector.parallel
        idxs, dists = knn_parallel(model.tree, X, detector.k, true)
        return _lof(idxs, dists, model.ndists, model.lrds, detector.k)
    else
        idxs, dists = NN.knn(model.tree, X, detector.k, true)
        return _lof(idxs, dists, model.ndists, model.lrds, detector.k)
    end
end

function _lof(idxs:: AbstractVector, dists:: AbstractVector, model_dists::AbstractArray,
    model_lrds::AbstractVector, k::Int)::Scores
    lrds = _calculate_lrd(idxs, dists, model_dists, k)
    # calculate the local outlier factor from the lrds
    map((idx, lrd) -> mean(model_lrds[idx]) / lrd, idxs, lrds)
end

function _lof_from_lrd(idxs::AbstractVector, lrds::AbstractVector)::Scores
    # Directly calculate the local outlier factor from idxs with corresponding lrds.
    map((idx, lrd) -> mean(lrds[idx]) / lrd, idxs, lrds)
end

function _calculate_lrd(idxs::AbstractVector, dists::AbstractVector, ndists::AbstractArray, k::Int)::AbstractVector
    # The LRD of a sample is the inverse of the average reachability distance of its k-nearest neighbors. Epsilon is
    # added in case that there are more than k duplicates.
    map((is, ds) -> 1 / (mean(_max!(ndists[k, is], ds)) + 1e-10), idxs, dists)
end

function _max!(ar1:: AbstractVector, ar2:: AbstractVector)::AbstractVector
    # Calculate the element wise maximum between two vectors.
    for i in eachindex(ar1)
        @inbounds ar1[i] = ar1[i] > ar2[i] ? ar1[i] : ar2[i]
    end
    ar1
end


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

mutable struct DiskSku <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    tier::Any # spec type: Union{ Nothing, String } # spec name: tier

    function DiskSku(;name=nothing, tier=nothing)
        o = new()
        validate_property(DiskSku, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(DiskSku, Symbol("tier"), tier)
        setfield!(o, Symbol("tier"), tier)
        o
    end
end # type DiskSku

const _property_map_DiskSku = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("tier")=>Symbol("tier"))
const _property_types_DiskSku = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("tier")=>"String")
Base.propertynames(::Type{ DiskSku }) = collect(keys(_property_map_DiskSku))
Swagger.property_type(::Type{ DiskSku }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DiskSku[name]))}
Swagger.field_name(::Type{ DiskSku }, property_name::Symbol) =  _property_map_DiskSku[property_name]

const _allowed_DiskSku_name = ["Standard_LRS", "Premium_LRS", "StandardSSD_LRS", "UltraSSD_LRS"]

function check_required(o::DiskSku)
    true
end

function validate_property(::Type{ DiskSku }, name::Symbol, val)
    if name === Symbol("name")
        Swagger.validate_param(name, "DiskSku", :enum, val, _allowed_DiskSku_name)
    end
end


text:
# Autogenerated wrapper script for Cgl_jll for i686-linux-musl-cxx11
export libCgl

using Clp_jll
using Osi_jll
using CoinUtils_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Cgl")
JLLWrappers.@declare_library_product(libCgl, "libCgl.so.1")
function __init__()
    JLLWrappers.@generate_init_header(Clp_jll, Osi_jll, CoinUtils_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libCgl,
        "lib/libCgl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function dematel(comparisonMat::DataFrame; threshold=nothing)

    n, _ = size(comparisonMat)

    csums =  colsums(comparisonMat)
    rsums =  rowsums(comparisonMat)

    largest = maximum(vcat(csums, rsums))

    K = convert(Array{Float64,2}, comparisonMat)
    ND = K ./ largest

    T = ND * inv(I(n) - ND)

    c = rowsums(T)
    r = colsums(T)

    c_plus_r  = c .+ r
    c_minus_r = c .- r

    if threshold === nothing
        threshold = sum(T) / length(T)
    end

    influence_matrix = zeros(Float64, n, n)
    for i in 1:n
        for j in 1:n
            if T[i, j] > threshold
                influence_matrix[i, j] = 1
            end
        end
    end

    w_row = sqrt.((c_plus_r.^2.0) .+ (c_minus_r.^2.0))
    w = w_row ./ sum(w_row)

    result = DematelResult(
        comparisonMat,
        ND,
        T,
        c,
        r,
        c_plus_r,
        c_minus_r,
        threshold,
        influence_matrix,
        w
    )

    return result
end

max_stars_repo_name: ersagunkuruca/JMcDM

max_issues_repo_name: ersagunkuruca/JMcDM

max_forks_repo_name: ersagunkuruca/JMcDM


text:
#enum cusparseStatus_t
#error messages from CUSPARSE

typealias cusparseStatus_t UInt32
const CUSPARSE_STATUS_SUCCESS                   = 0
const CUSPARSE_STATUS_NOT_INITIALIZED           = 1
const CUSPARSE_STATUS_ALLOC_FAILED              = 2
const CUSPARSE_STATUS_INVALID_VALUE             = 3
const CUSPARSE_STATUS_ARCH_MISMATCH             = 4
const CUSPARSE_STATUS_MAPPING_ERROR             = 5
const CUSPARSE_STATUS_EXECUTION_FAILED          = 6
const CUSPARSE_STATUS_INTERNAL_ERROR            = 7
const CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED = 8

#enum cusparseAction_t
# perform operation only indices only or
# on both data and indices
typealias cusparseAction_t UInt32
const CUSPARSE_ACTION_SYMBOLIC = 0
const CUSPARSE_ACTION_NUMERIC  = 1

#enum cusparseDirection_t
# parse dense matrix by rows or cols
# to compute its number of non-zeros
typealias cusparseDirection_t UInt32
const CUSPARSE_DIRECTION_ROW = 0
const CUSPARSE_DIRECTION_COL = 1

#enum cusparseHybPartition_t
# how to partition the HYB matrix
typealias cusparseHybPartition_t UInt32
const CUSPARSE_HYB_PARTITION_AUTO = 0
const CUSPARSE_HYB_PARTITION_USER = 1
const CUSPARSE_HYB_PARTITION_MAX  = 2

#enum cusparseFillMode_t
# filling for HE,TR,SY matrices
typealias cusparseFillMode_t UInt32
const CUSPARSE_FILL_MODE_LOWER = 0
const CUSPARSE_FILL_MODE_UPPER = 1

#enum cusparseDiagType_t
# is the diagonal all ones
typealias cusparseDiagType_t UInt32
const CUSPARSE_DIAG_TYPE_NON_UNIT = 0
const CUSPARSE_DIAG_TYPE_UNIT     = 1

#enum cusparsePointerMode_t
# are scalars on the host or device
typealias cusparsePointerMode_t UInt32
const CUSPARSE_POINTER_MODE_HOST   = 0
const CUSPARSE_POINTER_MODE_DEVICE = 1

#enum cusparseOperation_t
# no op, transpose, conj transpose
typealias cusparseOperation_t UInt32
const CUSPARSE_OPERATION_NON_TRANSPOSE       = 0
const CUSPARSE_OPERATION_TRANSPOSE           = 1
const CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE = 2

#enum cusparseMatrixType_t
# marker to flag HE,TR,SY matrices
typealias cusparseMatrixType_t UInt32
const CUSPARSE_MATRIX_TYPE_GENERAL    = 0
const CUSPARSE_MATRIX_TYPE_SYMMETRIC  = 1
const CUSPARSE_MATRIX_TYPE_HERMITIAN  = 2
const CUSPARSE_MATRIX_TYPE_TRIANGULAR = 3

#enum cusparseSolvePolicy_t
# keep level info in solvers
typealias cusparseSolvePolicy_t UInt32
const CUSPARSE_SOLVE_POLICY_NO_LEVEL  = 0
const CUSPARSE_SOLVE_POLICY_USE_LEVEL = 1

#enum cusparseIndexBase_t
# 0- or 1-based indexing
typealias cusparseIndexBase_t UInt32
const CUSPARSE_INDEX_BASE_ZERO = 0
const CUSPARSE_INDEX_BASE_ONE  = 1

#struct cusparseMatDescr_t
# Describes shape and properties
# of a CUSPARSE matrix
type cusparseMatDescr_t
    MatrixType::cusparseMatrixType_t
    FillMode::cusparseFillMode_t
    DiagType::cusparseDiagType_t
    IndexBase::cusparseIndexBase_t
    function cusparseMatDescr_t(MatrixType,FillMode,DiagType,IndexBase)
        new(MatrixType,FillMode,DiagType,IndexBase)
    end
end

typealias cusparseSolveAnalysisInfo_t Ptr{Void}
typealias bsrsm2Info_t Ptr{Void}
typealias bsrsv2Info_t Ptr{Void}
typealias csrsv2Info_t Ptr{Void}
typealias csric02Info_t Ptr{Void}
typealias csrilu02Info_t Ptr{Void}
typealias bsric02Info_t Ptr{Void}
typealias bsrilu02Info_t Ptr{Void}

typealias cusparseContext Void
typealias cusparseHandle_t Ptr{cusparseContext}

#complex numbers

typealias cuComplex Complex{Float32}
typealias cuDoubleComplex Complex{Float64}

typealias CusparseFloat Union{Float64,Float32,Complex128,Complex64}
typealias CusparseReal Union{Float64,Float32}
typealias CusparseComplex Union{Complex128,Complex64}


text:
@testset "Initialization" begin

    @testset "Split Variables" begin
        st = build_scen_tree()
        js = create_model(PH.scid(0))
        vid_name_map = PH.report_variable_info(js, st)
        (branch_ids, leaf_ids) = PH._split_variables(st, collect(keys(vid_name_map)))

        for br_vid in branch_ids
            vname = vid_name_map[br_vid].name
            @test (vname == "y") || occursin("x", vname)
            if vname == "y"
                @test br_vid.stage == PH.stid(2)
            else
                @test br_vid.stage == PH.stid(1)
            end
        end

        for lf_vid in leaf_ids
            vname = vid_name_map[lf_vid].name
            @test occursin("z", vname)
            @test lf_vid.stage == PH.stid(3)
        end
    end

    @testset "Build Subproblems" begin
        st = build_scen_tree()

        wr = PH.WorkerRecord(myid())
        rc = RemoteChannel(()->Channel(2*length(PH.scenarios(st))), myid())

        PH._build_subproblems(rc, wr, PH.scenarios(st), st,
                              create_model, (), NamedTuple(), SubproblemCallback[]
                              )

        for (scid, sub) in pairs(wr.subproblems)
            @test length(JuMP.all_variables(sub.problem.model)) == 6
            ncons = 0
            for (ftype,stype) in JuMP.list_of_constraint_types(sub.problem.model)
                ncons += JuMP.num_constraints(sub.problem.model, ftype, stype)
            end
            @test ncons == 7
        end

        count = 0
        my_var_map = Dict{PH.VariableID,PH.VariableInfo}()
        while isready(rc)
            count += 1
            msg = take!(rc)
            @test typeof(msg) <: PH.VariableMap
            merge!(my_var_map, msg.var_info)
            for (vid,vinfo) in pairs(msg.var_info)
                @test vid.scenario == msg.scen
            end
        end
        @test length(my_var_map) == 24
    end

    @testset "Initialize Subproblems (Scalar)" begin
        st = build_scen_tree()

        sp_map = Dict(1 => Set([PH.ScenarioID(0), PH.ScenarioID(1), PH.ScenarioID(2)]),
                      2 => Set([PH.ScenarioID(3)])
                      )

        worker_input_queues = Dict(1 => RemoteChannel(()->Channel{PH.Message}(10), myid()),
                                   2 => RemoteChannel(()->Channel{PH.Message}(10), myid())
                                   )
        worker_output_queue = RemoteChannel(()->Channel{PH.Message}(10), myid())
        futures = Dict(1 => remotecall(PH.worker_loop,
                                       myid(),
                                       1,
                                       worker_input_queues[1],
                                       worker_output_queue),
                       2 => remotecall(PH.worker_loop,
                                       myid(),
                                       2,
                                       worker_input_queues[2],
                                       worker_output_queue),
                       )
        wi = PH.WorkerInf(worker_input_queues, worker_output_queue, futures)

        my_task = @async begin
            PH._initialize_subproblems(sp_map,
                                       wi,
                                       st,
                                       create_model,
                                       (),
                                       PH.ScalarPenaltyParameter(1.0),
                                       false,
                                       SubproblemCallback[])
        end

        if timeout_wait(my_task, 90) && !istaskfailed(my_task)

            var_map = fetch(my_task)

            for scen in PH.scenarios(st)
                for (vid, vinfo) in var_map[scen]
                    vname = vinfo.name
                    @test vid.scenario == scen
                    if vid.stage == PH.stid(1)
                        @test occursin("x", vname)
                    elseif vid.stage == PH.stid(2)
                        @test vname == "y"
                    elseif vid.stage == PH.stid(3)
                        @test occursin("z", vname)
                    else
                        stage_int = PH.value(vid.stage)
                        error("Stage $(stage_int)! There is no stage $(stage_int)!!")
                    end
                end
            end

        elseif istaskfailed(my_task)

            throw(my_task.exception)

        else

            error("Test timed out")

        end

        phd = PH.PHData(PH.ScalarPenaltyParameter(1.0),
                        st,
                        sp_map,
                        var_map,
                        TimerOutputs.TimerOutput()
                        )

        for (xhid, xhat) in pairs(phd.xhat)
            vname = phd.variable_data[first(xhat.vars)].name
            for vid in xhat.vars
                @test xhid == phd.variable_data[vid].xhat_id
                @test vname == phd.variable_data[vid].name
            end
        end

        for (scen, sinfo) in phd.scenario_map
            for vid in keys(sinfo.branch_vars)
                if occursin("x", phd.variable_data[vid].name)
                    @test vid.stage == PH.stid(1)
                else
                    @test vid.stage == PH.stid(2)
                end
            end

            for vid in keys(sinfo.leaf_vars)
                @test vid.stage == PH.stid(3)
            end
        end

        my_task = @async begin
            count = 0
            messages = Vector{PH.Message}()
            while count < length(PH.scenarios(st))
                msg = PH._retrieve_message_type(wi, PH.ReportBranch)
                push!(messages, msg)
                count += 1
            end
            return messages
        end

        if timeout_wait(my_task, 90) && !istaskfailed(my_task)

            msgs = fetch(my_task)
            @test length(msgs) == length(PH.scenarios(st))
            msg_from_scen = Set{PH.ScenarioID}()

            for msg in msgs
                scen = msg.scen
                push!(msg_from_scen, scen)

                sinfo = phd.scenario_map[scen]
                @test keys(sinfo.branch_vars) == keys(msg.vals)
                @test (msg.sts == MOI.OPTIMAL ||
                       msg.sts == MOI.LOCALLY_SOLVED ||
                       msg.sts == MOI.ALMOST_LOCALLY_SOLVED
                       )

                put!(wi.output, msg)
            end

            @test msg_from_scen == PH.scenarios(st)

        elseif istaskfailed(my_task)

            throw(my_task.exception)

        else

            error("Test timed out")

        end

        my_task = @async PH._set_initial_values(phd, wi)

        if timeout_wait(my_task, 90) && !istaskfailed(my_task)
            @test fetch(my_task) >= 0.0
        elseif istaskfailed(my_task)
            throw(my_task.exception)
        else
            error("Test timedout")
        end

        my_task = @async PH._set_penalty_parameter(phd, wi)
        if timeout_wait(my_task, 90) && !istaskfailed(my_task)
            @test fetch(my_task) >= 0.0
        elseif istaskfailed(my_task)
            throw(my_task.exception)
        else
            error("Test timedout")
        end

        PH._shutdown(wi)
        PH._wait_for_shutdown(wi)
    end

    @testset "Initialize Subproblems (Proportional)" begin
        st = build_scen_tree()

        sp_map = Dict(1 => Set([PH.ScenarioID(0), PH.ScenarioID(1), PH.ScenarioID(2)]),
                      2 => Set([PH.ScenarioID(3)])
                      )

        worker_input_queues = Dict(1 => RemoteChannel(()->Channel{PH.Message}(10), myid()),
                                   2 => RemoteChannel(()->Channel{PH.Message}(10), myid())
                                   )
        worker_output_queue = RemoteChannel(()->Channel{PH.Message}(10), myid())
        futures = Dict(1 => remotecall(PH.worker_loop,
                                       myid(),
                                       1,
                                       worker_input_queues[1],
                                       worker_output_queue),
                       2 => remotecall(PH.worker_loop,
                                       myid(),
                                       2,
                                       worker_input_queues[2],
                                       worker_output_queue),
                       )
        wi = PH.WorkerInf(worker_input_queues, worker_output_queue, futures)

        my_task = @async begin
            PH._initialize_subproblems(sp_map,
                                       wi,
                                       st,
                                       create_model,
                                       (),
                                       PH.ProportionalPenaltyParameter(1.0),
                                       false,
                                       SubproblemCallback[])
        end

        if timeout_wait(my_task, 90) && !istaskfailed(my_task)

            var_map = fetch(my_task)

            for scen in PH.scenarios(st)
                for (vid, vinfo) in var_map[scen]
                    vname = vinfo.name
                    @test vid.scenario == scen
                    if vid.stage == PH.stid(1)
                        @test occursin("x", vname)
                    elseif vid.stage == PH.stid(2)
                        @test vname == "y"
                    elseif vid.stage == PH.stid(3)
                        @test occursin("z", vname)
                    else
                        stage_int = PH.value(vid.stage)
                        error("Stage $(stage_int)! There is no stage $(stage_int)!!")
                    end
                end
            end

        elseif istaskfailed(my_task)

            throw(my_task.exception)

        else

            error("Test timed out")

        end

        phd = PH.PHData(PH.ProportionalPenaltyParameter(1.0),
                        st,
                        sp_map,
                        var_map,
                        TimerOutputs.TimerOutput()
                        )

        for (xhid, xhat) in pairs(phd.xhat)
            vname = phd.variable_data[first(xhat.vars)].name
            for vid in xhat.vars
                @test xhid == phd.variable_data[vid].xhat_id
                @test vname == phd.variable_data[vid].name
            end
        end

        for (scen, sinfo) in phd.scenario_map
            for vid in keys(sinfo.branch_vars)
                if occursin("x", phd.variable_data[vid].name)
                    @test vid.stage == PH.stid(1)
                else
                    @test vid.stage == PH.stid(2)
                end
            end

            for vid in keys(sinfo.leaf_vars)
                @test vid.stage == PH.stid(3)
            end
        end

        my_task = @async begin
            messages = Vector{PH.Message}()
            count = 0
            while count < 2*length(PH.scenarios(st))
                msg = PH._retrieve_message_type(wi, Union{PH.ReportBranch,PH.PenaltyInfo})
                push!(messages, msg)
                count += 1
            end
            return messages
        end

        if timeout_wait(my_task, 90) && !istaskfailed(my_task)

            msgs = fetch(my_task)
            @test length(msgs) == 2 * length(PH.scenarios(st))

            iv_msg_from_scen = Set{PH.ScenarioID}()
            pp_msg_from_scen = Set{PH.ScenarioID}()

            for msg in msgs

                scen = msg.scen
                sinfo = phd.scenario_map[scen]

                if typeof(msg) <: PH.ReportBranch
                    push!(iv_msg_from_scen, scen)
                    @test (msg.sts == MOI.OPTIMAL ||
                           msg.sts == MOI.LOCALLY_SOLVED ||
                           msg.sts == MOI.ALMOST_LOCALLY_SOLVED
                           )
                    @test keys(sinfo.branch_vars) == keys(msg.vals)
                elseif typeof(msg) <: PH.PenaltyInfo
                    push!(pp_msg_from_scen, scen)
                    @test keys(sinfo.branch_vars) == keys(msg.penalty)
                else
                    error("Unexpected message of type $(typeof(msg))")
                end

                put!(wi.output, msg)

            end

            @test iv_msg_from_scen == PH.scenarios(st)
            @test pp_msg_from_scen == PH.scenarios(st)

        elseif istaskfailed(my_task)

            throw(my_task.exception)

        else

            error("Test timed out")

        end

        my_task = @async PH._set_initial_values(phd, wi)
        if timeout_wait(my_task, 90) && !istaskfailed(my_task)
            @test fetch(my_task) >= 0.0
        elseif istaskfailed(my_task)
            throw(my_task.exception)
        else
            error("Test timedout")
        end

        my_task = @async PH._set_penalty_parameter(phd, wi)
        if timeout_wait(my_task, 90) && !istaskfailed(my_task)
            @test fetch(my_task) >= 0.0
        elseif istaskfailed(my_task)
            throw(my_task.exception)
        else
            error("Test timedout")
        end

        PH._shutdown(wi)
        PH._wait_for_shutdown(wi)
    end

    @testset "Initialize Subproblems (SEP)" begin
        st = build_scen_tree()

        sp_map = Dict(1 => Set([PH.ScenarioID(0), PH.ScenarioID(1), PH.ScenarioID(2)]),
                      2 => Set([PH.ScenarioID(3)])
                      )

        worker_input_queues = Dict(1 => RemoteChannel(()->Channel{PH.Message}(10), myid()),
                                   2 => RemoteChannel(()->Channel{PH.Message}(10), myid())
                                   )
        worker_output_queue = RemoteChannel(()->Channel{PH.Message}(10), myid())
        futures = Dict(1 => remotecall(PH.worker_loop,
                                       myid(),
                                       1,
                                       worker_input_queues[1],
                                       worker_output_queue),
                       2 => remotecall(PH.worker_loop,
                                       myid(),
                                       2,
                                       worker_input_queues[2],
                                       worker_output_queue),
                       )
        wi = PH.WorkerInf(worker_input_queues, worker_output_queue, futures)

        my_task = @async begin
            PH._initialize_subproblems(sp_map,
                                       wi,
                                       st,
                                       create_model,
                                       (),
                                       PH.SEPPenaltyParameter(),
                                       false,
                                       SubproblemCallback[])
        end

        if timeout_wait(my_task, 90) && !istaskfailed(my_task)

            var_map = fetch(my_task)

            for scen in PH.scenarios(st)
                for (vid, vinfo) in var_map[scen]
                    vname = vinfo.name
                    @test vid.scenario == scen
                    if vid.stage == PH.stid(1)
                        @test occursin("x", vname)
                    elseif vid.stage == PH.stid(2)
                        @test vname == "y"
                    elseif vid.stage == PH.stid(3)
                        @test occursin("z", vname)
                    else
                        stage_int = PH.value(vid.stage)
                        error("Stage $(stage_int)! There is no stage $(stage_int)!!")
                    end
                end
            end

        elseif istaskfailed(my_task)

            throw(my_task.exception)

        else

            error("Test timed out")

        end

        phd = PH.PHData(PH.SEPPenaltyParameter(),
                        st,
                        sp_map,
                        var_map,
                        TimerOutputs.TimerOutput()
                        )

        for (xhid, xhat) in pairs(phd.xhat)
            vname = phd.variable_data[first(xhat.vars)].name
            for vid in xhat.vars
                @test xhid == phd.variable_data[vid].xhat_id
                @test vname == phd.variable_data[vid].name
            end
        end

        for (scen, sinfo) in phd.scenario_map
            for vid in keys(sinfo.branch_vars)
                if occursin("x", phd.variable_data[vid].name)
                    @test vid.stage == PH.stid(1)
                else
                    @test vid.stage == PH.stid(2)
                end
            end

            for vid in keys(sinfo.leaf_vars)
                @test vid.stage == PH.stid(3)
            end
        end

        my_task = @async begin
            messages = Vector{PH.Message}()
            count = 0
            while count < 2*length(PH.scenarios(st))
                msg = PH._retrieve_message_type(wi, Union{PH.ReportBranch,PH.PenaltyInfo})
                push!(messages, msg)
                count += 1
            end
            return messages
        end

        if timeout_wait(my_task, 90) && !istaskfailed(my_task)

            msgs = fetch(my_task)
            @test length(msgs) == 2 * length(PH.scenarios(st))

            iv_msg_from_scen = Set{PH.ScenarioID}()
            pp_msg_from_scen = Set{PH.ScenarioID}()

            for msg in msgs

                scen = msg.scen
                sinfo = phd.scenario_map[scen]

                if typeof(msg) <: PH.ReportBranch
                    push!(iv_msg_from_scen, scen)
                    @test (msg.sts == MOI.OPTIMAL ||
                           msg.sts == MOI.LOCALLY_SOLVED ||
                           msg.sts == MOI.ALMOST_LOCALLY_SOLVED
                           )
                    @test keys(sinfo.branch_vars) == keys(msg.vals)
                elseif typeof(msg) <: PH.PenaltyInfo
                    push!(pp_msg_from_scen, scen)
                    @test keys(sinfo.branch_vars) == keys(msg.penalty)
                else
                    error("Unexpected message of type $(typeof(msg))")
                end

                put!(wi.output, msg)

            end

            @test iv_msg_from_scen == PH.scenarios(st)
            @test pp_msg_from_scen == PH.scenarios(st)

        elseif istaskfailed(my_task)

            throw(my_task.exception)

        else

            error("Test timed out")

        end

        my_task = @async PH._set_initial_values(phd, wi)
        if timeout_wait(my_task, 90) && !istaskfailed(my_task)
            @test fetch(my_task) >= 0.0
        elseif istaskfailed(my_task)
            throw(my_task.exception)
        else
            error("Test timedout")
        end

        my_task = @async PH._set_penalty_parameter(phd, wi)
        if timeout_wait(my_task, 90) && !istaskfailed(my_task)
            @test fetch(my_task) >= 0.0
        elseif istaskfailed(my_task)
            throw(my_task.exception)
        else
            error("Test timedout")
        end

        PH._shutdown(wi)
        PH._wait_for_shutdown(wi)
    end

    @testset "Lower Bound Initialization" begin

        st = build_scen_tree()
        sp_map = PH.assign_scenarios_to_procs(st)

        worker_input_queues = Dict(1 => RemoteChannel(()->Channel{PH.Message}(10), myid()))
        worker_output_queue = RemoteChannel(()->Channel{PH.Message}(10), myid())
        futures = Dict(1 => remotecall(PH.worker_loop,
                                       myid(),
                                       1,
                                       worker_input_queues[1],
                                       worker_output_queue)
                       )
        wi = PH.WorkerInf(worker_input_queues, worker_output_queue, futures)

        PH._initialize_lb_subproblems(sp_map, wi, st, create_model, (), false)
        PH._send_message(wi, 1, PH.DumpState())
        my_task = @async PH._retrieve_message(wi)
        if timeout_wait(my_task, 90) && !istaskfailed(my_task)
            msg = fetch(my_task)
        elseif istaskfailed(my_task)
            throw(my_task.exception)
        else
            error("Test timedout")
        end

        @test typeof(msg) <: PH.WorkerState{PH.WorkerRecord}
        @test length(msg.state.lb_subproblems) == length(scenarios(st))

        for (scen,sub) in pairs(msg.state.lb_subproblems)
            @test scen in [PH.scid(0), PH.scid(1), PH.scid(2), PH.scid(3)]
            @test length(sub.branch_vars) == 4
            @test length(sub.leaf_vars) == 2
        end

        PH._shutdown(wi)
        PH._wait_for_shutdown(wi)

    end
end


text:
using LinearAlgebra: svd, svdvals, cond, pinv, I, diagm
using Statistics
using Distributions
using Random
using JLD

include("../../../Potentials.jl/src/Potentials.jl")
include("solver_utils.jl")

# Parameters
kb              = 8.617e-5 # eV/K
Tc              = 120.0    # K
μ               = 39.948   # u
ϵ               = Tc*kb    # 0.0103404 eV
σ               = 3.4      # Angstroms Å

# Use reduced units
ϵ               = 1.0
σ               = 1.0
μ               = 1.0
Tc_kb           = 1.0
Tc              = 1.0

lj = Potentials.LennardJones(ϵ, σ)
c = Potentials.Configuration("examples/Argon7/DATA_training/DATA_start", atom_names = [:Ar], radii = [3.5], weights = [1.0])
# c = init_position(c, 2)
# c = init_velocity(c, 2, 0.25*Tc)
for Temp = [0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4]
    println("Temp $Temp")
    c = Potentials.Configuration("examples/Argon7/DATA_training/DATA_start", atom_names = [:Ar], radii = [3.5], weights = [1.0])
    Nt = 1000000
    dT = 1000
    dt = 0.01
    path = "examples/Argon7/DATA_training/TEMP_$(Int(100*Temp))/"
    # mkdir(path)
    for seed = 1:100
        if seed % 10 == 0
            println("seed $seed")
        end
        path = "examples/Argon7/DATA_training/TEMP_$(Int(100*Temp))/"
        path = path * "seed_$seed/"
        # mkdir(path)
        @time rdf, pe, msd = Potentials.run_md(c, lj, Nt, path; dt = dt, dim = 2, Temp = Temp, dT = dT, seed = seed)
    end
end

# Nt = 1000000
# dT = 1000
# dt = 0.01
# Temp = 0.30
# seed = 2
# T = dT:dT:Nt
# num_configs = length(T)
# r = Vector{PotentialUQ.Potentials.Configuration}(undef, num_configs)
# for (i,t) in enumerate(T)
#     file = "DATA.$t"
#     c_temp = PotentialUQ.Potentials.Configuration("examples/Argon7/DATA_training/TEMP_30/seed_$seed/"*file; atom_names = [:Ar],
#                     radii = [3.5], weights = [1.0])
#     r[i] = c_temp
# end

# for seed = 2:100
#     path = "examples/Argon7/DATA_md/LJ/seed_$seed/"
#     mkdir(path)
#     Nt = 50000
#     dT = 100
#     println("Seed $seed")
#     @time rdf, pe, msd = Potentials.run_md_lj(c, Nt, path; dim = 2, Temp = 0.75, dT = dT, seed = seed)
# end


text:
using BenchmarkTools

C_code = raw"""
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int compare(const void* l, const void* r){
    char* cl = *(char**)l+8;
    char* cr = *(char**)r+8;
    return strcmp(cl, cr);
}

void str_qsort(char **strings, size_t len) {  /* you need to pass len here */

    /* qsort(strings, len, sizeof(char*), (int (*)(void*, void*))strcmp); */
    qsort(strings, len, sizeof(char*), compare);

}
"""

const Clib = "strqsort" # tempname()   # make a temporary file

# compile to a shared library by piping C_code to gcc
# (works only if you have gcc installed):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

# define a Julia function that calls the C function:
str_qsort(X::Array{String}) = ccall(("str_qsort", Clib), Void, (Ptr{UInt64}, Cint), reinterpret(Ptr{UInt64}, pointer(X)), length(X))

a = [String([Char(i)]) for i in 'Z':-1:'A']
println("before sort: $a")
@btime str_qsort($a)
println("after sort: $a")

const M=100_000_000; const K=100
srand(1)
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M)
l = length(svec1)
println("before sort: length=$l")
# for i in 1:length(svec1)-1
#    if svec1[i]>svec1[i+1]
#       println("it is ok now to have unordered strings! svec1[$i]>svec1[$(i+1)] \"$(svec1[i])\">\"$(svec1[i+1])\"")
#       break
#    end
# end

@time str_qsort(svec1)
issorted(svec1)
# l = length(svec1)
# println("after sort length=$l")
# global sorted = true
# for i in 1:length(svec1)-1
#    if svec1[i]>svec1[i+1]
#       println("error svec1[$i]>svec1[$(i+1)] $(svec1[i])>$(svec1[i+1])")
#       sorted = false
#       break
#    end
# end
# println(sorted ? "sorted! :)" : "unsorted! :(")

const M=100_000_000; const K=100
using FastGroupBy
srand(1)
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M)
@time FastGroupBy.radixsort!(svec1) # 29 seconds
issorted(svec1)


text:
module TestFirstOrder

using ModiaLang
using DifferentialEquations
@usingModiaPlot

# using RuntimeGeneratedFunctions
# RuntimeGeneratedFunctions.init(@__MODULE__)

inputSignal(t) = sin(t)

FirstOrder = Model(
    T = 0.2,
    x = Var(init=0.3),
    equations = :[u = inputSignal(time/u"s"),
                  T * der(x) + x = u,
                  y = 2*x]
)

firstOrder = @instantiateModel(FirstOrder, logCode=false)

simulate!(firstOrder, Tsit5(), stopTime = 10, log=false, requiredFinalStates = [-0.3617373025974107])

plot(firstOrder, ["u", "x", "der(x)", "y"])

end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
# :PySet_New => (PyPtr,) => PyPtr,
# :PyFrozenSet_New => (PyPtr,) => PyPtr,
# :PySet_Add => (PyPtr, PyPtr) => Cint,

pyset_add(set::Py, x) = (errcheck(@autopy x C.PySet_Add(getptr(set), getptr(x_))); set)

function pyset_update_fromiter(set::Py, xs)
    for x in xs
        pyset_add(set, x)
    end
    return set
end
pyset_fromiter(xs) = pyset_update_fromiter(pyset(), xs)
pyfrozenset_fromiter(xs) = pyset_update_fromiter(pyfrozenset(), xs)

"""
    pyset(x=())

Convert `x` to a Python `set`.

If `x` is a Python object, this is equivalent to `set(x)` in Python.
Otherwise `x` must be iterable.
"""
pyset() = pynew(errcheck(C.PySet_New(C.PyNULL)))
pyset(x) = ispy(x) ? pybuiltins.set(x) : pyset_fromiter(x)
export pyset

"""
    pyfrozenset(x=())

Convert `x` to a Python `frozenset`.

If `x` is a Python object, this is equivalent to `frozenset(x)` in Python.
Otherwise `x` must be iterable.
"""
pyfrozenset() = pynew(errcheck(C.PyFrozenSet_New(C.PyNULL)))
pyfrozenset(x) = ispy(x) ? pybuiltins.frozenset(x) : pyfrozenset_fromiter(x)
export pyfrozenset


text:
using Test
using CompilerPluginTools

const GLOBAL_CONST = 2.0

struct Foo
    x::Int
end

@testset "inline_const" begin
    ir = @make_ircode begin
        Expr(:call, Core.Intrinsics.abs_float, 1.0)::Float64
        GlobalRef(Main, :GLOBAL_CONST)::Float64
        ReturnNode(SSAValue(1))::Float64
    end

    ir = inline_const!(ir)
    @test ir.stmts[1][:inst] == 1.0
    @test ir.stmts[1][:type] == Const(1.0)
    @test ir.stmts[2][:inst] == 2.0
    @test ir.stmts[2][:type] == Const(2.0)

    ir = @make_ircode begin
        Expr(:call, Core.Intrinsics.abs_float, 1.0)::Float64
        Expr(:new, Foo, 2)::Foo
        ReturnNode(SSAValue(1))::Float64
    end

    ir = inline_const!(ir)

    @test ir.stmts[2][:inst] == QuoteNode(Foo(2))
    @test ir.stmts[2][:type] == Const(Foo(2))

    ir = @make_ircode begin
        QuoteNode(1)::Int
        Expr(:new, Foo, SSAValue(1))::Foo
        ReturnNode(SSAValue(1))::Int
    end
    ir = inline_const!(ir)
    @test ir.stmts[2][:inst] == QuoteNode(Foo(1))
    @test ir.stmts[2][:type] == Const(Foo(1))

    ir = @make_ircode begin
        Expr(:call, Core.tuple, 1, 2, 3)::Tuple{Int, Int, Int}
        ReturnNode(SSAValue(1))::Tuple{Int, Int, Int}
    end
    ir = inline_const!(ir)
    @test ir.stmts[1][:inst] == (1, 2, 3)
    @test ir.stmts[1][:type] == Const((1, 2, 3))
end

ir = @make_ircode begin
    Expr(:call, Core.Intrinsics.abs_float, 1.0)::Float64
    Expr(:new, Foo, 2)::Foo
    ReturnNode(SSAValue(1))::Float64
end

ir = inline_const!(ir)

@testset "permute_stmts!" begin
    ir = @make_ircode begin
        QuoteNode(1.0)::Const(1.0)
        Expr(:call, sin, SSAValue(1))::Float64
        QuoteNode(2.0)::Const(2.0)
        Expr(:call, sin, SSAValue(3))::Float64
        Expr(:call, <, SSAValue(4), QuoteNode(1.0))::Bool
        GotoIfNot(SSAValue(5), 7)
        Expr(:call, +, SSAValue(2), SSAValue(4))::Float64
        ReturnNode(SSAValue(6))::Float64
    end

    ir = permute_stmts!(ir, [1, 3, 2, 4, 5, 6, 7, 8])
    @test ir.stmts[1][:inst] == QuoteNode(1.0)
    @test ir.stmts[2][:inst] == QuoteNode(2.0)
    @test ir.stmts[3][:inst] == Expr(:call, sin, SSAValue(1))
    @test ir.stmts[4][:inst] == Expr(:call, sin, SSAValue(2))
    @test ir.stmts[5][:inst] == Expr(:call, <, SSAValue(4), QuoteNode(1.0))
    @test ir.stmts[6][:inst] == GotoIfNot(SSAValue(5), 2)
    @test ir.stmts[7][:inst] == Expr(:call, +, SSAValue(3), SSAValue(4))
    @test ir.stmts[8][:inst] == ReturnNode(SSAValue(6))
end

function foo(x)
    2x
end

@testset "const_invoke!" begin
    mi = method_instances(foo, (Float64, ))[1]

    ir = @make_ircode begin
        Expr(:invoke, mi, GlobalRef(Main, :foo), 2.2)::Float64
        Expr(:invoke, mi, GlobalRef(Main, :foo), Argument(2))::Float64
        ReturnNode(SSAValue(1))::Float64
    end

    ir = const_invoke!(ir, GlobalRef(Main, :foo)) do x
        3x
    end

    @test ir.stmts[1][:inst] ≈ 3 * 2.2
    @test ir.stmts[2][:inst] == Expr(:invoke, mi, GlobalRef(Main, :foo), Argument(2))
end