text:
#####################################################################
#                                                                   #
# This file provides functions for interfacing with pti .raw files  #
#                                                                   #
#####################################################################

"""
    get_pti_sections()

Returns `Array` of the names of the sections, in the order that they
appear in a PTI file, v33+
"""
function get_pti_sections()::Array
    return ["CASE IDENTIFICATION", "BUS", "LOAD", "FIXED SHUNT", "GENERATOR", "BRANCH", "TRANSFORMER",
            "AREA INTERCHANGE", "TWO-TERMINAL DC", "VOLTAGE SOURCE CONVERTER", "IMPEDANCE CORRECTION",
            "MULTI-TERMINAL DC", "MULTI-SECTION LINE", "ZONE", "INTER-AREA TRANSFER", "OWNER",
            "FACTS CONTROL DEVICE", "SWITCHED SHUNT", "GNE DEVICE", "INDUCTION MACHINE"]
end

"""
    get_pti_dtypes(field_name)

Returns `OrderedDict` of data types for PTI file section given by `field_name`,
as enumerated by PSS/E Program Operation Manual
"""
function get_pti_dtypes(field_name::AbstractString)::Array
    transaction_dtypes = [("IC", Int64), ("SBASE", Float64), ("REV", Int64),
                          ("XFRRAT", Float64), ("NXFRAT", Float64), ("BASFRQ", Float64)]

    bus_dtypes = [("I", Int64), ("NAME", String), ("BASKV", Float64),
                  ("IDE", Int64), ("AREA", Int64), ("ZONE", Int64), ("OWNER", Int64),
                  ("VM", Float64), ("VA", Float64), ("NVHI", Float64), ("NVLO", Float64),
                  ("EVHI", Float64), ("EVLO", Float64)]

    load_dtypes = [("I", Int64), ("ID", String), ("STATUS", Int64), ("AREA", Int64),
                   ("ZONE", Int64), ("PL", Float64), ("QL", Float64), ("IP", Float64),
                   ("IQ", Float64), ("YP", Float64), ("YQ", Float64), ("OWNER", Int64),
                   ("SCALE", Int64), ("INTRPT", Int64)]

    fixded_shunt_dtypes = [("I", Int64), ("ID", String), ("STATUS", Int64), ("GL", Float64),
                           ("BL", Float64)]

    generator_dtypes = [("I", Int64), ("ID", String), ("PG", Float64), ("QG", Float64),
                        ("QT", Float64), ("QB", Float64), ("VS", Float64), ("IREG", Int64),
                        ("MBASE", Float64), ("ZR", Float64), ("ZX", Float64), ("RT", Float64),
                        ("XT", Float64), ("GTAP", Float64), ("STAT", Int64), ("RMPCT", Float64),
                        ("PT", Float64), ("PB", Float64), ("O1", Int64), ("F1", Float64),
                        ("O2", Int64), ("F2", Float64), ("O3", Int64), ("F3", Float64),
                        ("O4", Int64), ("F4", Float64), ("WMOD", Int64), ("WPF", Float64)]

    branch_dtypes = [("I", Int64), ("J", Int64), ("CKT", String), ("R", Float64), ("X", Float64),
                     ("B", Float64), ("RATEA", Float64), ("RATEB", Float64), ("RATEC", Float64),
                     ("GI", Float64), ("BI", Float64), ("GJ", Float64), ("BJ", Float64),
                     ("ST", Int64), ("MET", Int64), ("LEN", Float64), ("O1", Int64), ("F1", Float64),
                     ("O2", Int64), ("F2", Float64), ("O3", Int64), ("F3", Float64),
                     ("O4", Int64), ("F4", Float64)]

    transformer_3_dtypes = [("I", Int64), ("J", Int64), ("K", Int64), ("CKT", String),
                            ("CW", Int64), ("CZ", Int64), ("CM", Int64), ("MAG1", Float64),
                            ("MAG2", Float64), ("NMETR", Int64), ("NAME", String),
                            ("STAT", Int64),
                            ("O1", Int64), ("F1", Float64),
                            ("O2", Int64), ("F2", Float64),
                            ("O3", Int64), ("F3", Float64),
                            ("O4", Int64), ("F4", Float64),
                            ("VECGRP", String),

                            ("R1-2", Float64), ("X1-2", Float64), ("SBASE1-2", Float64),
                            ("R2-3", Float64), ("X2-3", Float64), ("SBASE2-3", Float64),
                            ("R3-1", Float64), ("X3-1", Float64), ("SBASE3-1", Float64),
                            ("VMSTAR", Float64), ("ANSTAR", Float64),

                            ("WINDV1", Float64), ("NOMV1", Float64), ("ANG1", Float64),
                            ("RATA1", Float64), ("RATB1", Float64), ("RATC1", Float64),
                            ("COD1", Int64), ("CONT1", Int64), ("RMA1", Float64), ("RMI1", Float64),
                            ("VMA1", Float64), ("VMI1", Float64), ("NTP1", Float64), ("TAB1", Int64),
                            ("CR1", Float64), ("CX1", Float64), ("CNXA1", Float64),

                            ("WINDV2", Float64), ("NOMV2", Float64), ("ANG2", Float64),
                            ("RATA2", Float64), ("RATB2", Float64), ("RATC2", Float64),
                            ("COD2", Int64), ("CONT2", Int64), ("RMA2", Float64), ("RMI2", Float64),
                            ("VMA2", Float64), ("VMI2", Float64), ("NTP2", Float64), ("TAB2", Int64),
                            ("CR2", Float64), ("CX2", Float64), ("CNXA2", Float64),

                            ("WINDV3", Float64), ("NOMV3", Float64), ("ANG3", Float64),
                            ("RATA3", Float64), ("RATB3", Float64), ("RATC3", Float64),
                            ("COD3", Int64), ("CONT3", Int64), ("RMA3", Float64), ("RMI3", Float64),
                            ("VMA3", Float64), ("VMI3", Float64), ("NTP3", Float64), ("TAB3", Int64),
                            ("CR3", Float64), ("CX3", Float64), ("CNXA3", Float64)]

    transformer_2_dtypes = [("I", Int64), ("J", Int64), ("K", Int64), ("CKT", String),
                            ("CW", Int64), ("CZ", Int64), ("CM", Int64), ("MAG1", Float64),
                            ("MAG2", Float64), ("NMETR", Int64), ("NAME", String),
                            ("STAT", Int64),
                            ("O1", Int64), ("F1", Float64),
                            ("O2", Int64), ("F2", Float64),
                            ("O3", Int64), ("F3", Float64),
                            ("O4", Int64), ("F4", Float64),
                            ("VECGRP", String),

                            ("R1-2", Float64), ("X1-2", Float64), ("SBASE1-2", Float64),

                            ("WINDV1", Float64), ("NOMV1", Float64), ("ANG1", Float64),
                            ("RATA1", Float64), ("RATB1", Float64), ("RATC1", Float64),
                            ("COD1", Int64), ("CONT1", Int64), ("RMA1", Float64), ("RMI1", Float64),
                            ("VMA1", Float64), ("VMI1", Float64), ("NTP1", Float64), ("TAB1", Int64),
                            ("CR1", Float64), ("CX1", Float64), ("CNXA1", Float64),

                            ("WINDV2", Float64), ("NOMV2", Float64)]

    area_interchange_dtypes = [("I", Int64), ("ISW", Int64),
                               ("PDES", Float64), ("PTOL", Float64),
                               ("ARNAME", String)]

    two_terminal_line_dtypes = [("NAME", String), ("MDC", Int64), ("RDC", Float64),
                                ("SETVL", Float64), ("VSCHD", Float64), ("VCMOD", Float64),
                                ("RCOMP", Float64), ("DELTI", Float64), ("METER", String),
                                ("DCVMIN", Float64), ("CCCITMX", Int64), ("CCCACC", Float64),
                                ("IPR", Int64), ("NBR", Int64), ("ANMXR", Float64),
                                ("ANMNR", Float64), ("RCR", Float64), ("XCR", Float64),
                                ("EBASR", Float64), ("TRR", Float64), ("TAPR", Float64),
                                ("TMXR", Float64), ("TMNR", Float64), ("STPR", Float64),
                                ("ICR", Int64), ("IFR", Int64), ("ITR", Int64), ("IDR", String),
                                ("XCAPR", Float64), ("IPI", Int64), ("NBI", Int64),
                                ("ANMXI", Float64), ("ANMNI", Float64), ("RCI", Float64),
                                ("XCI", Float64), ("EBASI", Float64), ("TRI", Float64),
                                ("TAPI", Float64), ("TMXI", Float64), ("TMNI", Float64),
                                ("STPI", Float64), ("ICI", Int64), ("IFI", Int64),
                                ("ITI", Int64), ("IDI", String), ("XCAPI", Float64)]

    vsc_line_dtypes = [("NAME", String), ("MDC", Int64), ("RDC", Float64),
                       ("O1", Int64), ("F1", Float64),
                       ("O2", Int64), ("F2", Float64),
                       ("O3", Int64), ("F3", Float64),
                       ("O4", Int64), ("F4", Float64)]

    vsc_subline_dtypes = [("IBUS", Int64), ("TYPE", Int64), ("MODE", Int64),
                          ("DCSET", Float64), ("ACSET", Float64),
                          ("ALOSS", Float64), ("BLOSS", Float64), ("MINLOSS", Float64),
                          ("SMAX", Float64), ("IMAX", Float64), ("PWF", Float64),
                          ("MAXQ", Float64), ("MINQ", Float64),
                          ("REMOT", Int64), ("RMPCT", Float64)]

    impedance_correction_dtypes = [("I", Int64), ("T1", Float64), ("F1", Float64),
                                   ("T2", Float64), ("F2", Float64), ("T3", Float64), ("F3", Float64),
                                   ("T4", Float64), ("F4", Float64), ("T5", Float64), ("F5", Float64),
                                   ("T6", Float64), ("F6", Float64), ("T7", Float64), ("F7", Float64),
                                   ("T8", Float64), ("F8", Float64), ("T9", Float64), ("F9", Float64),
                                   ("T10", Float64), ("F10", Float64), ("T11", Float64), ("F11", Float64)]

    multi_term_main_dtypes = [("NAME", String), ("NCONV", Int64), ("NDCBS", Int64), ("NDCLN", Int64),
                              ("MDC", Int64), ("VCONV", Int64), ("VCMOD", Float64), ("VCONVN", Float64)]

    multi_term_nconv_dtypes = [("IB", Int64), ("N", Int64), ("ANGMX", Float64), ("ANGMN", Float64),
                               ("RC", Float64), ("XC", Float64), ("EBAS", Float64), ("TR", Float64),
                               ("TAP", Float64), ("TPMX", Float64), ("TPMN", Float64), ("TSTP", Float64),
                               ("SETVL", Float64), ("DCPF", Float64), ("MARG", Float64), ("CNVCOD", Int64)]

    multi_term_ndcbs_dtypes = [("IDC", Int64), ("IB", Int64), ("AREA", Int64), ("ZONE", Int64),
                               ("DCNAME", String), ("IDC2", Int64), ("RGRND", Float64), ("OWNER", Int64)]

    multi_term_ndcln_dtypes = [("IDC", Int64), ("JDC", Int64), ("DCCKT", String), ("MET", Int64),
                               ("RDC", Float64), ("LDC", Float64)]

    multi_section_dtypes = [("I", Int64), ("J", Int64), ("ID", String), ("MET", Int64), ("DUM1", Int64),
                            ("DUM2", Int64), ("DUM3", Int64), ("DUM4", Int64), ("DUM5", Int64),
                            ("DUM6", Int64), ("DUM7", Int64), ("DUM8", Int64), ("DUM9", Int64)]

    zone_dtypes = [("I", Int64), ("ZONAME", String)]

    interarea_dtypes = [("ARFROM", Int64), ("ARTO", Int64), ("TRID", String), ("PTRAN", Float64)]

    owner_dtypes = [("I", Int64), ("OWNAME", String)]

    FACTS_dtypes = [("NAME", String), ("I", Int64), ("J", Int64), ("MODE", Int64), ("PDES", Float64),
                    ("QDES", Float64), ("VSET", Float64), ("SHMX", Float64), ("TRMX", Float64),
                    ("VTMN", Float64), ("VTMX", Float64), ("VSMX", Float64), ("IMX", Float64),
                    ("LINX", Float64), ("RMPCT", Float64), ("OWNER", Int64), ("SET1", Float64),
                    ("SET2", Float64), ("VSREF", Int64), ("REMOT", Int64), ("MNAME", String)]

    switched_shunt_dtypes = [("I", Int64), ("MODSW", Int64), ("ADJM", Int64), ("STAT", Int64),
                             ("VSWHI", Float64), ("VSWLO", Float64), ("SWREM", Int64),
                             ("RMPCT", Float64), ("RMIDNT", String), ("BINIT", Float64),
                             ("N1", Int64), ("B1", Float64), ("N2", Int64), ("B2", Float64),
                             ("N3", Int64), ("B3", Float64), ("N4", Int64), ("B4", Float64),
                             ("N5", Int64), ("B5", Float64), ("N6", Int64), ("B6", Float64),
                             ("N7", Int64), ("B7", Float64), ("N8", Int64), ("B8", Float64)]

    # TODO: Account for multiple lines in GNE Device entries
    gne_device_dtypes = [("NAME", String), ("MODEL", String), ("NTERM", Int64), ("BUSi", Int64),
                         ("NREAL", Int64), ("NINTG", Int64), ("NCHAR", Int64),

                         ("STATUS", Int64), ("OWNER", Int64), ("NMETR", Int64),

                         ("REALi", Float64),
                         ("INTGi", Int64),
                         ("CHARi", String)]

    induction_machine_dtypes = [("I", Int64), ("ID", String), ("STAT", Int64), ("SCODE", Int64),
                                ("DCODE", Int64), ("AREA", Int64), ("ZONE", Int64), ("OWNER", Int64),
                                ("TCODE", Int64), ("BCODE", Int64), ("MBASE", Float64),
                                ("RATEKV", Float64), ("PCODE", Int64), ("PSET", Float64),
                                ("H", Float64), ("A", Float64), ("B", Float64), ("D", Float64),
                                ("E", Float64), ("RA", Float64), ("XA", Float64), ("XM", Float64),
                                ("R1", Float64), ("X1", Float64), ("R2", Float64), ("X2", Float64),
                                ("X3", Float64), ("E1", Float64), ("SE1", Float64), ("E2", Float64),
                                ("SE2", Float64), ("IA1", Float64), ("IA2", Float64), ("XAMULT", Float64)]

    dtypes = Dict{String,Array}("BUS" => bus_dtypes,
                                "LOAD" => load_dtypes,
                                "FIXED SHUNT" => fixded_shunt_dtypes,
                                "GENERATOR" => generator_dtypes,
                                "BRANCH" => branch_dtypes,
                                "TRANSFORMER TWO WINDING" => transformer_2_dtypes,
                                "TRANSFORMER THREE WINDING" => transformer_3_dtypes,
                                "AREA INTERCHANGE" => area_interchange_dtypes,
                                "TWO-TERMINAL DC" => two_terminal_line_dtypes,
                                "VOLTAGE SOURCE CONVERTER" => vsc_line_dtypes,
                                "VOLTAGE SOURCE CONVERTER SUBLINES" => vsc_subline_dtypes,
                                "IMPEDANCE CORRECTION" => impedance_correction_dtypes,
                                "MULTI-TERMINAL DC" => multi_term_main_dtypes,
                                "MULTI-TERMINAL DC NCONV" => multi_term_nconv_dtypes,
                                "MULTI-TERMINAL DC NDCBS" => multi_term_ndcbs_dtypes,
                                "MULTI-TERMINAL DC NDCLN" => multi_term_ndcln_dtypes,
                                "MULTI-SECTION LINE" => multi_section_dtypes,
                                "ZONE" => zone_dtypes,
                                "INTER-AREA TRANSFER" => interarea_dtypes,
                                "OWNER" => owner_dtypes,
                                "FACTS CONTROL DEVICE" => FACTS_dtypes,
                                "SWITCHED SHUNT" => switched_shunt_dtypes,
                                "CASE IDENTIFICATION" => transaction_dtypes,
                                "GNE DEVICE" => gne_device_dtypes,
                                "INDUCTION MACHINE" => induction_machine_dtypes)

    return dtypes[field_name]
end

"""
    parse_line_element!(data, elements, section)

Parses a single "line" of data elements from a PTI file, as given by `elements`
which is an array of the line, typically split at `,`. Elements are parsed into
data types given by `section` and saved into `data::Dict`

"""
function parse_line_element!(data::Dict, elements::Array, section::AbstractString)
    missing = []
    for (field, dtype) in get_pti_dtypes(section)
        try
            element = popfirst!(elements)
        catch message
            if isa(message, ArgumentError)
                @debug "Have run out of elements in $section at $field"
                push!(missing, field)
                continue
            end
        end

        if startswith(element, "'") && endswith(element, "'")
            dtype = String
            element = chop(reverse(chop(reverse(element))))
        end

        try
            if dtype != String
                data[field] = parse(dtype, element)
            else
                data[field] = element
            end

        catch message
            if isa(message, Meta.ParseError)
                data[field] = element
            else
                error("value '$element' for $field in section $section is not of type $dtype.")
            end
        end
    end

    if length(missing) > 0
        missing_str = join(missing, ", ")
        if !(section == "SWITCHED SHUNT" && startswith(missing_str, "N")) &&
            !(section == "MULTI-SECTION LINE" && startswith(missing_str, "DUM")) &&
            !(section == "IMPEDANCE CORRECTION" && startswith(missing_str, "T"))
            @info("The following fields in $section are missing: $missing_str")
        end
    end
end

"""
    add_section_data!(pti_data, section_data, section)

Adds `section_data::Dict`, which contains all parsed elements of a PTI file
section given by `section`, into the parent `pti_data::Dict`
"""
function add_section_data!(pti_data::Dict, section_data::Dict, section::AbstractString)
    try
        pti_data[section] = append!(pti_data[section], [deepcopy(section_data)])
    catch message
        if isa(message, KeyError)
            pti_data[section] = [deepcopy(section_data)]
        else
            @error(sprint(showerror, message))
        end
    end
end

"""
    get_line_elements(line)

Uses regular expressions to extract all separate data elements from a line of
a PTI file and populate them into an `Array{String}`. Comments, typically
indicated at the end of a line with a `'/'` character, are also extracted
separately, and `Array{Array{String}, String}` is returned.
"""
function get_line_elements(line::AbstractString)::Array
    match_string = r"(-*\d*\.*\d+[eE]*[+-]*\d*)|(\'.{12}\')|(\'[^\']*?\')|(\"[^\"]*?\")|(\w+)|\,(\s+)?\,|(\/.*)"
    matches = collect((m.match for m = eachmatch(match_string, line, overlap=false)))
    #matches = matchall(match_string, line)

    @debug "$line"
    @debug "$matches"

    elements = []
    comment = ""
    for item in matches
        if startswith(item, ',') && endswith(item, ',') && length(item) == 2
            elements = append!(elements, [""])
        elseif startswith(item, '/')
            comment = item
        else
            elements = append!(elements, [item])
        end
    end

    return [elements, comment]
end

"""
    parse_pti_data(data_string, sections)

Parse a PTI raw file into a `Dict`, given the `data_string` of the file and a
list of the `sections` in the PTI file (typically given by default by
`get_pti_sections()`.
"""
function parse_pti_data(data_io::IO, sections::Array)
    data_lines = readlines(data_io)
    skip_lines = 0
    skip_sublines = 0
    subsection = ""

    pti_data = Dict{String,Array{Dict}}()

    section = popfirst!(sections)
    section_data = Dict{String,Any}()

    for (line_number, line) in enumerate(data_lines)
        @debug "$line_number: $line"

        (elements, comment) = get_line_elements(line)

        if length(elements) != 0 && elements[1] == "Q" && line_number > 3
            break

        elseif length(elements) != 0 && elements[1] == "0" && line_number > 3
            if line_number == 4
                section = popfirst!(sections)
            end

            match_string = r"\s*END OF ([\w\s-]+) DATA(?:, BEGIN ([\w\s-]+) DATA)?"
            @debug "$comment"
            matches = match(match_string, comment)

            if !isa(matches, Nothing)
                guess_section = matches.captures[1]
            else
                guess_section = ""
            end

            if guess_section == section && length(sections) > 0
                section = popfirst!(sections)
                continue
            else
                if length(elements) > 1
                    @info("At line $line_number, new section started with '0', but additional non-comment data is present. Pattern '^\\s*0\\s*[/]*.*' is reserved for section start/end.")
                elseif length(comment) > 0
                    @info("At line $line_number, unexpected section: expected: $section, comment specified: $(guess_section)")
                end
                if !isempty(sections)
                    section = popfirst!(sections)
                end

                continue
            end
        else
            if line_number == 4
                section = popfirst!(sections)
                section_data = Dict{String,Any}()
            end

            if skip_lines > 0
                skip_lines -= 1
                continue
            end

            @debug join(["Section:", section], " ")
            if !(section in ["CASE IDENTIFICATION","TRANSFORMER","VOLTAGE SOURCE CONVERTER","MULTI-TERMINAL DC","TWO-TERMINAL DC","GNE DEVICE"])
                section_data = Dict{String,Any}()
                try
                    parse_line_element!(section_data, elements, section)
                catch message
                    throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                end

            elseif section == "CASE IDENTIFICATION"
                if line_number == 1
                    try
                        parse_line_element!(section_data, elements, section)
                    catch message
                        throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                    end
                    try
                        if section_data["REV"] < 33
                            @info("Version $(section_data["REV"]) of PTI format is unsupported, parser may not function correctly.")
                        end
                    catch message
                        if isa(message, KeyError)
                            @error("This file is unrecognized and cannot be parsed")
                        end
                    end
                else
                    section_data["Comment_Line_$(line_number - 1)"] = line
                end

                if line_number < 3
                    continue
                end

            elseif section == "TRANSFORMER"
                section_data = Dict{String,Any}()
                if length(get_line_elements(data_lines[line_number + 1])[1]) == 3  && parse(Int64, get_line_elements(line)[1][3]) == 0 # two winding transformer
                    temp_section = "TRANSFORMER TWO WINDING"
                    (elements, comment) = get_line_elements(join(data_lines[line_number:line_number + 3], ','))
                    skip_lines = 3
                elseif length(get_line_elements(data_lines[line_number + 1])[1]) == 11 && parse(Int64, get_line_elements(line)[1][3]) != 0 # three winding transformer
                    temp_section = "TRANSFORMER THREE WINDING"
                    (elements, comment) = get_line_elements(join(data_lines[line_number:line_number + 4], ','))
                    skip_lines = 4
                else
                    @error("Cannot detect type of Transformer")
                end

                try
                    parse_line_element!(section_data, elements, temp_section)
                catch message
                    throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                end

            elseif section == "VOLTAGE SOURCE CONVERTER"
                if length(get_line_elements(line)[1]) == 11
                    section_data = Dict{String,Any}()
                    try
                        parse_line_element!(section_data, elements, section)
                    catch message
                        throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                    end
                    skip_sublines = 2
                    continue

                elseif skip_sublines > 0
                    skip_sublines -= 1

                    (elements, comment) = get_line_elements(line)
                    subsection_data = Dict{String,Any}()

                    for (field, dtype) in get_pti_dtypes("$section SUBLINES")
                        element = popfirst!(elements)
                        subsection_data[field] = parse(dtype, element)
                    end

                    try
                        section_data["CONVERTER BUSES"] = append!(section_data["CONVERTER BUSES"], [deepcopy(subsection_data)])
                    catch message
                        if isa(message, KeyError)
                            section_data["CONVERTER BUSES"] = [deepcopy(subsection_data)]
                            continue
                        else
                            @error(message)
                        end
                    end
                end

            elseif section == "TWO-TERMINAL DC"
                section_data = Dict{String,Any}()
                if length(get_line_elements(line)[1]) == 12
                    (elements, comment) = get_line_elements(join(data_lines[line_number:line_number + 2], ','))
                    skip_lines = 2
                end

                try
                    parse_line_element!(section_data, elements, section)
                catch message
                    throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                end

            elseif section == "MULTI-TERMINAL DC"
                if skip_sublines == 0
                    section_data = Dict{String,Any}()
                    try
                        parse_line_element!(section_data, elements, section)
                    catch message
                        throw(@error("Parsing failed at line $line_number: $(sprint(showerror, message))"))
                    end

                    if section_data["NCONV"] > 0
                        skip_sublines = section_data["NCONV"]
                        subsection = "NCONV"
                        continue
                    elseif section_data["NDCBS"] > 0
                        skip_sublines = section_data["NDCBS"]
                        subsection = "NDCBS"
                        continue
                    elseif section_data["NDCLN"] > 0
                        skip_sublines = section_data["NDCLN"]
                        subsection = "NDCLN"
                        continue
                    end
                end

                if skip_sublines > 0
                    skip_sublines -= 1

                    subsection_data = Dict{String,Any}()

                    for (field, dtype) in get_pti_dtypes("$section $subsection")
                        element = popfirst!(elements)
                        if startswith(element, "'") && endswith(element, "'")
                            subsection_data[field] = element[2:end-1]
                        else
                            subsection_data[field] = parse(dtype, element)
                        end
                    end

                    try
                        section_data["$(subsection[2:end])"] = append!(section_data["$(subsection[2:end])"], [deepcopy(subsection_data)])
                        if skip_sublines > 0 && subsection != "NDCLN"
                            continue
                        end
                    catch message
                        if isa(message, KeyError)
                            section_data["$(subsection[2:end])"] = [deepcopy(subsection_data)]
                            if skip_sublines > 0 && subsection != "NDCLN"
                                continue
                            end
                        else
                            @error(sprint(showerror, message))
                        end
                    end

                    if skip_sublines == 0 && subsection != "NDCLN"
                        if subsection == "NDCBS"
                            skip_sublines = section_data["NDCLN"]
                            subsection = "NDCLN"
                            continue
                        elseif subsection == "NCONV"
                            skip_sublines = section_data["NDCBS"]
                            subsection = "NDCBS"
                            continue
                        end
                    elseif skip_sublines == 0 && subsection == "NDCLN"
                        subsection = ""
                    else
                        continue
                    end
                end

            elseif section == "GNE DEVICE"
                # TODO: handle multiple lines of GNE Device
                @info("GNE DEVICE parsing is not supported.")
            end
        end
        if subsection != ""
            @debug "appending data"
        end
        add_section_data!(pti_data, section_data, section)
    end

    return pti_data
end

"""
    parse_pti(filename::String)

Open PTI raw file given by `filename`, returning a `Dict` of the data parsed
into the proper types.
"""
function parse_pti(filename::String)::Dict
    pti_data = open(filename) do f
        parse_pti(f)
    end

    return pti_data
end

"""
    parse_pti(io::IO)

Reads PTI data in `io::IO`, returning a `Dict` of the data parsed into the
proper types.
"""
function parse_pti(io::IO)::Dict
    pti_data = parse_pti_data(io, get_pti_sections())
    pti_data["CASE IDENTIFICATION"][1]["NAME"] = match(r"^\<file\s[\/\\]*(?:.*[\/\\])*(.*)\.raw\>$", lowercase(io.name)).captures[1]

    return pti_data
end


text:
using Test
using SimpleTensorNetworks

@testset "rm multi edges" begin
    t1 = randn(2,3,4,5)
    l1 = ['a', 'b', 'c', 'd']
    t2 = randn(6,4,7,2,3)
    l2 = ['e', 'c', 'f', 'a', 'b']
    t1_, t2_ = SimpleTensorNetworks.rm_multiedge(LabeledTensor(t1, l1), LabeledTensor(t2, l2))
    @test t1_.array == reshape(permutedims(t1, (4,1,2,3)), 5, :)
    @test t2_.array == reshape(permutedims(t2, (1,3,4,5,2)), 6,7,:)
    @test t1_.labels == ['d', 'a']
    @test t2_.labels == ['e', 'f', 'a']
end

@testset "rm single nodes" begin
    tn = TensorNetwork([
        LabeledTensor(randn(2), ["a"]),
        LabeledTensor(randn(2,2,2,2,2), ["g", "c", "a", "k", "l"]),
        LabeledTensor(randn(2), ["c"]),
        LabeledTensor(randn(2,2,2), ["g", "k", "l"]),
        ])
    r1 = contract(tn, [[[1,2], 3], 4])
    factor, tn = rm_degree12(tn)
    @test length(tn) == 2
    @test tn.tensors[1].labels == ["g", "k", "l"]
    @test tn.tensors[2].labels == ["g", "k", "l"]
    r2 = contract(tn, [1,2])
    @test r1 ≈ r2
    @test factor == 1.0

    tn = TensorNetwork([
        LabeledTensor(randn(2), ["a"]),
        LabeledTensor(randn(2,2,2,2), ["g", "c", "a", "k"]),
        LabeledTensor(randn(2), ["c"]),
        LabeledTensor(randn(2,2), ["g", "k"]),
        ])
    r1 = contract(tn, [[[1,2], 3], 4])
    factor, tn = rm_degree12(tn)
    @test length(tn) == 0
    @test r1.array[] ≈ factor
end

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

types = Any[
    Bool,
    Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Float32, Float64,
    Rational{Int8}, Rational{UInt8}, Rational{Int16}, Rational{UInt16},
    Rational{Int32}, Rational{UInt32}, Rational{Int64}, Rational{UInt64}
]
vals = vcat(
    typemin(Int64),
    -Int64(maxintfloat(Float64)) .+ Int64[-4:1;],
    typemin(Int32),
    -Integer(maxintfloat(Float32)) .+ (-4:1),
    -2:2,
    Integer(maxintfloat(Float32)) .+ (-1:4),
    typemax(Int32),
    Int64(maxintfloat(Float64)) .+ Int64[-1:4;],
    typemax(Int64),
)

function coerce(T::Type, x)
    if T<:Rational
        convert(T, coerce(typeof(numerator(zero(T))), x))
    elseif !(T<:Integer)
        convert(T, x)
    else
        x % T
    end
end

for T = types[2:end],
    x = vals,
    a = coerce(T, x)
    @test hash(a,zero(UInt)) == invoke(hash, Tuple{Real, UInt}, a, zero(UInt))
end

for T = types,
    S = types,
    x = vals,
    a = coerce(T, x),
    b = coerce(S, x)
    #println("$(typeof(a)) $a")
    #println("$(typeof(b)) $b")
    @test isequal(a,b) == (hash(a)==hash(b))
    # for y=vals
    #     println("T=$T; S=$S; x=$x; y=$y")
    #     c = convert(T,x//y)
    #     d = convert(S,x//y)
    #     @test !isequal(a,b) || hash(a)==hash(b)
    # end
end

# issue #8619
@test hash(nextfloat(2.0^63)) == hash(UInt64(nextfloat(2.0^63)))
@test hash(prevfloat(2.0^64)) == hash(UInt64(prevfloat(2.0^64)))

# issue #9264
@test hash(1//6,zero(UInt)) == invoke(hash, Tuple{Real, UInt}, 1//6, zero(UInt))
@test hash(1//6) == hash(big(1)//big(6))
@test hash(1//6) == hash(0x01//0x06)

# hashing collections (e.g. issue #6870)
vals = Any[
    [1,2,3,4], [1 3;2 4], Any[1,2,3,4], [1,3,2,4],
    [1.0, 2.0, 3.0, 4.0], BigInt[1, 2, 3, 4],
    [1,0], [true,false], BitArray([true,false]),
    # Irrationals
    Any[1, pi], [1, pi], [pi, pi], Any[pi, pi],
    # Overflow with Int8
    Any[Int8(127), Int8(-128), -383], 127:-255:-383,
    # Loss of precision with Float64
    Any[-2^53-1, 0.0, 2^53+1], [-2^53-1, 0, 2^53+1], (-2^53-1):2^53+1:(2^53+1),
    Set([1,2,3,4]),
    Set([1:10;]),                # these lead to different key orders
    Set([7,9,4,10,2,3,5,8,6,1]), #
    Dict(42 => 101, 77 => 93), Dict{Any,Any}(42 => 101, 77 => 93),
    (1,2,3,4), (1.0,2.0,3.0,4.0), (1,3,2,4),
    ("a","b"), (SubString("a",1,1), SubString("b",1,1)),
    # issue #6900
    Dict(x => x for x in 1:10),
    Dict(7=>7,9=>9,4=>4,10=>10,2=>2,3=>3,8=>8,5=>5,6=>6,1=>1),
    [], [1], [2], [1, 1], [1, 2], [1, 3], [2, 2], [1, 2, 2], [1, 3, 3],
    zeros(2, 2), spzeros(2, 2), Matrix(1.0I, 2, 2), sparse(1.0I, 2, 2),
    sparse(ones(2, 2)), ones(2, 2), sparse([0 0; 1 0]), [0 0; 1 0],
    [-0. 0; -0. 0.], SparseMatrixCSC(2, 2, [1, 3, 3], [1, 2], [-0., -0.]),
    # issue #16364
    1:4, 1:1:4, 1:-1:0, 1.0:4.0, 1.0:1.0:4.0, linspace(1, 4, 4),
    'a':'e', ['a', 'b', 'c', 'd', 'e'],
    # check that hash is still consistent with heterogeneous arrays for which - is defined
    # for some pairs and not others
    ["a", "b", 1, 2], ["a", 1, 2], ["a", "b", 2, 2], ["a", "a", 1, 2], ["a", "b", 2, 3]
]

for a in vals, b in vals
    @test isequal(a,b) == (hash(a)==hash(b))
    if a isa AbstractArray
        @test hash(a) == hash(collect(a)) == hash(collect(Any, a))
    end
end

vals = Any[
    Int[], Float64[],
    [0], [1], [2],
    # test vectors starting with ranges
    [1, 2], [1, 2, 3], [1, 2, 3, 4], [1, 2, 3, 4, 5], [1, 2, 3, 4, 5, 6],
    [2, 1], [3, 2, 1], [4, 3, 2, 1], [5, 4, 3, 2, 1], [5, 4, 3, 2, 1, 0, -1],
    # test vectors starting with ranges which trigger overflow with Int8
    [124, 125, 126, 127], [124, 125, 126, 127, -128], [-128, 127, -128],
    # test vectors including ranges
    [2, 1, 2, 3], [2, 3, 2, 1], [2, 1, 2, 3, 2], [2, 3, 2, 1, 2],
    # test various sparsity patterns
    [0, 0], [0, 0, 0], [0, 1], [1, 0],
    [0, 0, 1], [0, 1, 0], [1, 0, 0], [0, 1, 2],
    [0 0; 0 0], [1 0; 0 0], [0 1; 0 0], [0 0; 1 0], [0 0; 0 1],
    [5 1; 0 0], [1 1; 0 1], [0 2; 3 0], [0 2; 4 6], [4 0; 0 1],
    [0 0 0; 0 0 0], [1 0 0; 0 0 1], [0 0 2; 3 0 0], [0 0 7; 6 1 2],
    [4 0 0; 3 0 1], [0 2 4; 6 0 0],
    # various stored zeros patterns
    sparse([1], [1], [0]), sparse([1], [1], [-0.0]),
    sparse([1, 2], [1, 1], [-0.0, 0.0]), sparse([1, 2], [1, 1], [0.0, -0.0]),
    sparse([1, 2], [1, 1], [-0.0, 0.0], 3, 1), sparse([1, 2], [1, 1], [0.0, -0.0], 3, 1),
    sparse([1, 3], [1, 1], [-0.0, 0.0], 3, 1), sparse([1, 3], [1, 1], [0.0, -0.0], 3, 1),
    sparse([1, 2, 3], [1, 1, 1], [-1, 0, 1], 3, 1), sparse([1, 2, 3], [1, 1, 1], [-1.0, -0.0, 1.0], 3, 1),
    sparse([1, 3], [1, 1], [-1, 0], 3, 1), sparse([1, 2], [1, 1], [-1, 0], 3, 1)
]

for a in vals
    b = Array(a)
    @test hash(convert(Array{Any}, a)) == hash(b)
    @test hash(convert(Array{supertype(eltype(a))}, a)) == hash(b)
    @test hash(convert(Array{Float64}, a)) == hash(b)
    @test hash(sparse(a)) == hash(b)
    if !any(x -> isequal(x, -0.0), a)
        @test hash(convert(Array{Int}, a)) == hash(b)
        if all(x -> typemin(Int8) <= x <= typemax(Int8), a)
            @test hash(convert(Array{Int8}, a)) == hash(b)
        end
    end
end

# Test that overflow does not give inconsistent hashes with heterogeneous arrays
@test hash(Any[Int8(1), Int8(2), 255]) == hash([1, 2, 255])
@test hash(Any[Int8(127), Int8(-128), 129, 130]) ==
    hash([127, -128, 129, 130]) != hash([127,  128, 129, 130])

# Test hashing sparse matrix with type which does not support -
struct CustomHashReal
    x::Float64
end
Base.hash(x::CustomHashReal, h::UInt) = hash(x.x, h)
Base.:(==)(x::CustomHashReal, y::Number) = x.x == y
Base.:(==)(x::Number, y::CustomHashReal) = x == y.x
Base.zero(::Type{CustomHashReal}) = CustomHashReal(0.0)

let a = sparse([CustomHashReal(0), CustomHashReal(3), CustomHashReal(3)])
    @test hash(a) == hash(Array(a))
end

vals = Any[
    0.0:0.1:0.3, 0.3:-0.1:0.0,
    0:-1:1, 0.0:-1.0:1.0, 0.0:1.1:10.0, -4:10,
    'a':'e', 'b':'a',
    linspace(1, 1, 1), linspace(0.3, 1.0, 3),  linspace(1, 1.1, 20)
]

for a in vals
    @test hash(collect(a)) == hash(a)
end

@test hash(SubString("--hello--",3,7)) == hash("hello")
@test hash(:(X.x)) == hash(:(X.x))
@test hash(:(X.x)) != hash(:(X.y))

@test hash([1,2]) == hash(view([1,2,3,4],1:2))

let a = QuoteNode(1), b = QuoteNode(1.0)
    @test (hash(a)==hash(b)) == (a==b)
end

let a = Expr(:block, TypedSlot(1, Any)),
    b = Expr(:block, TypedSlot(1, Any)),
    c = Expr(:block, TypedSlot(3, Any))
    @test a == b && hash(a) == hash(b)
    @test a != c && hash(a) != hash(c)
    @test b != c && hash(b) != hash(c)
end

@test hash(Dict(),hash(Set())) != hash(Set(),hash(Dict()))

# issue 15659
for prec in [3, 11, 15, 16, 31, 32, 33, 63, 64, 65, 254, 255, 256, 257, 258, 1023, 1024, 1025],
    v in Any[-0.0, 0, 1, -1, 1//10, 2//10, 3//10, 1//2, pi]
    setprecision(prec) do
        x = convert(BigFloat, v)
        @test precision(x) == prec
        num, pow, den = Base.decompose(x)
        y = num*big(2.0)^pow/den
        @test precision(y) == prec
        @test isequal(x, y)
    end
end

# issue #20744
@test hash(:c, hash(:b, hash(:a))) != hash(:a, hash(:b, hash(:c)))

# issue #5849, object_id of types
@test Vector === (Array{T,1} where T)
@test (Pair{A,B} where A where B) !== (Pair{A,B} where B where A)
let vals_expr = :(Any[Vector, (Array{T,1} where T), 1, 2, Union{Int, String}, Union{String, Int},
                      (Union{String, T} where T), Ref{Ref{T} where T}, (Ref{Ref{T}} where T),
                      (Vector{T} where T<:Real), (Vector{T} where T<:Integer),
                      (Vector{T} where T>:Integer),
                      (Pair{A,B} where A where B), (Pair{A,B} where B where A)])
    vals_a = eval(vals_expr)
    vals_b = eval(vals_expr)
    for (i, a) in enumerate(vals_a), (j, b) in enumerate(vals_b)
        @test i != j || (a === b)
        @test (a === b) == (object_id(a) == object_id(b))
    end
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
"""
```
pseudo_transient_relaxation(stategrid::StateGrid, value_functions::NTuple{N, T}, Gs::NTuple{N, T},
    μ::AbstractArray{S}, Σ²::AbstractArray{S}, Δ::S;
    uniform::Bool = false, bc::Vector{Tuple{S, S}} = Vector{Tuple{T, T}}(undef, 0),
    Q = nothing, banded::Bool = false) where {S <: Real, T <: AbstractArray{<: Real}, N}

pseudo_transient_relaxation!(stategrid::StateGrid, new_value_functions::NTuple{N, T},
    value_functions::NTuple{N, T}, Gs::NTuple{N, T},
    μ::AbstractArray{S}, Σ²::AbstractArray{S}, Δ::S;
    uniform::Bool = false, bc::Vector{Tuple{S, S}} = Vector{Tuple{S, S}}(undef, 0),
    L₁ = nothing, L₂ = nothing, Q = nothing,
    banded::Bool = false) where {S <: Real, T <: AbstractArray{<: Real}, N}
```

The homogeneous diffusion must take the form

```
dXₜ / Xₜ = μₜ dt + Σₜ dW_t,
```

where `Σₜ` is a matrix and every component of the n-dimensional Brownian motion `Wₜ`
is assumed to be independent of the other components.

For a broad class of problems, the HJB can be represented in the form

```
0 = μᵥ - G(v)
```

where `G(v)` collects all the terms unrelated to `μᵥ` and thus contains the HJB's nonlinearities.
By Ito's lemma, this equation can be rewritten as

```
v G(v) = ∇v⋅(μₜXₜ) + ½ Tr[(ΣₜXₜ)ᵀ * (H∘v) * (ΣₜXₜ)] + ∂v / ∂t.
```

Because this PDE is derived from an HJB,
it is solved backward in time. The pseudo-transient relaxation method solves this problem
progressively by treating the nonlinear term `G(v)` as constant and solving the linear part
via an implicit Euler step.

Let `A` be the matrix approximation of the infinitessimal generator of an Ito diffusion. Then
the implicit Euler step is
```
-∂v / ∂t ≈ (vₙ₊₁ - vₙ) / dt + G(vₙ) vₙ₊₁ = A vₙ₊₁
(I / dt + diag(G(vₙ)) - A) vₙ₊₁   = vₙ / dt
(I + dt * (diag(G(vₙ)) - A)) vₙ₊₁ = vₙ
```
In general, `A` may be affine such that `A v = L v + b`, hence

```
(I + dt * (diag(G(vₙ)) - A)) vₙ₊₁ = (vₙ + b).
```

It is also convenient to transform `dt` from `(0, ∞)` to `(0, 1)` by
`Δ = dt / (1 + dt) ⇒ dt = Δ / (1 - Δ)`, hence

```
(I + Δ / (1 - Δ) * (diag(G(vₙ)) - A)) vₙ₊₁  = (vₙ + b)
(diag(1 - Δ)  + Δ * (diag(G(vₙ)) - A)) vₙ₊₁ = (1 - Δ) (vₙ + b)
```
"""
# TODO: Fix Σ² to take in a matrix and to not assume it has been squared already
function pseudo_transient_relaxation(stategrid::StateGrid,
                                     value_functions::NTuple{N, T}, Gs::NTuple{N, T},
                                     μ::AbstractArray{S}, Σ²::AbstractArray{S}, Δ::S;
                                     uniform::Bool = false, bc::Vector{Tuple{S, S}} = Vector{Tuple{T, T}}(undef, 0),
                                     Q = nothing, banded::Bool = false) where {S <: Real, T <: AbstractArray{<: Real}, N}
    new_value_functions = map(x -> similar(x), value_functions)
    _, _, err = pseudo_transient_relaxation!(stategrid, new_value_functions, value_functions, Gs, μ, Σ², Δ,
                                             uniform = uniform, bc = bc, Q = Q, banded = banded)
    return new_value_functions, err
end
function pseudo_transient_relaxation!(stategrid::StateGrid, new_value_functions::NTuple{N, T},
                                      value_functions::NTuple{N, T}, Gs::NTuple{N, T},
                                      μ::AbstractArray{S}, Σ²::AbstractArray{S}, Δ::S;
                                      uniform::Bool = false, bc::Vector{Tuple{S, S}} = Vector{Tuple{S, S}}(undef, 0),
                                      L₁ = nothing, L₂ = nothing, Q = nothing,
                                      banded::Bool = false) where {S <: Real, T <: AbstractArray{<: Real}, N}
    @assert (ndims(μ) == ndims(Σ²) && length(μ) == length(Σ²))  "The dimensions of `drift` and `volatility` must match."
    @assert length(value_functions) == length(Gs) "For each value function, we need the associated nonlinear component."
    if ndims(μ) == 1
        return _ptr_1D!(stategrid, new_value_functions, value_functions, Gs, μ, Σ², Δ;
                        uniform = uniform, bc = bc, L₁ = L₁, L₂ = L₂, Q = Q, banded = banded)
    else
        error("not implemented yet")
    end
end

function _ptr_1D!(stategrid::StateGrid, new_value_functions::NTuple{N, T},
                  value_functions::NTuple{N, T}, Gs::NTuple{N, T},
                  μ::AbstractVector{S}, Σ²::AbstractVector{S}, Δ::S; uniform::Bool = false,
                  bc::Vector{Tuple{S, S}} = Vector{Tuple{S, S}}(undef, 0), L₁ = nothing, L₂ = nothing,
                  Q = nothing, banded::Bool = false) where {S <: Real, T <: AbstractVector{<: Real}, N}

    if banded
        @warn "The BandedMatrix concretization for DiffEqOperators does not work for GhostDerivativeOperator types yet. The output is an Array."
    end

    # Set up finite differences and construct the infinitessimal generator
    x = values(stategrid.x)[1]
    n = length(stategrid)
    dx = uniform ?  x[2] - x[1] : nonuniform_grid_spacing(stategrid, bc)

    # If no DiffEqOperator objects are passed, then we assume the user wants a BandedMatrix w/first-order upwind first derivatives
    # and second order central difference second derivatives. Reflecting boundaries are assumed.
    if isnothing(Q) && isnothing(L₁) && isnothing(L₂)
        # Set up the implicit time step
        for (nvf, vf, G) in zip(new_value_functions, value_functions, Gs)
            _default_ptr_1D!(nvf, vf, G, x, dx, n, μ, Σ², Δ)
        end
    else
        if banded
            # Calculate number of diagonals for BandedMatrix representation
            L₁_n_offdiags = if isnothing(L₁)
                # first order upwinded by default
                1
            elseif typeof(L₁.parameters)[3] || L₁.approximation_order > 2
                # If upwinded, need diagonals for stencil length (minus 0th diagonal) on both sides of 0th diagonal
                # If CenteredDifference, approximation order > 2, need whole stencil length for boundaries,
                # and subtract 1 to not double count 0th diagonal
                L₁.stencil_length - 1
            elseif L₁.approximation_order == 2 # CenteredDifference, approximation order = 2
                1 # just need 1st upper and lower diagonals, plus main diagonal b/c handling boundary values
            end
            L₂_n_offdiags = if isnothing(L₂) # see CenteredDifference cases of L₁_ndiags
                1
            elseif L₂.approximation_order == 2
                1
            else
                L₂.stencil_length - 1
            end
        end

        if isnothing(Q) # Assume reflecting boundaries,
            # typical case for models where the volatility vanishes at boundaries and drift point "inward"
            Q = RobinBC((0., 1., 0.), (0., 1., 0.), dx)
        end
        if isnothing(L₁)
            # Filling coefficients inside this function can cause undefineds to be used in the concretization, reason unknown
            L₁ = UpwindDifference(1, 1, dx, n, ones(n)) # So just use 1s first
        end
        L₁.coefficients .= μ .* x
        L₁ = sparse(L₁ * Q)
        L₂ = sparse((Σ² .* x.^2 ./ 2.) * (isnothing(L₂) ? CenteredDifference(2, 2, dx, n) : L₂) * Q)

        # Set up the implicit time step
        if banded
            # Currently, BandedMatrix concretization doesn't work for GhostDerivativeOperator,
            # so form the BandedMatrix from the sparse concretization. The hope is that time is
            # saved during the left divide step.
            min_n_offdiags = min(L₁_n_offdiags, L₂_n_offdiags)
            diags_dict = Dict(i => diag(L₁[1], i) + diag(L₂[1], i) for i in 1:min_n_offdiags)
            for i in -min_n_offdiags:-1
                diags_dict[i] = -Δ * (diag(L₁[1], i) + diag(L₂[1], i)) # multiplied by -Δ here to obtain further speed up
            end
            orig_0th_diag = -Δ .* (Array(diag(L₁[1], 0)) + Array(diag(L₂[1], 0))) # Handle 0th diagonal separately b/c will add terms
            diags_dict[0] = similar(orig_0th_diag)                                # Also convert them to Array rather than SparseVector

            for (nvf, vf, G) in zip(new_value_functions, value_functions, Gs)
                diags_dict[0] .= (1 - Δ) .+ Δ .* G + orig_0th_diag # add nonlinear terms to main diagonal
                A = BandedMatrix(diags_dict...)

                # nvf .= A \ (1. - Δ) .* (vf + b)
                ldiv!(nvf, qr(A), (1. - Δ) .* (vf + b)) # Typically faster b/c fewer allocations, qr also faster than factorize/cholesky
            end
        else
            A = L₁[1] + L₂[1]
            b = L₁[2] + L₂[2]

            for (nvf, vf, G) in zip(new_value_functions, value_functions, Gs)
                nvf .= (Diagonal((1 - Δ) .+ Δ .* G) - (Δ .* A)) \ ((1 - Δ) .* (vf + b))
            end
        end
    end

    return new_value_functions, value_functions
end

function _default_ptr_1D!(nvf::T, vf::T, G::T, x::AbstractVector{S},
                          dx::AbstractVector{S}, n::Int, μ::AbstractVector{S}, Σ²::AbstractVector{S},
                          Δ::S) where {N <: Int, T <: AbstractVector{<: Real}, S <: Real}

    # Process Σ² matrix
    Σ²in = similar(Σ²)
    Σ²in[2:n - 1] .= Σ²[2:n - 1] .* x[2:n - 1].^2 ./ 2. # In non-uniform case, we utilize Fornberg weights, so just calculate
    Σ²in[1] = 0.                                        # the raw coeficients.
    Σ²in[n] = 0.

    # Populate diagonals and construct BandedMatrix
    DD, D0, DU = _centered_difference_reflecting_bc_weights(2, x, -Δ .* Σ²in) # Initialize diagonals for centered difference operator
    μx         = μ .* x
    DU       .-= Δ .* (max.( μx[1:n - 1], 0.) ./ dx[1:n - 1]) # up diagonal
    DD       .-= Δ .* (max.(-μx[2:n],     0.) ./ dx[2:n])     # down diagonal, note FD scheme makes DD negative, multiplied by negative drift ⇒ positive, then subtracted ⇒ negative

    D0          .+= (1 - Δ) .+ Δ .* G
    D0[1:n - 1] .+= Δ .* (max.( μx[1:n - 1], 0.) ./ dx[1:n - 1]) # Also have to add to the main diagonal
    D0[2:n]     .+= Δ .* (max.(-μx[2:n],     0.) ./ dx[2:n])     # for the first-order FD operator

    A             = BandedMatrix(0 => D0, 1 => DU, -1 => DD)

    # Solve linear system
    # nvf .= A \ ((1. - Δ) .* vf)
    ldiv!(nvf, qr(A), (1. - Δ) .* vf) # Typically faster b/c fewer allocations, qr also faster than factorize/cholesky
end

function _default_ptr_1D!(nvf::T, vf::T, G::T, x::AbstractVector{S},
                          dx::S, n::Int, μ::AbstractVector{S}, Σ²::AbstractVector{S},
                          Δ::S) where {N <: Int, T <: AbstractVector{<: Real}, S <: Real}

    # Process Σ² matrix
    Σ²in = similar(Σ²)                                           # In non-uniform case, for 2:n - 1, want to divide by the
    Σ²in[2:n - 1] .= Σ²[2:n - 1] .* x[2:n - 1].^2 ./ (2 .* dx^2) # forward and backward difference, but for the boundaries,
    Σ²in[1] = 0.                                                 # we impose reflecting boundaries.
    Σ²in[n] = 0.                                                 # Note that dx is length (n + 1).

    # Populate diagonals and construct BandedMatrix
    μx = μ .* x
    DU = -Δ .* (max.( μx[1:n - 1], 0.) ./ dx + Σ²in[1:n - 1]) # up diagonal
    DD = -Δ .* (max.(-μx[2:n],     0.) ./ dx + Σ²in[2:n])     # down diagonal, note FD scheme makes DD negative, multiplied by negative drift ⇒ positive, then subtracted ⇒ negative

    D0            = (1 - Δ) .+ Δ .* G
    D0[1:n - 1] .-= DU
    D0[2:n]     .-= DD

    A = BandedMatrix(0 => D0, 1 => DU, -1 => DD)

    # Solve linear system
    # nvf .= A \ ((1. - Δ) .* vf)
    ldiv!(nvf, qr(A), (1. - Δ) .* vf) # Typically faster b/c fewer allocations, qr also faster than factorize/cholesky
end


text:
module Optim

"""
`argmax(v, ℓ)`

Find `ω` such that `v(ω) != ⊥` and ℓ(ω) is maximised
"""
function argmax end

"""
`solution(v)`

Find `ω` that satisfies constraints of `v` (`v(ω) != ⊥`)
"""
function solution end

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using CVChannel
using Test
using DelimitedFiles

function noisyWHnonMult()
    print("\n\nThis script explores the Werner-Holevo channel's non-multiplicativity's robustness to noise.")
    print("\nTo do this we just take a linear combination of the Werner-Holevo channel with the maximally")
    print("depolarizing channel.")

    println("First we initialize the solver (One moment please...)")
    identChan(X) = X
    eaCVDual(choi(identChan,2,2),2,2)

    println("Now we actually get our data.")
    println("We consider d=3 and λ in [0,0.025,...,0.3] and mixing parameter α in [0,0.5,...,0.2].")
    println("NOTE: This will take a few minutes.")
    dim = 3
    depolChan(X) = 1/dim * [1 0 0 ; 0 1 0 ; 0 0 1]
    depol_choi = Choi(depolChan,dim,dim)
    λ_vals = 0:0.025:0.3
    α_vals = 0:0.05:0.2
    len, wid = length(λ_vals), length(α_vals)
    mult_table, cv_table = zeros(len,wid), zeros(len,wid)
    λ_ctr, α_ctr = 0, 1
    for λ in λ_vals
        println("Now evaluating for λ = ", λ)
        wern_choi = Choi(dim*wernerState(dim,λ),dim,dim)
        λ_ctr = λ_ctr + 1
        for α_amt in α_vals
            if α_amt == α_vals[Int(round(length(α_vals)/2))]
                println("Now halfway through this value of p.")
            end
            mix_choi = Choi((1-α_amt)*wern_choi.JN + α_amt*depol_choi.JN,dim,dim) #This uses linearity of Choi mapping
            cv_1, cv_2, par_cv, mult_val = pptCVMultiplicativity(mix_choi, mix_choi)
            mult_table[λ_ctr,α_ctr] = mult_val
            cv_table[λ_ctr,α_ctr] = cv_1
            α_ctr = α_ctr + 1
        end
        α_ctr = 1
    end

    info_vec = Vector{Union{Nothing,String}}(nothing, 2*len+3)
    info_vec[1] = "INFO:"
    info_vec[2] = "Generated by werner-holevo-non-multiplicativity-noise-robustness.jl"
    info_vec[3] = "dim = " * string(dim)
    info_vec[4] = "lambda_vals range: " * string(λ_vals)
    info_vec[5] = "alpha_vals range: " * string(α_vals)
    info_vec[6] = "Fixed lambda each row"
    info_vec[7] = "Fixed alpha each column"
    info_vec[8] = "First table is multiplicativity"
    info_vec[9] = "Second table is single copy cvPPT values"
    λ_label = vcat("", [λ_vals;])
    data_to_save_1 = hcat(λ_label,vcat([α_vals;]',mult_table))
    data_to_save_2 = hcat(λ_label,vcat([α_vals;]',cv_table))
    data_to_save = hcat(vcat(vcat(data_to_save_1,zeros(1,wid+1)),data_to_save_2),info_vec)
    println("Please name the file you'd like to write the results to: \n")
    file_name = readline()
    file_to_open = string(file_name,".csv")
    writedlm(file_to_open, data_to_save, ',')
    if all(mult_table -> mult_table < 1e-4 , mult_table)
        println("Something went awry. Please check your results.")
        return false
    else
        #We build printable multiplicativity table
        top_row = hcat(["λ↓  α:"],[α_vals;]')
        results = vcat(top_row,hcat([λ_vals;],mult_table))
        println("This is the table of multiplicativity values:")
        show(stdout, "text/plain", results)
        println("")
        println("Your results may also be found in your saved data file.")
        return true
    end
end

@testset "non-Multiplicativity with Noise" begin
    @test noisyWHnonMult()
end
println("Goodbye!")

max_forks_repo_name: ChitambarLab/cv-channel


text:
__precompile__(true)
module Peg

using HTTP, JSON2, Reexport

@reexport using DataStreams
@reexport using Dates

JSON2.write(io::IO, ::Missing) = write(io, "null")

const DATAPATH = string(@__DIR__, "/../app/data")

icontype(::Type{Union{T, Missing}}) where {T} = icontype(T)
icontype(::Type{T}) where {T <: AbstractString} = "TEXT"
icontype(::Type{T}) where {T <: Real} = "REAL"
icontype(::Type{T}) where {T <: Dates.TimeType} = "DATE"
icontype(::Type{Any}) = "OTHER"
icontype(x) = "OTHER"

mutable struct View
    index::Int
    name::String
    hide::Bool
    preCode::String
    sourceCode::String
    generatedQuery::String
    cacheSource::Bool
    sortSearch::String
    columnSearch::String
    columns::Vector{NamedTuple}
    limit::Union{Void,Int}
    offset::Union{Void,Int}
end

const DEFAULT_VIEW = View(0, "view_0", false, "", "(a=[1,2,3,1], b=[\"hey\", \"ho\", \"neighbor\", \"hi\"], c=[4.0, 5.0, 6.6, 6.0], d=[missing,0,0,0], new_col=[2,2,2,2], new_col2=[3,4,4,4])", "", true, "", "", [], 500, 0)

function generateQuery(view)
    columns = view.columns
    actions = :([])
    foreach(columns) do col
        if !get(col, :unselected, false)
            expr = :(())
            get(col, :group, false) && push!(expr.args, :(group=true))
            !isempty(get(col, :aggregate, "")) && push!(expr.args, :(aggregate=$(Meta.parse(col.aggregate))))
            if get(col, :sort, false)
                push!(expr.args, :(sort=true))
                push!(expr.args, :(sortasc=$(col.sortasc)))
                push!(expr.args, :(sortindex=$(col.sortidx)))
            end
            !isempty(get(col, :filter, "")) && push!(expr.args, :(filter=$(Meta.parse(col.filter))))
            !isempty(get(col, :having, "")) && push!(expr.args, :(having=$(Meta.parse(col.having))))
            if !isempty(get(col, :compute, ""))
                push!(expr.args, :(name=$(QuoteNode(Symbol(col.name)))))
                push!(expr.args, :(computeargs=$(Tuple(col.computeargs))))
                if col.aggregated
                    push!(expr.args, :(computeaggregate=$(Meta.parse(col.computeaggregate))))
                else
                    push!(expr.args, :(compute=$(Meta.parse(col.compute))))
                end
            else
                push!(expr.args, :(col=$(col.col)))
                push!(expr.args, :(name=$(QuoteNode(Symbol(col.name)))))
            end
            push!(actions.args, expr)
        end
    end
    source = view.cacheSource && haskey(SOURCES, view.index) ? SOURCES[view.index] : Meta.parse(view.sourceCode)
    q = :(Data.query($source, $actions, NamedTuple; limit=$(view.limit), offset=$(view.offset)))
    @show q
    return q
end

function queryHandler(req, resp)
    # parse request body into structure
    view = JSON2.read(String(req), View)
    # generate query
    query = generateQuery(view)
    view.generatedQuery = string(query)
    # return result
    res = @eval Main $query
    sch = Data.schema(res)
    rows, cols = size(sch)
    names = Data.header(sch)
    types = Data.types(sch)
    columnNames = Dict(x.name=>i for (i, x) in enumerate(view.columns))
    foreach(1:cols) do i
        col = columnNames[names[i]]
        view.columns[col] = merge(view.columns[col], (name=names[i], T=types[i], icon=icontype(types[i]), data=res[i]))
    end
    json = JSON2.write(view)
    resp = HTTP.Response(200, json)
    resp.headers["Access-Control-Allow-Origin"] = "*"
    open(joinpath(DATAPATH, view.name), "w+") do f
        write(f, json)
    end
    return resp
end

function preSourceHandler(req, resp)
    view = JSON2.read(String(req), View)
    code = Meta.parse(view.preCode)
    @show code
    @eval Main $code
    resp = HTTP.Response(200)
    resp.headers["Access-Control-Allow-Origin"] = "*"
    return resp
end

const SOURCES = Dict()

function sourceHandler(req, resp)
    view = JSON2.read(String(req), View)
    source = "Data.stream!($(view.sourceCode), NamedTuple)"
    code = Meta.parse(source)
    view.generatedQuery = view.sourceCode
    res = @eval Main $code
    sch = Data.schema(res)
    rows, cols = size(sch)
    names = Data.header(sch)
    types = Data.types(sch)
    empty!(view.columns)
    foreach(1:cols) do i
        push!(view.columns, (col=i, origname=names[i], name=names[i], T=types[i], icon=icontype(types[i]), data=res[i], filter=""))
    end
    json = JSON2.write(view)
    resp = HTTP.Response(200, json)
    resp.headers["Access-Control-Allow-Origin"] = "*"
    if view.cacheSource
        SOURCES[view.index] = res
    end
    open(joinpath(DATAPATH, view.name), "w+") do f
        write(f, json)
    end
    return resp
end

function initialize()
    # walk /data dir, JSON2.read each view, sort by index
    views = [JSON2.read(read(joinpath(DATAPATH, x), String), View) for x in readdir(DATAPATH)]
    sort!(views, by=x->x.index)
    if isempty(views)
        push!(views, DEFAULT_VIEW)
    end
    return views
end

function initialRequest(req, resp)
    resp = HTTP.Response(200, JSON2.write(initialize()))
    resp.headers["Access-Control-Allow-Origin"] = "*"
    return resp
end

function run()
    r = HTTP.Router()
    HTTP.register!(r, HTTP.GET, "/initialize", initialRequest)
    HTTP.register!(r, HTTP.POST, "/presource", preSourceHandler)
    HTTP.register!(r, HTTP.POST, "/source", sourceHandler)
    HTTP.register!(r, HTTP.POST, "/query", queryHandler)
    HTTP.serve(port=8083, handler=r)
end

end # module


text:
include("examples/meNumerics/discontinuities.jl")

d_fv = SFFM.SFMDistribution{SFFM.FVMesh}(x1_FV,model,fvmesh)
d_me = SFFM.SFMDistribution{SFFM.FRAPMesh}(x1_ME,model,frapmesh)
d_dg = SFFM.SFMDistribution{SFFM.DGMesh}(x1_DG,model,dgmesh)

plot(layout = SFFM.NPhases(model))
for i in SFFM.phases(model)
    f_fv(x) = SFFM.cdf(d_fv,model,x,i)
    f_me(x) = SFFM.cdf(d_me,model,x,i)
    f_dg(x) = SFFM.cdf(d_dg,model,x,i)
    plot!(f_fv,-1,13; subplot = i)
    plot!(f_me,-1,13; subplot = i)
    plot!(f_dg,-1,13; subplot = i)
end
# scatter!([f_fv(0)],[6],subplot=1)
# scatter!([f_me(0)],[6],subplot=1)
# scatter!([f_dg(0)],[6],subplot=1)
plot!()

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
newcolfunc(x) = length(x) == 1 ? @inbounds(x[1]) : rowtable(x)
newcolfunc(x::AbstractArray) = x

process_other!(table, plans::SelectionResult{<:SelectionTerm{Symbol}}) = plans

function process_other!(tab, plans)
    isempty(plans) && (return plans)
    colset = Set(Symbol[])
    other_idxs = Int[]
    other_keyfunc, other_valfunc = nothing, nothing
    for (i, term) in enumerate(plans)
        if colname(term) isa OtherSelection
            push!(other_idxs, i)
            other_keyfunc = extend(other_keyfunc, keyfunc(term))
            other_valfunc = extend(other_valfunc, valfunc(term))
        else
            push!(colset, colname(term))
        end
    end
    other_colnames = isempty(colset) ? colnames(tab) : setdiff(colnames(tab), colset)
    if !isempty(other_colnames)
        if length(other_idxs) > 1
            deleteat!(plans.s, @inbounds(other_idxs[2:end]))
        end
        splice!(
            plans.s,
            @inbounds(other_idxs[1]),
            SelectionTerm.(other_colnames, Ref(other_keyfunc), Ref(other_valfunc))
        )
    else
        deleteat!(plans.s, other_idxs)
    end
end

function rename_colnames(column_names, renamings)
    out::Vector{Symbol} = apply_rename.(renamings, column_names)
    seen = Set()
    to_report = Dict()

    for (renamed, renaming) in zip(out, renamings)
        isnothing(renaming) && (push!(seen, renamed))
    end

    for (i, (column_name, renamed, renaming)) in enumerate(zip(column_names, out, renamings))
        isnothing(renaming) && continue
        new_renamed = renamed
        index = 0

        while new_renamed in seen
            index += 1
            new_renamed = Symbol(string(renamed, index))
        end

        if index > 0
            to_report[column_name] = renamed => new_renamed
            renamed = new_renamed
        end

        push!(seen, renamed)
        out[i] = renamed
    end

    !isempty(to_report) && @warn "Following columns' renaming was modified to preserve uniqueness" to_report
    Tuple(out)
end

function transform_columns(tab, column_names, transforms)
    Tuple(apply_trans_nested(f, c, tab) for (c, f) in zip(column_names, transforms))
end

function select(tab, args...; kwargs...)
    has_args = !isempty(args)
    has_kwargs = !isempty(kwargs)

    if has_args
        # Selection queries -> triplets of selection, renamings and transforms that are generic,
        # can be applied to any Table. Multiple `args` are chained with an `|`.
        # The translate_types function also tries to guess the correct meaning of inputs,
        # i.e. that `:a => :A` is a `selection(:a) => renaming(:A)`.
        queries = mapfoldl(translate_types, OrSelection, ensure_tuple(args))
        # Selection results -> triplets of column names, renamings and transforms that are fitted
        # to this particular table, the generic selections are replaces with actual column names.
        # If multiple selections were overlapping, their renamings and transformations were combined.
        plans = resolve_nested(tab, queries)::SelectionResult
        # the other_cols() are resolved if present
        process_other!(tab, plans)
        # Prepare renamings -- produces unique output names to be applied later
        outnames = rename_colnames(colnames(plans), keyfuncs(plans))
        # Apply transformations to tables columns
        outcols = transform_columns(tab, colnames(plans), valfuncs(plans))
        # At this point, results are in a form of a NamedTuple
        nt = NamedTuple{outnames}(outcols)
    else
        nt = columntable(tab)
    end
    if has_kwargs
        for (new_column_name, _args) in pairs(kwargs)
            # Produce source_columns => transformation pairs.
            plan = resolve_flat(nt, _args)
            # Create the new column
            new_column = apply_trans_flat(last(plan), first(plan), new_column_name, nt)
            # Add the result to our table (NamedTuple), so the next column definition can use them
            nt = merge(nt, (; new_column_name => newcolfunc(new_column)))
        end
    end
    materializer(tab)(nt) # Return a Table of the same "kind" as the input.
end

function rename(tab, args...) end
function transform(tab, args...; kwargs...) end
function select_colnames(tab, args...) end
function select_renames(tab, args...) end

# outnames = Tuple{Vararg{Symbol,length(plans)}}(apply_rename(f, x) for (f, x) in zip(keyfuncs(plans), colnames(plans)))
# outcols = Tuple(apply_trans_nested(f, x, tab) for (f, x) in zip(valfuncs(plans), colnames(plans)))


text:
function riemann_zeta(s::Float64)

	degree = 1E6
	z = 0
	for n in 1:degree
		z += n^-s
	end

	return z::Float64

end

function object(d_object::Float64, q::Array{Float64, 1}, r::Float64)

	c = r^(2*(3-d_object))

	intensity = min.(c .* (r.*q).^-(d_object), c .* (r.*q).^-4)

	return intensity::Array{Float64, 1}

end

function fractal(d_fractal::Float64, d_object::Float64, r_monomer::Float64, form_factor)

	intensity = 0 .* form_factor(1.0)
	p = (d_fractal+3*d_object-9)/(3-d_object)
	z = riemann_zeta(-(p+1-(4/(3-d_object)-1)))
	for n in 1:1E6
		r = r_monomer * n^(1/(3-d_object))
		N = (n^p)/z
		intensity .+= N .* form_factor(r)
	end
	return intensity::Array{Float64, 1}

end

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

using Plots

q = 10 .^ (-1:0.1:1) # [nm^-1]

gr()
plt = plot()
plot!(q[1:6], 1E2.*q[1:6].^-2.5;
	seriestype = :line,
	xaxis = :log,
	yaxis = :log,
	minorgrid = true,
	xticks = 10.0 .^ (-1:1),
	yticks = 10.0 .^ (-5:6),
	linewidth = 1.0,
	xlabel = "Q-value [nm\u207B\u00B9]",
	ylabel = "Intensity [a.u.]",
	foreground_color_legend = nothing,
	background_color_legend = nothing,
	framestyle = :box,
	legend = :topright,
	legendfontsize = 7,
	guidefontsize = 7,
	tickfontsize = 7,
	linestyle = :solid,
	linecolor = RGB(0, 1, 0),
	linealpha = 1.0,
	lab = "\u221D q\u207B\u00B2\u22C5\u2075")

plot!([NaN], [NaN];
	linealpha = 0,
	lab = " ")

plot!([NaN], [NaN];
	linestyle = :dot,
	linecolor = RGB(0, 0, 0),
	linealpha = 0.5,
	lab = "Normal")

plot!([NaN], [NaN];
	linestyle = :solid,
	linecolor = RGB(0, 0, 0),
	linealpha = 0.5,
	lab = "Fractal-like")

plot!([NaN], [NaN];
	linealpha = 0,
	lab = "(d = 2.5)")

plot!([NaN], [NaN];
	linealpha = 0,
	lab = " ")

plot!([1], [1];
	seriestype = :scatter,
	markershape = :circle,
	markersize = 1,
	markercolor = RGB(1, 0, 1),
	markerstrokecolor = RGB(1, 0, 1),
	markeralpha = 0.5,
	lab = "Cylindrish")

plot!([1], [1];
	seriestype = :scatter,
	markershape = :circle,
	markersize = 1,
	markercolor = RGB(0, 0, 1),
	markerstrokecolor = RGB(0, 0, 1),
	markeralpha = 0.5,
	lab = "Spherish")

plot!(q, 1E1.*object(1.0, q, 1.0);
	linestyle = :dot,
	linecolor = RGB(1, 0, 1),
	linealpha = 0.5,
	lab = nothing)

plot!(q, 1E0.*object(0.0, q, 1.0);
	linestyle = :dot,
	linecolor = RGB(0, 0, 1),
	linealpha = 0.5,
	lab = nothing)

plot!(q, 1E1.*fractal(2.5, 1.0, 1.0, r -> object(1.0, q, r));
	linestyle = :solid,
	linecolor = RGB(1, 0, 1),
	linealpha = 0.5,
	lab = nothing)

plot!(q, 1E0.*fractal(2.5, 0.0, 1.0, r -> object(0.0, q, r));
	linestyle = :solid,
	linecolor = RGB(0, 0, 1),
	linealpha = 0.5,
	lab = nothing)

savefig(plt, "figure.pdf")

max_stars_repo_name: JonPallbo/Scatt

max_issues_repo_name: JonPallbo/Scatt

max_forks_repo_name: JonPallbo/Scatt


text:
#! /usr/bin/julia

# Rosetta Code, Case-sensitivity of identifiers

dog = "Benjamin"
Dog = "Samba"
DOG = "Bernie"

if is(dog, Dog)
    println("There is only one dog, ", DOG)
else
    println("The three dogs are:  ", dog, ", ", Dog, " and ", DOG)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module test_ctriangulate
using Triangulate

function test()
    nodes=Matrix{Cdouble}([1.0 0.0 ; 0.0 1.0 ; -1.0 0.0 ; 0.0 -1.0]')
    faces=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 1 ]')
    faceregions=Matrix{Cint}([1 2 3 4]')
    regionpoints=Matrix{Cdouble}([0.5 0.5 1 0.01;]')
    triin=Triangulate.CTriangulateIO()
    triout=Triangulate.CTriangulateIO()
    vorout=Triangulate.CTriangulateIO()

    triin.numberofpoints=Cint(size(nodes,2))
    triin.pointlist=pointer(nodes)
    triin.numberofsegments=size(faces,2)
    triin.segmentlist=pointer(faces)
    triin.segmentmarkerlist=pointer(faceregions)
    triin.numberofregions=size(regionpoints,2)
    triin.regionlist=pointer(regionpoints)

    triangulate("paAqQ",triin,triout,vorout)

    triout.numberofpoints==177 && triout.numberoftriangles==319 && triout.numberofsegments==33
end
end


text:
using FvCFD
include("../test/shockTubeFunctions.jl")

### Select a mesh ###
# meshPath = "test/OFshockTube_100"
meshPath = "test/OFshockTube_400"

# Load it
OFmesh = OpenFOAMMesh(meshPath)
nCells = size(OFmesh.cells, 1)

### Initialize non-uniform initial conditions ###
_, cellPrimitives = initializeShockTube3DFVM(nCells...)

### Run ###
boundaryConditions = [ zeroGradientBoundary, [], emptyBoundary, [] ]
cellPrimitives = solve(OFmesh, meshPath, cellPrimitives, boundaryConditions, ShuOsher, initDt=0.0001, endTime=0.14267, outputInterval=0.14267, targetCFL=0.5, silent=false, createVTKOutput=true)

max_stars_repo_name: henrystoldt/fvCFD


text:
using RungeKutta: name, order, nstages, coefficients, weights, nodes

@testset "$(rpad("Partitioned Tableaus",80))" begin

    @test typeof(PartitionedTableau(:PRK4, TableauRK4(), TableauRK4())) <: PartitionedTableau
    @test PartitionedTableau(:PRK4, TableauRK4(), TableauRK4()) == PartitionedTableau(:PRK4, TableauRK4())

    @test typeof(PartitionedTableauGauss(1)) <: PartitionedTableau
    @test order(PartitionedTableauGauss(1)) == 2
    @test nstages(PartitionedTableauGauss(1)) == 1
    @test issymplectic(PartitionedTableauGauss(1))
    @test issymmetric(PartitionedTableauGauss(1))
    @test all(all.(check_symmetry(PartitionedTableauGauss(1))))

    @test !isexplicit(PartitionedTableauGauss(1))
    @test  isimplicit(PartitionedTableauGauss(1))
    @test !isdiagnonallyimplicit(PartitionedTableauGauss(1))
    @test  isfullyimplicit(PartitionedTableauGauss(1))

    @test typeof(PartitionedTableauGauss(2)) <: PartitionedTableau
    @test order(PartitionedTableauGauss(2)) == 4
    @test nstages(PartitionedTableauGauss(2)) == 2
    @test issymplectic(PartitionedTableauGauss(2))
    @test issymmetric(PartitionedTableauGauss(2))
    @test all(all.(check_symmetry(PartitionedTableauGauss(2))))

    @test !isexplicit(PartitionedTableauGauss(2))
    @test  isimplicit(PartitionedTableauGauss(2))
    @test !isdiagnonallyimplicit(PartitionedTableauGauss(2))
    @test  isfullyimplicit(PartitionedTableauGauss(2))

    @test typeof(TableauLobattoIIIAIIIB(2)) <: PartitionedTableau
    @test order(TableauLobattoIIIAIIIB(2)) == 2
    @test nstages(TableauLobattoIIIAIIIB(2)) == 2
    @test issymplectic(TableauLobattoIIIAIIIB(2))
    @test issymmetric(TableauLobattoIIIAIIIB(2))

    @test  isexplicit(TableauLobattoIIIAIIIB(2))
    @test !isimplicit(TableauLobattoIIIAIIIB(2))
    @test !isdiagnonallyimplicit(TableauLobattoIIIAIIIB(2))
    @test !isfullyimplicit(TableauLobattoIIIAIIIB(2))

    @test typeof(TableauLobattoIIIAIIIB(3)) <: PartitionedTableau
    @test order(TableauLobattoIIIAIIIB(3)) == 4
    @test nstages(TableauLobattoIIIAIIIB(3)) == 3
    @test issymplectic(TableauLobattoIIIAIIIB(3))
    @test issymmetric(TableauLobattoIIIAIIIB(3))

    @test !isexplicit(TableauLobattoIIIAIIIB(3))
    @test  isimplicit(TableauLobattoIIIAIIIB(3))
    @test !isdiagnonallyimplicit(TableauLobattoIIIAIIIB(3))
    @test  isfullyimplicit(TableauLobattoIIIAIIIB(3))

    @test typeof(TableauLobattoIIIBIIIA(2)) <: PartitionedTableau
    @test order(TableauLobattoIIIBIIIA(2)) == 2
    @test nstages(TableauLobattoIIIBIIIA(2)) == 2
    @test issymplectic(TableauLobattoIIIBIIIA(2))
    @test issymmetric(TableauLobattoIIIBIIIA(2))

    @test  isexplicit(TableauLobattoIIIBIIIA(2))
    @test !isimplicit(TableauLobattoIIIBIIIA(2))
    @test !isdiagnonallyimplicit(TableauLobattoIIIBIIIA(2))
    @test !isfullyimplicit(TableauLobattoIIIBIIIA(2))

    @test typeof(TableauLobattoIIIBIIIA(3)) <: PartitionedTableau
    @test order(TableauLobattoIIIBIIIA(3)) == 4
    @test nstages(TableauLobattoIIIBIIIA(3)) == 3
    @test issymplectic(TableauLobattoIIIBIIIA(3))
    @test issymmetric(TableauLobattoIIIBIIIA(3))

    @test !isexplicit(TableauLobattoIIIBIIIA(3))
    @test  isimplicit(TableauLobattoIIIBIIIA(3))
    @test !isdiagnonallyimplicit(TableauLobattoIIIBIIIA(3))
    @test  isfullyimplicit(TableauLobattoIIIBIIIA(3))

    @test typeof(TableauLobattoIIICIIIC̄(2)) <: PartitionedTableau
    @test order(TableauLobattoIIICIIIC̄(2)) == 2
    @test nstages(TableauLobattoIIICIIIC̄(2)) == 2
    @test  issymplectic(TableauLobattoIIICIIIC̄(2))
    @test !issymmetric(TableauLobattoIIICIIIC̄(2))

    @test !isexplicit(TableauLobattoIIICIIIC̄(2))
    @test  isimplicit(TableauLobattoIIICIIIC̄(2))
    @test !isdiagnonallyimplicit(TableauLobattoIIICIIIC̄(2))
    @test  isfullyimplicit(TableauLobattoIIICIIIC̄(2))

    @test typeof(TableauLobattoIIICIIIC̄(3)) <: PartitionedTableau
    @test order(TableauLobattoIIICIIIC̄(3)) == 4
    @test nstages(TableauLobattoIIICIIIC̄(3)) == 3
    @test  issymplectic(TableauLobattoIIICIIIC̄(3))
    @test !issymmetric(TableauLobattoIIICIIIC̄(3))

    @test !isexplicit(TableauLobattoIIICIIIC̄(3))
    @test  isimplicit(TableauLobattoIIICIIIC̄(3))
    @test !isdiagnonallyimplicit(TableauLobattoIIICIIIC̄(3))
    @test  isfullyimplicit(TableauLobattoIIICIIIC̄(3))

    @test typeof(TableauLobattoIIIC̄IIIC(2)) <: PartitionedTableau
    @test order(TableauLobattoIIIC̄IIIC(2)) == 2
    @test nstages(TableauLobattoIIIC̄IIIC(2)) == 2
    @test  issymplectic(TableauLobattoIIIC̄IIIC(2))
    @test !issymmetric(TableauLobattoIIIC̄IIIC(2))

    @test !isexplicit(TableauLobattoIIIC̄IIIC(2))
    @test  isimplicit(TableauLobattoIIIC̄IIIC(2))
    @test !isdiagnonallyimplicit(TableauLobattoIIIC̄IIIC(2))
    @test  isfullyimplicit(TableauLobattoIIIC̄IIIC(2))

    @test typeof(TableauLobattoIIIC̄IIIC(3)) <: PartitionedTableau
    @test order(TableauLobattoIIIC̄IIIC(3)) == 4
    @test nstages(TableauLobattoIIIC̄IIIC(3)) == 3
    @test  issymplectic(TableauLobattoIIIC̄IIIC(3))
    @test !issymmetric(TableauLobattoIIIC̄IIIC(3))

    @test !isexplicit(TableauLobattoIIIC̄IIIC(3))
    @test  isimplicit(TableauLobattoIIIC̄IIIC(3))
    @test !isdiagnonallyimplicit(TableauLobattoIIIC̄IIIC(3))
    @test  isfullyimplicit(TableauLobattoIIIC̄IIIC(3))

end


text:
module test_julia_inferred

using Test

function f1(a::Bool)
    Any[a ? 1 : 1.0]
end

function f2(a::Bool)::Vector{Any}
       [a ? 1 : 1.0]
end

function f3(a::Bool)
       [a ? 1 : 1.0]
end

                       @inferred(f1(true))
                       @inferred(f2(true))
@test_throws Exception @inferred(f3(true))

end # module test_julia_inferred


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Statistics, DataFrames, CSV
	using Random, Distributions, Plots
	using HypothesisTests
	Random.seed!(0)
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.03"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data = CSV.read(datadir("machine1.csv"), DataFrame; delim=",")[:, 1]
end;

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	xBar, n = mean(data), length(data)
	sig = 1.2
	alpha = 0.1
	z = quantile(Normal(),1-alpha/2)

	Text("Calculating formula: $((xBar - z*sig/sqrt(n), xBar + z*sig/sqrt(n)))")
end

# ╔═╡ c69d5eae-3b0c-11eb-38aa-b578cce81150
Text("Using confint() function: $(confint(OneSampleZTest(xBar,sig,n),alpha))")

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.03"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═c69d5eae-3b0c-11eb-38aa-b578cce81150
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using Base.Test
# using Lint

# lintpkg("NGSIM")
using NGSIM


text:
using FixedPointAcceleration
function funcfunc(x::Array{Float64,1})
    # The first coordinate convergences to 4.0 by 1 unit per iterate.
    output = Array{Float64,1}(x)
    if abs(x[1] - 4.0) <= 1.0
        output[1] = 4.0
    elseif x[1] > 4.0
        output[1] = x[1] - 1.0
    else
        output[1] = x[1] + 1.0
    end
    # The second does aitken convergence to 2.3
    output[2] = x[2] + (2.3-x[2])/2.0
    return output
end
Inputs = [19.0,10.0]
fp_anderson = fixed_point(funcfunc, Inputs; Algorithm = :Anderson, PrintReports = false, ReplaceInvalids = :ReplaceElements)
fp_anderson.TerminationCondition_ == :ReachedConvergenceThreshold
fp_aitken = fixed_point(funcfunc, Inputs; Algorithm = :Aitken, PrintReports = false, ReplaceInvalids = :ReplaceElements)
fp_aitken.TerminationCondition_ == :ReachedConvergenceThreshold
fp_newton = fixed_point(funcfunc, Inputs; Algorithm = :Newton, PrintReports = false, ReplaceInvalids = :ReplaceElements)
fp_newton.TerminationCondition_ == :ReachedConvergenceThreshold
fp_simple = fixed_point(funcfunc, Inputs; Algorithm = :Simple, PrintReports = false, ReplaceInvalids = :ReplaceElements, ReportingSigFig =  Integer(10))
fp_simple.TerminationCondition_ == :ReachedConvergenceThreshold
fp_SEA = fixed_point(funcfunc, Inputs; Algorithm = :SEA, PrintReports = false, ReplaceInvalids = :ReplaceElements)
fp_SEA.TerminationCondition_ == :ReachedConvergenceThreshold
fp_VEA = fixed_point(funcfunc, Inputs; Algorithm = :VEA, PrintReports = false, ReplaceInvalids = :ReplaceElements)
fp_VEA.TerminationCondition_ == :ReachedConvergenceThreshold
fp_RRE = fixed_point(funcfunc, Inputs; Algorithm = :RRE, PrintReports = false, ReplaceInvalids = :ReplaceElements)
fp_RRE.TerminationCondition_ == :ReachedMaxIter # This one fails because it keeps proposing bad ideas.
fp_MPE = fixed_point(funcfunc, Inputs; Algorithm = :MPE, PrintReports = false, ReplaceInvalids = :ReplaceElements)
fp_MPE.TerminationCondition_ == :ReachedConvergenceThreshold


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

# parameters that govern counterfactual

remove_shock!(parameters, :global_sectoral_shock_njs)
remove_shock!(parameters, :idiosyncratic_shock_njs)


text:
using Test
using Corpuscles
using Unitful

const DATA_DIR = joinpath(@__DIR__, "../data")

@testset "conversions" begin

    # Particle Identites
    # PDGID <-> Geant3ID
    geant_id = Geant3ID(1)
    pdg_id = convert(PDGID, geant_id)
    @test isequal(22, pdg_id.value)
    geant_id = Geant3ID(3)
    pdg_id = convert(PDGID, geant_id)
    @test isequal(11, pdg_id.value)
    geant_id = Geant3ID(3)
    pdg_id = convert(PDGID, geant_id)
    @test !isequal(13, pdg_id.value)
    pdg_id = PDGID(22)
    geant_id = convert(Geant3ID, pdg_id)
    @test isequal(1, geant_id.value)
    pdg_id = PDGID(11)
    geant_id = convert(Geant3ID, pdg_id)
    @test isequal(3, geant_id.value)
    pdg_id = PDGID(13)
    geant_id = convert(Geant3ID, pdg_id)
    @test !isequal(1, geant_id.value)

    try
        convert(PDGID, Geant3ID(5000))
        @test false
    catch e
        @test e isa Corpuscles.IDException
        msg = sprint(showerror, e)
        @test msg == "ParticleID Error: No corresponding PDGID for Geant3ID(5000) found!"
    end

    try
        convert(Geant3ID, PDGID(2222212))
        @test false
    catch e
        @test e isa Corpuscles.IDException
        msg = sprint(showerror, e)
        @test msg == "ParticleID Error: No corresponding Geant3ID for PDGID(2222212) found!"
    end

    for particle in particles()
        if particle.pdgid in [PDGID(-9000321), PDGID(9000321), PDGID(9020113),
                              PDGID(-9000311), PDGID(9020213), PDGID(9000311),
                              PDGID(-9020213)]
            continue
        end
        @test convert(Geant3ID, particle.pdgid) !== nothing
    end

end

@testset "catalogs" begin
    catalog_files = map(basename, Corpuscles.available_catalog_files())
    @test "particle2008.csv" in catalog_files
    @test "particle2018.csv" in catalog_files
    @test "particle2019.csv" in catalog_files

    Corpuscles.use_catalog_file(joinpath(DATA_DIR, "particle2008.csv"))
    @test 0.054u"MeV" == Particle(553).width.value
    Corpuscles.use_catalog_file(joinpath(DATA_DIR, "particle2019.csv"))
    @test 0.05402u"MeV" == Particle(553).width.value

end

@testset "comparison" begin
    a = Corpuscles.MeasuredValue(1u"m", 10u"cm", 10u"cm")
    b = Corpuscles.MeasuredValue(1u"m", 10u"cm", 10u"cm")
    c = Corpuscles.MeasuredValue(1u"m", 5u"cm", 5u"cm")
    d = Corpuscles.MeasuredValue(0.8u"m", 5u"cm", 5u"cm")
    e = Corpuscles.MeasuredValue(0.9u"m", 5u"cm", 5u"cm")
    f = Corpuscles.MeasuredValue(90u"s", 5u"s", 5u"s")
    g = 91u"cm"

    @test a == b
    @test !(a == c)
    @test d < a
    @test !(e < a)
    @test !(e < c)
    @test g < c
    @test isapprox(g, e)
    @test isapprox(g, a)
    @test isapprox(a, b)
    @test isapprox(a, c)
    @test isapprox(a, e)
    @test !isapprox(a, d)

    try
        e < f
        @test false
    catch e
        @test e isa Unitful.DimensionError
    end
end

@testset "show and print" begin
    # TODO: use an `IOBuffer` instead of `devnull` and check output
    show(devnull, Particle(1))
    print(devnull, Particle(1))
end


text:
using MPI, Knet, ParticleSwarmOptimization
include("mnist.jl")
include("Args.jl")

function minibatch(x, y, batchsize)
    data = Any[]
    for i=1:batchsize:size(x,2)-batchsize+1
        j=i+batchsize-1
        push!(data, (x[:,i:j], y[:,i:j]))
    end
    return data
end

@knet function mnist2layer(x; winit=Gaussian(0,.1))
    h = wbf(x; out=64, f=:relu, winit=winit)
    return wbf(h; out=10, f=:soft, winit=winit)
end

function test(f, data, loss)
    sumloss = numloss = 0
    for (x,ygold) in data
        ypred = forw(f, x)
        sumloss += loss(ypred, ygold)
        numloss += 1
    end
    sumloss / numloss
end

function main()
    MPI.Init()
    args = parse_commandline()
    gpu(args["gpu"])
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    world = MPI.Comm_size(comm)

    if rank == 0
        println("Parsed args:")

        for (arg,val) in args
            println("  $arg  =>  $val")
        end
    end

    dtrn = minibatch(identity(MNIST.xtrn), identity(MNIST.ytrn), args["batchsize"])
    dtst = minibatch(identity(MNIST.xtst), identity(MNIST.ytst), args["batchsize"])

    problem = NNProblem(:mnist2layer, args["l"], args["h"], MiniBatch(dtrn), softloss)

    #last parameter for the parallelization
    gbest = pso(problem, args["num"], args["epoch"], args["phi1"], args["phi2"], args["w"], args["vmin"], args["vmax"], true)

    if rank == 0
        println("\nBest: $(gbest.pbest_fitness)\n")
        println("Acc: $(1 - test(gbest.x, dtst, zeroone))")
    end

    MPI.Finalize()
end

main()


text:
using Plots

function diffusion(t)
    x = range(-0.12,0.12,length = 2^14)
    Ayz = pi*(13e-6)^2
    M = 1e23
    D = 2.14e-6
    return M/(Ayz*sqrt.(4pi*D.*t)).*exp.(-x.^2.0./(4*D.*t))
end

plot(diffusion(50))
#
# function test(t)
#     x = range(-0.12,0.12,length = 2^14)
#     dx = x[2] - x[1]
#     k = 2pi.*fftshift(fftfreq(length(x), 1/dx))
#     Ayz = pi*(13e-6)^2
#     M = 1e23
#     t0 = 0.1
#     D = 2.14e-6
#     C = 0.04
#     u0 = diffsion0(x,t0, M, D, Ayz, 0)
#     uw0 = fftshift(fft(u0))
#     uw = @.  uw0*exp(-k^2*D*t)
#     ux = real.(ifft(ifftshift(uw)))
#     ux
#     prob = ODEProblem(func, uw0, ts, (D, Progress(floor(Int, Float64(0.24*1000)),2)))
#     sol = solve(prob, Vern7(), reltol=1e-4, abstol=1e-10)
#     res = zeros(Complex{Float64}, size(sol))
#     for i in eachindex(sol.t)
#         res[:,i] = ifft(sol[:,i])
#     end
#     res
# end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
import Pkg
Pkg.activate(mktempdir())
Pkg.develop(Pkg.PackageSpec(path=ARGS[1]))
import Pluto
s = Pluto.ServerSession()
urls = [
    "https://raw.githubusercontent.com/fonsp/Pluto.jl/v0.12.16/sample/Basic.jl",
    "https://raw.githubusercontent.com/fonsp/Pluto.jl/v0.12.16/sample/Basic.jl",
    "https://gist.githubusercontent.com/fonsp/4e164a262a60fc4bdd638e124e629d64/raw/8ffe93c680e539056068456a62dea7bf6b8eb622/basic_pkg_notebook.jl",
]
success = all(urls) do url
    @show url
    nb = Pluto.SessionActions.open_url(s, url; run_async=false)
    @show [c.output.body for c in nb.cells]
    nb_success = !any(c.errored for c in nb.cells)
    @info "Done" url nb_success
    nb_success
end
exit(success ? 0 : 1)


text:
# This is a comment.

# çŻ

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia