text:
#-------------------------------------------------------------------------------
# Representation of raw acquisition data

using Dates

"""
Data stored in the header preceding each channel of a raw acquisition
"""
struct ChannelHeader
    meas_uid      ::UInt32
    scan_counter  ::UInt32
    sequence_time ::UInt32
    channel_id    ::UInt16
end

struct LoopCounters <: FieldVector{14, UInt16}
    line::UInt16
    acquisition::UInt16
    slice::UInt16
    partition::UInt16
    echo::UInt16
    phase::UInt16
    repetition::UInt16
    set::UInt16
    seg::UInt16
    ida::UInt16
    idb::UInt16
    idc::UInt16
    idd::UInt16
    ide::UInt16
end

loop_counter_index(name) = findfirst(isequal(name), fieldnames(MagneticResonanceSignals.LoopCounters))::Int

"""
Data from a single acquisition by a Siemens MR scanner (software version VD?)
"""
struct Acquisition
    # Measurement Data Header (MDH) section.
    meas_uid                    ::UInt32
    scan_counter                ::UInt32
    time_stamp                  ::UInt32  # Mars time stamp?
    pmu_time_stamp              ::UInt32  # phase measurement unit time stamp (on Rx??)
    system_type                 ::UInt16
    # ptab?
    ptab_pos_delay              ::UInt16
    ptab_pos_x                  ::UInt32
    ptab_pos_y                  ::UInt32
    ptab_pos_z                  ::UInt32
    # flags (bitpacked in raw MDH data)
    # TODO: Consider leaving these packed & using `getproperty()`?
    acq_end                     ::Bool    # Set on placeholder packet after sequence end.
    rt_feedback                 ::Bool    # Real time feedback (ICE will process acquisition ASAP)
    hp_feedback                 ::Bool    # High priority feedback (?)
    sync_data                   ::Bool    # Data following will be seq->ICE raw data transfer (not an acquisition)
    raw_data_correction         ::Bool
    ref_phase_stab_scan         ::Bool
    phase_stab_scan             ::Bool
    sign_rev                    ::Bool
    phase_correction            ::Bool
    pat_ref_scan                ::Bool
    pat_ref_ima_scan            ::Bool
    reflect                     ::Bool
    noise_adj_scan              ::Bool
    # raw data size from ADC
    num_samples                 ::UInt16
    num_channels                ::UInt16
    # Hypercube indexing system for comms with reconstruction code
    loop_counters               ::LoopCounters
    # Info about precise timing of digital samples
    cutoff_pre                  ::UInt16
    cutoff_post                 ::UInt16
    kspace_centre_column        ::UInt16
    coil_select                 ::UInt16
    readout_offcentre           ::Float32
    time_since_rf               ::UInt32
    kspace_centre_line_num      ::UInt16
    kspace_centre_partition_num ::UInt16
    # 3D coordinage system
    slice_position              ::SVector{3,Float32} # offset (saggital,coronal,transverse) ?
    slice_rotation_quat         ::SVector{4,Float32}
    # per-acquisition sequence-dependent parameters as passed to ICE program.
    ice_program_params          ::SVector{24,UInt16}
    reserved_params             ::SVector{4,UInt16}
    application_counter         ::UInt16
    application_mask            ::UInt16
    # Measured magnetization
    channel_info::Vector{ChannelHeader} # num_channels
    data::Matrix{ComplexF32}           # num_samples × num_channels
end

#-------------------------------------------------------------------------------
# Metadata

"""
Data quality control flags for raw file parsing.

We generally use these flags rather than returning hard errors because it's
useful to be able to parse apparently incomplete or partially corrupt raw data.

Flags include:

* `AcquisitionsIncomplete` - the acquisition list appeared to be truncated
* `MeasHeaderEmpty` - the measurement was empty (likely due to a cancelled measurement)
"""
@enum DataStatus AcquisitionsIncomplete MeasHeaderEmpty

"""
    RxCoilElementData(
        element, coil_id, unique_key, coil_copy, coil_type,
        rx_channel_connected, adc_channel_connected, mux_channel_connected,
        insertion_time, element_selected
    )

Siemens MR receiver coil element data.

## Logical connectivity
adc_channel_connected corresponds to the `channel_id` in the channel
acquisition header, though off by 1 (?!)

## Physical connectivity:
`rx_channel_connected` is duplicated in Siemens "dual-density" systems - two
coils are transmitted on a single wire with different carrier frequencies.
These are (presumably) split before ADC, and ADC channels come in pairs.
"""
struct RxCoilElementData
    # CoilElementID
    element::String
    coil_id::String
    unique_key::UInt32
    coil_copy::Int
    # CoilProperties
    coil_type::Int
    # Physical / logical connectivity
    rx_channel_connected::Int
    adc_channel_connected::Int
    mux_channel_connected::Int
    insertion_time::Int
    element_selected::Int
end

#-------------------------------------------------------------------------------
"""
    load_twix(twix_file_name)
    load_twix(io)

    MRExperiment(metadata, coils, quality_control, acq_data)

Container for data from a generic magnetic resonance experiment: a series of
`Acquisition`s, each of which records the coil response due to induced nuclear
magnetization. This is currently specialized for the type of raw data produced
by a Siemens scanner but could be generalized.

Note that the input pulse sequence is not recorded in Siemens raw format, so
it's not available here. The only way to know this in full detail is to
simulate the sequence with the same input parameters using the Siemens `poet`
tool from the Siemens proprietary IDEA development environment.
"""
struct MRExperiment
    metadata # TODO: rename to yaps, as we're only parsing yaps data.
    coils::Vector{RxCoilElementData}
    quality_control::Set{DataStatus}
    data::Vector{Acquisition}
end

"""
Get labels for loop counters (valid for VE11C)
"""
counter_labels(::MRExperiment) =
    ["line", "acquisition", "slice", "partition", "echo",
     "phase", "repetition", "set", "seg",
     "ida", "idb", "idc", "idd", "ide"]

"""
Extract the date and time of the localizer images used as reference
"""
function ref_epoch(expt::MRExperiment)
    # On the IDEA forum, Eddie Auerbach suggests:
    #   "tReferenceImage0,1,2 are the unique IDs for the 3 localizer images
    #   used for slice prescription. In YAPS you have tFrameOfReference,
    #   which I'm guessing is an ID for the last acquired shim phasemap."
    dts = []
    for (k,str) in meta_search(expt, "tReferenceImage")
        try
            push!(dts, DateTime(last(split(str, '.'))[1:14], dateformat"yyyymmddHHMMSS"))
        catch
            # There seems to commonly (always?) be an invalid datetime in the
            # year 3000 here - just ignore that one.
        end
    end
    if !isempty(dts)
        return minimum(dts)
    end
    # Some scans may not have tReferenceImage however it seems there can be a
    # tFrameOfReference in the Yaps section of the Meas header. (For whatever
    # reason this doesn't make it into the MeasYaps header??)
    try
        m = match(r"\.(2\d{13})", expt.metadata["Meas.tFrameOfReference"])
        return DateTime(m[1], dateformat"yyyymmddHHMMSS")
    catch exc
        @warn "Could not find date of scan (reference epoch) in twix metadata" exception=exc
        return missing
    end
end

function ref_epoch(expt::MRExperiment, search_key)
    dts = []
    for (k,str) in meta_search(expt, search_key)
        try
            push!(dts, DateTime(last(split(str, '.'))[1:14], dateformat"yyyymmddHHMMSS"))
        catch
            # There seems to commonly (always?) be an invalid datetime in the
            # year 3000 here - just ignore that one.
        end
    end
    isempty(dts) ? missing : minimum(dts)
end

"""
Extract scanner software version from experiment metadata
"""
function software_version(expt::MRExperiment)
    m = expt.metadata
    get(m, "sProtConsistencyInfo.tBaselineString") do
        get(m, "sProtConsistencyInfo.tMeasuredBaselineString") do
            # On Numaris NX, detailed software version info seems missing from
            # twix, or at least it's not clear how to get it :-(
            get(m, "Dicom.SoftwareVersions", missing)
        end
    end
end

function standard_metadata(expt::MRExperiment)
    protname = get(expt.metadata, "tProtocolName", missing)
    seqname = get(expt.metadata, "tSequenceFileName", missing)
    frequency = get(expt.metadata, "sTXSPEC.asNucleusInfo[0].lFrequency",missing)*u"Hz"
    device_id = ("Siemens",
                 get(expt.metadata, "Dicom.ManufacturersModelName", missing),
                 get(expt.metadata, "Dicom.DeviceSerialNumber", missing))
    epoch = ref_epoch(expt)
    version = software_version(expt)
    MRMetadata(protname, seqname, version, device_id, epoch, frequency)
end

"""
    scanner_time(acq | expt)

Get internal scanner time stamp from an acquisition or sequence of timestamps
from an experiment.
"""
scanner_time(expt::MRExperiment) = scanner_time.(expt.data)
# 2.5 ms per count... is this reliable?
scanner_time(acq::Acquisition) = 0.0025u"s" * acq.time_stamp

dwell_time(expt::MRExperiment) = expt.metadata["sRXSPEC.alDwellTime[0]"]*u"ns"

function Base.show(io::IO, expt::MRExperiment)
    println(io,
         "MRExperiment metadata:")
    meta = standard_metadata(expt)
    print(io, """
            Protocol           = $(meta.protocol_name)
            Sequence File Name = $(meta.sequence_name)
            Software Version   = $(meta.software_version)
            Device Id          = $(meta.device_id)
            Reference Date     = $(meta.ref_epoch)
            Frequency          = $(meta.frequency)
            Coils              = $(unique([c.coil_id for c in expt.coils]))
          """)
    if isempty(expt.data)
        return
    end
    tmin,tmax = extrema(scanner_time(expt))
    print(io, """
          Acquisition summary:
            Number   = $(length(expt.data))
            Duration = $(uconvert(u"s", 1.0*(tmax-tmin)))
          """)
    counters = [a.loop_counters for a in expt.data]
    for (i,label) in enumerate(counter_labels(expt))
        cntmin,cntmax = extrema(c[i] for c in counters)
        if cntmin != 0 || cntmax != 0
            println(io, "  Loop index $i in [$cntmin,$cntmax]   ($label)")
        end
    end
end

#-------------------------------------------------------------------------------
# Functions for loading Siemens TWIX data

function check_twix_type(io)
    # Detect whether this is a twix file from VB or VD software version,
    # using a magic number heuristic from suspect.py
    m1,m2 = read(io, SVector{2,UInt32})
    seek(io, 0)
    if m1 == 0 && m2 < 64
        return :vd
    else
        error("""
              Unknown TWIX type (TWIX VB?  Or not a TWIX file?).
              You could try siemens_to_ismrmd or suspect.py as alternatives.""")
    end
end

"""
    load_twix(filename; header_only=false, acquisition_filter=(acq)->true,
              meas_selector=last)

Load raw Siemens twix ".dat" format, producing an `MRExperiment` containing a
sequence of acqisitions.

For large files, acquisitions may be filtered out during file loading by
providing a predicate `acquisition_filter(acq)` which returns true when `acq`
should be kept.  By default, all acquisitions are retained.

N4 VD-version twix may have data from more than one sequence (for example, the
tune-up data may be included). Normally the data you're looking for is in the
`meas_selector=last` chunk, but meas_selector is provided for the cases where
you want to select other measurements.
"""
function load_twix(io::IO; header_only=false, acquisition_filter=(acq)->true,
                   meas_selector=last)
    # TWIX is little endian binary data, with ascii header
    check_twix_type(io)
    header_sections,quality_control,acquisitions =
        load_twix_vd(io, header_only, acquisition_filter, meas_selector)
    # For now parse the MeasYaps section as it's is the easiest to parse and
    # contains parameters relevant to downstream interpretation by the ICE
    # program (...I think?)
    @debug "Header section names" keys(header_sections)
    metadata = Dict{String,Any}()
    try
        phoenix_ascconv = match(r"(### ASCCONV BEGIN.*### ASCCONV END ###)"s,
                                header_sections["Phoenix"]).captures[1]
        yaps_meta = parse_header_yaps(header_sections["MeasYaps"])
        phoenix_ascconv_meta = parse_header_yaps(phoenix_ascconv)

        # Check if the same key has different value in yaps and phoenix to warn user
        differences = Dict()
        for k in intersect(Set(keys(yaps_meta)), Set(keys(phoenix_ascconv_meta)))
            if yaps_meta[k] != phoenix_ascconv_meta[k]
                differences[k] = (MeasYapsASCCONV=yaps_meta[k],
                                  PhoenixASCCONV=phoenix_ascconv_meta[k])
            end
        end
        if length(differences) > 0
            @warn """
                  Different values found in MeasYaps and Phoenix ASCCONV section. As per
                  our testing in getting the correct frequency, MeasYaps sometimes is not
                  consistent with other sections, so this metadata is extracted from Phoenix.
                  """ differences
        end
        dicom_meta = match_xprot_header(header_sections["Dicom"], "Dicom.",
                                        ["SoftwareVersions", "DeviceSerialNumber", "InstitutionName", "Manufacturer", "ManufacturersModelName"])
        meas_meta  = match_xprot_header(header_sections["Meas"], "Meas.",
                                        ["tReferenceImage0", "tReferenceImage1", "tReferenceImage2",
                                         "tFrameOfReference"])
        metadata = merge(metadata, yaps_meta, phoenix_ascconv_meta, dicom_meta, meas_meta)
    catch exc
        @error "Could not parse header metadata" exception=(exc,catch_backtrace())
    end
    coils = RxCoilElementData[]
    try
        coils = parse_yaps_rx_coil_selection(metadata)
    catch exc
        @error "Could not parse coil metadata" exception=(exc,catch_backtrace())
    end
    MRExperiment(metadata, coils, quality_control, acquisitions)
end

load_twix(filename::String; kws...) = open(io->load_twix(io; kws...), filename)

"""
    dump_twix_headers(filename, dump_dir; meas_selector=last)

Dump twix header sections verbatim into files `dump_dir/sect_name`, one for
each section.  Mostly useful for debugging.
"""
function dump_twix_headers(filename, dump_dir; meas_selector=last)
    if !isdir(dump_dir)
        mkdir(dump_dir)
    end
    # TWIX is little endian binary data, with ascii header
    open(filename) do io
        check_twix_type(io)
        header_sections,_,_ = load_twix_vd(io, true, (a)->true, meas_selector)
        for (sect_name,sect_data) in header_sections
            open(joinpath(dump_dir, sect_name), "w") do out
                write(out, sect_data)
            end
        end
    end
end

function read_meas_headers(io)
    twix_id, num_measurements = read(io, SVector{2,Int32})
    # vd file can contain multiple measurements, but we only want the MRS.
    # Assume that the MRS is the last measurement.
    @debug "Opened TWIX VD" twix_id num_measurements

    measurement_data = []
    for i=1:num_measurements
        meas_uid, file_id = read(io, SVector{2,Int32})
        meas_offset, meas_length = read(io, SVector{2,Int64})
        patient_name = read(io, FixedString{64})
        protocol_name = read(io, FixedString{64})
        @debug "Reading TWIX VD Header"   #=
            =# meas_uid file_id           #=
            =# meas_offset meas_length    #=
            =# patient_name protocol_name
        push!(measurement_data,
              (meas_uid=meas_uid, file_id=file_id,
               meas_offset=meas_offset, meas_length=meas_length,
               patient_name=patient_name, protocol_name=protocol_name))
    end
    return twix_id, num_measurements, measurement_data
end

# The following is inspired by the twix reader in suspect.py.
function load_twix_vd(io, header_only, acquisition_filter, meas_selector)
    twix_id, num_measurements, measurement_data = read_meas_headers(io)
    meas_header = meas_selector(measurement_data)

    seek(io, meas_header.meas_offset)

    acquisitions = Acquisition[]
    quality_control = Set{DataStatus}()

    header_size = read(io, UInt32)
    if header_size < sizeof(UInt32)
        push!(quality_control, MeasHeaderEmpty)
        meas_remaining_bytes = read(io, meas_header.meas_length - sizeof(header_size))
        @warn """
              Unexpected empty measurement header sections; measurement is
              probably corrupt.
              """ header_size length(meas_remaining_bytes) iszero(meas_remaining_bytes)
        return Dict{String,String}(), quality_control, acquisitions
    end
    header      = read(io, header_size - sizeof(UInt32))
    header_sections = parse_twix_header_sections(IOBuffer(header))

    if header_only
        return header_sections, quality_control, acquisitions
    end

    # read each scan until we hit the acq_end flag
    while !eof(io)
        # the first four bytes contain some composite information
        temp = read(io, UInt32)
        DMA_length = temp & (2^26 - 1)
        pack_flag = Bool((temp >> 25) & 1)
        PCI_rx = temp >> 26
        #@show DMA_length pack_flag PCI_rx
        if DMA_length < sizeof(UInt32)
            # From observing various partially complete twix files, this seems
            # to happen when the scan is cancelled.
            @warn "Unexpected empty meas packet before acq_end flag - assuming truncated acquisition" DMA_length
            push!(quality_control, AcquisitionsIncomplete)
            break
        end

        acq_len = DMA_length-sizeof(UInt32)
        acq_buf = read(io, acq_len)
        if length(acq_buf) < acq_len
            @warn "Twix acquisition truncated at position $(position(io))" io
            push!(quality_control, AcquisitionsIncomplete)
            break
        end
        iob = IOBuffer(acq_buf)

        meas_uid       = read(iob, UInt32)
        scan_counter   = read(iob, UInt32)
        time_stamp     = read(iob, UInt32)
        pmu_time_stamp = read(iob, UInt32)

        # NB: Following block is not in VB format data
        system_type    = read(iob, UInt16)
        ptab_pos_delay = read(iob, UInt16)
        ptab_pos_x     = read(iob, UInt32)
        ptab_pos_y     = read(iob, UInt32)
        ptab_pos_z     = read(iob, UInt32)
        reserved       = read(iob, UInt32)

        # more composite information
        eval_info_mask = read(iob, UInt64)

        acq_end             = Bool((eval_info_mask      ) & 1)
        rt_feedback         = Bool((eval_info_mask >> 1 ) & 1)
        hp_feedback         = Bool((eval_info_mask >> 2 ) & 1)
        sync_data           = Bool((eval_info_mask >> 5 ) & 1)
        raw_data_correction = Bool((eval_info_mask >> 10) & 1)
        ref_phase_stab_scan = Bool((eval_info_mask >> 14) & 1)
        phase_stab_scan     = Bool((eval_info_mask >> 15) & 1)
        sign_rev            = Bool((eval_info_mask >> 17) & 1)
        phase_correction    = Bool((eval_info_mask >> 21) & 1)
        pat_ref_scan        = Bool((eval_info_mask >> 22) & 1)
        pat_ref_ima_scan    = Bool((eval_info_mask >> 23) & 1)
        reflect             = Bool((eval_info_mask >> 24) & 1)
        noise_adj_scan      = Bool((eval_info_mask >> 25) & 1)
        # if acq_end is set, there is no more data
        if acq_end
            break
        end

        # Real time data comes in several flavours:
        #   * MR data (coil potentials) as recorded by the ADC
        #   * Physiological data as recorded by the PMU (physiological measurement unit)
        #   * Real time sync data from sequence program to ICE
        #
        # For now, we ignore non-MR data which has a different binary layout.
        if sync_data || hp_feedback
            @debug "Ignoring unknown data packet" sync_data hp_feedback maxlog=10
            continue
        end
        #=
        if rt_feedback || hp_feedback || phase_correction || noise_adj_scan || sync_data
            continue
        end
        =#

        num_samples   = read(iob, UInt16)
        num_channels  = read(iob, UInt16)
        # TODO: Ordering of loop counters have changed between software versions!!
        #
        # read_meas_dat_mdh_binary__alt.m Lists the ordering of the loop counters as
        # line,acquisition,slice,partition,echo,phase,repetition,set,seg,ida,idb,idc,idd,ide.
        #
        # It looks like this was right as of software version D11, but changed
        # in a later version (at least set & repetition were swapped) (FIXME!)
        loop_counters = read(iob, SVector{14,UInt16})

        #@show Int(num_samples) Int(num_channels)
        #@show Int.(loop_counters)

        cutoff_pre                  = read(iob, UInt16)
        cutoff_post                 = read(iob, UInt16)
        kspace_centre_column        = read(iob, UInt16)
        coil_select                 = read(iob, UInt16)   # Dummy in VB ?
        readout_offcentre           = read(iob, Float32)
        time_since_rf               = read(iob, UInt32)
        kspace_centre_line_num      = read(iob, UInt16)
        kspace_centre_partition_num = read(iob, UInt16)
        # TWIX VB Would have following here instead of below; note different length too.
        # ice_program_params          = read(iob, SVector{4,UInt16})
        # reserved_params             = read(iob, SVector{4,UInt16})
        slice_position              = read(iob, SVector{3,Float32})
        slice_rotation_quat         = read(iob, SVector{4,Float32})
        ice_program_params          = read(iob, SVector{24,UInt16})
        reserved_params             = read(iob, SVector{4,UInt16})
        # TWIX VB doesn't have these. Instead
        # channel_id::UInt16
        # ptab_pos_neg::UInt16
        application_counter         = read(iob, UInt16)
        application_mask            = read(iob, UInt16)
        crc                         = read(iob, UInt32)

        #fid_start_offset = Int(ice_program_params[5])
        #num_dummy_points = Int(reserved_params[1])
        #fid_start = fid_start_offset + num_dummy_points
        # Suspect slices the raw data as follows. Why!?  (Is this specific to
        # Siemens SVS program?)
        #np = 2^floor(Int, log2(num_samples - fid_start))

        # read the data for each channel in turn
        channel_info = ChannelHeader[]
        data = zeros(ComplexF32, (num_samples, num_channels))
        for channel_index = 1:num_channels
            dma_length    = read(iob, UInt32)
            meas_uid      = read(iob, UInt32)
            scan_counter  = read(iob, UInt32)
                            read(iob, UInt32) # skip (padding? observed to be zeroed)
            sequence_time = read(iob, UInt32)
                            read(iob, UInt32) # skip
            channel_id    = read(iob, UInt16)
                            read(iob, UInt16) # skip
                            read(iob, UInt32) # skip
            push!(channel_info, ChannelHeader(meas_uid, scan_counter, sequence_time, channel_id))
            raw_data = read!(iob, Vector{ComplexF32}(undef,num_samples))
            # NB: The quadrature convetion used in twix is such that the
            # positive frequencies after a simple fft to correspond to a
            # *negative* offset from the spectrometer reference frequency.
            data[:,channel_index] .= raw_data
        end
        acq = Acquisition(
            meas_uid, scan_counter, time_stamp, pmu_time_stamp,
            system_type, ptab_pos_delay, ptab_pos_x, ptab_pos_y, ptab_pos_z,
            acq_end, rt_feedback, hp_feedback, sync_data, raw_data_correction,
            ref_phase_stab_scan, phase_stab_scan, sign_rev, phase_correction,
            pat_ref_scan, pat_ref_ima_scan, reflect, noise_adj_scan,
            num_samples, num_channels, loop_counters,
            cutoff_pre, cutoff_post, kspace_centre_column, coil_select, readout_offcentre,
            time_since_rf, kspace_centre_line_num, kspace_centre_partition_num, slice_position,
            slice_rotation_quat, ice_program_params, reserved_params, application_counter,
            application_mask, channel_info, data)
        if acquisition_filter(acq)
            push!(acquisitions, acq)
        end
    end
    header_sections, quality_control, acquisitions
end

function parse_twix_header_sections(io)
    # The VD twix header is a list of header sections containing ascii data but
    # connected with binary section sizes.
    #
    # Here's some (typical?) section names:
    # Config, Dicom, Meas, MeasYaps, Phoenix, Spice
    num_sec = read(io, UInt32)
    sections = Dict{String,String}()
    for sec_index = 1:num_sec
        sec_name    = String(readuntil(io, '\0', keep=false))
        sec_size    = read(io, UInt32)
        sec_content = read(io, sec_size)
        if sec_content[end] == 0
            sec_content = sec_content[1:end-1]
        end
        sections[sec_name] = String(sec_content)
    end
    sections
end

function match_xprot_header(xprot_text, key_prefix, xprot_keys)
    # Very heuristic "parsing" of XProtocol header section.
    # Ideally we'd have a proper XProtocol parser (see xprotocol.jl), but that
    # seems like a lot of work.
    metadata = Dict{String,String}()
    for key in xprot_keys
        m = match(Regex("""<ParamString\\."$key">\\s*\\{\\s*"([^}]*)"\\s*\\}""", "s"), xprot_text)
        if m !== nothing
            metadata[key_prefix*key] = m[1]
        end
    end
    metadata
end

function parse_header_yaps(yaps)
    metadata = Dict{String,Any}()
    for line in split(yaps, '\n')
        if startswith(line, "### ASCCONV BEGIN") || startswith(line, "### ASCCONV END") || isempty(line)
            continue
        end
        m = match(r"^([^#\s]+)\s*=\s*(.*?)\s*$", line)
        if m !== nothing
            if occursin(r"^\".*\"$", m[2])
                metadata[m[1]] = m[2][2:end-1]
            elseif occursin(r"^[+-]?\d+\.\d*([eE][-+]?\d+)?$", m[2])
                metadata[m[1]] = parse(Float64, m[2])
            else # assume integer (though may be a float)
                metadata[m[1]] = parse(Int, m[2])
            end
        else
            @warn "could not match YAPS header line \"$(Vector{UInt8}(line))\""
        end
    end
    metadata
end

function parse_yaps_rx_coil_selection(yaps)
    nucleus_index = 0
    prefix = "sCoilSelectMeas.aRxCoilSelectData[$nucleus_index]"
    nuc_name = yaps[prefix*".tNucleus"]
    if nuc_name != "1H"
        @warn "Expected nucleus 0 to be 1H - you will get coil data for $nuc_name instead"
    end
    ncoil = get(yaps, prefix*".asList.__attribute__.size", 0)
    if ncoil == 0
        @warn "No receive coils found in yaps data" yaps
        return
    end
    coils = RxCoilElementData[]
    for i = 0:ncoil-1
        try
            push!(coils, RxCoilElementData(
                yaps[prefix*".asList[$i].sCoilElementID.tElement"],
                yaps[prefix*".asList[$i].sCoilElementID.tCoilID"],
                reinterpret(UInt32, Int32(yaps[prefix*".asList[$i].sCoilElementID.ulUniqueKey"])),
                yaps[prefix*".asList[$i].sCoilElementID.lCoilCopy"],
                yaps[prefix*".asList[$i].sCoilProperties.eCoilType"],
                yaps[prefix*".asList[$i].lRxChannelConnected"],
                yaps[prefix*".asList[$i].lADCChannelConnected"],
                get(yaps, prefix*".asList[$i].lMuxChannelConnected", -1),
                get(yaps, prefix*".asList[$i].uiInsertionTime", -1),
                yaps[prefix*".asList[$i].lElementSelected"],
           ))
        catch exc
            exc isa KeyError || rethrow()
            # Unclear why this happens, but it seems that `ncoil` is not an
            # accurate estimation of the number of coils (is it an upper bound!?)
            # May be common, in which case we should downgrade this to a debug message.
            @warn "Could not find some metadata for coil $i" #=
                =# exc.key
        end
    end
    coils
end

"""
    meta_search(metadata, pattern)

Search through MR experiment `metadata` for a given regular expression,
`pattern` or for a case insensitive string `pattern`.
"""
function meta_search(expt::MRExperiment, pattern)
    Dict(k=>v for (k,v) in expt.metadata if occursin(pattern, k))
end
function meta_search(expt::MRExperiment, pattern::String)
    Dict(k=>v for (k,v) in expt.metadata if occursin(Regex(pattern,"i"), k))
end

#-------------------------------------------------------------------------------
# Raw data extraction and light processing specific to Siemens twix format.
#
# Ideally the functions here don't do signal processing, but defer that to the
# generic signal processing code instead.
#

function sampledata(expt, index; downsample=1)
    acq = expt.data[index]
    # The siemens sequence SVS_SE provides a few additional samples before and
    # after the desired ones. They comment that this is mainly to allow some
    # samples to be cut off after downsampling, (presumably to remove some of
    # the ringing artifacts of doing this with a simple FFT).
    cutpre  = acq.cutoff_pre
    cutpost = acq.cutoff_post
    data = acq.data
    # Adjust time so that the t=0 occurs in the first retained sample
    t = ((0:size(data,1)-1) .- cutpre)*dwell_time(expt)
    t,z = downsample_and_truncate(t, data, cutpre, cutpost, downsample)
    # We want the positive frequencies after a simple fft to correspond to a
    # positive offset from the spectrometer reference frequency. But in twix
    # they are negative so directly doing an fft of the raw data would flip the
    # frequency axis.  We add a conj to standardize the convention.
    z = conj.(z)
    coilsyms = if isempty(expt.coils)
        [Symbol("C$i") for i in 1:length(acq.channel_info)]
    else
        # Match channels.
        #
        # This coil data appears to connect to the measurement data via the
        # channel header channel_id field, when the relation
        #
        #     channel_id-1 == adc_channel_connected
        #
        # holds. Why the -1 is here is a mystery - perhaps the channel_id field
        # uses 1-based indexing.
        [Symbol(expt.coils[findfirst(e->e.adc_channel_connected-1 == c.channel_id, expt.coils)].element)
         for c in acq.channel_info]
    end
    channel_ids = [Int(c.channel_id) for c in acq.channel_info]
    AxisArray(z, Axis{:time}(t), Axis{:channel}(coilsyms))
end

"""
    count_cycles(twix::MRExperiment)

Get phase cycle length for Siemens twix.
"""
function count_cycles(twix::MRExperiment)
    num_averages = twix.metadata["lAverages"]
    phase_cycle_types = OrderedDict(0x01 => 1, # PHASE_CYCLING_NONE
                                    0x04 => 2, # PHASE_CYCLING_TWOSTEP
                                    0x08 => 8, # PHASE_CYCLING_EIGHTSTEP
                                    0x10 => 4, # PHASE_CYCLING_EXORCYCLE
                                    0x20 => 16)# PHASE_CYCLING_SIXTEENSTEP_EXOR
    if haskey(phase_cycle_types, twix.metadata["sSpecPara.lPhaseCyclingType"])
        n_cycles = phase_cycle_types[twix.metadata["sSpecPara.lPhaseCyclingType"]]
    elseif twix.metadata["sSpecPara.lPhaseCyclingType"] == 0x02 # PHASE_CYCLING_AUTO (Default)
        # As per the source code, this mode sets the cycling automatically, which can be computed
        # by longest phase cycle possible from the other phase cycling types, starting from 2 steps
        n_cycles = maximum(c for c in values(phase_cycle_types) if num_averages % c == 0)
    end
    n_cycles
end

"""
    mr_load(twix::MRExperiment, repair=false)

If `repair=true`, attempt to repair the experiment by discarding a subset of
broken acquisitions where possible.

Recognizes spectro experiments in Siemens twix format, producing one of:
  * PRESS (Siemens product sequence SVS_SE)
  * L-COSY (Customer sequences srcosy, svs_lcosy)
"""
function mr_load(twix::MRExperiment; repair=false)
    @debug "Recognizing twix input" twix

    meta = standard_metadata(twix)
    @debug "Found standard metadata" meta

    AUTO_REF_SCAN_OFF = 1
    is_svs_lcosy = occursin("svs_lcosy", meta.sequence_name)
    is_srcosy    = occursin("srcosy", meta.sequence_name) ||
                   occursin("sr_cosy", meta.sequence_name)
    is_svs_cosy  = occursin("svs_cosy", meta.sequence_name)
    is_press     = "%SiemensSeq%\\svs_se" == meta.sequence_name
    num_averages = twix.metadata["lAverages"]
    if is_svs_lcosy || is_srcosy || is_svs_cosy
        release_versions = ["%CustomerSeq%\\srcosy",
                            "%CustomerSeq%\\sr_cosy",
                            "%CustomerSeq%\\svs_cosy",
                            "%CustomerSeq%\\svs_lcosy-1"]
        if !(meta.sequence_name in release_versions)
            @warn """
                  Sequence `$(meta.sequence_name)` looks like an L-COSY sequence,
                  but is not one of the versions I recognize. We will continue
                  trying to read the data but if it's a development version of
                  the sequence you might get strange results.
                  """  meta.sequence_name twix
        end

        t1_inc_loop_idx = 0
        avg_loop_idx = 0

        if is_srcosy
            old_srcosy = occursin("sr_cosy", meta.sequence_name) ||
                         occursin("N4_VD", meta.software_version)
            # Uuugh. Number of t1 increments must be inferred from repetitions
            # loop counter.
            t1_inc_loop_idx = loop_counter_index(:repetition)
            avg_loop_idx = loop_counter_index(:acquisition)
            nsamp_t1 = Int(maximum(d.loop_counters.repetition for d in twix.data)) + 1
            dt1_key = "sWipMemBlock.adFree[1]"
            if !haskey(twix.metadata, dt1_key)
                # Guuuuh! It gets worse. It seems that likely that for this
                # very old sr_cosy sequence data we have no `dt1` value in the
                # WIP mem block. Probably it was just hardcoded so we will have
                # to guess.
                # TODO: Check this!!!
                dt1 = 0.8u"ms"
                #if !old_srcosy # TODO.
                    @warn """
                          No T1 increment found in twix metadata. We are guessing a value of `$dt1`.
                          """ meta
                #end
            else
                dt1 = twix.metadata[dt1_key]*u"ms"
            end
        elseif is_svs_cosy
            # Uuugh. Number of t1 increments must be inferred from acquisitions
            # loop counter.
            t1_inc_loop_idx = loop_counter_index(:partition)
            avg_loop_idx = loop_counter_index(:acquisition)
            nsamp_t1 = Int(maximum(d.loop_counters.partition for d in twix.data)) + 1
            dt1_key = "sWipMemBlock.adFree[1]"
            dt1 = twix.metadata[dt1_key]*u"ms"
        elseif is_svs_lcosy
            nsamp_t1 = twix.metadata["sWipMemBlock.alFree[3]"]
            t1_inc_loop_idx = loop_counter_index(:partition)
            avg_loop_idx = loop_counter_index(:acquisition)
            # TODO - stash legacy_timing somewhere
            legacy_timing = get(twix.metadata, "sWipMemBlock.alFree[5]", 0) == 1
            dt1           = twix.metadata["sWipMemBlock.adFree[1]"]*u"ms"
        end

        TE = get(twix.metadata,"alTE[0]", 0)u"μs"
        if TE != 30000u"μs"
            @warn "Found unexpected TE=$TE in L-COSY sequence"
        end

        has_refscans = get(twix.metadata, "sSpecPara.lAutoRefScanMode", 1) != AUTO_REF_SCAN_OFF

        ref_scans = Int[]
        lcosy_scans = zeros(Int, num_averages, nsamp_t1)

        did_loop_warning = false
        for (i,acq) in enumerate(twix.data)
            if !has_refscans && acq.loop_counters.phase != 0
                @warn """
                      Phase loop is nonzero but AutoRefScanMode is OFF.
                      I don't know what to do with this acquisition!
                      """
                continue
            end
            is_refscan = has_refscans && acq.loop_counters.phase == 0
            if is_refscan
                push!(ref_scans, i)
            end
            t1_index = acq.loop_counters[t1_inc_loop_idx] + 1
            avg_index = acq.loop_counters[avg_loop_idx] + 1
            if avg_index == 1 && acq.loop_counters.set > 0
                avg_index = acq.loop_counters.set + 1
                if !did_loop_warning
                    # There's two possibilities here.
                    # 1) VD13A and older data have a different ABI for twix,
                    #    with the loop order switched.
                    # 2) `set` was used for spectro averaging in the older software.
                    # TODO: Figure out which one of these it was, and remove this warning!
                    @warn """
                          We're guessing that the averaging loop is using the
                          `set` loop counter for this twix version.
                          """ metadata=meta
                    did_loop_warning = true
                end
            end
            lcosy_scans[avg_index, t1_index] = i
        end

        if any(lcosy_scans .== 0)
            bad_t1_incs = vec(any(lcosy_scans .== 0, dims=1))
            if repair && AcquisitionsIncomplete in twix.quality_control
                first_bad_avg = findfirst(bad_t1_incs)
                nsamp_t1 = first_bad_avg-1
                @warn """Attempting to repair LCOSY by discarding trailing T1 increment(s):
                          lcosy_scans[:, $first_bad_avg:$(lastindex(lcosy_scans,2))]
                      """
                lcosy_scans = lcosy_scans[:,1:first_bad_avg-1]
            else
                nmissing = length(findall(lcosy_scans .== 0))
                error("Missing increments or averages in L-COSY scan ($nmissing of $(length(lcosy_scans)))")
            end
        end

        t1 = (0:nsamp_t1-1)*dt1
        return LCOSY(t1, meta, ref_scans, lcosy_scans, twix)
    elseif is_press
        n_cycles = count_cycles(twix)
        TE = get(twix.metadata,"alTE[0]", missing)u"μs"

        if ismissing(TE)
            throw(MissingError("TE not found in PRESS sequence"))
        end
        # Extract ref_scans and press_scans
        has_refscans = get(twix.metadata, "sSpecPara.lAutoRefScanMode", 1) != AUTO_REF_SCAN_OFF
        if !has_refscans
            @warn """
                  The auto ref scan mode is turned off.
                  This feature has not been tested, please use with caution
                  """
        end

        ref_scans = Int[]
        press_scans = zeros(Int, Int(num_averages/n_cycles), n_cycles)
        navigator = Int[]

        # Loop through twix data and check the phase loop counter
        # if phase loop counter = 0, it's a ref scan
        # else if phase loop counter = 1, it's water suppressed scan
        for (i,acq) in enumerate(twix.data)
            if !has_refscans && acq.loop_counters.phase != 0
                @warn """
                      Phase loop is nonzero but AutoRefScanMode is OFF.
                      I don't know what to do with this acquisition!
                      """
                continue
            end
            is_refscan = has_refscans && acq.loop_counters.phase == 0
            if acq.rt_feedback
                push!(navigator, i)
                continue
            elseif is_refscan
                push!(ref_scans, i)
                continue
            end

            avg_index           = acq.loop_counters.acquisition % n_cycles + 1
            cycle_no            = Int(floor(acq.loop_counters.acquisition/n_cycles) + 1)
            press_scans[cycle_no, avg_index] = i
        end
        return PRESS(meta, ref_scans, press_scans, navigator, twix)
    end

    error("Siemens raw data with sequence \"$(meta.sequence_name)\" is unrecognized")
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
function xmlASCII2file(xml_string::String, file)
    xdoc = parse_string(xml_string)
    save_file(xdoc, file)
end

ext: jl


text:
#Getting lifetimes of 87Rb from Branden et al. (2009)

using JLD

nns = 28:45;
nnp = 34:44;
nnd = 29:44;

sexpt = [15.2, 14.9, 18.9, 17.9, 18.1, 20.1, 23.8, 24.4, 25.7, 26.8, 28.6, 34.3, 34.2, 38.1, 40.9, 40.1, 47, 50.8];
sexpterr = [0.7, 0.8, 0.8, 0.8, 0.5, 0.8, 1.7, 1.0, 1.4, 1.0, 1.1, 1.5, 0.9, 2.1, 1.2, 1.6, 1.7, 2.4];
pexpt = [33.2, 33.7, 39.8, 43.5, 45.1, 45.9, 52.6, 58.7, 60.3, 60.2, 64.2];
pexpterr = [1.4, 1.5, 1.3, 1.9, 2.0, 1.3, 1.3, 3.0, 1.7, 2.8, 2.6];
dexpt = [13.5, 15.8, 19.5, 19.9, 20.5, 24.1, 25.4, 28.6, 31.5, 32.3, 36.4, 35.0, 36.9, 38.7, 38.9, 42];
dexpterr = [0.7, 0.5, 1.5, 1.2, 2.3, 0.6, 1.8, 1.6, 0.8, 1.2, 1.4, 1.0, 1.8, 2.4, 3.0, 3.2];

save("Lifetimes.jld","sexpt",sexpt,"sexpterr",sexpterr,"pexpt",pexpt,"pexpterr",pexpterr,"dexpt",dexpt,"dexpterr",dexpterr);

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
@testset "Wellbore object" begin

md_bad = [-1.,2.,3.,4.]
md_good = [1.,2.,3.,4.]
inc = [1.,2.,3.,4.]
tvd_bad = [-1.,2.,3.,4.]
tvd_good = [1.,2.,3.,4.]
id = [1.,1.,1.,1.]

#implicit test for adding the leading 0,0 survey point:
w = Wellbore(md_good, inc, tvd_good, id)

@test w.id[1] == w.id[2]

#implicit test for allowing negatives:
Wellbore(md_bad, inc, tvd_bad, id, true)

try
    Wellbore(md_bad, inc, tvd_good, id)
catch e
    @test e isa Exception
end

try
    Wellbore(md_good, inc, tvd_bad, id)
catch e
    @test e isa Exception
end

end #testset for Wellbore object

@testset "Wellbore with valves" begin

    md = [1.,2,3,4]
    tvd = [2.,4,5,6]
    inc = [0.,2,3,4]
    id = [1.,1,1,1]
    starting_length = length(md)

    valve_md = [1.5,3.5]
    PTRO = [1000, 900]
    R = [0.07,0.07]
    ports = [16,16]
    valves = GasliftValves(valve_md, PTRO, R, ports)

    well = Wellbore(md, inc, tvd, id, valves)

    @test all(length.([md,tvd,inc,id]) .== starting_length)
    @test well.md == [0,1,1.5,2,3,3.5,4]
    @test well.tvd == [0,2.,3,4,5,5.5,6]
    @test well.inc == [0,0.,1,2,3,3.5,4]
    @test well.id == [1,1.,1,1,1,1,1]

end


text:
println("This is an example.")

x = 3
y = x^2

max_stars_repo_name: dpsanders/RepoBrowser

max_issues_repo_name: dpsanders/RepoBrowser

max_forks_repo_name: dpsanders/RepoBrowser


text:
"""
    growth_rate(x,q,[args])

Computes the growth rate of the system, i.e. the highest imaginary part of the physical eigenvales.
Empties the tables during the run.
"""
function growth_rate(x::Float64=0.0, q::Float64=1.0,
    tabOmega=tabOmega_serial, tabKappaSqOverTwoOmega=tabKappaSqOverTwoOmega_serial,
        tabAln=tabAln_serial, tabBln=tabBln_serial, tabCln=tabCln_serial,
        tabDln=tabDln_serial, tabFln=tabFln_serial, tabGln=tabGln_serial,
        tabHln=tabHln_serial, tabTruncMln=tabTruncMln_serial,
        tabEigValsMln=tabEigValsMln_serial)
    table_function_fill!(x,q,tabOmega,tabKappaSqOverTwoOmega)
    matrix_fill!(x,q,tabAln,tabBln,tabCln,tabDln,tabFln,tabGln,tabHln,tabOmega,tabKappaSqOverTwoOmega)
    tabTruncMln!(tabTruncMln,tabAln,tabBln,tabCln,tabDln,tabFln,tabGln,tabHln)
    tabEigValsMln!(tabTruncMln,tabEigValsMln)
    physical_eig = getPhysicalEigvals(tabEigValsMln)
    kEgv, egv = get_max(physical_eig)

    # clear the temp table
    tabEigValsMln_clear!(tabEigValsMln)
    tabTruncMln_clear!(tabTruncMln)
    matrix_clear!(tabAln,tabBln,tabCln,tabDln,tabFln,tabGln,tabHln)
    table_function_clear!(tabOmega,tabKappaSqOverTwoOmega)

    # return maximum growth rate
    return imag(egv)
end

"""
    growth_rate_with_rotation(x,q,[args])

Computes the physical eigenvalue with the highest imaginary part.
Empties the tables during the run.
"""
function growth_rate_with_rotation(x::Float64=0.0, q::Float64=1.0,
    tabOmega=tabOmega_serial, tabKappaSqOverTwoOmega=tabKappaSqOverTwoOmega_serial,
        tabAln=tabAln_serial, tabBln=tabBln_serial, tabCln=tabCln_serial,
        tabDln=tabDln_serial, tabFln=tabFln_serial, tabGln=tabGln_serial,
        tabHln=tabHln_serial, tabTruncMln=tabTruncMln_serial,
        tabEigValsMln=tabEigValsMln_serial)
    table_function_fill!(x,q,tabOmega,tabKappaSqOverTwoOmega)
    matrix_fill!(x,q,tabAln,tabBln,tabCln,tabDln,tabFln,tabGln,tabHln,tabOmega,tabKappaSqOverTwoOmega)
    tabTruncMln!(tabTruncMln,tabAln,tabBln,tabCln,tabDln,tabFln,tabGln,tabHln)
    tabEigValsMln!(tabTruncMln,tabEigValsMln)
    physical_eig = getPhysicalEigvals(tabEigValsMln)
    kEgv, egv = get_max(physical_eig)

    # clear the temp table
    tabEigValsMln_clear!(tabEigValsMln)
    tabTruncMln_clear!(tabTruncMln)
    matrix_clear!(tabAln,tabBln,tabCln,tabDln,tabFln,tabGln,tabHln)
    table_function_clear!(tabOmega,tabKappaSqOverTwoOmega)

    # return maximum growth rate
    return real(egv), imag(egv)
end

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
using NaNStatistics
using Test, StatsBase, Statistics, LinearAlgebra

@testset "ArrayStats" begin include("testArrayStats.jl") end
@testset "Covariance and Correlation" begin include("testCovCor.jl") end
@testset "Sorting" begin include("testSorting.jl") end

@testset "Histograms" begin include("testHistograms.jl") end
@testset "Binning" begin include("testBinning.jl") end


text:
using FactCheck
include(Pkg.dir("PseudoGenomes", "src","pseudogenome-snp.jl"))

ref_fasta = joinpath( Pkg.dir(),"PseudoGenomes","testdata","test.fasta")
vcf_path = joinpath( Pkg.dir(),"PseudoGenomes","testdata", "vcf-file2.vcf")
vcf_noheader_path = joinpath( Pkg.dir(),"PseudoGenomes","testdata", "vcf-file2.noheader.vcf")
outpath = joinpath( Pkg.dir(),"PseudoGenomes","testdata","out.fasta")
cpgbedfile1 = joinpath( Pkg.dir(),"PseudoGenomes","testdata","cpgs1.bedgraph")
#cpgbedfile2 = joinpath( Pkg.dir(),"PseudoGenomes","testdata","cpgs2.bedgraph")
# - test read_vcf_columns
read_snp_positions_from_columns(vcf_path, gzip=false )

# - test substitute_
expected_path = joinpath( Pkg.dir(),"PseudoGenomes","testdata","expected.fasta")
substitute_snps_in_fasta_file(ref_fasta,vcf_path,outpath)

expected =  readall(expected_path)
got = readall(outpath)

facts("fasta-check") do
    @fact got => expected
end

# -  test cpg
a = ['A','C','C','C','G','T','G','C','G','A']
cgpos = find_cpgs(a)
# default zero based
expected_cgpos = [3,7]
facts("cpg-check") do
    @fact cgpos => expected_cgpos
end

# - make cpg bedfile
#make_cpg_bedfile_type1(outpath,cpgbedfile1)
make_cg_bedfile(outpath,cpgbedfile1)

# - test check variant alleles
# - switch logging to file and check file is correct
check_variant_alleles_in_fasta_file(outpath,vcf_path)

# -
# -
check_variant_alleles_in_fasta_file(outpath,vcf_noheader_path)


text:
# DO NOT MODIFY THIS CODE

function make_g1()
  g = fill(0, 63, 63)
  for i in 1:21
    for j in 22:63
      g[i, j] = 1
      g[j, i] = 1
    end
  end
  return g
end

function make_g2()
  g = fill(0, 63, 63)
  g[1,62] = 1
  g[62,1] = 1
  for i in 1:62
    g[i,63] = 1
    g[63,i] = 1
    if i - 1 > 0
      g[i, i-1] = 1
      g[i - 1, i] = 1
    end
  end
  return g
end

function make_g3()
  g = fill(0, 63, 63)
  for i in 1:31
    g[i, 2 * i] = 1
    g[2 * i, i] = 1
    g[i, 2 * i + 1] = 1
    g[2 * i + 1, i] = 1
  end
  return g
end

g1 = make_g1()
g2 = make_g2()
g3 = make_g3()

function fit(g, x)
  if length(x) == 1 # clever! to save time, return 1 if length of x is 1 as it's obviously an independent set so we don't have to perform the checks below
    return 1
  end

  #traverse the array that represents our solution and check on the adjacency matrix if any of the nodes are connected
  #if so, our fitness is zero as the set is not a solution
  for i in 1:length(x) - 1
    for j in i:length(x)
      if g[x[i], x[j]] == 1 || g[x[j], x[i]] == 1
        return 0
      end
    end
  end

  return length(x) #if execution reached this point, array representing the set is an independent set, so return the length of the array as fitness
end

function gen()
  return [rand(1:63)] # sweet and elegant function generating an array with a single randomly generated vertex
end

function mut(x, s)
  rate = Int(round(5 * s)) # mutation rate goes down proportionally to the number of generations executed, approximating 0 and defining how many vertices will be added to the array that represents a solution set
  new = fill(0, rate)
  filled = 0
  while filled < rate
    r = rand(1:63) #gets random node
    if r in new || r in x # if node already exists in the received set or in the new one, ignore node and continue
      continue
    end
    new[filled + 1] = r #if node is not in the set, add to it
    filled += 1 #variable that controls the loop
  end
  return cat(1, x, new) #after defining vertices that will be added, concatenate with original set and return
end

# END OF UNMODIFIABLE SECTION

# You can add more functions, but you can only change the ev method
# Do NOT modify p or g!

function ev(graph)
  p = 10 # DO NOT CHANGE!
  g = 100 # DO NOT CHANGE!
  initialSelectionSize = 5 # number of individuals that will be selected to be mutated in each generation (initial selection)
                          # initial selection strategy = random
  finalSelectionSize = 5  # number of individuals that will be selected to be compared against a mutated individual

  # initializing population with empty items
  pop = fill([0], p)

  # populate with randomly generated individuals: for this algorithm, we've defined that our initial individuals will be
  # sets with only one vertex randomly chosen (rand(1:63)) -- trivial graph
  for i in 1:p
    pop[i] = gen()
  end

  # run for "g" generations
  for generation in 1:g

    # select a predefined number of individuals to be mutated in each generation (initial selection)
    # for each mutated individual, get a random opponent, make them fight and
    for count in 1:initialSelectionSize
      x = rand(1:p) # random element to be mutated
      x_mutant = mut(pop[x], ( (g - generation + 1.0) / g ) ) # random element is mutated based on a mutation rate

      for count2 in 1:finalSelectionSize
        y = rand(1:p) # another random individual that will fight the mutated individual
        parent1 = rand(1:p) # randomly selected parent 1
        parent2 = rand(1:p) # randomly selected parent 2

        child = cross(pop[parent1], pop[parent2]) # naughty parents have sex and make a baby

        # compare mutant against the other unlucky randomly selected individual
        # if mutation did not produce a fitter individual than this other, mutant is going to be discarded and the
        # other random individual lives.
        if fit(graph, x_mutant) > fit(graph, pop[y])
          pop[y] = x_mutant
          break
        end

        # after crossover, replaces parent with resulting child if its fitter
        if fit(graph, child) > fit(graph, pop[parent1])
          pop[parent1] = child
        elseif fit(graph, child) > fit(graph, pop[parent2])
          pop[parent2] = child
        end
      end #end for

    end #end for

  end #end ev function

  return pop # return a population the same size as the initial one, but beautiful, fit, healthy and hopefully converging to the optimal solution =]

end

# reproduction: get half of each set and combine
function cross(parent1, parent2)
  halfParent1Index = Int(round(length(parent1)/2))
  halfParent2Index = Int(round(length(parent2)/2))

  if rand(0:1) === 0
    return vcat(parent1[1:halfParent1Index], parent2[halfParent2Index:length(parent2)])
  else
    return vcat(parent2[1:halfParent2Index], parent1[halfParent1Index:length(parent1)])
  end
end

# DO NOT MODIFY!

function do_and_print(graph, which)
  answer = ev(graph) # run evolutionary algorithm for g1,g2 or g3 depending on what we're receiving
  best = answer[1] # answer = evolved population. get the first individual and put it in the "best" variable
  best_fit = fit(graph, answer[1]) # calculate the fit of the individual above and store it in the "best_fit" variable
  @printf("Results for graph %d\n", which) # Results for graph 1, 2 or 3

  # for each individual in the returned population that has been put through evolution, print
  # the number of the answer corresponding to its position in the population, the array itself and its fit
  # e.g.: 1: [1,23,5,8], fitness: 4
  # find the best fit in the population and assign the individual to "best" and the fit to "best_fit"
  for i in 1:length(answer)
    @printf("%d: %s, fitness %d\n", i, answer[i], fit(graph, answer[i]))
    if (best_fit < fit(graph, answer[i]))
      best = answer[i]
      best_fit = fit(graph, answer[i])
    end
  end

  # print best answer alongside its fitness
  # e.g.: Best: 1, fitness: 16
  @printf("Best: %s, fitness %d\n", best, best_fit)

  # print "End of results for graph 1, 2 or 3"
  @printf("End of results for graph %d\n\n", which)
end

# execution starts here
#do_and_print(g1, 1)
#do_and_print(g2, 2)
#do_and_print(g3, 3)

# END OF UNMODIFIABLE SECTION

max_stars_repo_name: kathleen-saz/Independent-set-EA

max_issues_repo_name: kathleen-saz/Independent-set-EA

max_forks_repo_name: kathleen-saz/Independent-set-EA


text:
using Random
using MLJ
using DataFrames
using LinearAlgebra

"""
Split a given Dataframe into 2 sub-dataframe allowing
for training/testing dataset split
"""
function train_test_split(df::DF; split=0.8, seed=42, shuffled=true) where {DF <: Union{DataFrame, SubDataFrame}}
  Random.seed!(seed)
  (nr, nc) = size(df)
  nrp = round(Int, nr * split)
  row_ixes = shuffled ? shuffle(1:nr) : collect(1:nr)

  df_train = view(df[row_ixes, :], 1:nrp, 1:nc)
  df_test = view(df[row_ixes, :], nrp+1:nr, 1:nc)

  (df_train, df_test)
end

"""
From a Dataframe to matrix
from week2
"""
function get_data(df::DF, features, output;
                  colname=:constant) where {DF <: Union{DataFrame, SubDataFrame}}
  df[:, colname] .= 1.0
  features = [colname, features...]
  X_matrix = convert(Matrix, select(df, features)) # to get a matrix
  y = df[!, output]                                # => to get a vector
  return (X_matrix, y)
end

"""
Assume feature_matrix is a matrix containing the features as columns
and weights is a corresponding array

from week2
"""
function predict_output(X::Matrix{T}, weights::Vector{T}) where {T <: Real}
    X * weights
end

"""
Calculate RSS given weights vector

depends on predict_output
from week2

y::SubArray{Float64,1, Array{Float64,1}, Tuple{UnitRange{Int64}}, true}
"""
function calc_rss(X::Matrix{T}, y, weights::Vector{T}) where {T <: Real}
  preds = predict_output(X, weights)
  sum((preds .- y) .^ 2)
end

"""
Calculate RSS given a (MLJ) machine which is a wrapper combining a model
and its data.

depends on MLJ predict()

mach::MLJBase.Machine{MLJLinearModels.RidgeRegressor,true}
"""
function get_rss(mach::Machine, X::DF, y::Vector{T}) where {DF <: Union{DataFrame, SubDataFrame}, T <: Real}
  ŷ = predict(mach, X)     # First get the predictions
  diff = y .- ŷ            # Then compute the residuals/errors
  rss = sum(diff .* diff)  # Then square and add them up
  return rss
end

function normalize_features(X::Matrix)
  n = size(X)[2]
  norms = zeros(eltype(X), n)'
  for ix ∈ 1:n
    norms[ix] = norm(X[:, ix])
  end
  (X ./ norms, norms)
end

"""
Select a random sample percentage of an index range
"""
function sampling(range::Integer, perc::Float64; seed=42)
  up_lim = ceil(Integer, perc * range)
  randperm(MersenneTwister(seed), range)[1:up_lim]
end

"""
Given a DF and the features to encode, do on hot encoding of each
feature

on-place transformation
"""
function hot_encode!(df::DF;
                     features=Features) where {DF <: Union{DataFrame, SubDataFrame}}
  for f ∈ features
    ## get unique value over column feature
    fval = reduce(vcat, df[!, f]) |> unique

    ## apply hotencoding inplace
    ## cf. https://stackoverflow.com/questions/64565276/julia-dataframes-how-to-do-one-hot-encoding
    transform!(df,
      f .=> [ByRow(v -> occursin(x, v) ? 1 : 0) for x ∈ fval] .=> Symbol.(f, "_", fval))

    ## get rid of initial feature (column)
    select!(df, Not(f))
  end
end

"""
Condition the data, adjusting to right MLG scientific type

features, target of type F_ST (Feature SciType)
"""
function data_prep_Xy(df::DF, features, target) where {DF <: Union{DataFrame, SubDataFrame}}

  y_target = target.names[1]

  X = select(df, features.names)
  coerce!(X, features.st_map...)

  y = select(df, y_target)
  coerce!(y, target.st_map...)
  y = y[:, y_target]   ## cast to vector

  return (X, y)
end

function data_prep(df::DF, features) where {DF <: Union{DataFrame, SubDataFrame}}
  coerce(df, features.st_map...)
end


text:
# TODO: handling singularity

const quon_atol = 1e-8

function change_direction(α::Complex)
    is_zero(α) && return InfZero
    is_pi(α) && return InfPi
    return to_quon_const(log(-tanh(α/2)))
end

"""
    yang_baxter_param(α1, β1, γ1)

Return α2, β2, γ2 according to the Yang-Baxter equation:

    \\    /    |       |    \\    /
     \\  /     |       |     \\  /
      α1      |       |      α2
     |  \\    /         \\    /  |
     |   \\  /           \\  /   |
     |    β1      =      β2    |
     |   /  \\           /  \\   |
     |  /    \\         /    \\  |
      γ1      |       |      γ1
     /  \\     |       |     /  \\
    /    \\    |       |    /    \\
"""
function yang_baxter_param(α1::Complex, β1::Complex, γ1::Complex)
    try
        if all(p -> (p isa QuonConst), to_quon_const.([α1, β1, γ1]))
            return yang_baxter_param(to_quon_const(α1), to_quon_const(β1), to_quon_const(γ1))
        end
        s1 = (α1+γ1)/2
        m1 = (α1-γ1)/2
        if !isapprox(cosh(s1), 0; atol = quon_atol)
            tanh_s2 = tanh(β1/2)*cosh(m1)/cosh(s1)
            @assert !(tanh_s2 ≈ 1)
            @assert !(tanh_s2 ≈ -1)
            s2 = atanh(tanh_s2)
        elseif !isapprox(cosh(m1)*sinh(β1), 0; atol = quon_atol)
            s2 = π*im/2
        else
            error("This case should not be handled here")
        end
        if !isapprox(sinh(s1), 0; atol = quon_atol)
            tanh_m2 = -tanh(β1/2)*sinh(m1)/sinh(s1)
            @assert !(tanh_m2 ≈ 1)
            @assert !(tanh_m2 ≈ -1)
            m2 = atanh(tanh_m2)
        elseif !isapprox(sinh(m1)*sinh(β1), 0; atol = quon_atol)
            m2 = π*im/2
        else
            error("This case should not be handled here")
        end

        α2 = s2+m2
        γ2 = s2-m2
        if !isapprox(cosh(m2), 0; atol = quon_atol) && !isapprox(cosh(s1), 0; atol = quon_atol)
            β2 = 2*atanh(tanh(s1)*cosh(s2)/cosh(m2))
        elseif !isapprox(sinh(m2), 0; atol = quon_atol) && !isapprox(cosh(m1), 0; atol = quon_atol)
            β2 = 2*atanh(tanh(m1)*sinh(s2)/sinh(-m2))
        else
            error("This case should not be handled here")
        end

        return to_quon_const(α2), to_quon_const(β2), to_quon_const(γ2)
    catch e
        error("$((α1, β1, γ1)) is a singular point of Yang-Baxter equation")
    end
end

function yang_baxter_param_inv(α2, β2, γ2)
    a2 = change_direction(α2)
    b2 = change_direction(β2)
    α1, c1, b1 = yang_baxter_param(b2, a2, γ2)
    β1 = change_direction(b1)
    γ1 = change_direction(c1)
    return α1, β1, γ1
end


text:
module EnergySystemModeling

include("model.jl")
include("io.jl")
export EnergySystemModel,
    Specs,
    Params,
    JuMPObj,
    JuMPVar,
    retrieve_data,
    Expressions,
    equivalent_annual_cost,
    save_json,
    load_json,
    create_nodedata,
    replace_nans,
    getdispatch,
    read_clusters,
    unpack_clusters_features

include("plotting.jl")
export plot_objective_values,
    plot_generation_dispatch,
    plot_generation_capacities,
    plot_generation_capacities_stacked,
    plot_transmission_flow,
    plot_transmission_capacities,
    plot_transmission_bars,
    plot_storage_level,
    plot_storage_capacities,
    plot_loss_of_load,
    plot_box,
    plot_box_all,
    plot_dispatch_bars

include("aggreg.jl")
export SeriesInstance,
    ClustInstance,
    AggregInstance,
    DistUpdate,
    load_series_instance,
    load_clust_instance,
    aggreg1D,
    cdad,
    search_min_dist,
    compute_dist,
    update_marker,
    replace_lines,
    update_clust!,
    update_k!,
    find_clusters!
end # module


text:
### HssMatrices.jl module
# A simple Julia package that allows working with HSS matrices
# The package aims to be simple, intuitive and efficient
# Written by Boris Bonev, Nov. 2020
__precompile__()
module HssMatrices

  # dependencies, trying to keep this list to a minimum if possible
  using LinearAlgebra
  using SparseArrays # introduce custom constructors from sparse matrices
  using AbstractTrees
  using Plots # in the future, move to RecipesBase instead
  using LowRankApprox # optional for now , to replace prrqr with an efficient implementation

  # load BLAS/LAPACK routines only used within the library# load efficient BLAS and LAPACK routines for factorizations
  import LinearAlgebra.LAPACK.geqlf!
  import LinearAlgebra.LAPACK.gelqf!
  import LinearAlgebra.LAPACK.ormql!
  import LinearAlgebra.LAPACK.ormlq!
  import LinearAlgebra.BLAS.ger!
  import LinearAlgebra.BLAS.trsm
  import LinearAlgebra.ishermitian

  # using InvertedIndices, DataStructures
  import Base: +, -, *, \, /, ^, ==, !=, copy, size, getindex, adjoint, transpose, convert, Matrix
  import LinearAlgebra: ldiv!, rdiv!, mul!

  # HssMatrices.jl
  export HssOptions
  # binarytree.jl
  export BinaryNode, leftchild, rightchild, isleaf, isbranch, depth, descendants, compatible, nleaves, prune_leaves!
  # clustertree.jl
  export bisection_cluster, cluster
  # linearmap.jl
  export LinearMap, HermitianLinearMap
  # hssmatrix.jl
  export HssMatrix, hss, isleaf, isbranch, ishss, hssrank, full, checkdims, prune_leaves!, blkdiagm
  # compression.jl
  export compress, randcompress, randcompress_adaptive, recompress!, hss_blkdiag
  # generators.jl
  export generators, orthonormalize_generators!
  # matmul.jl
  # ulvfactor.jl
  export ulvfactsolve
  # ulvdivide.jl
  export ldiv!, rdiv!
  # constructors.jl
  export lowrank2hss
  # visualization.jl
  export plotranks, pcolor

  mutable struct HssOptions
    atol::Float64
    rtol::Float64
    leafsize::Int
    noversampling::Int
    stepsize::Int
    verbose::Bool
  end

  # set default values
  atol = 1e-9
  rtol = 1e-9
  leafsize = 64
  noversampling = 10
  stepsize = 20
  verbose = false

  # struct to pass around with options
  function HssOptions(::Type{T}; args...) where T
    opts = HssOptions(atol, rtol, leafsize, noversampling, stepsize, verbose)
    for (key, value) in args
      setfield!(opts, key, value)
    end
    opts
  end
  HssOptions(; args...) = HssOptions(Float64; args...)

  function copy(opts::HssOptions; args...)
    opts_ = HssOptions()
    for field in fieldnames(typeof(opts))
      setfield!(opts_, field, getfield(opts, field))
    end
    for (key, value) in args
      setfield!(opts_, key, value)
    end
    opts_
  end

  function chkopts!(opts::HssOptions)
    opts.atol ≥ 0. || throw(ArgumentError("atol"))
    opts.rtol ≥ 0. || throw(ArgumentError("rtol"))
    opts.leafsize ≥ 1 || throw(ArgumentError("leafsize"))
    opts.stepsize ≥ 1 || throw(ArgumentError("stepsize"))
    opts.noversampling ≥ 1 || throw(ArgumentError("noversampling"))
  end

  function setopts(opts=HssOptions(Float64); args...)
    opts = copy(opts; args...)
    chkopts!(opts)
    global atol = opts.atol
    global rtol = opts.rtol
    global leafsize = opts.leafsize
    global noversampling = opts.noversampling
    global stepsize = opts.stepsize
    global verbose = opts.verbose
    return opts
  end

  include("binarytree.jl")
  include("clustertree.jl")
  include("linearmap.jl")
  include("hssmatrix.jl")
  include("compression.jl")
  include("generators.jl")
  include("matmul.jl")
  include("ulvfactor.jl")
  include("ulvdivide.jl")
  include("constructors.jl")
  include("visualization.jl")
end


text:
#
# Test script for convex optimization solver
#
# author: Atsushi Sakai
#

using Test
using ArgMin

function test_solve_linear_programming_with_simplex_method1()
    println("test_solve_linear_programming_with_simplex_method1")

	c = [2.0 3.0]'

	G = [1.0 2.0;
		 1.0 1.0;
		 3.0 1.0]

	h = [14.0 8.0 18.0]'

	x_hat, status = solve_linear_programming_with_simplex_method(
			c, G=G, h=h)
	# display(x_hat)
	@test x_hat[1] ≈ 2.0 atol=0.01
	@test x_hat[2] ≈ 6.0 atol=0.01
	@test status["objectives"] ≈ 22.0 atol=0.01
end

function test_solve_linear_programming_with_simplex_method2()
    println("test_solve_linear_programming_with_simplex_method2")
	# A Blending Problem — PuLP 1.6.0 documentation https://pythonhosted.org/PuLP/CaseStudies/a_blending_problem.html

	c = [2.0 3.0]'
	# println("c")
	# display(c)
	# println("")
	G = [1.0 2.0;
		 2.0 1.0]
	# println("A")
	# display(A)
	# println("")
	h = [10.0 8.0]'
	# println("b")
	# display(b)
	# println("")

	x_hat, status = solve_linear_programming_with_simplex_method(
			c, G=G, h=h)
	# display(x_hat)
	@test x_hat[1] ≈ 2.0 atol=0.01
	@test x_hat[2] ≈ 4.0 atol=0.01
end

function test_solve_linear_programming_with_simplex_method3()
    println("test_solve_linear_programming_with_simplex_method3")
	# A Blending Problem — PuLP 1.6.0 documentation
    # https://pythonhosted.org/PuLP/CaseStudies/a_blending_problem.html

	c = [0.013 0.008 0.01 0.002 0.005 0.001]'
	# println("c")
	# display(c)
	# println("")
	G = [-0.1  -0.2 -0.15  0.0    -0.04 0.0;
         -0.08 -0.1 -0.11 -0.01   -0.01 0.0;
          0.001 0.005 0.003 0.1   0.15  0.0;
          0.002 0.005 0.007 0.002 0.008 0.0]
	# println("A")
	# display(A)
	# println("")
	h = [8.0 6.0 2.0 0.4]'
	# println("b")
	# display(b)
	# println("")

    A = [1.0 1.0 1.0 1.0 1.0 1.0]

    b = [100]

	x_hat, status = solve_linear_programming_with_simplex_method(
			c, G=G, h=h, A=A, b=b)

    display(x_hat)
	@test x_hat[1] ≈ 0.0 atol=0.01
	@test x_hat[2] ≈ 60.0 atol=0.01
	@test x_hat[3] ≈ 0.0 atol=0.01
	@test x_hat[4] ≈ 0.0 atol=0.01
	@test x_hat[5] ≈ 0.0 atol=0.01
	@test x_hat[6] ≈ 40.0 atol=0.01
end

function test_solve_quadratic_programming1()
    println("test_solve_quadratic_programming1")

	P = [1.0 0.0;
		 0.0 2.0]
	# display(P)
	q = [3.0,4.0]'

	x_hat = solve_quadratic_programming(P, q)
	# display(x_hat)
	@test x_hat[1] ≈ -3.0 atol=0.01
	@test x_hat[2] ≈ -2.0 atol=0.01
end

function test_solve_quadratic_programming2()
    println("test_solve_quadratic_programming2")

	P = [1.0 0.0;
		 0.0 0.0]
	# display(P)
	q = [3.0 4.0]
	# display(q)
	A = [1.0 1.0]
	# display(A)
	b = [1.0]
	# display(b)

	x_hat = solve_quadratic_programming(P, q, A, b)
	# display(x_hat)
	@test x_hat[1] ≈ 1.0 atol=0.01
	@test x_hat[2] ≈ 0.0 atol=0.01
end

# test_solve_linear_programming_with_simplex_method1()
# test_solve_linear_programming_with_simplex_method2()
# test_solve_linear_programming_with_simplex_method3()
test_solve_quadratic_programming1()
test_solve_quadratic_programming2()


text:
# KeysightConstants.jl

"""
Compatibility Options for function `open` in `SD_Module.jl`
- Legacy : Channel enumeration starts with CH0 --> COMPATIBILITY_LEGACY = 0
- Keysight : Channel enumeration starts with CH1 --> COMPATIBILITY_KEYSIGHT = 1
"""
const COMPATIBILITY_LEGACY              = Cint(0)
const COMPATIBILITY_KEYSIGHT            = Cint(1)

"""
Options for the input impedance setting for function `channelInputConfig` in
`SD_AIN.jl`
- High Impedance : Input impedance is high --> AIN_IMPEDANCE_HZ = 0
- 50Ω : Input impedance is 50Ω  --> AIN_IMPEDANCE_50 = 1
"""
const AIN_IMPEDANCE_HZ                  = Cint(0)
const AIN_IMPEDANCE_50                  = Cint(1)

"""
Options for the input coupling setting for function `channelInputConfig` in
`SD_AIN.jl`
- DC : DC coupling --> AIN_COUPLING_DC = 0
- AC : AC coupling --> AIN_COUPLING_AC = 1
"""
const AIN_COUPLING_DC                   = Cint(0)
const AIN_COUPLING_AC                   = Cint(1)

"""
Options for the analog trigger (parameter `analogTriggerMode` in function
`channelTriggerConfig` in 'SD_AIN.jl'
- Rising Edge : Trigger is generated when the input signal is rising and crosses
the threshold --> AIN_RISING_EDGE = 0
- Falling Edge : Trigger is generated when the input signal is falling and
crosses the threshold --> AIN_FALLING_EDGE = 1
- Both Edges : Trigger is generated when the input signal crosses the threshold,
no matter if it is rising or falling --> AIN_BOTH_EDGES = 3
"""
const AIN_RISING_EDGE                   = Cint(1)
const AIN_FALLING_EDGE                  = Cint(2)
const AIN_BOTH_EDGES                    = Cint(3)

"""
Options for trigger modes
### DAQ trigger modes
parameter `triggerMode` in function `DAQconfig` in `SD_AIN.jl`
- Auto (Immediate) : The acquisition starts automatically after a call to
function DAQstart --> AUTOTRIG = 0
- Software / HVI : Software trigger. The acquisition is triggered by the
function DAQtrigger, provided that the DAQ is running. DAQtrigger can be
executed from the user application (VI) or from an HVI --> SWHVITRIG = 1
- Hardware Digital Trigger : Hardware trigger. The DAQ waits for an external
digital trigger --> HWDIGTRIG = 2
- Hardware Analog Trigger : Hardware trigger. The DAQ waits for an external
analog trigger (only products with analog inputs) --> HWANATRIG = 3
### Trigger methods for the waveforms queued in an AWG
parameter `triggerMode` in function `AWG`, or in function `AWGqueueWaveform` in
`SD_AOU.jl`.
- Auto : The waveform is launched automatically after function AWGstart, or
when the previous waveform in the queue finishes --> AUTOTRIG = 0
- Software / HVI : Software trigger. The AWG is triggered by the function
AWGtrigger, provided that the AWG is running. AWGtrigger can be executed from
the user application (VI) or from an HVI --> SWHVITRIG = 1
- Software / HVI (per cycle) : Software trigger. Identical to the previous
option, but the trigger is required per each waveform cycle
--> SWHVITRIG_CYCLE = 5
- External Trigger : Hardware trigger. The AWG waits for an external trigger
--> EXTTRIG = 2
- External Trigger (per cycle) : Hardware trigger. Identical to the previous
option, but the trigger is required per each waveform cycle --> EXTTRIG_CYCLE = 6
"""
const AUTOTRIG                      = Cint(0)
const SWHVITRIG                     = Cint(1)
const HWDIGTRIG                     = Cint(2)
const HWANATRIG                     = Cint(3)
const EXTTRIG                       = Cint(2)
const SWHVITRIG_CYCLE               = Cint(5)
const EXTTRIG_CYCLE                 = Cint(6)

"""
External trigger source
### External trigger source for the DAQ
parameter `externalSource` in function `DAQdigitalTriggerConfig`
- External I/O Trigger : The DAQ trigger is a TRG connector/line of the product.
PXI form factor only: this trigger can be synchronized to CLK10
--> TRIG_EXTERNAL = 0
- PXI Trigger : PXI form factor only. The DAQ trigger is a PXI trigger line
and it is synchronized to CLK10 --> TRIG_PXI = 1
### External trigger source for the AWGs
parameter `externalSource` in function `AWGtriggerExternalConfig`
- External I/O Trigger : The AWG trigger is a TRG connector/line of the module
I/O Triggers. PXI form factor only: this trigger can be synchronized to CLK10
--> TRIG_EXTERNAL = 0
- PXI Trigger [0..n] : PXI form factor only. The AWG external trigger is a PXI
trigger line and it is synchronized to CLK10 --> TRIG_PXI_AWG(=4000) + Trigger No.
"""
const TRIG_EXTERNAL                 = Cint(0)
const TRIG_PXI                      = Cint(1)
const TRIG_PXI_AWG                  = Cint(4000)

"""
Hardware digital trigger behavior for the DAQ (parameter `triggerBehavior`
in function `DAQdigitalTriggerConfig` in 'SD_AIN.jl') or External trigger
behavior for the AWGs (parameter `triggerBehavior` in function
`AWGtriggerExternalConfig` in 'SD_AOU.jl')
- Active High : Trigger is active when it is at level high --> TRIG_HIGH = 1
- Active Low : Trigger is active when it is at level low --> TRIG_LOW = 2
- Rising Edge : Trigger is active on the rising edge --> TRIG_RISE = 3
- Falling Edge : Trigger is active on the falling edge --> TRIG_FALL = 4
"""
const TRIG_HIGH                     = Cint(1)
const TRIG_LOW                      = Cint(2)
const TRIG_RISE                     = Cint(3)
const TRIG_FALL                     = Cint(4)

"""
### Trigger I/O Options
M31/M33XXA Digitizers & M32/M33XXA AWGs TRG connectors configuration (parameter
`direction` in function `triggerIOconfig` in both 'SD_AIN.jl' and 'SD_AOU.jl'
- Trigger Output (readable) : TRG operates as a general purpose digital output
signal, which can be written by the user software --> AOU_TRIG_OUT = 0
- Trigger Input : TRG operates as a trigger input, or as general purpose digital
input signal, which can be read by the user software --> AOU_TRIG_IN = 1
"""
const AOU_TRG_OUT                   = Cint(0)
const AOU_TRG_IN                    = Cint(1)

"""
### Trigger Synchronization/Sampling Options
M31/M33XXA Digitizers & M32/M33XXA AWGs TRG connectors configuration (parameter
`syncMode` in function `triggerIOconfig` in both 'SD_AIN.jl' and 'SD_AOU.jl'
- Non-synchronized mode : The trigger is sampled with an internal 100 MHz clock
--> SYNC_NONE = 0
- Synchronized mode : (PXI form factor only) The trigger is sampled using CLK10
--> SYNC_CLK_0 = 1
"""
const SYNC_NONE                     = Cint(0)
const SYNC_CLK_0                    = Cint(1)

"""
### Clock set frequency mode
parameter `mode` in function 'clockSetFrequency' in both 'SD_AIN.jl' and
'SD_AOU.jl'
- Low Jitter Mode : The clock system is set to achieve the lowest jitter,
sacrificing tuning speed --> CLK_LOW_JITTER = 0
- Fast Tuning Mode : The clock system is set to achieve the lowest tuning time,
sacrificing jitter performance --> CLK_FAST_TUNE = 1
"""
const CLK_LOW_JITTER                = Cint(0)
const CLK_FAST_TUNE                 = Cint(1)

"""
### Window types used in FFT function
M31/M33XXA digitizer FFT window types. Parameter `windowType` in function `FFT`
in `SD_AIN.jl`
- Rectangular : Simplest B-spine window --> WINDOW_RECTANGULAR = 0 (default)
- Bartlett : Hybrid window --> WINDOW_BARTLETT = 1
- Hanning : Side-lobes roll off about 18 dB per octave --> WINDOW_HANNING = 2
- Hamming : Optimized to minimize the maximum nearest side lobe
--> WINDOW_HAMMING = 3
- Blackman : Higher-order generalized cosine windows for applications that
require windowing by the convolution in the frequency-domain
--> WINDOW_BLACKMAN = 4
- Kaiser : Adjustable window maximizing energy concentration in the main lobe
--> WINDOW_KAISER = 5
- Gauss : Adjustable window (can be used for quadratic interpolation in
frequency estimation) --> WINDOW_GAUSS = 6
"""
const WINDOW_RECTANGULAR            = Cint(0)
const WINDOW_BARTLETT               = Cint(1)
const WINDOW_HANNING                = Cint(2)
const WINDOW_HAMMING                = Cint(3)
const WINDOW_BLACKMAN               = Cint(4)
const WINDOW_KAISER                 = Cint(5)
const WINDOW_GAUSS                  = Cint(6)

"""
### Output signal selection
parameter `waveShape` in function `channelWaveShape` in `SD_AOU.jl`
- HIZ : The output signal is set to HIZ (no output signal is provided, only
available for M3202A) AOU_HIZ = -1
- No Signal : The output signal is set to 0. All other channel settings are
maintained --> AOU_OFF = 0 (default)
- Sinusoidal : Generated by the Function Generator --> AOU_SINUSOIDAL = 1
- Triangular : Generated by the Function Generator --> AOU_TRIANGULAR = 2
- Square : Generated by the Function Generator --> AOU_SQUARE = 4
- DC Voltage : The output DC voltage is set by the channel amplitude setting
--> AOU_DC = 5
- Arbitrary Waveform : Generated by the Arbitrary Waveform Generator
--> AOU_AWG = 6
- Partner Channel : Only for odd channels. It is the output of the previous
channel (to create differential signals, etc.) --> AOU_PARTNER = 8
"""
const AOU_HIZ                       = Cint(-1)
const AOU_OFF                       = Cint(0)
const AOU_SINUSOIDAL                = Cint(1)
const AOU_TRIANGULAR                = Cint(2)
const AOU_SQUARE                    = Cint(4)
const AOU_DC                        = Cint(5)
const AOU_AWG                       = Cint(6)
const AOU_PARTNER                   = Cint(8)

"""
## Modulation Options
### Angle Modulation Options
Parameter `modulationType` in function `modulationAngleConfig` in `SD_AOU.jl`
- No Modulation : Modulation is disabled --> AOU_MOD_OFF = 0 (default)
- Frequency Modulation : The AWG is used to modulate the channel frequency
--> AOU_MOD_FM = 1
- Frequency Modulation (32 bits): The AWG is used to modulate the channel
frequency. For models with -DM1 dual modulation option --> AOU_MOD_FM_32b = 1
- Phase Modulation : The AWG is used to modulate the channel phase
--> AOU_MOD_PM = 2
### Amplitude Modulation Options
- No Modulation : Modulation is disabled. The channel amplitude and offset are
only set by the main registers set by the main registers
--> AOU_MOD_OFF = 0 (default)
- Amplitude Modulation : The modulating signal is used to modulate the channel
amplitude --> AOU_MOD_AM = 1
- Offset Modulation : The modulating signal is used to modulate the channel
offset --> AOU_MOD_OFFSET = 2
"""
const AOU_MOD_OFF                   = Cint(0)
const AOU_MOD_FM                    = Cint(1)
const AOU_MOD_FM_32b                = Cint(1)
const AOU_MOD_PM                    = Cint(2)
const AOU_MOD_AM                    = Cint(1)
const AOU_MOD_OFFSET                = Cint(2)

"""
## CLK Output Options
M32/M33XXA AWGs output clock configuration. Parameter `clockConfig` in function
`clockIOconfig` in `SD_AOU.jl`
"""
const CLK_CONN_DISABLED             = Cint(0)
const CLK_REF_OUTPUT                = Cint(1)

"""
## AWG Waveform Types
Waveform types for the AWGs (parameter `waveformType` in the waveform file or in
function `new` in `SD_Wave.jl`
- Analog 16 Bits : Analog normalized waveforms (-1..1) defined with doubles
--> WAVE_ANALOG_16 = 0
- Analog 32 Bits : Analog normalized waveforms (-1..1) defined with doubles
--> WAVE_ANALOG_32 = 1
- Analog 16 Bits Dual : Analog normalized waveforms (-1..1) defined with
doubles, with two components (A and B) --> WAVE_ANALOG_16_DUAL = 4
- Analog 32 Bits Dual : Analog normalized waveforms (-1..1) defined with
doubles, with two components (A and B) --> WAVE_ANALOG_32_DUAL = 6
- IQ : Analog normalized waveforms (-1..1) defined with doubles, with two
components (I and Q) --> WAVE_IQ = 2
- IQ Polar : Analog waveforms (-1..1 module, -180..+180 phase) defined with
doubles, with two components (Module and Phase) --> WAVE_IQPOLAR = 3
- Digital : Digital waveforms defined with integers --> WAVE_DIGITAL = 5
"""
const WAVE_ANALOG_16                = Cint(0)
const WAVE_ANALOG_32                = Cint(1)
const WAVE_ANALOG_16_DUAL           = Cint(4)
const WAVE_ANALOG_32_DUAL           = Cint(6)
const WAVE_IQ                       = Cint(2)
const WAVE_IQPOLAR                  = Cint(3)
const WAVE_DIGITAL                  = Cint(5)


text:
"""
    GAE(enc[, σ])

Graph autoencoder.

# Arguments
- `enc`: encoder. It can be any graph convolutional layer.

Encoder is specified by user and decoder will be `InnerProductDecoder` layer.
"""
struct GAE{T,S}
    encoder::T
    decoder::S
end

GAE(enc, σ::Function=identity) = GAE(enc, InnerProductDecoder(σ))

@functor GAE

function (g::GAE)(X::AbstractMatrix)
    Z = g.encoder(X)
    A = g.decoder(Z)
    A
end

"""
    VGAE(enc[, σ])

Variational graph autoencoder.

# Arguments
- `enc`: encoder. It can be any graph convolutional layer.

Encoder is specified by user and decoder will be `InnerProductDecoder` layer.
"""
struct VGAE{T,S}
    encoder::T
    decoder::S
end

function VGAE(enc, h_dim::Integer, z_dim::Integer, σ::Function=identity)
    VGAE(VariationalEncoder(enc, h_dim, z_dim), InnerProductDecoder(σ))
end

@functor VGAE

function (g::VGAE)(fg::FeaturedGraph)
    fg_ = g.encoder(fg)
    fg_ = g.decoder(fg_)
    fg_
end

"""
    InnerProductDecoder(σ)

Inner-product decoder layer.

# Arguments
- `σ`: activation function.
"""
struct InnerProductDecoder
    σ
end

@functor InnerProductDecoder

(i::InnerProductDecoder)(Z::AbstractMatrix)::AbstractMatrix = i.σ.(Z'*Z)

function (i::InnerProductDecoder)(fg::FeaturedGraph)::FeaturedGraph
    Z = node_feature(fg)
    A = i(Z)
    FeaturedGraph(graph(fg), nf=A)
end

"""
    VariationalEncoder(nn, h_dim, z_dim)

Variational encoder layer.

# Arguments
- `nn`: neural network. It can be any graph convolutional layer.
- `h_dim`: dimension of hidden layer. This should fit the output dimension of `nn`.
- `z_dim`: dimension of latent variable layer. This will be parametrized into `μ` and `logσ`.

Encoder can be any graph convolutional layer.
"""
struct VariationalEncoder
    nn
    μ
    logσ
    z_dim::Integer
end

function VariationalEncoder(nn, h_dim::Integer, z_dim::Integer)
    VariationalEncoder(nn,
                       GCNConv(h_dim=>z_dim),
                       GCNConv(h_dim=>z_dim),
                       z_dim)
end

@functor VariationalEncoder

function (ve::VariationalEncoder)(fg::FeaturedGraph)::FeaturedGraph
    μ, logσ = summarize(ve, fg)
    Z = sample(μ, logσ)
    FeaturedGraph(graph(fg), nf=Z)
end

function summarize(ve::VariationalEncoder, fg::FeaturedGraph)
    fg_ = ve.nn(fg)
    fg_μ, fg_logσ = ve.μ(fg_), ve.logσ(fg_)
    node_feature(fg_μ), node_feature(fg_logσ)
end

sample(μ::AbstractArray{T}, logσ::AbstractArray{T}) where {T<:Real} =
    μ + exp.(logσ) .* randn(T, size(logσ))


text:
@testset "word vectors" begin
    wn = W2VNetwork(doc, 5, 100)
    wv = wn.wv
    a = cosine_similar_words(wn, "something", 2)
    println(a)

end


text:
using Revise
using AdFem
using LinearAlgebra
using PyPlot
using DelimitedFiles
using MAT

function f(x, y)
    10/(1+x^2) + x * y + 10*y^2
end

function calc_residual_and_jacobian(θ, u)
    E = θ[1:get_ngauss(mmesh)]
    nu = θ[get_ngauss(mmesh)+1:end]
    μ = E / (2*(1 + nu))
    λ = E*nu/((1 + nu)*(1 - 2*nu))

    ψ, J = neo_hookean(u, μ, λ, mmesh)
    B = compute_fem_source_term(zeros(get_ngauss(mmesh)), -0.5*ones(get_ngauss(mmesh)), mmesh)
    ψ = ψ - B
    J, ψ = impose_Dirichlet_boundary_conditions(J, ψ, bdnode, zeros(length(bdnode)))
    ψ, J
end

mmesh = Mesh(joinpath(PDATA, "twoholes_large.stl"))
mmesh = Mesh(mmesh.nodes * 10, mmesh.elems)

E = eval_f_on_gauss_pts(f, mmesh)
nu = 0.3 * ones(get_ngauss(mmesh))
θ = placeholder([E; nu])
bdnode = bcnode((x,y)->(x<1e-5) || (x>0.49-1e-5), mmesh)
bdnode = [bdnode; bdnode .+ mmesh.ndof]
nr = newton_raphson_with_grad(calc_residual_and_jacobian, zeros(2mmesh.ndof), θ)

ADCME.options.newton_raphson.verbose = true
sess = Session(); init(sess)
x = run(sess, nr)

matwrite("fenics/fwd2.mat", Dict("x"=>x, "E"=>E))

close("all")
figure(figsize = (15, 5))
subplot(121)
u = x[1:mmesh.nnode]
v = x[mmesh.ndof+1:mmesh.ndof + mmesh.nnode]
visualize_scalar_on_fem_points(sqrt.(u.^2+v.^2), mmesh)
subplot(122)
visualize_displacement(x, mmesh)
savefig("fenics/fwd2.png")

close("all")
visualize_scalar_on_gauss_points(E, mmesh)
savefig("fenics/E.png")


text:
@testset "Joint Life" begin

    @testset "Joint Last Survivor unknown status" begin

        @testset "ALMCR §9.4" begin
            ℓ₁ = [43302, 42854, 42081, 41351, 40050]
            ℓ₂ = [47260, 47040, 46755, 46500, 46227]

            ps₁ = ℓ₁ ./ ℓ₁[1]
            qs₁ = [1 - ps₁[t] / ps₁[t-1] for t = 2:5]

            ps₂ = ℓ₂ ./ ℓ₂[1]
            qs₂ = [1 - ps₂[t] / ps₂[t-1] for t = 2:5]

            m1 = UltimateMortality(qs₁, start_age = 65)
            m2 = UltimateMortality(qs₂, start_age = 60)

            @test decrement(m1, 65, 66) == 1 - ℓ₁[2] / ℓ₁[1]
            @test decrement(m2, 60, 61) == 1 - ℓ₂[2] / ℓ₂[1]

            l1 = SingleLife(mortality = m1, issue_age = 65)
            l2 = SingleLife(mortality = m2, issue_age = 60)

            jl = JointLife(lives = (l1, l2), contingency = LastSurvivor(), joint_assumption = Frasier())

            @test isapprox(survival(jl, 2), 0.9997, atol = 1e-4)

            ins = LifeContingency(jl, Yields.Constant(0.05))
            ins_l1 = LifeContingency(jl.lives[1], Yields.Constant(0.05))
            ins_l2 = LifeContingency(jl.lives[2], Yields.Constant(0.05))
            # problem 9.1.f
            @test isapprox(present_value(AnnuityDue(ins, 5)), 4.5437, atol = 1e-4)

            @test isapprox(present_value(AnnuityDue(ins)), 4.5437, atol = 1e-4)
        end

        @testset "CIA tables" begin
            m1 = MortalityTables.table("1986-92 CIA – Male Smoker, ANB")
            m2 = MortalityTables.table("1986-92 CIA – Female Nonsmoker, ANB")
            l1 = SingleLife(mortality = m1.ultimate, issue_age = 40)
            l2 = SingleLife(mortality = m2.ultimate, issue_age = 37)

            jl = JointLife(lives = (l1, l2), contingency = LastSurvivor(), joint_assumption = Frasier())

            @testset "independent lives" begin
                for time = 1:40
                    tpx = survival(l1, time)
                    tpy = survival(l2, time)

                    @test survival(jl, time) == tpx + tpy - tpx * tpy
                end
            end

            q_annual = [0.00000141120, 0.00000478349, 0.00000921100, 0.00001508953, 0.00002255325, 0.00003179413, 0.00004327035, 0.00005782855, 0.00007524178, 0.00009695039, 0.00012370408, 0.00015606419, 0.00019612090, 0.00024529620, 0.00030502659, 0.00037768543, 0.00046587280, 0.00057325907, 0.00070238876, 0.00085952824, 0.00104786439, 0.00127420204, 0.00154640464, 0.00187374100, 0.00226186101, 0.00272457482, 0.00327193743, 0.00391984214, 0.00468084563, 0.00557778027, 0.00662135901, 0.00783684326, 0.00924475551, 0.01086696400, 0.01272977389, 0.01485150445, 0.01727277853, 0.02000566587, 0.02308982834, 0.02654184882, 0.03039819991, 0.03469299691, 0.03944737905, 0.04469931521, 0.05048294677, 0.05683522818, 0.06379854263, 0.07142629169, 0.07976664920, 0.08888579420, 0.09884411368, 0.10971505011, 0.12159455247, 0.13456410066, 0.14873755448, 0.16421406941, 0.18111392905, 0.19954137714, 0.21962070363, 0.24144744481, 0.26515968560, 0.29099967654, 0.31901974549, 0.35111309929, 0.39733960046, 0.47004540325, 0.58255000000, 0.74852030000, 1.00000000000]
            q_cumulative = [0.00000141120, 0.00000619469, 0.00001540563, 0.00003049492, 0.00005304748, 0.00008483992, 0.00012810660, 0.00018592774, 0.00026115553, 0.00035808060, 0.00048174038, 0.00063772939, 0.00083372522, 0.00107881691, 0.00138351443, 0.00176067732, 0.00222572986, 0.00279771301, 0.00349813669, 0.00435465818, 0.00539795948, 0.00666528343, 0.00820138084, 0.01005975458, 0.01229886182, 0.01498992747, 0.01821281880, 0.02206126957, 0.02663884979, 0.03206804441, 0.03847706939, 0.04601237389, 0.05483175625, 0.06510286552, 0.07700389465, 0.09071177542, 0.10641770955, 0.12429441828, 0.14451430984, 0.16722048170, 0.19253547997, 0.22054884408, 0.25129614928, 0.28476269870, 0.32086998531, 0.35946849466, 0.40033347120, 0.44316542761, 0.48758225561, 0.53312891378, 0.57927637250, 0.62543608637, 0.67098101782, 0.71525516126, 0.75760741222, 0.79741168545, 0.83410325108, 0.86720651682, 0.89637071504, 0.92139174110, 0.94223548231, 0.95904493828, 0.97211041164, 0.98190281145, 0.98909354112, 0.99422007198, 0.99758716905, 0.99939322200, 1.00000000000]
            @testset "precalced vectors" begin
                for time = 1:40
                    @test isapprox(decrement(jl, time), q_cumulative[time], atol = 1e-6)
                end

                for time = 1:40
                    q′ = 1 - survival(jl, time) / survival(jl, time - 1)

                    @test isapprox(q′, q_annual[time], atol = 1e-6)
                end

                for time = 1:40
                    @test isapprox(survival(jl, time), 1 - q_cumulative[time], atol = 1e-6)
                end
            end
        end

    end
end


text:
include("diophantine_solution.jl")
include("diophantine_linear.jl")
include("diophantine_quadratic.jl")

function test_diophantine_all()
    print(rpad("Math.NumberTheory.Diophantine...", 50, ' '))

    test_diophantine_solutions()
    test_diophantine_solutions_nonex_noney()
    test_diophantine_solutions_anyx_anyy()
    test_diophantine_solutions_onex_anyy()
    test_diophantine_solutions_anyx_oney()
    test_diophantine_solutions_linearx_lineary()
    test_diophantine_solutions_quadraticx_quadraticy()

    test_diophantine_solve_diophantine_linear()
    test_diophantine_solve_diophantine_quadratic()

    println("PASS")
end


text:
using LinearAlgebra
"""
    lu!(A)
Compute LU factorization of the matrix `A`.

# Example
```julia-repl
julia> lu!([1 2; 3 4])
```
"""
function lu!(A::AbstractMatrix)
    n,n = size(A)
    p = collect(1:n)
    for k in 1:n-1
        m = argmax(abs.(A[k:n,k]))
        m = m + k - 1
        if A[m,k] != 0
            if m != k
                p[k], p[m] = p[m], p[k]
                for j in 1:n
                    A[k,j], A[m,j] = A[m,j], A[k,j]
                end
            end
            for i in k+1:n
                A[i,k] /= A[k,k]
                for j in k+1:n
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    L = tril(A,-1) + Matrix{Float64}(I, n, n)
    U = triu(A)
    return L, U, p
end

"""
    chol!(A)
Compute cholesky factorization of the spd matrix `A`.

# Argument
- `A:AbstractMatrix`: symmetric positive definite matrix

# Example
```julia-repl
julia> chol([1 2; 2 5])
```
"""
function chol(A::AbstractMatrix)
    n,n = size(A)
    L = similar(A)
    @inbounds begin
         for i in 1:n
             tmp = 0
             for k in 1:i-1
                tmp += L[i,k].^2
            end
            L[i,i] = sqrt(A[i,i] - tmp)
            for j in i+1:n
                L[j,i] = A[i,j] / L[i,i]
                for k in 1:i-1
                    L[j,i] -=  L[i,k].*L[j,k] / L[i,i]
                end
            end
        end
    end
    return L
end

"""
    gepp!(A,b)
Solve the linear system ``Ax=b`` using Gaussian elimination with partial pivot

# Example
```julia-repl
julia> A = rand(3,3)
julia> b = rand(3)
julia> x = gepp!(A, b)
```
"""
function gepp!(A::AbstractMatrix, b::AbstractVector)
    n,n = size(A)
    for k in 1:n-1
        # pivot
        m = argmax(abs.(A[k:n,k]))
        m = m + k - 1
        if A[m,k] != 0
            if m != k
                b[k], b[m] = b[m], b[k]
                for j in 1:n
                    A[k,j], A[m,j] = A[m,j], A[k,j]
                end
            end
            for i in k+1:n
                A[i,k] /= A[k,k]
                b[i] -= A[i,k]*b[k]
                for j in k+1:n
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    # back substituion
    x = zeros(n)
    for k in n:-1:1
        for j = n:-1:k+1
            b[k] -= A[k,j]*x[j]
        end
        x[k] = b[k] / A[k,k]
    end
    return x
end

"""
    jacobi(A, b, x0)
    jacobi(A, b, x0, tol)

Solve the linear system using Jacobi method.

# Example
```julia-repl
julia> A = rand(3,3)
julia> b = rand(3)
julia> x0 = rand(3)
julia> x = jacobi(A, b, x0)
```
"""
function jacobi(A::AbstractMatrix, b::AbstractVector, x::AbstractVector, tol::Float64=1e-10)
    xc = similar(b)
    n, = size(b)
    err = 1.0
    it = 0
    while err > tol
        @inbounds for i in 1:n
            xc[i] = b[i]
            @inbounds for k in 1:n
                if k != i
                    xc[i] -= A[i,k]*x[k]
                end
            end
            xc[i] /= A[i,i]
        end
        err = norm(x-xc)
        x = copy(xc)
        it += 1
    end
    return xc, it
end

"""
    gaussseidel(A, b, x0)
    gaussseidel(A, b, x0, tol)

Solve the linear system using Gauss-Seidel method.

# Example
```julia-repl
julia> A = rand(3,3)
julia> b = rand(3)
julia> x0 = rand(3)
julia> x = gaussseidel(A, b, x0)
```
"""
function gaussseidel(A::AbstractMatrix, b::AbstractVector, x::AbstractVector, tol::Float64=1e-10)
    n, = size(b)
    err = 1.0
    it = 0
    xc = copy(x)
    while err > tol
        @inbounds for i in 1:n
            xc[i] = b[i]
            @inbounds for k in 1:n
                if k != i
                    xc[i] -= A[i,k]*xc[k]
                end
            end
            xc[i] /= A[i,i]
        end
        err = norm(x-xc)
        x = copy(xc)
        it += 1
    end
    return xc, it
end


text:
module Syntax

using ..OmegaCore.Util: mapf
using ..OmegaCore.Var: pw
using Distributions
export @joint

"""
Random variable over named tuple using varnames as keys

`@joint randvar1 randvar2 ...`

```julia
using Distributions
a = 1 ~ Normal(0, 1)
b = 2 ~ Normal(0, 1)
c = a >=ₚ b
randsample(@joint a b c)
=> (Ashoots = false, Bshoots = false)
```
"""
macro joint(args::Symbol...)
  esc(:(ω -> NamedTuple{$args}(OmegaCore.Util.mapf(ω, tuple($(args...))))))
end

export ==ₚ, >=ₚ, <=ₚ, >ₚ, <ₚ, !ₚ, &ₚ, |ₚ, ifelseₚ, +ₚ, -ₚ, *ₚ, /ₚ
@inline x ==ₚ y = pw(==, x, y)
@inline x >=ₚ y = pw(>=, x, y)
@inline x >ₚ y = pw(>, x, y)
@inline x <ₚ y = pw(<, x, y)
@inline x <=ₚ y = pw(<=, x, y)
@inline x +ₚ y = pw(+, x, y)
@inline x -ₚ y = pw(-, x, y)
@inline x *ₚ y = pw(*, x, y)
@inline x /ₚ y = pw(/, x, y)

@inline x |ₚ y = pw(|, x, y)
@inline x &ₚ y = pw(&, x, y)
@inline !ₚ(x) = pw(!, x)
@inline ifelseₚ(a, b, c) = pw(ifelse, a, b, c)

#FIXme generalize this
# Normalₚ(args...) = pw(Distributions.Normal, args...)
# Uniformₚ(args...) = pw(Distributions.Uniform, args...)
# Gammaₚ(args...) = pw(Distributions.Gamma, args...)
# DiscreteUniformₚ(args...) = pw(Distributions.DiscreteUniform, args...)
# Poissonₚ(args...) = pw(Distributions.Poisson, args...)
# NegativeBinomialₚ(args...) = pw(Distributions.NegativeBinomial, args...)

# export Normalₚ,
#        Uniformₚ,
#        Gammaₚ,
#        DiscreteUniformₚ,
#        Poissonₚ,
#        NegativeBinomialₚ

end


text:
using FileIO, FixedPointNumbers, Colors, ImageCore, Test, Statistics

workdir = joinpath(tempdir(), "ImagesNRRD")
isdir(workdir) || mkdir(workdir)
const teembase = "http://teem.sourceforge.net/nrrd/files"
const base3D   = "http://www.tc18.org/code_data_set/3D_greyscale"

function getfile(name, base=teembase)
    file = joinpath(workdir, name)
    if !isfile(file)
        file = download(joinpath(base, name), file)
    end
    file
end

@testset "readremote" begin
    fool_file = getfile("fool.nrrd")
    imgu = load(fool_file)
    @test eltype(imgu) == UInt8
    @test size(imgu) == (128,128)
    @test imgu[120,5] < imgu[5,5] < imgu[120,120]  # a simple orientation test
    foolf_file = getfile("foolf.nrrd")
    imgf = load(foolf_file)
    @test eltype(imgf) == Float32
    @test size(imgf) == (128,128)
    @test imgf[120,5] < imgf[5,5] < imgf[120,120]
    @test maximum(abs.(imgf-imgu))/maximum(imgu) <= 1/200
    foolc_file = getfile("foolc.nrrd")
    imgc = load(foolc_file, RGB)
    @test eltype(imgc) == RGB{N0f8}
    @test size(imgc) == (128,128)
    imgg = convert(Array{Gray}, imgc)
    @test mean(abs.(Int.(rawview(channelview(imgg))) - Int.(imgu))) < 12
#     aneurism_rawfile = getfile("aneurism.raw.gz", base3D)
#     aneurism_file = joinpath(workdir, "aneurism.nhdr")
#     aneurism_rawfile_path = joinpath(workdir, "aneurism.raw.gz")
#     open(aneurism_file, "w") do io
#         write(io, """
# NRRD0001
# content: aneurism
# # Courtesy of Philips Research, Hamburg, Germany
# dimension: 3
# type: uchar
# sizes: 256 256 256
# spacings: 1 1 1
# data file: $aneurism_rawfile_path
# encoding: gzip
# """)
#     end
#     imga = load(aneurism_file)
#     @test eltype(imga) == UInt8
#     @test size(imga) == (256,256,256)
#     @test sdims(imga) == 3
#     @test pixelspacing(imga) == (1.0,1.0,1.0)
end

GC.gc()  # to close any mmapped files
rm(workdir, recursive=true)

nothing


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

process predefined datasets into txt files for fast importing via readdlm
=#

using LinearAlgebra
import DataFrames
import CSV
import DelimitedFiles

# naics5811 dataset (save as naics5811.csv)
# see https://arxiv.org/pdf/1509.08165v1.pdf (example 3)
# data obtained from http://www.nber.org/data/nbprod2005.html
function get_naics5811_data()
    df = CSV.read(joinpath(@__DIR__, "naics5811.csv"), copycols = true)
    DataFrames.deleterows!(df, 157) # outlier
    # number of non production employees
    df[!, :prode] .= df[!, :emp] - df[!, :prode]
    # group by industry codes
    df_aggr = DataFrames.aggregate(DataFrames.dropmissing(df), :naics, sum)

    # four covariates: non production employees, production worker hours, production workers, total capital stock
    # use the log transform of covariates
    X = log.(convert(Matrix{Float64}, df_aggr[!, [:prode_sum, :prodh_sum, :prodw_sum, :cap_sum]])) # n = 4
    # value of shipment
    y = convert(Vector{Float64}, df_aggr[!, :vship_sum])
    # mean center
    X .-= sum(X, dims = 1) ./ size(X, 1)
    y .-= sum(y) / length(y)
    # normalize to unit norm
    X ./= norm.(eachcol(X))'
    y /= norm(y)

    return hcat(X, y) # concatenate columns of X and y
end

open(joinpath(@__DIR__, "naics5811.txt"), "w") do io
    DelimitedFiles.writedlm(io, get_naics5811_data())
end


text:
#Modifying impute function from https://github.com/invenia/Impute.jl/blob/master/src/imputors/interp.jl

"""
Uses linear interpolation to fill for missing values.
# Example:
```
x = rand(1:10.0, 10)
x[5] = missing
linear_interpolation(x)
```
"""
function linear_interpolation(timeseries)
    if  count(i->(ismissing(i)), timeseries) > length(timeseries) *1/2
        return zero(1:length(timeseries))
    end
    data = copy(timeseries)
    i = findfirst(!ismissing, data) + 1
    while i < length(data)
        if ismissing(data[i])
            prev_idx = i - 1
            next_idx = findnext(!ismissing, data, i + 1)
            if next_idx !== nothing
                gap_sz = (next_idx - prev_idx) - 1

                diff = data[next_idx] - data[prev_idx]
                incr = diff / (gap_sz + 1)
                val = data[prev_idx] + incr
                for j in i:(next_idx - 1)
                    data[j] = val
                    val += incr
                end
                i = next_idx
            else
                break
            end
        end
        i += 1
    end
    if count(i->(ismissing(i)), data)>0
        data[1:findfirst(!ismissing, data)] .= data[findfirst(!ismissing, data)]
        data = reverse(data)
        data[1:findfirst(!ismissing, data)] .= data[findfirst(!ismissing, data)]
        data = reverse(data)
    end
    return data
end