text:
#!/usr/bin/env julia

using ArgParse.ArgParseSettings
using ArgParse.@add_arg_table
using ArgParse.parse_args
using DataFrames
using CSV
using StatsBase

NEW_PATH=join(split(Base.source_path(),"/")[1:end-3],"/")
push!(LOAD_PATH, NEW_PATH)
using FidelityUMI

function main(args)
    arg_parse_settings = ArgParseSettings(description="Calculates N of variants based on UMIs groups and filtering sequencing errors.
                                                      UMI_tools should be used for tagging readname and/or grouping.")
    @add_arg_table arg_parse_settings begin

        "--bam", "-b"
            arg_type = String
            help = "name of bam file."
            required = true
            dest_name = "bam"

        "--output", "-o"
            arg_type = String
            help = "Output csv file"
            default = "output.csv"
            dest_name = "out"

        "--reference", "-r"
            arg_type = String
            help = "Reference fasta."
            required = true
            dest_name = "ref"

        "--alignment-quality", "-q"
            arg_type = Int64
            help = "Minimum alignment quality of a read."
            default = 0
            dest_name = "alignQ"

        "--cluster-size", "-c"
            arg_type = Int64
            help = "How many read pairs in cluster should be considered."
            default = 0
            dest_name = "clust_size"

        "--mut-per-read", "-m"
            arg_type = Int64
            help = "Allowed max number of mutations in read."
            default = 3
            dest_name = "mut_n_read"

        "--consider-as-ref-error", "-e"
            arg_type = Float64
            help = "Consider mutations at position as reference error if fraction is >= x."
            default = 1.0
            dest_name = "ref_error"

        "--base-quality", "-Q"
            arg_type = Int64
            help = "Base quality [>=]."
            default = 30
            dest_name = "q_limit"

        "--umi-string", "-s"
            arg_type = String
            help = "Umi place holder string."
            default = "_"
            dest_name = "umi_s"

        "--umi-place", "-p"
            arg_type = Int64
            help = "Umi place after spliting by 'umi-string'. 1 based array."
            default = 2
            dest_name = "umi_p"

        "--verbose", "-v"
            help = "Print log messages to STDOUT."
            action = :store_true
            dest_name = "verb"

        "--use-grouped-reads", "-g"
            help = ("Force to use already grouped reads by UMI (using `umi_tools group`).
                    If so, provide which tag to use with --group-tag option.")
            action = :store_true
            dest_name = "use_grouped"

        "--group-tag", "-t"
            arg_type = String
            help = "Which tag should be used if --use-grouped-reads option is selected."
            default = "BX"
            dest_name = "group_tag"

        "--report-variants", "-R"
            help = "If to report sets of mutations."
            action = :store_true
            dest_name = "rep_vars"

        "--gene-conde", "-C"
            arg_type = Int64
            help = "Gene code from BioSequences.ncbi_trans_table"
            default = 11
            dest_name = "gene_code"

        "--translate", "-T"
            help = ("Set to translte genes. --genes-csv is needed.")
            action = :store_true
            dest_name = "translate"

        "--gene-csv", "-G"
            arg_type = String
            help = "Gene csv file if translate. Header: Chromosome,Start,End,Gene,Strand."
            dest_name = "gene_csv"

        "--algorithm", "-a"
            arg_type = String
            help = (
                "Algorithm for calculating fidelity. Available options are:
                umi - group reads to families/clusters by unique molecular identifiers.
                position - group reads to families/clusters by mapping positions.
                paired - do not group but determine confident error only if found in both paired end reads.
                    This means that insert of a library should be small to get as many overlaping/read-through as possible.
                naive - do not group and do not filter by anything.
                ")
            default = "umi"
            dest_name = "algo"

        "--frequency-in-cluster", "-F"
            arg_type = Float64
            help = "Mutation frequency in cluster. If larger, mutation passes."
            default = 1.0
            dest_name = "freq_in_clust"

        "--min-cluster-counts", "-L"
            arg_type = Int64
            help = "Minimum count of clusters for the mutation to be reported to a single_mutaitons file."
            default = 1
            dest_name = "min_cluster_counts"
        "--min-mutation-coverage", "-M"
            arg_type = Int64
            help = "Minimum confident coverage for the mutation to be reported to a single_mutaitons file."
            default = 1
            dest_name = "min_mutation_coverage"
    end

    parsed_args = parse_args(arg_parse_settings)

    file = parsed_args["bam"]
    alignQ = parsed_args["alignQ"]
    mut_n_read = parsed_args["mut_n_read"]
    q_limit = parsed_args["q_limit"]
    umi_s = parsed_args["umi_s"]
    umi_p = parsed_args["umi_p"]
    use_grouped = parsed_args["use_grouped"]
    group_tag = parsed_args["group_tag"]
    verbose = parsed_args["verb"]
    reference = parsed_args["ref"]
    clust_size = parsed_args["clust_size"]
    ref_error = parsed_args["ref_error"]
    output = parsed_args["out"]
    algo = parsed_args["algo"]
    freq_in_clust = parsed_args["freq_in_clust"]
    min_cluster_counts = parsed_args["min_cluster_counts"]
    min_mutation_coverage = parsed_args["min_mutation_coverage"]

    if !(algo in ["umi", "position", "paired", "naive"])
        println(STDERR, "ERROR: Selected --algorithm does not exist: $algo")
        println(STDERR, "ERROR: Options are: umi position paired naive")
        exit(1)
    end

    referenceSruct = ReferenceSequences(reference)
    cvCoverage = Coverage(referenceSruct)
    coverage_counter!(cvCoverage, file; mappingQuality=alignQ)

    clust_coverages = Coverage(referenceSruct)

    if algo == "paired" || algo == "naive"
        clust_size = 1
    end

    (vars,
    mut_coverages,
    raw_bases,
    total_reads,
    filtered_reads,
    filtered_R1,
    filtered_R2,
    error_distr_per_read
    ) = collect_variants(file, clust_coverages, alignQ, mut_n_read, q_limit,
                         umi_s, umi_p, use_grouped, group_tag; verbose=verbose, algo=algo)

    (passing_vars,
    err_distr_per_cluster,
    clust_size_distr,
    found_mutations,
    found_insertions,
    found_deletions,
    found_ns,
    total_clusters,
    conf_bases,
    clusts_with_vars,
    passing_mutations,
    total_bases,
    conf_coverage
    ) = filter_variants(vars, mut_coverages,
                           cvCoverage, referenceSruct,
                           clust_size, ref_error;
                           verbose=verbose, algo=algo,
                           freq_in_clust=freq_in_clust)

    df = DataFrame(
        Sample=String[],
        TotalErrors=Int64[],
        RawBases=Int64[],
        TotalClusterBases=Int64[],
        ConfidentBases=Int64[],
        Frequency=Float64[],
        MeanClustSize=Float64[],
        MinClusterSize=Int64[],
        MaxClusterSize=Int64[],
        SumClusterSize=Int64[],
        Mutations=Int64[],
        Insertions=Int64[],
        Deletions=Int64[],
        Ns=Int64[],
        TotalClusters=Int64[],
        ClustersWithVars=Int64[],
        TotalReads=Int64[],
        FilteredReads=Int64[],
        FilteredR1=Int64[],
        FilteredR2=Int64[]
    )

    sample_name = split(join(split(file,".")[1:end-1], "_"),"/")[end]
    push!(df, [
            sample_name,
            passing_vars,
            raw_bases,
            total_bases,
            conf_bases,
            passing_vars/conf_bases,
            mean(clust_size_distr),
            minimum(clust_size_distr),
            maximum(clust_size_distr),
            sum(clust_size_distr),
            found_mutations,
            found_insertions,
            found_deletions,
            found_ns,
            total_clusters,
            clusts_with_vars,
            total_reads,
            filtered_reads,
            filtered_R1,
            filtered_R2
        ]
    )

    CSV.write(output, df, delim=',')
    outfile = replace(output, ".csv", "")

    df = DataFrame(
        Sample=String[],
        Chromosome=String[],
        ClustOfVariants=String[],
        Count=Int64[],
        Coverage=Float64[],
        ConfidentCoverage=Float64[]
    )

    write_array(outfile*"_distributions.csv", err_distr_per_cluster)
    write_array(outfile*"_size_distribution.csv", clust_size_distr)
    write_array(outfile*"_error_distributions_per_read.csv", error_distr_per_read; new_line=false)

    transitions, transversions, types = transvers(passing_mutations)

    write_array(outfile*"_transitversions.csv", ["Sample","Transitions","Transversions"]; new_line=false)
    write_array(outfile*"_transitversions.csv", [sample_name, transitions, transversions]; new_line=false, meth="a")
    header = collect(keys(types))
    values = [string(types[i]) for i in header]
    write_array(outfile*"_transvers_types.csv", unshift!(header, "Sample"); new_line=false)
    write_array(outfile*"_transvers_types.csv", unshift!(values, sample_name); new_line=false, meth="a")

    if parsed_args["rep_vars"]
        mapped = mapp_mutations(passing_mutations)

        for mutations in keys(mapped)
            # average coverage
            cov = 0
            conf_cov = 0
            ct = 0
            clust_of_variants = ""
            chr = ""

            for mutation in mutations
                chr = mutation.chrom
                try
                    if mutation.mStart != mutation.mEnd
                        println(mutation)
                    end
                    conf_cov += mean(conf_coverage.coverage[chr][mutation.mStart:mutation.mEnd])
                    cov += mean(clust_coverages.coverage[chr][mutation.mStart:mutation.mEnd])
                    ct += 1
                catch error
                    println(STDERR, "ERROR: $error")
                    continue
                end

                clust_of_variants = clust_of_variants*";"*stringify(mutation)
            end

            push!(df, [
                    sample_name,
                    chr,
                    clust_of_variants[2:end],
                    mapped[mutations],
                    round(cov/ct, 0),
                    round(conf_cov/ct, 0),
                ]
            )
        end

        CSV.write(outfile*"_passing_mutations.csv", df, delim=',')

        df = DataFrame(
            Sample=String[],
            Chromosome=String[],
            ClustOfVariants=String[],
            Count=Int64[],
            Coverage=Float64[],
            ConfidentCoverage=Float64[]
        )
        mapped = mapp_mutations(passing_mutations; cluster=false)

        for mutation in keys(mapped)
            # average coverage
            cov = 0
            conf_cov = 0
            chr = mutation.chrom

            try
                conf_cov += sum(conf_coverage.coverage[chr][mutation.mStart:mutation.mEnd])
                cov += sum(clust_coverages.coverage[chr][mutation.mStart:mutation.mEnd])
            catch error
                println(STDERR, "ERROR: $error")
                continue
            end

            if mapped[mutation] >= min_cluster_counts && conf_cov >= min_mutation_coverage
                push!(df, [
                        sample_name,
                        chr,
                        stringify(mutation),
                        mapped[mutation],
                        cov,
                        conf_cov,
                    ]
                )
            end
        end

        CSV.write(outfile*"_passing_single_mutations.csv", df, delim=',')

    end

    # Tranlate mutations

    if parsed_args["translate"]
        if verbose
            println(STDOUT, "INFO: translating passed mutations. Using ", parsed_args["gene_csv"])
        end

        if parsed_args["gene_csv"] == ""
            println(STDERR, "ERROR: gene bed file was not provided.")
            println(STDERR, "ERROR: Provide gene bed file if you want to translate mutations.")
            exit(1)
        end
        aa_substitutions, frame_shifts = translate_mutations(
            referenceSruct, parsed_args["gene_csv"],
            passing_mutations;
            genetic_code=parsed_args["gene_code"],
            verbose=verbose
            )

        df = DataFrame(
            Sample=String[],
            Chromosome=String[],
            ClustOfVariants=String[],
            Count=Int64[],
            Coverage=Float64[],
            ConfidentCoverage=Float64[]
        )

        for chr in keys(aa_substitutions)
            mapped = countmap(aa_substitutions[chr])
            for mutations in keys(mapped)
                # average coverage
                conf_cov = 0
                cov = 0
                ct = 0
                clust_of_variants = ""

                for mutation in mutations
                    position = mutation[2:end-1]
                    position = parse(Int64, String(position))
                    position = position * 3
                    try
                        conf_cov += sum(conf_coverage.coverage[chr][position-2:position])
                        cov += sum(clust_coverages.coverage[chr][position-2:position])
                        ct += 3
                    catch error
                        println(STDERR, "ERROR: $error")
                        continue
                    end

                    clust_of_variants = clust_of_variants*";"*mutation
                end

                if mapped[mutations] >= min_cluster_counts && round(conf_cov/ct, 0) >= min_mutation_coverage
                    push!(df, [
                            sample_name,
                            chr,
                            clust_of_variants[2:end],
                            mapped[mutations],
                            round(cov/ct, 0),
                            round(conf_cov/ct, 0)
                        ]
                    )
                end
            end
        end

        CSV.write(outfile*"_passing_aa_mutations.csv", df, delim=',')
        df = DataFrame(Sample=String[], FrameShifts=Int64[])
        push!(df, [sample_name frame_shifts])
        CSV.write(outfile*"_passing_frameshifts.csv", df, delim=',')

        df = DataFrame(
            Sample=String[],
            Chromosome=String[],
            Variants=String[],
            Count=Int64[],
            Coverage=Float64[],
            ConfidentCoverage=Float64[]
        )

        for chr in keys(aa_substitutions)
            mapped = mapp_single_aa_subs(aa_substitutions[chr])
            for mutation in keys(mapped)
                # average coverage
                conf_cov = 0
                cov = 0
                ct = 0
                clust_of_variants = ""

                position = mutation[2:end-1]
                position = parse(Int64, String(position))
                position = position * 3
                try
                    conf_cov += sum(conf_coverage.coverage[chr][position-2:position])
                    cov += sum(clust_coverages.coverage[chr][position-2:position])
                    ct += 3
                catch error
                    println(STDERR, "ERROR: $error")
                    continue
                end

                if mapped[mutation] >= min_cluster_counts && round(conf_cov/ct, 0) >= min_mutation_coverage
                    push!(df, [
                            sample_name,
                            chr,
                            mutation,
                            mapped[mutation],
                            round(cov/ct, 0),
                            round(conf_cov/ct, 0)
                        ]
                    )
                end
            end
        end
        CSV.write(outfile*"_passing_single_aa_mutations.csv", df, delim=',')
    end
end

function write_array(flname::String, arr; new_line::Bool=true, meth::String="w")

    out = open(flname, meth)

    if new_line
        for i in arr
            write(out, string(i)*"\n")
        end
    else
        ct = 0
        for i in arr
            ct += 1
            if ct == 1
                write(out, string(i))
            else
                write(out, ","*string(i))
            end
        end
        write(out, "\n")
    end

    close(out)
end

main(ARGS)

max_stars_repo_name: KarolisM/rnd-it-molbio-fidelity-umi

max_issues_repo_name: KarolisM/rnd-it-molbio-fidelity-umi

max_forks_repo_name: KarolisM/rnd-it-molbio-fidelity-umi


text:
export iauC2i00b
"""
**  Form the celestial-to-intermediate matrix for a given date using the
**  IAU 2000B precession-nutation model.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards Of Fundamental Astronomy) software collection.
**
**  Status:  support function.
**
**  Given:
**     date1,date2 double       TT as a 2-part Julian Date (Note 1)
**
**  Returned:
**     rc2i        double[3][3] celestial-to-intermediate matrix (Note 2)
**
**  Notes:
**
**  1) The TT date date1+date2 is a Julian Date, apportioned in any
**     convenient way between the two arguments.  For example,
**     JD(TT)=2450123.7 could be expressed in any of these ways,
**     among others:
**
**            date1          date2
**
**         2450123.7           0.0       (JD method)
**         2451545.0       -1421.3       (J2000 method)
**         2400000.5       50123.2       (MJD method)
**         2450123.5           0.2       (date & time method)
**
**     The JD method is the most natural and convenient to use in
**     cases where the loss of several decimal digits of resolution
**     is acceptable.  The J2000 method is best matched to the way
**     the argument is handled internally and will deliver the
**     optimum resolution.  The MJD method and the date & time methods
**     are both good compromises between resolution and convenience.
**
**  2) The matrix rc2i is the first stage in the transformation from
**     celestial to terrestrial coordinates:
**
**        [TRS]  =  RPOM * R_3(ERA) * rc2i * [CRS]
**
**               =  rc2t * [CRS]
**
**     where [CRS] is a vector in the Geocentric Celestial Reference
**     System and [TRS] is a vector in the International Terrestrial
**     Reference System (see IERS Conventions 2003), ERA is the Earth
**     Rotation Angle and RPOM is the polar motion matrix.
**
**  3) The present function is faster, but slightly less accurate (about
**     1 mas), than the iauC2i00a function.
**
**  Called:
**     iauPnm00b    classical NPB matrix, IAU 2000B
**     iauC2ibpn    celestial-to-intermediate matrix, given NPB matrix
**
**  References:
**
**     "Expressions for the Celestial Intermediate Pole and Celestial
**     Ephemeris Origin consistent with the IAU 2000A precession-
**     nutation model", Astron.Astrophys. 400, 1145-1154
**     (2003)
**
**     n.b. The celestial ephemeris origin (CEO) was renamed "celestial
**          intermediate origin" (CIO) by IAU 2006 Resolution 2.
**
**     McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
**     IERS Technical Note No. 32, BKG (2004)
**
**  This revision:  2013 June 18
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauC2i00b(double date1, double date2, double rc2i[3][3])

function iauC2i00b(date1::Real, date2::Real)
   # Allocate return value
   rc2i = zeros(Float64, 3, 3)

   ccall((:iauC2i00b, libsofa_c), Cvoid,
         (Cdouble, Cdouble,
         Ptr{Cdouble}),
         convert(Float64, date1), convert(Float64, date2),
         rc2i)

   return SMatrix{3,3}(rc2i')
end

ext: jl


text:
# Automa.jl generated readrecord! and readmetainfo! functions
# ========================================

import Automa
import Automa.RegExp: @re_str
import Automa.Stream: @mark, @markpos, @relpos, @abspos

# file   = header . body
# header = metainfo*
# body   = record*
const sam_machine_metainfo, sam_machine_record, sam_machine_header, sam_machine_body, sam_machine = (function ()

    isinteractive() && @info "compiling SAM"

    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    alt = Automa.RegExp.alt
    opt = Automa.RegExp.opt
    any = Automa.RegExp.any

    metainfo = let
        tag = re"[A-Z][A-Z]" \ cat("CO")
        tag.actions[:enter] = [:pos1]
        tag.actions[:exit]  = [:metainfo_tag]

        dict = let
            key = re"[A-Za-z][A-Za-z0-9]"
            key.actions[:enter] = [:pos2]
            key.actions[:exit]  = [:metainfo_dict_key]
            val = re"[ -~]+"
            val.actions[:enter] = [:pos2]
            val.actions[:exit]  = [:metainfo_dict_val]
            keyval = cat(key, ':', val)

            cat(keyval, rep(cat('\t', keyval)))
        end
        dict.actions[:enter] = [:pos1]
        dict.actions[:exit]  = [:metainfo_val]

        co = cat("CO")
        co.actions[:enter] = [:pos1]
        co.actions[:exit]  = [:metainfo_tag]

        comment = re"[^\r\n]*" # Note: Only single line comments are allowed.
        comment.actions[:enter] = [:pos1]
        comment.actions[:exit]  = [:metainfo_val]

        cat('@', alt(cat(tag, '\t', dict), cat(co, '\t', comment)))
    end
    metainfo.actions[:enter] = [:mark]
    metainfo.actions[:exit]  = [:metainfo]

    record = let
        qname = re"[!-?A-~]+"
        qname.actions[:enter] = [:pos]
        qname.actions[:exit]  = [:record_qname]

        flag = re"[0-9]+"
        flag.actions[:enter] = [:pos]
        flag.actions[:exit]  = [:record_flag]

        rname = re"\*|[!-()+-<>-~][!-~]*"
        rname.actions[:enter] = [:pos]
        rname.actions[:exit]  = [:record_rname]

        pos = re"[0-9]+"
        pos.actions[:enter] = [:pos]
        pos.actions[:exit]  = [:record_pos]

        mapq = re"[0-9]+"
        mapq.actions[:enter] = [:pos]
        mapq.actions[:exit]  = [:record_mapq]

        cigar = re"\*|([0-9]+[MIDNSHPX=])+"
        cigar.actions[:enter] = [:pos]
        cigar.actions[:exit]  = [:record_cigar]

        rnext = re"\*|=|[!-()+-<>-~][!-~]*"
        rnext.actions[:enter] = [:pos]
        rnext.actions[:exit]  = [:record_rnext]

        pnext = re"[0-9]+"
        pnext.actions[:enter] = [:pos]
        pnext.actions[:exit]  = [:record_pnext]

        tlen = re"[-+]?[0-9]+"
        tlen.actions[:enter] = [:pos]
        tlen.actions[:exit]  = [:record_tlen]

        seq = re"\*|[A-Za-z=.]+"
        seq.actions[:enter] = [:pos]
        seq.actions[:exit]  = [:record_seq]

        qual = re"[!-~]+"
        qual.actions[:enter] = [:pos]
        qual.actions[:exit]  = [:record_qual]

        field = let
            tag = re"[A-Za-z][A-Za-z0-9]"
            val = alt(
                re"A:[!-~]",
                re"i:[-+]?[0-9]+",
                re"f:[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?",
                re"Z:[ !-~]*",
                re"H:([0-9A-F][0-9A-F])*",
                re"B:[cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+")

            cat(tag, ':', val)
        end
        field.actions[:enter] = [:pos]
        field.actions[:exit]  = [:record_field]

        cat(
            qname, '\t',
            flag,  '\t',
            rname, '\t',
            pos,   '\t',
            mapq,  '\t',
            cigar, '\t',
            rnext, '\t',
            pnext, '\t',
            tlen,  '\t',
            seq,   '\t',
            qual,
            rep(cat('\t', field)))
    end
    record.actions[:enter] = [:mark]
    record.actions[:exit]  = [:record]

    newline = let
        lf = re"\n"
        lf.actions[:enter] = [:countline]

        cat(re"\r?", lf)
    end

    header = rep(cat(metainfo, newline))
    header.actions[:exit] = [:header]

    body = record * rep(newline * record) * opt(newline)
    body.actions[:exit]  = [:body]

    sam = cat(header, body)

    return map(Automa.compile, (metainfo, record, header, body, sam))
end)()

# write("sam_machine_metainfo.dot", Automa.machine2dot(sam_machine_metainfo))
# run(`dot -Tsvg -o sam_machine_metainfo.svg sam_machine_metainfo.dot`)
#
# write("sam_machine_record.dot", Automa.machine2dot(sam_machine_record))
# run(`dot -Tsvg -o sam_machine_record.svg sam_machine_record.dot`)
#
# write("sam_machine_header.dot", Automa.machine2dot(sam_machine_header))
# run(`dot -Tsvg -o sam_machine_header.svg sam_machine_header.dot`)
#
# write("sam_machine_body.dot", Automa.machine2dot(sam_machine_body))
# run(`dot -Tsvg -o sam_machine_body.svg sam_machine_body.dot`)
#
# write("sam_machine.dot", Automa.machine2dot(sam_machine))
# run(`dot -Tsvg -o sam_machine.svg sam_machine.dot`)

function appendfrom!(dst, dpos, src, spos, n)
    if length(dst) < dpos + n - 1
        resize!(dst, dpos + n - 1)
    end
    unsafe_copyto!(dst, dpos, src, spos, n)
    return dst
end

const action_metainfo = quote
    appendfrom!(metainfo.data, 1, data, @markpos, p-@markpos)
    metainfo.filled = 1:(p-@markpos)

    found_metainfo = true
end

const sam_actions_metainfo = Dict(
    :mark => :(@mark),
    :pos1  => :(pos1 = @relpos(p)),
    :pos2  => :(pos2 = @relpos(p)),
    :metainfo_tag => :(metainfo.tag = pos1:@relpos(p-1)),
    :metainfo_val => :(metainfo.val = pos1:@relpos(p-1)),
    :metainfo_dict_key => :(push!(metainfo.dictkey, pos2:@relpos(p-1))),
    :metainfo_dict_val => :(push!(metainfo.dictval, pos2:@relpos(p-1))),
    :metainfo => action_metainfo
)

const sam_actions_header = merge(
    sam_actions_metainfo,
    Dict(
        :countline => :(linenum += 1),
        :metainfo => quote
            $(action_metainfo)
            push!(header, metainfo)
            metainfo = MetaInfo()
        end,
        :header => quote

            finish_header = true

            @escape
        end
    )
)

const sam_actions_record = Dict(
    :mark => :(@mark),
    :pos => :(pos = @relpos(p)),
    :record_qname => :(record.qname = pos:@relpos(p-1)),
    :record_flag  => :(record.flag  = pos:@relpos(p-1)),
    :record_rname => :(record.rname = pos:@relpos(p-1)),
    :record_pos   => :(record.pos   = pos:@relpos(p-1)),
    :record_mapq  => :(record.mapq  = pos:@relpos(p-1)),
    :record_cigar => :(record.cigar = pos:@relpos(p-1)),
    :record_rnext => :(record.rnext = pos:@relpos(p-1)),
    :record_pnext => :(record.pnext = pos:@relpos(p-1)),
    :record_tlen  => :(record.tlen  = pos:@relpos(p-1)),
    :record_seq   => :(record.seq   = pos:@relpos(p-1)),
    :record_qual  => :(record.qual  = pos:@relpos(p-1)),
    :record_field => :(push!(record.fields, pos:@relpos(p-1))),
    :record       => quote
        appendfrom!(record.data, 1, data, @markpos, p-@markpos)
        record.filled = 1:(p-@markpos)

        found_record = true
        @escape
    end
)

const sam_actions_body = merge(
    sam_actions_record,
    Dict(
        :countline => :(linenum += 1),
        :body => quote
            finish_body = true
            @escape
        end
    )
)

# const sam_actions = merge(
#     sam_actions_header,
#     sam_actions_body
# )

const sam_context = Automa.CodeGenContext(
    generator = :goto,
    checkbounds = false,
    loopunroll = 0
)

const sam_initcode_metainfo = quote
    pos1 = 0
    pos2 = 0
    found_metainfo = false
end

const sam_initcode_record = quote
    pos = 0
    found_record = false
end

const sam_initcode_header = quote
    $(sam_initcode_metainfo)
    metainfo = MetaInfo()
    finish_header = false
    cs, linenum = state
end

const sam_initcode_body = quote
    $(sam_initcode_record)
    finish_body = false
    cs, linenum = state
end

const sam_loopcode_metainfo = quote

    if cs < 0
        throw(ArgumentError("malformed metainfo at pos $(p)"))
    end

    if found_metainfo
        @goto __return__
    end
end

const sam_returncode_metainfo = quote
    return found_metainfo
end

Automa.Stream.generate_reader(
    :index!,
    sam_machine_metainfo,
    arguments = (:(metainfo::MetaInfo),),
    actions = sam_actions_metainfo,
    context = sam_context,
    initcode = sam_initcode_metainfo,
    loopcode = sam_loopcode_metainfo,
    returncode = sam_returncode_metainfo
) |> eval

const sam_loopcode_header = quote

    if finish_header
        @goto __return__
    end
end

const sam_returncode_header = quote
    return cs, linenum, finish_header
end

Automa.Stream.generate_reader(
    :readheader!,
    sam_machine_header,
    arguments = (:(header::SAM.Header), :(state::Tuple{Int,Int})),
    actions = sam_actions_header,
    context = sam_context,
    initcode = sam_initcode_header,
    loopcode = sam_loopcode_header,
    returncode = sam_returncode_header
) |> eval

const sam_loopcode_record = quote

    if cs < 0
        throw(ArgumentError("malformed SAM record at position $(p), line $(linenum)"))
    end

    # # if cs != 0
    # #     throw(ArgumentError(string("failed to index ", $(record_type), " ~>", repr(String(data[p:min(p+7,p_end)])))))
    # # end

    if found_record
        @goto __return__
    end

end

const sam_loopcode_body = quote

    $(sam_loopcode_record)

    if finish_body
        @goto __return__
    end
end

const sam_returncode_record = quote
    return found_record
end

Automa.Stream.generate_reader(
    :index!,
    sam_machine_record,
    arguments = (:(record::Record),),
    actions = sam_actions_record,
    context = sam_context,
    initcode = sam_initcode_record,
    loopcode = sam_loopcode_record,
    returncode = sam_returncode_record
) |> eval

const sam_returncode_body = quote
    return cs, linenum, found_record
end

Automa.Stream.generate_reader(
    :readrecord!,
    sam_machine_body,
    arguments = (:(record::Record), :(state::Tuple{Int,Int})),
    actions = sam_actions_body,
    context = sam_context,
    initcode = sam_initcode_body,
    loopcode = sam_loopcode_body,
    returncode = sam_returncode_body
) |> eval


text:
"""
	colorbar(cmd0::String="", arg1=nothing; kwargs...)

Plots gray scales or color scales on maps.

Full option list at [`psscale`]($(GMTdoc)psscale.html)

- **D** | **pos** | **position** :: [Type => Str]

    Defines the reference point on the map for the color scale using one of four coordinate systems.
    ($(GMTdoc)psscale.html#d)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **F** | **box** :: [Type => Str]

    Draws a rectangular border around the scale.
    ($(GMTdoc)psscale.html#f)
- **G** | **truncate** :: [Type => Str]

    Truncate the incoming CPT so that the lowest and highest z-levels are to zlo and zhi.
    ($(GMTdoc)psscale.html#g)
- **I** | **shade** :: [Type => Number | Str]

    Add illumination effects.
    ($(GMTdoc)psscale.html#i)
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- **L** | **equal** | **equal_size** :: [Type => Str | Bool]		`Arg = [i][gap]`

    Gives equal-sized color rectangles. Default scales rectangles according to the z-range in the CPT.
    ($(GMTdoc)psscale.html#l)
- **M** | **monochrome** :: [Type => Bool]

    Force conversion to monochrome image using the (television) YIQ transformation.
    ($(GMTdoc)psscale.html#m)
- **N** | **dpi** :: [Type => Str | Number]

    Controls how the color scale is represented by the PostScript language.
    ($(GMTdoc)psscale.html#n)
- **Q** | **log** :: [Type => Str]

    Selects a logarithmic interpolation scheme [Default is linear].
    ($(GMTdoc)psscale.html#q)
- $(GMT.opt_R)
- **S** | **nolines** :: [Type => Bool | []]

    Do not separate different color intervals with black grid lines.
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** | **zscale** :: [Type => Number]

    Multiply all z-values in the CPT by the provided scale.
    ($(GMTdoc)psscale.html#w)
- **Z** | **zfile** :: [Type => Str]

    File with colorbar-width per color entry.
    ($(GMTdoc)psscale.html#z)
"""
function colorbar(cmd0::String="", arg1=nothing; first=true, kwargs...)

	gmt_proggy = (IamModern[1]) ? "colorbar "  : "psscale "
	(length(kwargs) == 0 && cmd0 != "") && return monolitic(gmt_proggy, cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	cmd, opt_B, = parse_BJR(d, "", "", O, "")
	cmd  = parse_JZ(d, cmd)[1]
	cmd, = parse_common_opts(d, cmd, [:F :UVXY :params :c :p :t], first)
	cmd  = parse_these_opts(cmd, d, [[:G :truncate], [:I :shade], [:M :monochrome], [:N :dpi],
	                                [:Q :log], [:S :nolines], [:W :zscale], [:Z :zfile]])
	#cmd = parse_type_anchor(d, cmd, [:D :pos :position])
	cmd = parse_type_anchor(d, cmd, [:D :pos :position],
							(map=("-g", arg2str, 1), outside=("J", nothing, 1), inside=("j", nothing, 1), norm=("-n", arg2str, 1), paper=("-x", arg2str, 1), anchor=("", arg2str, 2), length=("+w", arg2str), size=("+w", arg2str), justify="+j", triangles="+e", horizontal="_+h", move_annot="+m", neon="_+mc", nan="+n", offset=("+o", arg2str)), 'J')

	cmd, arg1, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', 0, arg1)
	if (!isa(arg1, GMTcpt) && !occursin("-C", cmd))	# If given no CPT, try to see if we have a current one stored in global
		if (!isempty(current_cpt[1]))
			cmd *= " -C";	arg1 = current_cpt[1]
		end
	end

	cmd = add_opt(d, cmd, 'L', [:L :equal :equal_size], (range="i", gap=""))	# Aditive
	(!occursin(" -D", cmd)) && (cmd *= " -DJMR")			#  So that we can call it with just a CPT

	r = finish_PS_module(d, gmt_proggy * cmd, "", K, O, true, arg1)
	gmt("destroy")      # Probably because of the rasters in cpt
	return r
end

# ---------------------------------------------------------------------------------------------------
colorbar(arg1; kw...) = colorbar("", arg1; first=true, kw...)
colorbar!(arg1; kw...) = colorbar("", arg1; first=false, kw...)
colorbar!(cmd0::String="", arg1=nothing; kw...) = colorbar(cmd0, arg1; first=false, kw...)

const psscale  = colorbar         # Alias
const psscale! = colorbar!        # Alias
const colorscale  = colorbar      # Alias
const colorscale! = colorbar!     # Alias


text:
using DiscreteFiltering
using LinearAlgebra
using Plots
using Symbolics
using Latexify

## Plot theme
theme(:juno)

## Parameters
# Domain
a = 0.0
b = 1.0
# domain = PeriodicIntervalDomain(a, b)
domain = ClosedIntervalDomain(a, b)

# Time
T = 1.00

## Symbolics
@variables x t

# Exact solution (heat equation, Borggaard test case)
u = t + sin(2π * x) + sin(8π * x)
u_int = t * x - 1 / 2π * cos(2π * x) - 1 / 8π * cos(8π * x)

# Exact solution (heat equation, more complicated test case)
# u = 1 + sin(t) * (1 - 8 / 10 * x^2) + exp(-t) / 15 * sin(20π * x) + 1 / 5 * sin(10x)
# u_int =
#     x + sin(t) * (x - 8 / 30 * x^3) - exp(-t) / 15 / 20π * cos(20π * x) - 1 / 50 * cos(10x)

# Compute deduced quantities
dₜ = Differential(t)
dₓₓ = Differential(x)^2
f = expand_derivatives(dₜ(u) - dₓₓ(u))
g_b = substitute(u, Dict(x => b))
g_a = substitute(u, Dict(x => a))

# for sym ∈ [:u, :u_int, :f, :g_b, :g_a]
#     open("output/$sym.tex", "w") do io
#         @eval write($io, latexify($sym))
#     end
# end

u = eval(build_function(u, x, t))
u_int = eval(build_function(u_int, x, t))
f = eval(build_function(f, x, t))
g_a = eval(build_function(g_a, t))
g_b = eval(build_function(g_b, t))

## Ode solver tolerances
# tols = (;)
# tols = (; abstol = 1e-6, reltol = 1e-4)
tols = (; abstol = 1e-9, reltol = 1e-8)
degmax = 50
λ = 1e-9

# Number of mesh points
nrefine = 5
NN = floor.(Int, 10 .^ LinRange(2, 4, nrefine))
MM = round.(Int, 9 // 10 .* NN)

# Errors
err = zeros(nrefine)
err_bar = zeros(nrefine)
err_adbc = zeros(nrefine)

## Solve
@time for i = 1:nrefine
    M = MM[i]
    N = NN[i]

    println("Solving for M = $M and N = $N")

    # Discretization
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    Δx = (b - a) / M

    # h(x) = Δx / 2
    # filter = TopHatFilter(h)

    h₀ = 3.1Δx
    h(x) = h₀ # * (1 - 1 / 2 * cos(x))
    σ = Δx / 2
    filter = GaussianFilter(h, σ)

    # Equations
    equation = DiffusionEquation(domain, IdentityFilter(), f, g_a, g_b)
    equation_filtered = DiffusionEquation(domain, filter, f, g_a, g_b)

    u_use = ξ -> u(ξ, T)

    # Exact filtered solution
    ū = apply_filter(u_use, filter, domain)

    # Exact extended-filtered solution (for ADBC)
    ū_ext = apply_filter_extend(u_use, filter, domain)

    # Solve discretized problem
    # sol = solve(
    #     equation,
    #     ξ -> u(ξ, 0.0),
    #     (0.0, T),
    #     M,
    #     N;
    #     method = "discretizefirst",
    #     boundary_conditions = "derivative",
    #     tols...,
    #     degmax,
    #     λ,
    # )

    # Solve discretized-then-filtered problem
    sol_bar = solve(
        equation_filtered,
        ξ -> u(ξ, 0.0),
        (0.0, T),
        M,
        N;
        method = "discretizefirst",
        boundary_conditions = "derivative",
        tols...,
        degmax,
        λ,
    )

    # Solve filtered-then-discretized problem with ADBC
    # ū_adbc = solve_adbc(equation_filtered, ξ -> u(ξ, 0.0), (0.0, T), M, T / 20_000_000)

    # Relative error
    u_exact = u.(ξ, T)
    ū_exact = ū.(x)
    ū_ext_exact = ū_ext.(x)
    # err[i] = norm(sol(T) - u_exact) / norm(u_ exact)
    err_bar[i] = norm(sol_bar.u[end] - ū_exact) / norm(ū_exact)
    # err_adbc[i] = norm(ū_adbc - ū_ext_exact) / norm(ū_ext_exact)
end

##
plot(x, [ū_exact, sol_bar.u[end]])

## Set GR backend for fast plotting
gr()

## Set PGFPlotsX plotting backend to export to tikz
pgfplotsx()

## Plot exact solution
ξ = LinRange(a, b, 101)
p = plot(
    xlabel = raw"$x$",
    # size = (400, 300),
    legend = :topright,
)
for t ∈ LinRange(0.0, T, 5)
    plot!(p, ξ, u.(ξ, t), label = "t = $t")
end
display(p)
savefig(p, "output/diffusion/solution.tikz")

## Plot convergence
p = plot(
    xaxis = :log10,
    yaxis = :log10,
    # size = (400, 300),
    minorgrid = true,
    legend = :topright,
    # xlims = (NN[1], NN[end]),
    ylims = (1e-7, 1e0),
    xticks = 10 .^ (1:4),
)
# plot!(p, NN, err, label = "Discretized")
plot!(p, NN, err_adbc, marker = :c, label = "Filtered-then-discretized with ADBC")
plot!(p, NN, err_bar, marker = :d, label = "Discretized-then-filtered")
plot!(p, NN, 10 ./ NN .^ 2, linestyle = :dash, label = "\$10 / N^2\$")
# plot!(p, NN, 20 ./ NN .^ 2, linestyle = :dash, label = raw"$20 N^{-2}$")
# plot!(p, NN, 10 ./ NN .^ 1.5, linestyle = :dash, label = raw"$10 N^{-3/2}}$")
xlabel!(p, raw"$N$")
# title!(p, raw"Heat equation, $h(x) = \Delta x / 2$")
display(p)
# savefig(p, "output/diffusion/convergence.tikz")


text:
function get_component_to_pm(ix::Int, bus::PSY.Bus)
    number = PSY.get_number(bus)
    zone = PSY.get_load_zone(bus)
    PM_bus = Dict{String, Any}(
        "zone" => isnothing(zone) ? "1" : PSY.get_name(zone),
        "bus_i" => number,
        "bus_type" => PM_BUSTYPES[PSY.get_bustype(bus)],
        "vmax" => PSY.get_voltage_limits(bus).max,
        "area" => 1,
        "vmin" => PSY.get_voltage_limits(bus).min,
        "index" => number,
        "va" => PSY.get_angle(bus),
        "vm" => PSY.get_magnitude(bus),
        "base_kv" => PSY.get_base_voltage(bus),
        "name" => PSY.get_name(bus),
        "source_id" => ["bus", number],
    )

    return PM_bus
end

function get_components_to_pm(sys::PSY.System, ::Type{PSY.Bus})
    buses = PSY.get_components(PSY.Bus, sys)
    PM_buses = Dict{String, Any}()

    for (ix, bus) in enumerate(buses)
        PM_buses[string(PSY.get_number(bus))] = get_component_to_pm(ix, bus)
    end
    return PM_buses
end

function get_pm_map(sys::PSY.System, ::Type{PSY.Bus})
    return Dict(string(PSY.get_number(b)) => b for b in PSY.get_components(PSY.Bus, sys))
end

max_stars_repo_name: NREL-SIIP/PowerModelsInterface

max_issues_repo_name: NREL-SIIP/PowerModelsInterface

max_forks_repo_name: NREL-SIIP/PowerModelsInterface


text:
using LinearAlgebra: normalize!

export Soliton, SolitonIdeal, SolitonRobust, get_probabilities

abstract type Soliton end

#=

Derivation of the ideal soliton distribution:

(See [mackay] Exercise 50.2)

It's easy to see that at the t-th iteration, when t out of K packets have been recovered, the
expected number of packets of degree d reduced to (d - 1) is

$$
h_t(d) (d / (K - t))
$$

So in order to have $\forall i \in \{ 0, ..., K - 1\}, h_t(1) = 1$, we start with $h_0(1) = 1$ and
$h_0(2) = K / 2$. In general,

$$
h_t(d) = (K - t) / ((d - 1) * d)
$$

This can be proved using mathematical induction, that is, because of

$$
h_{t + 1}(d) = h_t(d) * (1 - d / (K - t))   + h_t(d + 1) * ((d + 1) / (K - t))
$$

we can verify the equation holds for d = 2, and by

$$
(K - t) / ((d - 1) * d) * (1 - d / (K - t)) + h_t(d + 1) * ((d + 1) / (K - t)) = (K - (t + 1)) / ((d - 1) * d)
$$

which is simplified to

$$
h_t(d + 1) = (K - t) / (d * (d + 1))
$$

so it holds for d = 2, 3, ..., K.

Also, it's easy to see the soliton distribution sums to 1 by `1/(d * (d - 1)) = 1/(d - 1) - 1/d`.
=#

"The ideal soliton distribution."
struct SolitonIdeal <: Soliton
    probabilities::Vector{Float64}

    function SolitonIdeal(n::Int)
        probabilities = [1 / (k * (k - 1)) for k in 1:n]
        probabilities[1] = 1 / n

        new(probabilities)
    end
end

Base.length(soliton::SolitonIdeal) = length(soliton.probabilities)
Base.getindex(soliton::SolitonIdeal, i::Int) = soliton.probabilities[i]

function get_probabilities(soliton::SolitonIdeal)
    soliton.probabilities
end

const ROBUST_FAILURE_PROB = 0.01

#=

The derivation of the robust soliton distribution is similar to the ideal distribution, but it aims
to have the expected number of 1-degree packets be

$$
h_t(1) = 1 + S
$$

See [mackay] Exercise 50.4 for details.

=#

"The robust soliton distribution."
struct SolitonRobust <: Soliton
    probabilities::Vector{Float64}

    function SolitonRobust(n::Int, delta::Float64=0.01)
        # $S ≈ 2$ in [mackay] Eq. (50.4)
        m = n ÷ 2 + 1
        probabilities = zeros(Float64, n)
        probabilities[1:(m - 1)] = [1 / (i * m) for i in 1:(m - 1)]
        probabilities[m] = log(n / (m * delta)) / m

        ideal_prob = SolitonIdeal(n).probabilities
        probabilities .+= ideal_prob
        normalize!(probabilities, 1)

        new(probabilities)
    end
end

Base.length(soliton::SolitonRobust) = length(soliton.probabilities)
Base.getindex(soliton::SolitonRobust, i::Int) = soliton.probabilities[i]

function get_probabilities(soliton::SolitonRobust)
    soliton.probabilities
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
######################################## Model 1
struct MyBaseModel <: ModelName end
myparameter1 = (μ = Param(0.0, Normal()), σ = Param(10.0, Gamma()))
mymodel1 = ModelWrapper(MyBaseModel(), myparameter1)

data_uv = rand(_rng, Normal(mymodel1.val.μ, mymodel1.val.σ), N)
#Create objective for both μ and σ and define a target function for it
myobjective1 = Objective(mymodel1, data_uv, (:μ, :σ))
function (objective::Objective{<:ModelWrapper{MyBaseModel}})(θ::NamedTuple)
	@unpack data = objective
	lprior = Distributions.logpdf(Distributions.Normal(),θ.μ) + Distributions.logpdf(Distributions.Exponential(), θ.σ)
    llik = sum(Distributions.logpdf( Distributions.Normal(θ.μ, θ.σ), data[iter] ) for iter in eachindex(data))
	return lprior + llik
end
myobjective1(myobjective1.model.val)

function ModelWrappers.generate(_rng::Random.AbstractRNG, objective::Objective{<:ModelWrapper{MyBaseModel}})
    @unpack model, data = objective
    @unpack μ, σ = model.val
    return Float16(μ[1])
end

function ModelWrappers.predict(_rng::Random.AbstractRNG, objective::Objective{<:ModelWrapper{MyBaseModel}})
    @unpack model, data = objective
    @unpack μ, σ = model.val
	return rand(_rng, Normal(μ, σ))
end


text:
using NURBS
using Base.Test

@testset "Base Generator" begin

    @testset "basis" begin
        x = knot(4,4)
        @test typeof(basis(4,4,0.5,x)) == Array{Float64,1}
        @test isapprox(basis(4,4,0.5,x),[0.125, 0.375, 0.375, 0.125],atol=1e-5)
    end

    @testset "dbasis" begin
        x = knot(4,4)
        @test typeof(dbasis(4,4,0.5,x)) == Tuple{Array{Float64,1}, Array{Float64,1}, Array{Float64,1}}
        @test isapprox(dbasis(4,4,0.5,x)[1],[0.125, 0.375, 0.375, 0.125],atol=1e-5)
        @test isapprox(dbasis(4,4,0.5,x)[2],[-0.75, -0.75, 0.75, 0.75],atol =1e-5)
        @test isapprox(dbasis(4,4,0.5,x)[3],[3.0, -3.0, -3.0, 3.0],atol =1e-5)
    end

    @testset "dbasisu" begin
        x = knotu(5,3)
        @test typeof(dbasisu(5,3,0.5,x)) == Tuple{Array{Float64,1}, Array{Float64,1}, Array{Float64,1}}
        @test isapprox(dbasisu(5,3,0.5,x)[1],[0.125, 0.0, 0.0, 0.0, 0.0],atol=1e-5)
        @test isapprox(dbasisu(5,3,0.5,x)[2],[0.5, 0.0, 0.0, 0.0, 0.0],atol=1e-5)
        @test isapprox(dbasisu(5,3,0.5,x)[3],[1.0, 0.0, 0.0, 0.0, 0.0],atol=1e-5)
    end
end


text:
#=
Metadata for offline datasets included in this respository...
=#

data_manifest = OrderedDict(
"IDRC_2002_Pharma_Shootout" =>
    OrderedDict(
        "usage" => ["Calibration transfer", "Regression"],
        "keywords" => [ "API", "Pharmaceuticals", "IDRC", "shootout"],
        "property values" => ["weight","hardness","assay"],
        "measurements" => ["NIR"],
        "spectral range" => "600-1398nm",
        "spectral resolution" => "2nm",
        "observations" => "654",
        "format" => "csv",
        "description" => "\"Spectra of 654 pharmaceutical tablets from two Multitab spectrometers (FossNIRsystems, Silverspring, MD) will be the basis of the ShootOut at the next Chambersburg meeting (August 10-16, 2002). All tablets were measured in the transmittance mode.  Each individual tablet was subsequently analyzed for ASSAY value, tablet weight, and tablet hardness. All log (1/T) spectra cover the region from 600 to 1898 nm in 2 nm increments.\"... \"Data from each instrument has been split into a calibration file (155 spectra), and a validation file (459 spectra). In addition, a small test file (40 spectra) from each instrument is avaliable for testing calibrations. The challenge is to develop the best possible calibration for the ASSAY value on one instrument, and to provide the best means of transferring the calibration to the other instrument. The reference method has an error of +/- 1.3 mg on these tablets with a nominal value of 200 mg per tablet. The data set includes tablets with a wide ASSAY range, 152 to 239 mg, for developing calibration(s). No deliberate outliers are included in the data set, but there may be outliers due to experimental errors. The name of the active ingredient in these tablets can not be disclosed for proprietary reasons.\"",
        "URL" => "https://web.archive.org/web/20031204005143/http://www.idrc-chambersburg.org/shootout_2002.htm",
        "Date Accessed" => "7/3/2020",
        "references" => ["Gary E. Ritchie, Chairman NIRVWoG and Sr. Scientist, Pharmaceutical Analysis/New Technology, Purdue Pharma L.P., 444 Saw Mill River Road, Ardsley, NY 10502",
                        "http://www.foss.dk/c/p/solutions/products/showprodfamily.asp?prodfamilypkid=98"],
        "permission" => "Public domain. Permission to redistribute given to Dr. Casey Kneale via e-mail in 2018."
    ),

"Tecator_Meat_Data" =>
    OrderedDict(
        "usage" => ["Regression"],
        "keywords" => ["food", "fat", "meat"],
        "property values" => ["fat","water","protein"],
        "format" => "csv",
        "measurements" => ["NIR"],
        "spectral range" => "850-1050nm",
        "spectral resolution" => "",
        "observations" => "240",
        "description" => "Statement of permission from Tecator (the original data source).These data are recorded on a Tecator Infratec Food and Feed Analyzer working in the wavelength range 850 - 1050 nm by the Near Infrared Transmission (NIT) principle. Each sample contains finely chopped pure meat with different moisture, fat and protein contents.If results from these data are used in a publication we want you to mention the instrument and company name (Tecator) in the publication. In addition, please send a preprint of your article to Karin Thente, Tecator AB, Box 70, S-263 21 Hoganas, Sweden. The data are available in the public domain with no responsability from the original data source. The data can be redistributed as long as this permission note is attached. For more information about the instrument - call Perstorp Analytical's representative in your area.",
        "URL" => "https://www.openml.org/d/505",
        "Date Accessed" => "7/3/2020",
        "references" => [""],
        "permission" => "Public domain with statement of permission from Tecator."
    ),

"Cargill_Corn" =>
    OrderedDict(
        "usage" => ["calibration transfer", "regression"],
        "keywords" => ["food", "corn"],
        "property values" => ["moisture", "oil", "protein", "starch"],
        "format" => "csv",
        "measurements" => ["NIR"],
        "spectral range" => "1100-2498nm",
        "spectral resolution" => "2nm",
        "observations" => "80",
        "description" => "This dataset contains 80 NIR measures of corn taken from 3 different instruments( M5, MP5,and MP6 ).
        Measures of NBS glass standards (with replicates) on each instrument are also contained in seperate files.",
        "URL" => "https://eigenvector.com/resources/data-sets/",
        "Date Accessed" => "7/3/2020",
        "references" => [""],
        "permission" => "Public domain. Permission to redistribute given to Dr. Casey Kneale granted by Mike Blackburn 6/3/2020"
    ),

# "Southwest_Research_Institute_Diesel_Fuels" =>
#     OrderedDict(
#         "usage" => ["regression", "impution", "semisupervised learning"],
#         "keywords" => ["diesel", "fuel", "oil", "gas", "SWRI"],
#         "property values" => ["bp50", "CN", "d4052", "freeze", "total"],
#         "format" => "csv",
#         "measurements" => ["NIR"],
#         "spectral range" => "750-1550nm",
#         "spectral resolution" => "2nm",
#         "observations" => "785",
#         "description" => "Various diesel fuels were measured with a single NIR instrument
#         bp50 - boilingpoint (*C) at 50% recovery(ASTM D 86)
#         CN - Cetane Number (ASTM D 613)
#         d4052 - density, g/mL (15*C ASTM D 4052)
#         freeze - freezing temperature of the fuel, deg C
#         total - total aromatics (mass% ASTM D 5186)
#         visc - viscosity, cSt(40*C)",
#         "URL" => "https://eigenvector.com/wp-content/uploads/2019/06/SWRI_Diesel_NIR_CSV.zip",
#         "Date Accessed" => "7/3/2020",
#         "references" => [""],
#         "permission" => "Public Domain. Collection funded by the US Army."
#     ),

"Fruit_Puree_Adulteration" =>
    OrderedDict(
        "usage" => ["Classification", "One class classification"],
        "keywords" => ["food", "fruit", "adulteration", "interferent"],
        "format" => "csv",
        "measurements" => ["MIR"],
        "spectral range" => "899.327-1802.564 cm-1",
        "spectral resolution" => "3.86cm-1",
        "observations" => "983",
        "description" => "...\"a collection of 983 Mid-infrared spectra collected from
        different authenticated fruit purees in one of two classes:
        \"Strawberry\" (purees prepared from fresh whole fruits by the researchers)
        \"NON-Strawberry\" (diverse collection of other purees, including: strawberry adulterated with other fruits and sugar solutions; raspberry; apple; blackcurrant; blackberry; plum; cherry; apricot; grape juice and mixtures of these.\"
         ",
        "URL" => "https://csr.quadram.ac.uk/example-datasets-for-download/",
        "Date Accessed" => "7/3/2020",
        "references" => ["\"Use of Fourier transform infrared spectroscopy and partial least squares regression for the detection of adulteration of strawberry purees\" Holland JK, Kemsley EK, Wilson RH. (1998). Journal of the Science of Food and Agriculture, 76, 263-269"],
        "permission" => "BBSRC policy on data sharing. Free to analyse and redistribute the data; if the original source is cited."
    ),

"Fresh_Meats" =>
    OrderedDict(
        "usage" => ["Classification", "Authentication"],
        "keywords" => ["food", "chicken", "pork", "turkey"],
        "format" => "csv",
        "measurements" => ["MIR"],
        "spectral range" => "1005.3495-1867.864 cm-1",
        "spectral resolution" => "1.93cm-1",
        "observations" => "120",
        "description" => "...\"120 Mid-infrared (MIR) spectra of fresh minced meats – chicken, pork and turkey. Duplicate acquisitions from 60 independent samples. Raw data matrix size [448 x 120]. Obtained using Fourier transform infrared (FTIR) spectroscopy with attenuated total reflectance (ATR) sampling\"",
        "URL" => "https://csr.quadram.ac.uk/example-datasets-for-download/",
        "Date Accessed" => "7/3/2020",
        "references" => ["\"Mid-infrared spectroscopy and authenticity problems in selected meats: a feasibility study\" Al-Jowder O, Kemsley E K, Wilson R. H.(1997) Food Chemistry 59 195-20."],
        "permission" => "BBSRC policy on data sharing. Free to analyse and redistribute the data; if the original source is cited."
    ),

"Olive_Oils_Quadrum" =>
    OrderedDict(
        "usage" => ["Classification", "Provenance"],
        "keywords" => ["food", "olive", "oil"],
        "format" => "csv",
        "measurements" => ["FTIR"],
        "spectral range" => "798.892-1896.8085 cm-1",
        "spectral resolution" => "1.93 cm-1",
        "observations" => "120",
        "description" => "...\" contains a collection of 120 Mid-infrared spectra collected from 60 different authenticated extra virgin olive oils, supplied to the Institute of Food Research, UK, by the International Olive Oil Council.
        Spectra were acquired from each oil using attenuated total reflectance (ATR) sampling.  The acquisition order was randomized with respect to the country of origin code.
        Once all the samples had been examined once, a second acquisition session commenced, to produce a second spectrum from each sample.  Again the acquisition order was randomized with respect to country of origin.  thus, duplicate spectra were collected from all samples.
        \"",
        "URL" => "https://csr.quadram.ac.uk/example-datasets-for-download/",
        "Date Accessed" => "7/3/2020",
        "references" => ["\"FTIR spectroscopy and multivariate analysis can distinguish the geographic origin of extra virgin olive oils\" Tapp H.S. et al, J. Agric. Food Chem. 51 (21) 6110-5 (2003) "],
        "permission" => "BBSRC policy on data sharing. Free to analyse and redistribute the data; if the original source is cited."
    ),

"Instant_Coffee_Quadrum" =>
    OrderedDict(
        "usage" => ["Classification", "Provenance", "Fraud"],
        "keywords" => ["food", "instant", "coffee", "arabica", "robusta"],
        "format" => "csv",
        "measurements" => ["MIR", "DRIFT", "MIR-DRIFT"],
        "spectral range" => "810.548-1910.644 cm-1",
        "spectral resolution" => "3.86 cm-1",
        "observations" => "56",
        "description" => "...\"Mid-infrared (MIR) spectra of authenticated freeze-dried coffee samples (arabica and robusta species, respectively 29 and 27 of each). Raw data matrix size is [286 x 56]. Obtained by Fourier transform infrared (FTIR) spectroscopy with diffuse reflectance (DRIFT) sampling.\"",
        "URL" => "https://csr.quadram.ac.uk/example-datasets-for-download/",
        "Date Accessed" => "7/3/2020",
        "references" => [   "\"Discrimination of Arabica and Robusta in Instant Coffee by Fourier Transform Infrared Spectroscopy and Chemometrics\" Briandet R et al, J. Agric. Food Chem., 44 (1),170–174 (1996)",
                            "\"Near- and Mid-Infrared Spectroscopies in Food Authentication: Coffee Varietal Identification\" Downey G. et al, J. Agric. Food Chem. 45 (11) 4357-4361 (1997)."],
        "permission" => "BBSRC policy on data sharing. Free to analyse and redistribute the data; if the original source is cited."
    ),

"Raman_PUFA" =>
    OrderedDict(
        "usage" => ["Regression"],
        "keywords" => ["food", "polyunsaturated", "PUFA", "fat"],
        "property values" => ["PUFA%emul", "%PUFA"],
        "format" => "csv",
        "measurements" => ["Raman"],
        "spectral range" => "675-1697 cm-1",
        "spectral resolution" => "1.0 cm-1",
        "observations" => "69",
        "description" => " see citation for more details. ",
        "URL" => "https://data.mendeley.com/datasets/3ympjxywdm/1",
        "Date Accessed" => "7/4/2020",
        "references" => [ "Rimal, Raju; Almøy, Trygve; Sæbø, Solve (2019), \"Data for: Comparison of Multi-response Prediction Methods\", Mendeley Data, v1 http://dx.doi.org/10.17632/3ympjxywdm.1 "],
        "permission" => "CC BY NC 3.0."
    ),

"Dough_NIR" =>
    OrderedDict(
        "usage" => ["Regression"],
        "keywords" => ["food", "dough", "biscuit"],
        "property values" => ["fat", "sucrose", "flour", "water"],
        "format" => "csv",
        "measurements" => ["Raman"],
        "spectral range" => "1000-2498 nm",
        "spectral resolution" => "2.0 nm",
        "observations" => "73",
        "description" => "see citation for more details. ",
        "URL" => "https://data.mendeley.com/datasets/3ympjxywdm/1",
        "Date Accessed" => "7/4/2020",
        "references" => [ "Brown, P.J., Fearn, T. and Vannucci, M. (2001). \"Bayesian wavelet regression on curves with application to a spectroscopic calibration problem.\" Journal of the American Statistical Association 96, 398-408.",
                            "Rimal, Raju; Almøy, Trygve; Sæbø, Solve (2019), \"Data for: Comparison of Multi-response Prediction Methods\", Mendeley Data, v1 http://dx.doi.org/10.17632/3ympjxywdm.1 "],
        "permission" => "CC BY NC 3.0. Permission to redistribute given to Dr. Casey Kneale by Marina Vannucci.",
        "Usage Statement" => "Please acknowledge the following paper if utilizing \"Brown, P.J., Fearn, T. and Vannucci, M. (2001). \"Bayesian wavelet regression on curves with application to a spectroscopic calibration problem.\" Journal of the American Statistical Association 96, 398-408."
    ),

"Penicillin_Monitoring" =>
    OrderedDict(
        "usage" => ["Regression", "control", "monitoring", "ode"],
        "keywords" => ["fermentation", "process", "semisynthetic"],
        "property values" => ["penicillin"],
        "format" => "csv",
        "measurements" => ["Raman"],
        "spectral range" => "?-? cm-1",
        "spectral resolution" => "? cm-1",
        "observations" => "?",
        "description" => "100 batches of penicillin femernations monitored via Raman spectroscopy.
        See references for more details. ",
        "URL" => "https://data.mendeley.com/datasets/pdnjz7zz5x/2",
        "Date Accessed" => "7/4/2020",
        "references" => [ "Goldrick, Stephen (2019), \"Data for: Modern day monitoring and control challenges outlined on an industrial-scale benchmark fermentation process\", Mendeley Data, v2 http://dx.doi.org/10.17632/pdnjz7zz5x.2",
                        "Goldrick S., Duran-Villalobos C., K. Jankauskas, Lovett D., Farid S. S, Lennox B., (2019) Modern day control challenges for industrial-scale fermentation processes. Computers and Chemical Engineering."],
        "permission" => "CC BY NC 4.0. Permission given to Dr. Casey Kneale to redistribute via LinkedIn provided a statement is displayed on GitHub.",
        "Usage Statement" => "Please acknowledge the following paper if utilizing the spectral data which can be freely downloaded at www.industrialpenicillinsimulation.com. Goldrick S., Duran-Villalobos C., K. Jankauskas, Lovett D., Farid S. S, Lennox B., (2019) Modern day control challenges for industrial-scale fermentation processes. Computers and Chemical Engineering."
    )
)


text:
using Core.IR
const Compiler = Core.Compiler

let code = Any[
    Expr(:gotoifnot, SlotNumber(2), 4),
    LabelNode(2),
    Expr(:(=), SlotNumber(3), 2),
    LabelNode(4),
    # Test a SlotNumber as a value of a PhiNode
    Expr(:(=), SSAValue(0), PhiNode(Any[2,3], Any[1, SlotNumber(3)])),
    Expr(:return, SSAValue(0))
]

    ci = eval(Expr(:new, CodeInfo,
        code,
        nothing,
        Any[Any, Any, Any],
        Any[Any],
        UInt8[0, 0, 0],
        Any[Symbol("Self"), :arg, :slot],
        false, false, false, false
    ))

    Compiler.run_passes(ci, 1, Compiler.LineInfoNode[Compiler.NullLineInfo])
end

max_stars_repo_name: ederag/julia

max_issues_repo_name: ederag/julia

max_forks_repo_name: ederag/julia


text:
# Test special handling of scalars

import Devectorize
import Devectorize.@devec
import Devectorize.@inspect_devec
import Devectorize.dump_devec
import Devectorize.sqr

using Base.Test

a = [1., 2., 3., 4., 5., 6., 7., 8.]

r = zeros(1)

@devec v = 1.
@test isequal(v, 1.)

@devec r[1] = 2.
@test isequal(r, [2.])

@devec v = r[1]
@test isequal(v, 2.)

r[1] = 0.
@devec r[1] = v
@test isequal(r, [v])

r[1] = 0.
@devec r[1,1] = v
@test isequal(r, [v])

i = 2
@devec v = a[i]
@test isequal(v, a[2])

@devec v = -a[3]
@test isequal(v, -a[3])

@devec v = 2. + a[3]
@test isequal(v, 2. + a[3])

r[1] = 0.
@devec r[1] = scalar(v)
@test isequal(r, [v])

i = 2
@devec v = scalar(a[i])
@test isequal(v, a[i])

@devec v = -scalar(a[3])
@test isequal(v, -a[3])

@devec v = 2. + scalar(a[3])
@test isequal(v, 2. + a[3])


text:
module GenieCMS

greet() = print("Hello World!")

using Genie

function install(dest::String; force = false)
  src = abspath(normpath(joinpath(@__DIR__, "..", Genie.Plugins.FILES_FOLDER)))

  for f in readdir(src)
    isdir(f) || continue
    Genie.Plugins.install(joinpath(src, f), dest, force = force)
  end
end

end # module


text:
#import WavesAndEigenvalues.Meshutils: create_rotation_matrix_around_axis, get_rotated_index, get_reflected_index
using WavesAndEigenvalues.Helmholtz
## load mesh from file
case="NTNU_new"
mesh=Mesh("../NTNU/NTNU.msh",scale=1.0)
mesh=octosplit(mesh)
#create unit and full mesh from half-cell
#doms=[("Interior",:full),("Inlet",:full), ("Outlet_high",:full), ("Outlet_low",:full), ("Walls",:full),("Flame",:unit),("Bloch",:unit)]#
doms=[("Interior",:full),("Inlet",:full), ("Outlet_high",:full), ("Outlet_low",:full), ("Flame",:unit),]#

collect_lines!(mesh)
unit_mesh=extend_mesh(mesh,doms,unit=true)
full_mesh=extend_mesh(mesh,doms,unit=false)
D=Dict()
D["mesh"]=zeros(Float64,6044)
vtk_write(case*"_mesh",unit_mesh,D)
## meshing
#naxis=unit_mesh.dos.naxis
#nxbloch=unit_mesh.dos.nxbloch

## Helmholtz solver
#define speed of sound field
function speedofsound(x,y,z)
    if z<0.415
        return 347.0#m/s
    else
        return 850.0#m/s
    end
end

c=generate_field(unit_mesh,speedofsound,0)
C=generate_field(full_mesh,speedofsound,0)
#D=Dict()

#vtk_write("NTNU_c",full_mesh,D)
##describe model
dscrp=Dict()
dscrp["Interior"]=(:interior,())
dscrp["Outlet_high"]=(:admittance, (:Y_in,1E15))
dscrp["Outlet_low"]=(:admittance, (:Y_out,1E15))
##discretize models
#L=discretize(full_mesh, dscrp, C, order=:1)
l=discretize(unit_mesh, dscrp, c, order=:1,b=:b)
## solve model using beyn
#(type Γ by typing \Gamma and Enter)
Γ=[150.0+5.0im, 150.0-5.0im, 1000.0-5.0im, 1000.0+5.0im].*2*pi #corner points for the contour (in this case a rectangle)
#Ω, P = beyn(L,Γ,l=10,N=64, output=true)
l.params[:b]=1
#ω, p = beyn(l,Γ,l=10,N=4*128, output=true)
##
sol, n, flag = householder(l,914*2*pi,maxiter=16, tol=0*1E-10,output = true, n_eig_val=3)

##
# identify surface points
surface_points, tri_mask, tet_mask =get_surface_points(unit_mesh)
# compute normal vectors on surface triangles
normal_vectors=get_normal_vectors(unit_mesh)
## DA approach
blochify_surface_points!(unit_mesh, surface_points, tri_mask, tet_mask)
sens=discrete_adjoint_shape_sensitivity(unit_mesh,dscrp,c,surface_points,tri_mask,tet_mask,l,sol)

## correct sens for bloch formalism
#sens[:,unit_mesh.dos.naxis+1:unit_mesh.dos.naxis+unit_mesh.dos.nxbloch].+=sens[:,end-unit_mesh.dos.nxbloch+1:end]
sens[:,end-unit_mesh.dos.nxbloch+1:end]=sens[:,unit_mesh.dos.naxis+1:unit_mesh.dos.naxis+unit_mesh.dos.nxbloch]

##
# normalize point sensitivity with directed surface triangle area
normed_sens=normalize_sensitivity(surface_points,normal_vectors,tri_mask,sens)
# boundary-mass-based-normalizations
nsens = bound_mass_normalize(surface_points,normal_vectors,tri_mask,unit_mesh,sens)
# compute sensitivity in unit normal direction
normal_sens = normal_sensitivity(normal_vectors, normed_sens)

## save data
DD=Dict()
DD["real normed_sens1"]=real.(normed_sens[1,:])
DD["real normed_sens2"]=real.(normed_sens[2,:])
DD["real normed_sens3"]=real.(normed_sens[3,:])
DD["real normal_sens"]=real.(normal_sens)
DD["imag normed_sens1"]=imag.(normed_sens[1,:])
DD["imag normed_sens2"]=imag.(normed_sens[2,:])
DD["imag normed_sens3"]=imag.(normed_sens[3,:])
DD["imag normal_sens"]=imag.(normal_sens)
vtk_write_tri(case*"_tri",unit_mesh,DD)

##
mode="[$(string(round((sol.params[:ω]/2/pi),digits=2)))]Hz"
data=Dict()
data["real DA x"*mode]=real.(sens[1,:])
data["real DA y"*mode]=real.(sens[2,:])
data["real DA z"*mode]=real.(sens[3,:])
data["imag DA x"*mode]=imag.(sens[1,:])
data["imag DA y"*mode]=imag.(sens[2,:])
data["imag DA z"*mode]=imag.(sens[3,:])
data["real nDA x"*mode]=real.(nsens[1,:])
data["real nDA y"*mode]=real.(nsens[2,:])
data["real nDA z"*mode]=real.(nsens[3,:])
data["imag nDA x"*mode]=imag.(nsens[1,:])
data["imag nDA y"*mode]=imag.(nsens[2,:])
data["imag nDA z"*mode]=imag.(nsens[3,:])
#data["finite difference"*mode]=real.(sens_FD)
#data["central finite difference"]=real.(sens_CFD)
#data["diff DA -- FD"*mode]=abs.((sens.-sens_FD))
vtk_write(case, unit_mesh, data)
data=Dict()
v=bloch_expand(unit_mesh,sol,:b)
data["abs_p"*mode]=abs.(v)./maximum(abs.(v))
data["phase_p"*mode]=angle.(v)
#data["abs_p_adj"*mode]=abs.(sol.v_adj)./maximum(abs.(sol.v_adj))
#data["phase_p_adj"*mode]=angle.(sol.v_adj)
vtk_write(case*"_modes", full_mesh, data)

## FD approach

fd_sens=forward_finite_differences_shape_sensitivity(unit_mesh,dscrp,c,surface_points,tri_mask,tet_mask,l,sol)
fd_sens[:,end-unit_mesh.dos.nxbloch+1:end]=fd_sens[:,unit_mesh.dos.naxis+1:unit_mesh.dos.naxis+unit_mesh.dos.nxbloch]

## write output
mode="[$(string(round((sol.params[:ω]/2/pi),digits=2)))]Hz"
data=Dict()
data["abs_p"*mode]=abs.(sol.v)./maximum(abs.(sol.v))
data["phase_p"*mode]=angle.(sol.v)
data["abs_p_adj"*mode]=abs.(sol.v_adj)./maximum(abs.(sol.v_adj))
data["phase_p_adj"*mode]=angle.(sol.v_adj)
data["real DA x"*mode]=real.(fd_sens[1,:])
data["real DA y"*mode]=real.(fd_sens[2,:])
data["real DA z"*mode]=real.(fd_sens[3,:])
data["imag DA x"*mode]=imag.(fd_sens[1,:])
data["imag DA y"*mode]=imag.(fd_sens[2,:])
data["imag DA z"*mode]=imag.(fd_sens[3,:])
#data["finite difference"*mode]=real.(fd_sens_FD)
#data["central finite difference"]=real.(sens_CFD)
#data["diff DA -- FD"*mode]=abs.((sens.-sens_FD))
vtk_write(case*"_fd", unit_mesh, data)
##
findmax(abs.(sens-fd_sens))


text:
using ThreeBodyDecay

const mΛb = 5.61960
const mπ = 0.14
# intermediate resonances
const mΣb = 5.81065; const ΓΣb = 0.005
const mΣb_x = 5.83032; const ΓΣb_x = 0.0094

const mΛb2S = 6.3; # just a peak of the plot

tbs = ThreeBodySystem(mπ,mΛb,mπ; m0=mΛb2S,
    two_jps = ([0, 1, 0, 1], ['-','+','-','+']))  # 1/2+ 0- 0- 1/2+
dpp = randomPoint(tbs)

# lineshape
dc_Σb1 = DecayChainLS(1, σ->BW(σ,mΣb,ΓΣb); tbs=tbs, two_s = 1, parity='+')
dc_Σb3 = DecayChainLS(3, σ->BW(σ,mΣb,ΓΣb); tbs=tbs, two_s = 1, parity='+')
full_AΣb(σs,two_λs)   = sum(amplitude(σs, two_λs, ch) for ch in [dc_Σb1, dc_Σb3])

dc_Σb_x1 = DecayChainLS(1, σ->BW(σ,mΣb_x,ΓΣb_x); tbs=tbs, two_s = 3, parity='+')
dc_Σb_x3 = DecayChainLS(3, σ->BW(σ,mΣb_x,ΓΣb_x); tbs=tbs, two_s = 3, parity='+')
full_AΣb_x(σs,two_λs) = sum(amplitude(σs, two_λs, ch) for ch in [dc_Σb_x1, dc_Σb_x3])

dc_f0 = DecayChainLS(2, σ->1.0; tbs=tbs, two_s = 0, parity='+')
full_Af0(σs,two_λs)   = amplitude(σs, two_λs, dc_f0)

####

Φij = let
    As = [full_AΣb, full_AΣb_x, full_Af0]
    [let
        Aij(σs,two_λs)=Ai(σs,two_λs) * conj(Aj(σs,two_λs))
        I = summed_over_polarization(Aij, tbs.two_js)
        three_body_phase_space_integral(I, tbs)
        end for (Ai,Aj) in Iterators.product(As,As)]
end

using LinearAlgebra
Φij_n = let
    d = sqrt.(diag(Φij))
    Φij ./ d ./ d'
end

###
one_AΣb(σs,two_λs)   = amplitude(σs, two_λs, dc_Σb1)
one_AΣb_x(σs,two_λs) = amplitude(σs, two_λs, dc_Σb_x1)

Φ0ij = let
    As = [one_AΣb, one_AΣb_x, full_Af0]
    [let
        Aij(σs,two_λs)=Ai(σs,two_λs) * conj(Aj(σs,two_λs))
        I = summed_over_polarization(Aij, tbs.two_js)
        three_body_phase_space_integral(I, tbs)
        end for (Ai,Aj) in Iterators.product(As,As)]
end

Φ0ij_n = Φ0ij .* [√2; √2; 1] .* [√2 √2 1]
Φ0ij_n

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
struct samples
    samples :: Vector{Vector}
end

function num_samples(s::samples)
    length(s.samples)
end


text:
# This file was generated, do not modify it. # hide
r = report(m)

figure(figsize=(8,6))

res = r.plotting

vals_sf = res.parameter_values[:, 1]
vals_bf = res.parameter_values[:, 2]

tricontourf(vals_sf, vals_bf, res.measurements)
xticks(1:3:12, fontsize=12)
xlabel("Number of sub-features", fontsize=14)
yticks(0.4:0.2:1, fontsize=12)
ylabel("Bagging fraction", fontsize=14)

savefig(joinpath(@OUTPUT, "A-ensembles-2-heatmap.svg")) # hide


text:
using SparseArrays
using LinearAlgebra
using SparseMatricesCSR
using Test

include("SparseMatrixCSC.jl")
include("SparseMatrixCSR.jl")
include("SymSparseMatrixCSR.jl")


text:
module PsychoJulia
include("clocks.jl")
include("visual.jl")
include("logic.jl")
end#module


text:
## define structures
struct lPreconditioner{T}
    f::T
    function lPreconditioner(pfunc)
        f = pfunc
        return new{typeof(f)}(f)
    end
end

struct CholT
    ld::SparseMatrixCSC{Float64, Int64}
    ldT::SparseMatrixCSC{Float64, Int64}
    d::Vector{Float64}
    p::Vector{Int64}
    invp::Vector{Int64}
end

struct HierarchyLevel
    sd::Bool
    islast::Bool
    iterative::Bool
    A::SparseMatrixCSC{Float64, Int64}
    invD::Vector{Float64}
    cI::Vector{Int64}
    nc::Int64
    n::Int64
    nnz::Int64
    #chol::CholT
    chol::LDLFactorizations.LDLFactorization{Float64, Int64, Int64, Int64}
end

struct Hierarchy
    A::SparseMatrixCSC{Float64, Int64}
    invD::Vector{Float64}
    cI::Vector{Int64}
    #chol::CholT
    chol::LDLFactorizations.LDLFactorization{Float64, Int64, Int64, Int64}
end

mutable struct LevelAux
    fwd::Bool
    rc::Int32
    repeat::Int32
    #
    sd::Bool
    islast::Bool
    iterative::Bool
    n::Int64
    nc::Int64
end

struct Workspace
    x::Vector{Float64}
    b::Vector{Float64}
    tmp::Vector{Float64}
end

## Hierarchy Structures Initialization

function init_LevelAux(H::Vector{HierarchyLevel})
    n_levels = length(H)
    LevelAux_ = Vector{LevelAux}(undef, n_levels)

    @inbounds for j in 1:n_levels
        if j==1
            repeat = Int32(1)
        elseif j==n_levels
            repeat = Int32(0)
        else
            repeat = Int32(max(floor(nnz(H[j-1].A)/nnz(H[j].A) - 1), 1))
        end

        LevelAux_[j] = LevelAux(true,1, repeat, H[j].sd, H[j].islast, H[j].iterative, H[j].n, H[j].nc)
    end

    return LevelAux_
end

function init_Workspace(H::Vector{HierarchyLevel})
    n_levels = length(H)
    Workspace_ = Vector{Workspace}(undef, n_levels)

    @inbounds for j in 1:n_levels
        n_ = size(H[j].A, 1)
        if j==n_levels && !H[j].iterative
            n_ = n_ + 1
        end
        Workspace_[j] = Workspace(zeros(n_), zeros(n_), zeros(n_))
    end

    return Workspace_
end

function init_Hierarchy(H::Vector{HierarchyLevel})
    n_levels = length(H)
    Hierarchy_ = Vector{Hierarchy}(undef, n_levels)

    @inbounds for j in 1:n_levels
        Hierarchy_[j] = Hierarchy(H[j].A, H[j].invD, H[j].cI, H[j].chol)
    end

    return Hierarchy_
end

##

function cmg_preconditioner_lap(A::SparseMatrixCSC)
    cmg_!(A, A)
end

function cmg_preconditioner_adj(A::SparseMatrixCSC)
    msg = undef
    flag = undef
    A_lap = lap(A)

    if size(A_lap, 1) < 500 # handle small input
        throw(ArgumentError("Input Matrix is Small. Solve Ax=B with A\\b"))
    else # validate input
        flag, A_lap_ = validateInput!(A_lap)
        if flag == 1
            throw(ArgumentError("Input Matrix Must Be Symmetric!"))
        elseif flag == 2
            throw(ArgumentError("Current Version of CMG Does Not Support Positive Off-Diagonals!"))
        end
    end

    cmg_!(A_lap, A_lap_)
end

function cmg_!(A::T, A_::T) where {T<:SparseMatrixCSC}
    A_o = A
    flag = Int64(0)
    sflag = Int64(0)
    loop = true
    sd = true
    iterative = true
    j = Int64(0)
    h_nnz = Int64(0)
    n = Int64(0)
    H = Vector{HierarchyLevel}()

    if size(A_, 1) > size(A, 1)
        sflag = 1
        sd = true
    else
        sflag = 1
        sd = false
    end

    # build up H
    while loop
        iterative = true
        n = size(A_, 1)
        # direct method for small size
        if (n < 500)
            iterative = false
            break
        end
        dA_ = Array(diag(A_))
        (cI, ~) = steiner_group(A_, dA_)
        nc = maximum(cI)
        islast = false
        A = A_ # !
        invD = 1 ./ (2 * dA_) # !
        R = sparse(cI, 1:n, ones(n), nc, n) # added for efficiency

        if nc == 1
            islast = true
            iterative = true
            flag = 1
        end

        # check for hierarchy stagnation for potentially bad reasons
        h_nnz = h_nnz + nnz(A_)
        if (nc >= n - 1) || (h_nnz > 5 * nnz(A_o))
            islast = true
            iterative = true
            flag = 3 # indicates stagnation
            @warn "CMG convergence may be slow due to matrix density. Future versions of CMG will eliminate this problem."
            break
        end

        Rt = sparse(cI, 1:n, 1, nc, n) # ! take out double
        A_ = (Rt * A) * Rt'
        push!(H, HierarchyLevel(sd, islast, iterative, A, invD, cI, nc, n, nnz(A), ldl([1.0 0; 0 1.0])))

        if sflag == 1
            sd = true
            sflag = 0
        end

        if nc == 1
            break
        end
    end

    # code for last hierarchy level
    if flag == 0
        j += 1
        B = A_[1:n - 1, 1:n - 1]
        ldlt = ldl(B)
        push!(H, HierarchyLevel(true, true, false,B, Float64[], Int64[], 0, n, nnz(ldlt.L), ldlt))
    end

    X = init_LevelAux(H)
    W = init_Workspace(H)
    M = init_Hierarchy(H)

    # create precondition function
    pfunc = make_preconditioner(M, W, X)

    return (pfunc, H)
end

function findRowSumAndDominance(A::SparseMatrixCSC)
    colptr::Vector{Int64} = A.colptr
    nzval::Vector = A.nzval
    rowval::Vector = A.rowval
    n = length(nzval)
    s = Int64(1)
    flag = Int64(0)
    k = Int64(0)
    il = Vector{Int64}(undef, n)
    jl = Vector{Int64}(undef, n)
    vl = Vector{Float64}(undef, n)
    sumR = zeros(Float64, length(colptr)-1)

    @inbounds for i in 1:length(colptr)-1
        l = colptr[i+1] - colptr[i]
        t = s + l - 1
        sum = 0
        for j in s:t
            val = nzval[j]
            sum += val
            row = rowval[j]
            k += 1
            il[k] = row
            jl[k] = i
            vl[k] = val
            if row != i && flag == 0
                if val > 0
                    flag = 2
                end
            end
        end
        sumR[i] = sum
        s += l
    end

    return sumR, flag, il, jl, vl
end

function validateInput!(A::SparseMatrixCSC)
    flag = Int64(0)
    # check symmetry
    if !issymmetric(A)
        return 1, A
    end
    # detect strict dominance && positive off diagonals
    n = size(A, 1)
    sAp = Vector{Float64}(undef, n)
    sd = Vector{Int64}(undef, n)
    dA = diag(A)
    (sA, flag, i, j, v) = findRowSumAndDominance(A)

    if flag == 2
        return 2, A
    end

    @inbounds @simd for i in 1:length(sA)
        sAp[i] = (sA[i] + abs(sA[i]))/2
        sd[i] = (sAp[i]/dA[i]) > 1e-13 ? 1 : 0
    end

    # augment by extra coordinate if strictly dominant
    if maximum(sd) > 0.0
        ex_v = -sAp[sd]
        ex_v_sum = -sum(ex_v)
        exd = length(ex_v)
        exd_c = findall(!iszero, sd) #nonzeros(sd) # get coordinates
        i_ = ones(Int64, exd) * (n+1)
        i = vcat(i, i_, exd_c, n+1)
        j = vcat(j, exd_c, i_, n+1)
        v = vcat(v, ex_v, ex_v, ex_v_sum)
        A = sparse(i, j, v, n + 1, n + 1)
    end
    return 3, A
end

"""
    function(cI, nc) = steiner_group(A, dA_)
    Steiner groups
    # Arguments
    -`A`: Laplacian
"""

function steiner_group(A::SparseMatrixCSC, dA_::Vector{Float64})
    (C, M) = findMinSparse(A)
    split_forest_!(C)
    efd = abs.(M ./ dA_)
    if minimum(efd) < 1 / 8 # low effective degree nodes found
        C = update_groups_(A, C, dA_)
    end
    #return C, efd
    (cI, nc, ~) = forest_components_(C)
end

"""
    function C1 = update_groups_(A, C, dA_)
    update groups based on nodes with low effective degree
    # Arguments
    -`A`: Laplacian
"""
function update_groups_(A::SparseMatrixCSC, C::Vector{Int64}, dA_::Vector{Float64})
    n = length(C)
    B = zeros(Float64, n)
    # B[j] is the total tree weight incident to node j
    @inbounds for i in 1:n
        if C[i] != i
            B[i] = A[i, C[i]] + B[i]
            B[C[i]] = A[i,C[i]] + B[C[i]]
        end
    end
    ndx = findall(x->x > -0.125, B ./ dA_)
    @inbounds @simd for i in 1:length(ndx)
        C[ndx[i]] = Int32(ndx[i])
    end
    return C
end # update_groups_

"""
    function[cI, nc, csizes] = forest_components_(C)
    forest components, connected components in unimodal forest
    # Arguments
    -`C`: unimodal tree
"""

function forest_components_(C::Vector{Int64})
    n = length(C)
    cI = zeros(Int64, n)
    cSizes = zeros(Int64, n)
    buffer = zeros(Int64, 100)
    ccI = 1

    @inbounds for j in 1:n
        bufferI = 1
        jwalk = j
        # tree walk
        while cI[jwalk] == 0
            cI[jwalk] = ccI
            buffer[bufferI] = jwalk
            bufferI = bufferI + 1
            if bufferI == size(buffer, 1) # for memory efficiency
                buffer_ = zeros(Int64, min(2 * size(buffer, 1), n))
                buffer_[1:size(buffer, 1)] = buffer
                buffer = buffer_
            end
            jwalk = C[jwalk]
        end # while

        bufferI = bufferI - 1
        en = C[jwalk] # end node
        if cI[en] != ccI
            @simd for i in 1:bufferI
                cI[buffer[i]] = cI[en]
            end
        else
            ccI = ccI + 1
        end
        cSizes[en] = cSizes[en] + bufferI
    end # for
    if cSizes[ccI] == 0
        ccI = ccI - 1
    end

    cSizes = cSizes[1:ccI]
    return cI, ccI, cSizes
end

"""
   function pfun = make_preconditioner(H)
Make preconditioner
"""
function make_preconditioner(H::Vector{Hierarchy}, W::Vector{Workspace}, X::Vector{LevelAux})
    if !X[1].sd
        pfun = b->preconditioner_i(H, W, X, b)
        return pfun
    end

    if X[1].sd
        pfun = b->preconditioner_sd(b, H, X,  W)
        return pfun
    end
end

"""
   function x = preconditioner_sd(b,H)
preconditioner sd
"""
function preconditioner_sd(b::Array{Float64}, H::Array{Hierarchy}, X::Vector{LevelAux}, W::Vector{Workspace})
    n = length(b)
    push!(b, -sum(b))
    x = b->preconditioner_i(H, W, X, b)
    x = x[1:n] + x[n + 1]
    return x
end

"""
   function C = split_forest_(C1::Array{Int})
decompose unimodal forest into low conductance components
"""

function split_forest_!(C1::Vector{Int64})
    n = length(C1)
    C = C1
    new_front = Int64(0)
    removed_ancestors = Int64(0)
    k = Int64(0)
    jwalk = Int64(0)
    jwalka = Int64(0)
    jwalkb = Int64(0)
    cut_mode = false
    startwalk = false
    ancestors_in_path = Int64(0)
    ancestors = zeros(Int64, n)
    indegree = zeros(Int64, n + 2)
    visited = falses(n) # logical sparse array
    walkbuffer = zeros(Int64, 20)
    newancestorbuff = zeros(Int64, 20)

   # compute indegrees
    for i in 1:n
        indegree[C[i]] = indegree[C[i]] + 1
    end

   # partition into clusters of small diameter

    for j in 1:n
        jwalk = j
        startwalk = true

        while (startwalk && (indegree[jwalk] == 0 ) && !visited[jwalk])
            startwalk = false
            ancestors_in_path = 0 # change over C-CMG
            k = 1
            walkbuffer[k] = jwalk
            newancestorbuff[k] = 0
            while (k <= 6 || visited[jwalk]) # +1 for c-indexing adjust
                jwalk = C[jwalk]
                walkterminated = (jwalk == walkbuffer[k]) || ((k > 1) && (jwalk == walkbuffer[k - 1]))
                if walkterminated
                    break # while
                end
                k += 1
                walkbuffer[k] = jwalk
                if visited[jwalk]
                    newancestorbuff[k] = ancestors_in_path
                else
                    ancestors_in_path = ancestors_in_path + 1
                    newancestorbuff[k] = ancestors_in_path
                end
            end

            if k > 6 # large diameter - cut
                middlek = Int64(ceil(k / 2))
                C[walkbuffer[middlek]] = walkbuffer[middlek] # cut middle edge
                indegree[walkbuffer[middlek + 1]] = indegree[walkbuffer[middlek + 1]] - 1 # update indegree

                for ik in (middlek + 1):k
                    ancestors[walkbuffer[ik]] = ancestors[walkbuffer[ik]] - ancestors[walkbuffer[middlek]]
                end

            # update ancestors and visited flag
                for ik in 1:middlek
                    visited[walkbuffer[ik]] = true
                    ancestors[walkbuffer[ik]] = ancestors[walkbuffer[ik]] + newancestorbuff[ik]
                end

            # set first vertex in new walk
                jwalk = walkbuffer[middlek + 1]
                startwalk = true
            end # end cut procedure

         # commit walk changes
            if !startwalk
                for ik in 1:k
                    ancestors[walkbuffer[ik]] = ancestors[walkbuffer[ik]] + newancestorbuff[ik]
                    visited[walkbuffer[ik]] = true
                end
            end
        end # outer while
    end

   # tree partition into clusters of high conductance
    for j in 1:n
        jwalk = j
        startwalk = true

        while startwalk && (indegree[jwalk] == 0)
            startwalk = false
            jwalkb = jwalk
            cut_mode = false
            # initialize new_front
            new_front = 0
            removed_ancestors = 0

            while true
                jwalka = C[jwalk]
                walkterminated = (jwalka == jwalk) || (jwalka == jwalkb)
                if walkterminated
                    break; # while
                end

                if (!cut_mode && (ancestors[jwalk] > 2) && (ancestors[jwalka] - ancestors[jwalk] > 2)) # possibly low conductance - make cut
                    C[jwalk] = jwalk # cut edge
                    indegree[jwalka] = indegree[jwalka] - 1
                    removed_ancestors = ancestors[jwalk]
                    new_front = jwalka
                    cut_mode = true
                end # end making cut

                jwalkb = jwalk
                jwalk = jwalka
                if cut_mode
                    ancestors[jwalk] = ancestors[jwalk] - removed_ancestors
                end
            end
            if cut_mode
                startwalk = true
                jwalk = new_front
            end
        end
    end
end # split_forest_

function findMinSparse(A::SparseMatrixCSC)
    minval = Inf
    row_ = Int64(0)
    colptr = A.colptr
    nzval = A.nzval
    rowval = A.rowval
    (~, m) = size(A)
    n = length(nzval)
    min_cols = zeros(Int64, m)
    min_vals = zeros(Float64, m)

    for i in 1:m
        r_loc_1 = colptr[i]
        r_loc_2 = colptr[i+1]
        size_l = r_loc_2-r_loc_1

        if size_l == 0
            continue
        end
        minval = Inf
        row_ = 0

        for j in r_loc_1:r_loc_2-1
            val = nzval[j]
            row = rowval[j]
            if val < minval
                minval = val
                row_ = row
            end
        end

        min_cols[i] = row_
        min_vals[i] = minval
    end
    return (min_cols, min_vals)
end

function findSumSparse(A::SparseMatrixCSC)
    colptr = A.colptr
    nzval = A.nzval
    s = 1
    m = length(colptr)-1
    k = 0
    sumVec = Vector{Float64}(undef, m)

    @inbounds for i in 1:m
        l = colptr[i+1] - colptr[i]
        t = s + l - 1
        if l > 0
            k += 1
        end
        sum = 0.0
        @simd for j in s:t
            val = nzval[j]
            sum += val
        end
        sumVec[k] = sum
        s += l
    end
    return sumVec
end

@inline function interpolate!(x::Vector{Float64}, cI::Vector{Int64}, z::Vector{Float64})
    @inbounds @simd for i in 1:length(x)
        x[i] = 0
    end

    @inbounds @simd for i in 1:length(z)
        idx = cI[i]
        x[idx] += z[i]
    end
end

function preconditioner_i(H::Vector{Hierarchy},  W::Vector{Workspace}, X::Vector{LevelAux}, b::Vector{Float64})
    level = Int64(1)
    #@inbounds W[1].b = b
    BLAS.blascopy!(length(b),b,1,W[1].b,1)

    while level > 0
        x = W[level].x
        b = W[level].b
        tmp = W[level].tmp

        invD = H[level].invD
        A = H[level].A
        cI = H[level].cI

        if X[level].islast && !X[level].iterative
            @inbounds ldiv!(x, H[level].chol,b)
            level -= 1
        elseif X[level].islast && X[level].iterative
            W[level].x .= b .* invD
            level -= 1
        elseif !X[level].islast && X[level].fwd
            repeat = X[level].repeat   #number of 'recursive' calls

            if X[level].rc > repeat
                X[level].rc = 1
                level -= 1
            else
                if X[level].rc == 1
                    W[level].x .= b.* invD
                else
                    mul!(tmp, A, x)
                    tmp .-= b
                    tmp .*= -invD
                    W[level].x .+= tmp
                end

                mul!(tmp, A, x)
                tmp .-= b
                X[level].fwd = false
                interpolate!(W[level+1].b, cI, -tmp)
                level += 1
            end
        elseif !X[level].islast && !X[level].fwd
            z = W[level+1].x
            W[level].x .+= z[cI]
            mul!(tmp, A, x)
            tmp .-= b
            tmp .*= -invD
            W[level].x .+= tmp
            X[level].rc += 1
            X[level].fwd = true
        end
    end

    return W[1].x
end


text:
import RCall.sexp
import RCall.rcopy
import RCall: RClass, rcopytype

## Julia Tuple will be convert to S3 class JuliaTuple based on R list

function sexp(x ::Tuple)
    r = sexp(Array{Any}([x...]))
    setattrib!(r, :class, sexp("JuliaTuple"))
    r
end

function rcopy(::Type{Tuple}, x::Ptr{VecSxp})
    Tuple(rcopy(Array, x))
end

rcopytype(::Type{RClass{:JuliaTuple}}, x::Ptr{VecSxp}) = Tuple

max_stars_repo_name: KevCaz/JuliaCall

max_issues_repo_name: tinagui/JuliaCall

max_forks_repo_name: tinagui/JuliaCall


text:
push!(LOAD_PATH, string(pwd(),"/src"))

using Revise
using SpecialFunctions
using ModifiedStokesSolver
using Base.Test

@testset "NewPowerSeriesTests" begin

    z = 0.99
    #(T1Sp, T1Lp, T2Sp, T2Lp, T3Sp, T3Lp) =
    Tp = ModifiedStokesSolver.Ti_split_pow(z)
    #(T1Sd, T1Ld, T2Sd, T2Ld, T3Sd, T3Ld) =
    Td = ModifiedStokesSolver.Ti_split_direct(z)
    Ep = Tp .- Td
    @test maximum(abs.(Ep./Td)) < 9e-14

    z = 1.8
    #(T1Sp, T1Lp, T2Sp, T2Lp, T3Sp, T3Lp) =
    Tp = ModifiedStokesSolver.Ti_split_pow(z,true)
    #(T1Sd, T1Ld, T2Sd, T2Ld, T3Sd, T3Ld) =
    Td = ModifiedStokesSolver.Ti_split_direct(z,true)
    Ep = Tp .- Td
    @test maximum(abs.(Ep./Td)) < 9e-14

end

# @testset "PowerSeries" begin

# # Series evaluator
# function evalseries(fn, z, nstart)
#     y = 0.0
#     dy = 1.0
#     n = nstart
#     while abs(dy)>=eps(y)
#         dy = fn(z, n)*z^(2*n)
#         y += dy
#         n += 1
#     end
#     return y
# end
# # Definitions for modified Bessel functions
# i0(z,n) = 1/(4^n*factorial(n)^2)
# i1(z, n) = 1/(2*4^n*factorial(n)*factorial(n+1))
# k0(z,n) = ((digamma(n+1) + log(2) - log(z)) / (4^n * factorial(n)^2))
# k1(z, n) = (2*log(z) - 2*log(2) - digamma(n+1) - digamma(n+2)) /
#     (4^(n+1)*factorial(n)*factorial(n+1))
# # Splits
# k0S(z,n) = ((digamma(n+1) + log(2)) / (4^n * factorial(n)^2))
# k0L(z,n) = (-1) / (4^n * factorial(n)^2)
# k1S(z, n) = ( - 2*log(2) - digamma(n+1) - digamma(n+2)) /
#     (4^(n+1)*factorial(n)*factorial(n+1))
# k1L(z, n) = (2) / (4^(n+1)*factorial(n)*factorial(n+1))

# # Direct power series for split K0, K1
# # Returns: (K0S, K0L, K1S, K1L, K1C)
# function besselpowsplit(z)
#     # Outputs
#     K0S = 0
#     K0L = 0
#     K1S = 0
#     K1L = 0
#     K1C = 1
#     # Recursion variables
#     n = 0
#     z2 = z*z
#     z2n = 1; # z^(2*n)
#     digamma_np1 = -eulergamma # digamma(n+1)
#     digamma_np2 = -eulergamma + 1
#     ifact_n = 1 # 1/factorial(n)
#     ifact_np1 = 1
#     i4n = 1 # 1/4^n
#     i4np1 = 1/4
#     log2 = log(2)
#     # Loop
#     maxit = 20
#     while n<maxit
#         k0n_denom = i4n*ifact_n*ifact_n
#         k1n_denom = i4np1*ifact_n*ifact_np1
#         k0Ln = (-1)*k0n_denom
#         k0Sn = (digamma_np1 + log2)*k0n_denom
#         k1Ln = (2)*k1n_denom
#         k1Sn = (-2*log2-digamma_np1-digamma_np2)*k1n_denom
#         K0S += k0Sn*z2n
#         K0L += k0Ln*z2n
#         K1S += k1Sn*z2n*z
#         K1L += k1Ln*z2n*z
#         # Update recursions
#         n += 1
#         z2n *= z2;
#         inp1 = 1/(n+1)
#         (digamma_np1, digamma_np2) = (digamma_np2, digamma_np2+inp1)
#         (ifact_n, ifact_np1) = (ifact_np1, ifact_np1*inp1)
#         (i4n, i4np1) = (i4np1, i4np1*0.25)
#     end
#     return (K0S, K0L, K1S, K1L, K1C)
# end

# # Tests
# z = 1/10

# ref = besseli(0, z)
# I0test = evalseries(i0, z, 0)
# @test abs((ref-I0test)/ref) < 1e-14
# ref = besseli(1, z)
# I1test = z*evalseries(i1, z, 0)
# @test abs((ref-I1test)/ref) < 1e-14
# ref = besselk(0, z)
# K0test = evalseries(k0, z, 0)
# K0split = evalseries(k0S, z, 0) + log(z)*evalseries(k0L, z, 0)
# @test abs((ref-K0test)/ref) < 1e-14
# @test abs((ref-K0split)/ref) < 1e-14
# ref = besselk(1, z)
# K1test = 1/z + z*evalseries(k1, z, 0)
# K1split = 1/z + z*evalseries(k1S, z, 0) + z*log(z)*evalseries(k1L, z, 0)
# @test abs((ref-K1test)/ref) < 1e-14
# @test abs((ref-K1split)/ref) < 1e-14
# # Test direct split power series
# ref = (besselk(0, z), besselk(1, z))
# (K0S, K0L, K1S, K1L, K1C) = besselpowsplit(z)
# split = (log(z)*K0L+K0S, K1S+log(z)*K1L+K1C/z)
# @test maximum(@. abs((ref-split)/ref) ) < 1e-14

# # Evaluate T, high precision reference from Mathematica
# for n=1:2
#     if n==1
#         z = 0.001
#         T1exact = 0.00028941085628129279453492498775927
#         T2exact = -318.30984639510755467178629422316
#         T3exact = -0.00059871605412101857968611091116825
#     end
#     if n==2
#         z = 0.99
#         T1exact = 0.023658849757364503594228422585464
#         T2exact = -0.28952030153253490875314823065265
#         T3exact = -0.063320117528153811164290571635020
#     end
#     # z = 1
#     # (T1exact, T2exact, T3exact) =
#     #     (0.023620976507380996090769874829016,
#     #      -0.28607692796680641340443170447384,
#     #      -0.063358432123254121248207660793628)

#     # z = 2
#     # (T1exact, T2exact, T3exact) =
#     #     (0.016930005787967333888517613314308,
#     #      -0.11224071607388210664375507421808,
#     #      -0.057317125084941282339599571205834)

#     println("z = ", z)
#     println()

#     K0 = besselk(0, z)
#     K1 = besselk(1, z)
#     K0S = K0 + besseli(0, z)*log(z)
#     K1S = K1 - besseli(1, z)*log(z) - 1/z
#     I0 = besseli(0,z)
#     I1 = besseli(1,z)

#     ( T1Spow, T1Lpow,
#       T2Spow, T2Lpow, T2Cpow,
#       T3Spow, T3Lpow) = ModifiedStokesSolver.Ti_split_pow(z)

#     ## T1
#     T1 = -( z^2*(2*K0 + z*K1) + 4*(z*K1 - 1) ) / (2*pi*z*z^2)
#     t1(z, n) = 2*k0(z, n) + 4*k1(z, n) + k1(z, n-1)
#     t1S(z, n) = 2*k0S(z, n) + 4*k1S(z, n) + k1S(z, n-1)
#     t1L(z, n) = 2*k0L(z, n) + 4*k1L(z, n) + k1L(z, n-1)
#     T1series = -evalseries(t1,z,1)/(2*pi*z)
#     T1Sseries = -evalseries(t1S,z,1)/(2*pi*z)
#     T1Lseries = -evalseries(t1L,z,1)/(2*pi*z)
#     T1splitseries = (T1Sseries + T1Lseries*log(z))
#     T1pow = (T1Spow + T1Lpow*log(z))

#     #quit()

#     T1S = -(2*z*K0S + (z^2+4)*K1S + z) / (2*pi*z^2)
#     T1L = (2*z*I0 - (z^2+4)*I1) / (2*pi*z^2)
#     T1split = (T1S + T1L*log(z))

#     println("== T1")
#     println("exact  = ", T1exact)
#     println("series = ", T1series, "  relerr = ", (T1series-T1exact)/T1exact)
#     T1spliterr = (T1splitseries-T1exact)/T1exact
#     println("spl+ser= ", T1splitseries, "  relerr = ", T1spliterr)
#     #@test abs(T1spliterr)<1e-14
#     println("regular= ", T1, "  relerr = ", (T1-T1exact)/T1exact)
#     println("reg+spl= ", T1split, "  relerr = ", (T1split-T1exact)/T1exact)
#     #@test abs( (T1series-T1exact)/T1exact )<1e-14
#     T1powerr = (T1pow-T1exact)/T1exact
#     println("final  = ", T1pow, "  relerr = ", T1powerr)
#     #@test abs(T1powerr)<1e-14
#     println()

#     ## T2
#     T2 = (4*z^2*K0 + (z^2+8)*z*K1 - 8) / (pi*z^3)
#     t2(z, n) = 4*k0(z, n) + 8*k1(z, n) + k1(z, n-1)
#     T2series = -1/(pi*z) + evalseries(t2,z,1)/(pi*z)
#     println("== T2")
#     println("exact  = ", T2exact)
#     println("series = ", T2series, "  relerr = ", (T2series-T2exact)/T2exact)
#     println("regular= ", T2, "  relerr = ", (T2-T2exact)/T2exact)
#     @test abs( (T2series-T2exact)/T2exact )<1e-14
#     T2pow = (T2Spow + T2Lpow*log(z) + T2Cpow/z)
#     T2powerr = (T2pow-T2exact)/T2exact
#     println("final  = ", T2pow, "  relerr = ", T2powerr)
#     @test abs(T2powerr)<1e-14
#     println()

#     ## T3
#     T3 = (z*K1 - 1) / (2*pi*z)
#     t3(z, n) = k1(z, n)
#     T3series = z*evalseries(t3,z,0)/(2*pi)
#     println("== T3")
#     println("exact  = ", T3exact)
#     println("series = ", T3series, "  relerr = ", (T3series-T3exact)/T3exact)
#     println("regular= ", T3, "  relerr = ", (T3-T3exact)/T3exact)
#     @test abs( (T3series-T3exact)/T3exact )<1e-14
#     T3pow = (T3Spow + T3Lpow*log(z))
#     T3powerr = (T3pow-T3exact)/T3exact
#     println("final  = ", T3pow, "  relerr = ", T3powerr)
#     @test abs(T3powerr)<1e-14

#     println()
# end

#end # End testset


text:
#' ---
#' title: N-Queens
#' ---

#' **Originally Contributed by**: Matthew Helm

#' The N-Queens problem involves placing N queens on an N x N chessboard such that none of the queens attacks another. In chess, a
#' queen can move vertically, horizontally, and diagonally so there cannot be more than one queen on any given row, column, or
#' diagonal.

#'

#' *Note that none of the queens above are able to attack any other as a result of their careful placement.*

using GLPK
using JuMP
using LinearAlgebra

# N-Queens
N = 8

model = Model(GLPK.Optimizer);

#' Next, let's create an N x N chessboard of binary values. 0 will represent an empty space on the board and 1 will represent a
#' space occupied by one of our queens:

@variable(model, x[i=1:N, j=1:N], Bin)

#' Now we can add our constraints:

# There must be exactly one queen in a given row/column
for i=1:N
    @constraint(model, sum(x[i, :]) == 1)
    @constraint(model, sum(x[:, i]) == 1)
end

# There can only be one queen on any given diagonal
for i in -(N-1):(N-1)
    @constraint(model, sum(diag(x,i)) <= 1)
    @constraint(model, sum(diag(reverse(x,dims=1), i)) <=1)
end

#' That's it! We are ready to put our model to work and see if it is able to find a feasible solution:

optimize!(model)

#' We can now review the solution that our model found:

solution = convert.(Int,value.(x))

#'


text:
using ApproxFunOrthogonalPolynomials, ApproxFunBase, Test, SpecialFunctions, LinearAlgebra
import ApproxFunBase: testbandedbelowoperator, testbandedoperator, testspace, testtransforms, Vec,
                    maxspace, NoSpace, hasconversion, testfunctional,
                    reverseorientation, ReverseOrientation
import ApproxFunOrthogonalPolynomials: jacobip

@testset "Ultraspherical" begin
    @testset "Conversion" begin
        # Tests whether invalid/unimplemented arguments correctly throws ArgumentError
        @test_throws ArgumentError Conversion(Ultraspherical(2), Ultraspherical(1))
        @test_throws ArgumentError Conversion(Ultraspherical(3), Ultraspherical(1.9))
    end
end