text:
@testset "FASTQ" begin
    @test isa(FASTQ.Record("1", dna"AA", UInt8[10, 11]), FASTQ.Record)
    @test isa(FASTQ.Record("1", "desc.", dna"AA", UInt8[10, 11]), FASTQ.Record)
    @test_throws ArgumentError FASTQ.Record("1", dna"AA", UInt8[10])

    output = IOBuffer()
    writer = FASTQ.Writer(output)
    write(writer, FASTQ.Record("1", dna"AN", UInt8[11, 25]))
    write(writer, FASTQ.Record("2", "high quality", dna"TGA", UInt8[40, 41, 45]))
    flush(writer)
    @test String(take!(output)) == """
    @1
    AN
    +
    ,:
    @2 high quality
    TGA
    +
    IJN
    """

    output = IOBuffer()
    writer = FASTQ.Writer(output, quality_header=true)
    write(writer, FASTQ.Record("1", dna"AN", UInt8[11, 25]))
    write(writer, FASTQ.Record("2", "high quality", dna"TGA", UInt8[40, 41, 45]))
    flush(writer)
    @test String(take!(output)) == """
    @1
    AN
    +1
    ,:
    @2 high quality
    TGA
    +2 high quality
    IJN
    """

    @testset "Record" begin
        record = FASTQ.Record()
        @test !isfilled(record)

        record = FASTQ.Record("""
        @SRR1238088.1.1 HWI-ST499:111:D0G94ACXX:1:1101:1173:2105
        AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA
        +SRR1238088.1.1 HWI-ST499:111:D0G94ACXX:1:1101:1173:2105
        @BCFFFDFHHHHHJJJIJIJJIJJJJJJJJIJJJJIIIJJJIJJJ
        """)
        @test isfilled(record)
        @test FASTQ.hasidentifier(record) == hasseqname(record) == true
        @test FASTQ.identifier(record) == seqname(record) == "SRR1238088.1.1"
        @test FASTQ.hasdescription(record)
        @test FASTQ.description(record) == "HWI-ST499:111:D0G94ACXX:1:1101:1173:2105"
        @test FASTQ.hassequence(record) == hassequence(record) == true
        @test FASTQ.sequence(DNASequence, record) == dna"AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA"
        @test FASTQ.sequence(record) == sequence(record) == dna"AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA"
        @test FASTQ.sequence(String, record) == "AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA"
        @test FASTQ.hasquality(record)
        @test FASTQ.quality(record) == b"@BCFFFDFHHHHHJJJIJIJJIJJJJJJJJIJJJJIIIJJJIJJJ" .- 33

        record = FASTQ.Record("""
        @SRR1238088.1.1
        AAGCTCATGACCCGTCTTACCTACACCCTTGACGAGATCGAAGGA
        +
        @BCFFFDFHHHHHJJJIJIJJIJJJJJJJJIJJJJIIIJJJIJJJ
        """)
        @test isfilled(record)
        @test !FASTQ.hasdescription(record)
    end

    function test_records(rs1, rs2)
        if length(rs1) != length(rs2)
            return false
        end
        for (r1, r2) in zip(rs1, rs2)
            if FASTQ.identifier(r1) != FASTQ.identifier(r2) ||
               FASTQ.sequence(r1)   != FASTQ.sequence(r2)   ||
               FASTQ.quality(r1)    != FASTQ.quality(r2)
                return false
            end
        end
        return true
    end

    function test_fastq_parse(filename, valid)
        # Reading from a reader
        reader = open(FASTQ.Reader, filename)
        @test eltype(reader) == FASTQ.Record
        if valid
            for record in reader end
            @test true  # no error
            @test close(reader) === nothing
        else
            @test_throws Exception begin
                for record in reader end
            end
            return
        end

        # in-place parsing
        reader = open(FASTQ.Reader, filename)
        record = eltype(reader)()
        try
            while true
                read!(reader, record)
            end
        catch ex
            close(reader)
            if !isa(ex, EOFError)
                rethrow()
            end
        end

        # Check round trip
        output = IOBuffer()
        writer = FASTQ.Writer(output)
        expected_entries = FASTQ.Record[]
        for record in open(FASTQ.Reader, filename)
            write(writer, record)
            push!(expected_entries, record)
        end
        flush(writer)

        seekstart(output)
        read_entries = FASTQ.Record[]
        for record in FASTQ.Reader(output)
            push!(read_entries, record)
        end

        return test_records(expected_entries, read_entries)
    end

    get_bio_fmt_specimens()
    path = joinpath(dirname(@__FILE__), "..", "BioFmtSpecimens", "FASTQ")
    for specimen in YAML.load_file(joinpath(path, "index.yml"))
        tags = split(get(specimen, "tags", ""))
        valid = get(specimen, "valid", true)
        # currently unsupported features
        if any(t ∈ tags for t in ["gaps", "rna", "comments", "linewrap"])
            continue
        end
        filename = specimen["filename"]
        test_fastq_parse(joinpath(path, filename), valid)
    end

    @testset "invalid quality encoding" begin
        # Sanger full range (note escape characters before '$' and '\')
        record = FASTQ.Record("""
        @FAKE0001 Original version has PHRED scores from 0 to 93 inclusive (in that order)
        ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTAC
        +
        !"#\$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~
        """)

        # the range is not enough in these encodings
        for encoding in (:solexa, :illumina13, :illumina15)
            @test_throws ErrorException FASTQ.quality(record, encoding)
        end

        # the range is enough in these encodings
        for encoding in (:sanger, :illumina18)
            @test FASTQ.quality(record, encoding) == collect(0:93)
        end
    end

    @testset "fill ambiguous nucleotides" begin
        input = IOBuffer("""
        @seq1
        ACGTNRacgtnr
        +
        BBBB##AAAA##
        """)
        @test FASTQ.sequence(first(FASTQ.Reader(input, fill_ambiguous=nothing))) == dna"ACGTNRACGTNR"
        seekstart(input)
        @test FASTQ.sequence(first(FASTQ.Reader(input, fill_ambiguous=DNA_A)))   == dna"ACGTAAACGTAA"
        seekstart(input)
        @test FASTQ.sequence(first(FASTQ.Reader(input, fill_ambiguous=DNA_G)))   == dna"ACGTGGACGTGG"
        seekstart(input)
        @test FASTQ.sequence(first(FASTQ.Reader(input, fill_ambiguous=DNA_N)))   == dna"ACGTNNACGTNN"
        seekstart(input)
        @test FASTQ.sequence(BioSequence{DNAAlphabet{2}}, first(FASTQ.Reader(input, fill_ambiguous=DNA_A))) == dna"ACGTAAACGTAA"
    end
end

@testset "Quality scores" begin
    @testset "Decoding base quality scores" begin
        function test_decode(encoding, values, expected)
            result = Array{Int8}(undef, length(expected))
            FASTQ.decode_quality_string!(encoding, values, result, 1, length(result))
            @test result == expected
        end

        test_decode(FASTQ.SANGER_QUAL_ENCODING,
                    UInt8['!', '#', '$', '%', '&', 'I', '~'],
                    Int8[0, 2, 3, 4, 5, 40, 93])

        test_decode(FASTQ.SOLEXA_QUAL_ENCODING,
                    UInt8[';', 'B', 'C', 'D', 'E', 'h', '~'],
                    Int8[-5, 2, 3, 4, 5, 40, 62])

        test_decode(FASTQ.ILLUMINA13_QUAL_ENCODING,
                    UInt8['@', 'B', 'C', 'D', 'E', 'h', '~'],
                    Int8[0, 2, 3, 4, 5, 40, 62])

        test_decode(FASTQ.ILLUMINA15_QUAL_ENCODING,
                    UInt8['C', 'D', 'E', 'h', '~'],
                    Int8[3, 4, 5, 40, 62])

        test_decode(FASTQ.ILLUMINA18_QUAL_ENCODING,
                    UInt8['!', '#', '$', '%', '&', 'I', '~'],
                    Int8[0, 2, 3, 4, 5, 40, 93])
    end

    @testset "Encoding base quality scores" begin
        function test_encode(encoding, values, expected)
            result = Array{UInt8}(undef, length(expected))
            FASTQ.encode_quality_string!(encoding, values, result, 1, length(result))
            @test result == expected
        end

        test_encode(FASTQ.SANGER_QUAL_ENCODING,
                    Int8[0, 2, 3, 4, 5, 40, 93],
                    UInt8['!', '#', '$', '%', '&', 'I', '~'])

        test_encode(FASTQ.SOLEXA_QUAL_ENCODING,
                    Int8[-5, 2, 3, 4, 5, 40, 62],
                    UInt8[';', 'B', 'C', 'D', 'E', 'h', '~'])

        test_encode(FASTQ.ILLUMINA13_QUAL_ENCODING,
                    Int8[0, 2, 3, 4, 5, 40, 62],
                    UInt8['@', 'B', 'C', 'D', 'E', 'h', '~'])

        test_encode(FASTQ.ILLUMINA15_QUAL_ENCODING,
                    Int8[3, 4, 5, 40, 62],
                    UInt8['C', 'D', 'E', 'h', '~'])

        test_encode(FASTQ.ILLUMINA18_QUAL_ENCODING,
                    Int8[0, 2, 3, 4, 5, 40, 93],
                    UInt8['!', '#', '$', '%', '&', 'I', '~'])
    end
end


text:
doors = falses(100)
for a in 1:100, b in a:a:100
    doors[b] = !doors[b]
end
for a = 1:100
    println("Door $a is " * (doors[a] ? "open." : "closed."))
end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@inline add_fpfp_dd(x::T, y::T) where {T<:AbstractFloat} = two_sum(x, y)
@inline sub_fpfp_dd(x::T, y::T) where {T<:AbstractFloat} = two_diff(x, y)
@inline mul_fpfp_dd(x::T, y::T) where {T<:AbstractFloat} = two_prod(x, y)
@inline dvi_fpfp_dd(x::T, y::T) where {T<:AbstractFloat} = two_dvi(x, y)
# seems to help
@inline add_fpfp_dd(x::T, y::T) where {T<:IEEEFloat} = two_sum(x, y)
@inline sub_fpfp_dd(x::T, y::T) where {T<:IEEEFloat} = two_diff(x, y)
@inline mul_fpfp_dd(x::T, y::T) where {T<:IEEEFloat} = two_prod(x, y)
@inline dvi_fpfp_dd(x::T, y::T) where {T<:IEEEFloat} = two_dvi(x, y)


text:
using FEM

# Load the VTK exporter
FEM.vtkexportmod()
using FEM.VTKExportMod

# Generates a mesh of the shape known as the "Cook membrane" which is a
# quadraterial with corners [0.0, 0.0], [0.0, 44.0], [48.0, 60.0], [48.0, 44.0]
# Possible mesh elements are
# - GeoTrig for 3 node triangles,
# - GeoQTrig for 6 node triangles
# - GeoQuad for 4 node quadraterials
n_ele_x = 100
n_ele_y = 100
geomesh = gencook(n_ele_x, n_ele_y, GeoQTrig)

# We create two node sets, one on the right edge and one on the left.
# This is done by giving an anonymous function that is satisfied by
# the edges.
push!(geomesh, gennodeset(n->n.coords.x>47.999999, "right", geomesh.nodes))
push!(geomesh, gennodeset(n->n.coords.x<0.00001, "left", geomesh.nodes))

# We create an element set containing all the elements
push!(geomesh, ElementSet("all", collect(1:length(geomesh.elements))))

# We create a material section of a linear isotropic material
# and assigns that section to all the elements.
mat_section = MaterialSection(FEM.linearisotropicmod().LinearIsotropic(1, 0.3))
push!(mat_section, geomesh.element_sets["all"])

# We create an element section and assign it to all the elements.
ele_section = ElementSection(FEM.quadtrigmod().QuadTrig)
push!(ele_section, geomesh.element_sets["all"])

# Apply Dirichlet BC to the left side in dofs for x-displacement (Du)
# and dofs in y-displacement (Dv)
bcs = Any[DirichletBC("0.0", [FEM.Du], geomesh.node_sets["left"]),
          DirichletBC("0.0", [FEM.Dv], geomesh.node_sets["left"])]

# Since we currently don't have edge load we give a nodal load
# on the node set on the right edge.
loads = Any[NodeLoad("1/($n_ele_y+1)", [FEM.Dv], geomesh.node_sets["right"])]

# Create the fe problem
fp = FEM.create_feproblem("cook_example_quad", geomesh, [ele_section], [mat_section], bcs, loads)

vtkexp = VTKExporter()
# Output fields are added by pushing them into the exporter.
# We want to export the stress and strain tensor as well as
# the von Mises stress so we push them.
push!(vtkexp, Stress)
push!(vtkexp, Strain)
push!(vtkexp, VonMises)
# We set the output of the exporter to be in ascii
set_binary!(vtkexp, false)

solver = NRSolver(rel_tol = 1e-7, max_iters = 4)

# Solve the fe problem using the solver and using the vtk exporter.
solve(solver, fp, vtkexp)


text:
using CUDAnative
import CUDAdrv: CuArray
using StaticArrays

mutable struct conv2dWinograd{T}
    fil::AbstractArray{T,4}
    fil_trans::AbstractArray{T,4}
    img::AbstractArray{T,4}
    dest::AbstractArray{T,4}
    function conv2dWinograd(fil::AbstractArray{T,4}, img::AbstractArray{T,4}) where T
        fil_k, fil_w, fil_h, fil_c=size(fil)
        img_n, img_w, img_h, img_c =size(img)
        @assert img_c == fil_c
        fil_trans=similar(img,fil_k,4,4,fil_c)
        dest=similar(img,img_n,img_w,img_h ,fil_k)
        return new{T}(fil,fil_trans,img,dest)
    end
end

@inbounds function fil_trans_kernel(fil::AbstractArray{T,4}, fil_trans::AbstractArray{T,4}, channels) where T
    c= (blockIdx().x - 1) * blockDim().x + threadIdx().x
    k= blockIdx().y
    if c > channels
        return
    end
    fil_ck=zeros(MMatrix{4,4,T})
    for i in 1:3
        for j in 1:3
            fil_ck[j,i]=fil[k,j,i,c]
        end
    end
    mid = zeros(MMatrix{4,3,T})
    for i=1:3
        mid[1,i]= fil_ck[1,i]
        mid[2,i]= T(.5) * (fil_ck[1,i] + fil_ck[2,i] + fil_ck[3,i])
        mid[3,i]= mid[2,i] - fil_ck[2,i]
        mid[4,i]= fil_ck[3,i]
    end
    for i=1:4
        fil_ck[i,1]=mid[i,1]
        fil_ck[i,2]=T(.5) * (mid[i,1] + mid[i,2] + mid[i,3])
        fil_ck[i,3]=fil_ck[i,2] - mid[i,2]
        fil_ck[i,4]=mid[i,3]
    end
    for i in 1 % UInt32:4 % UInt32
        for j in 1 % UInt32:4 % UInt32
            fil_trans[k,j,i,c]=fil_ck[j,i]
        end
    end
    return
end

function forward(layer::conv2dWinograd, img::AbstractArray{T,4}) where T
    layer.img=img
    img_n,img_w,img_h,img_c=size(img)
    fil_k,fil_w,fil_h,fil_c=size(layer.fil)
    @cuda blocks=(cld(fil_c,32),fil_k) threads =32 fil_trans_kernel(layer.fil,layer.fil_trans,fil_c)
    synchronize()
    @cuda blocks=(cld(img_w,8) * cld(img_h,8),cld(img_n,32),fil_k) threads=(32,4,4) shmem=sizeof(T) * (32 * (10 * 10) + 4 * 4)  img_trans_kernel(layer.img,layer.fil_trans,layer.dest,img_w,img_h,img_c)
end
#
# function forward_new(layer::conv2dWinograd, img::AbstractArray{T,4}) where T
#     layer.img=img
#     img_c,img_w,img_h,img_n=size(img)
#     fil_c,fil_w,fil_h,fil_k=size(layer.fil)
#     @cuda blocks=(cld(fil_c,32),fil_k) threads =32 fil_trans_kernel(layer.fil,layer.fil_trans,fil_c)
#     synchronize()
#     @cuda blocks=(fld(img_w,8) * fld(img_h,8),fil_k,img_n) threads=(32,4,4) shmem=sizeof(T) * 32 * (12 * 12 + 4 * 4) img_trans_kernel(layer.img,layer.fil_trans,layer.dest,(fld(img_w,12),fld(img_h,12)),(img_w,img_h),img_c)
# end

# @inbounds function img_trans_kernel(img::AbstractArray{T,4}, fil_trans::AbstractArray{T,4}, dest::AbstractArray{T,4}, superTile_sz::NTuple{2,Int}, img_sz::NTuple{2,Int}, channels::Int) where T
#     temp=zeros(MMatrix{4,4,T})
#     csum=zeros(MMatrix{2,2,T})
#     mid=zeros(MMatrix{4,4,T})
#     superTile_h,superTile_w=superTile_sz
#     img_h,img_w=img_sz
#     n=blockIdx().z
#     k=blockIdx().y
#     z=threadIdx().z
#     y=threadIdx().y
#     tile_id_h=2 * (threadIdx().z - 1)
#     tile_id_w=2 * (threadIdx().y - 1)
#     sTile_id_h= 6 * div(blockIdx().x - 1,superTile_w) + tile_id_h
#     sTile_id_w= 6 * rem(blockIdx().x - 1,superTile_w) + tile_id_w
#     c=threadIdx().x
#     img_n_shared = @cuStaticSharedMem(T, (32, 8, 8))
#     fil_k_shared = @cuStaticSharedMem(T, (32, 4, 4))
#     for ch in c:32:channels
#         for i in 1:2, j in 1:2
#             x_h=sTile_id_h + i - 1
#             x_w=sTile_id_w + j - 1
#             if 0 < x_w <= img_w && 0 < x_h <= img_h
#                 img_n_shared[c,tile_id_w + j,tile_id_h + i]=img[ch,x_w,x_h,n]
#             else
#                 img_n_shared[c,tile_id_w + j,tile_id_h + i]=zero(T)
#             end
#         end
#         fil_k_shared[c,z,y]=fil_trans[ch,z,y]
#         sync_threads()
#         if y < 4 && z < 4
#             for i in 1:4, j in 1:4
#                 temp[j,i]=img_n_shared[c,tile_id_w + j,tile_id_h + i]
#             end
#             for i in 1:4
#                 mid[1,i]= temp[1,i] - temp[3,i]
#                 mid[2,i]= temp[2,i] + temp[3,i]
#                 mid[3,i]= temp[3,i] - temp[2,i]
#                 mid[4,i]= temp[2,i] - temp[4,i]
#             end
#             for i in 1:4
#                 temp[i,1]=(mid[i,1] - mid[i,3]) * fil_k_shared[c,i,1]
#                 temp[i,2]=(mid[i,2] + mid[i,3]) * fil_k_shared[c,i,2]
#                 temp[i,3]=(mid[i,3] - mid[i,2]) * fil_k_shared[c,i,3]
#                 temp[i,4]=(mid[i,2] - mid[i,4]) * fil_k_shared[c,i,4]
#             end
#             for i in 1:4
#                 mid[1,i]= temp[1,i] + temp[2,i] + temp[3,i]
#                 mid[2,i]= temp[2,i] - temp[3,i] - temp[4,i]
#             end
#             for i in 1:2
#                 csum[i,1] += mid[i,1] + mid[i,2] + mid[i,3]
#                 csum[i,2] += mid[i,2] - mid[i,3] - mid[i,4]
#             end
#         end
#     end
#     sync_threads()
#     for i in 1:2,j in 1:2
#         if z < 4 && y < 4
#             csum[j,i] = csum[j,i] + shfl_down(csum[j,i],16)
#             csum[j,i] = csum[j,i] + shfl_down(csum[j,i],8)
#             csum[j,i] = csum[j,i] + shfl_down(csum[j,i],4)
#             csum[j,i] = csum[j,i] + shfl_down(csum[j,i],2)
#             csum[j,i] = csum[j,i] + shfl_down(csum[j,i],1)
#             if c == 1
#                 x_w= sTile_id_w + j
#                 x_h= sTile_id_h + i
#                 if x_w <= img_w && x_h <= img_h
#                     dest[k,x_w,x_h,n]=csum[j,i]
#                 end
#             end
#         end
#     end
#     return
# end

# @inbounds function img_trans_kernel(img::AbstractArray{T,4}, fil_trans::AbstractArray{T,4}, dest::AbstractArray{T,4}, img_w::Int, img_h::Int, channels::Int, filters::Int) where T
#     tile_x, tile_y=divrem(blockIdx().x - 1, cld(img_w,2))
#     tile_x += 1
#     tile_y += 1
#     n=(blockIdx().y - 1) * blockDim().x + threadIdx().x
#     k=blockIdx().z
#     nx=threadIdx().x
#     fil_k=@cuDynamicSharedMem(T,(4,4,channels))
#     temp=zeros(MMatrix{4,4,T})
#     csum=zeros(MMatrix{2,2,T})
#     mid=zeros(MMatrix{4,4,T})
#     for c in 1:channels
#         if nx <= 16
#             fil_k[nx + (c - 1) * 16]=fil_trans[k + filters * (nx - 1 + (c - 1) * 16)]
#         end
#         for i in 1:4, j in 1:4
#             x_h=2 * tile_x + i - 1
#             x_w=2 * tile_y + j - 1
#             if 0 < x_h <= img_h && 0 < x_w <= img_w
#                 temp[j,i]=img[n,x_w,x_h,c]
#             else
#                 temp[j,i]=zero(T)
#             end
#         end
#         for i in 1:4
#             mid[1,i]= temp[1,i] - temp[3,i]
#             mid[2,i]= temp[2,i] + temp[3,i]
#             mid[3,i]= temp[3,i] - temp[2,i]
#             mid[4,i]= temp[2,i] - temp[4,i]
#         end
#         for i in 1:4
#             temp[i,1]=(mid[i,1] - mid[i,3]) * fil_k[i,1,c]
#             temp[i,2]=(mid[i,2] + mid[i,3]) * fil_k[i,2,c]
#             temp[i,3]=(mid[i,3] - mid[i,2]) * fil_k[i,3,c]
#             temp[i,4]=(mid[i,2] - mid[i,4]) * fil_k[i,4,c]
#         end
#         for i in 1:4
#             mid[1,i]= temp[1,i] + temp[2,i] + temp[3,i]
#             mid[2,i]= temp[2,i] - temp[3,i] - temp[4,i]
#         end
#         for i in 1:2
#             csum[i,1] += mid[i,1] + mid[i,2] + mid[i,3]
#             csum[i,2] += mid[i,2] - mid[i,3] - mid[i,4]
#         end
#     end
#     for i in 1:2,j in 1:2
#         x_h=2 * tile_x + i
#         x_w=2 * tile_y + j
#         if 0 < x_h <= img_h && 0 < x_w <= img_w
#             dest[n,x_w,x_h,k]=csum[j,i]
#         end
#     end
#     return
# end

@inbounds function img_trans_kernel(img::AbstractArray{T,4}, fil_trans::AbstractArray{T,4}, dest::AbstractArray{T,4}, stride_w::Int) where T
    _, img_w, img_h, _ =size(img)
    channels =size(fil_trans,4)
    n=(blockIdx().x - 1) * blockDim().x + threadIdx().x - 1
    nn=threadIdx().x - 1
    k=blockIdx().z - 1
    z=threadIdx().z - 1
    y=threadIdx().y - 1
    t_h,t_w= divrem(blockIdx().y - 1,stride_w)
    img_n= @cuStaticSharedMem(T, (32,12,12))
    fil_k = @cuStaticSharedMem(T,(4,4))
    csum=zeros(MMatrix{2,2,T})
    temp1=zeros(MMatrix{4,4,T})
    temp2=zeros(MMatrix{4,4,T})
    temp3=zeros(MMatrix{4,4,T})
    temp4=zeros(MMatrix{2,4,T})
    for c in 0:channels - 1
        for i in 0:2, j in 0:2
            x_h=8 * t_h + 3 * y + i - 1
            x_w=8 * t_w + 3 * z + j - 1
            img_n[nn + 1, 3 * z + j + 1 ,3 * y + i + 1]= 0 <= x_h < img_h && 0 <= x_w < img_w ? img[n + 1,x_w + 1,x_h + 1,c + 1] : zero(T)
        end
        if threadIdx().x - 1 == 0
            fil_k[z + 1,y + 1]=fil_trans[k + 1,z + 1,y + 1,c + 1]
        end
        sync_threads()
        for i in 1:4, j in 1:4
            temp1[j,i]=img_n[nn + 1,2 * z + j,2 * y + i]
        end
        for i in 1:4
            temp2[1,i]= temp1[1,i] - temp1[3,i]
            temp2[2,i]= temp1[2,i] + temp1[3,i]
            temp2[3,i]= temp1[3,i] - temp1[2,i]
            temp2[4,i]= temp1[2,i] - temp1[4,i]
        end
        for i in 1:4
            temp3[i,1]=(temp2[i,1] - temp2[i,3]) * fil_k[i,1]
            temp3[i,2]=(temp2[i,2] + temp2[i,3]) * fil_k[i,2]
            temp3[i,3]=(temp2[i,3] - temp2[i,2]) * fil_k[i,3]
            temp3[i,4]=(temp2[i,2] - temp2[i,4]) * fil_k[i,4]
        end
        for i in 1:4
            temp4[1,i]= temp3[1,i] + temp3[2,i] + temp3[3,i]
            temp4[2,i]= temp3[2,i] - temp3[3,i] - temp3[4,i]
        end
        for i in 1:2
            csum[i,1] += temp4[i,1] + temp4[i,2] + temp4[i,3]
            csum[i,2] += temp4[i,2] - temp4[i,3] - temp4[i,4]
        end
    end
    for i in 0:1,j in 0:1
        x_w=8 * t_w + 2 * z + j
        x_h=8 * t_h + 2 * y + i
        if x_w < img_w && x_h < img_h
            dest[n + 1,x_w + 1,x_h + 1,k + 1]=csum[j,i]
        end
    end
    return
end

max_stars_repo_name: JerryLingjieMei/WinogradConvolution

max_issues_repo_name: JerryLingjieMei/WinogradConvolution

max_forks_repo_name: JerryLingjieMei/WinogradConvolution


text:
module Scattering

using Parameters                #
using DocStringExtensions       # For simplifying docstring
using FastGaussQuadrature       # For fast Gauss-Legendre quadrature points
using JLD2                      # For saving and loading
using ProgressMeter             # For showing progress
using KernelAbstractions        # For heterogeneous (GPU+CPU) programming
using CUDA                      # For GPU programming
using Distributions             # Distributions from Julia
using ForwardDiff, DiffResults  # Automatic Differentiation tools
using LinearAlgebra
using StaticArrays
using StatsBase

using ..Architectures: device

include("Types/types.jl")                # All types used in this module
include("Mie/mie_helper_functions.jl") # Mie file-related functions
include("Legendre/legendre_functions.jl")   # Recursions for associated Legendre Polynomials
# include("mie_computations.jl")     # Functions for Mie calculations over size distribution
include("Truncation/phase_truncation.jl")     # Functions for truncation
include("Mie/wigner3j_recursive.jl")   # Recursive Wigner 3j calculations
include("Mie/mie_model.jl")
include("Mie/compute_NAI2.jl")
include("Mie/compute_PCW.jl")
include("Mie/phase_function_autodiff.jl")     # Auto-differentiation

# Export make functions/types
export make_univariate_aerosol, make_mie_model

# Export types
export NAI2, PCW, UnivariateAerosol, MieModel, Stokes_IQUV, Stokes_I, Stokes_IQU,
       δBGE, GreekCoefs, AerosolOptics, AbstractFourierDecompositionType

export compute_B, compute_ab, GreekCoefs, comp_ab, compute_mie_π_τ!,
       compute_wigner_values, save_wigner_values, load_wigner_values,
       compute_Sl, gausslegendre, compute_aerosol_optical_properties,
       compute_ref_aerosol_extinction,
       ConjugateTransposePairs, AbstractPolarizationType,
       AbstractAerosolType, AbstractAerosolType, MieModel,
       AbstractTruncationType

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 624310f3-37b1-49aa-9089-ad860f4690a8
begin
  using Pkg
  Pkg.activate(Base.Filesystem.homedir() * "/.config/julia/projects/oft")
end

# ╔═╡ 6a5292a2-c3ad-11eb-125a-3f1bbc0b7a2c
using Images#, FileIO

# ╔═╡ 4c03a12c-4128-4da6-a6e1-ef91afa02614
begin
  path_image = "./dog04.jpeg"
  image = load(path_image)
end

# ╔═╡ 541f8b1f-2773-442f-8104-9b430afef1b0
size(image)

# ╔═╡ e2eddba9-3c11-43f6-bfd8-722c66bfbf9d
image[1,1]

# ╔═╡ afbeff8b-f0f1-4ed2-a09e-3848db5c47fb
typeof(image)

# ╔═╡ 8e60dfbc-9a80-41e9-9091-67910940aeec
RGB(1,0,0)

# ╔═╡ 0e5ddc64-3943-45d5-ac8d-e00e7e10c1cf
typeof(RGB(1,0,0))

# ╔═╡ 0323a47c-3dcc-4162-b82b-d0504a7b029d
typeof(1)

# ╔═╡ 5446ce8a-5143-4209-a5ea-bd2b5bd20e1a
typeof(image[1,1])

# ╔═╡ 39595f13-7a31-401d-bc56-5a1106f567bc
n0f8(1)

# ╔═╡ 094184ad-325a-4b17-8e7a-89a29b5df473
image[1,1].r

# ╔═╡ a0b12e5b-31e8-4af9-9abb-6e676d13a4cf
n0f8(image[1,1].r)

# ╔═╡ f14e73f2-8d78-4d13-b7fd-3f7c2a729ebb
typeof(RGB(n0f8(0.5), 0, 0))

# ╔═╡ 9a793cd5-d7f2-496d-9a8e-2ec58716ec41
typeof(n0f8(0.5))

# ╔═╡ 948ceeda-6c71-41f7-a653-9ce961b586fd
typeof(RGB(n0f8(0.5), n0f8(0), n0f8(0)))

# ╔═╡ 8de562f7-b519-40a6-a46b-250fd0c365d3
function redify(color)
	RGB(color.r, 0, 0)
end

# ╔═╡ 2d24ddbe-58d1-4eb2-b176-08ac0842e083
R = redify.(image)

# ╔═╡ cf617485-1af4-4d0d-a9ef-1c3f986d413a
md"""
**(?)** Why the data type inside `RGB` is no more `N0f8` but `Float32`?
"""

# ╔═╡ 0fbffd3c-e548-4344-b9fd-013911219a27
typeof(redify.(image))

# ╔═╡ 2f6c7c75-0f53-42b1-8bcd-93fa489bf14b
G = (color -> RGB(0, color.g, 0)).(image)

# ╔═╡ 30050264-358d-46f6-b36e-ffc868034d06
typeof(G)

# ╔═╡ f9bbf516-c618-4ac8-a416-9716997d4358
B = (color -> RGB(0, 0, color.b)).(image)

# ╔═╡ c4d7b9a5-9aaf-4350-bd31-61fe1d999916
md"""
$(R+G+B) ``=`` $R ``\huge \bf +`` $G ``+`` $B
"""

# ╔═╡ 2b5f1501-430c-4ce3-bf89-827e906cc36e
begin
  step = 20
  md"""
$((R+G+B)[1:step:end, 1:step:end]) ``=`` $(R[1:step:end, 1:step:end]) ``+`` $(G[1:step:end, 1:step:end]) ``+`` $(B[1:step:end, 1:step:end])
  """
end

# ╔═╡ Cell order:
# ╠═624310f3-37b1-49aa-9089-ad860f4690a8
# ╠═6a5292a2-c3ad-11eb-125a-3f1bbc0b7a2c
# ╠═4c03a12c-4128-4da6-a6e1-ef91afa02614
# ╠═541f8b1f-2773-442f-8104-9b430afef1b0
# ╠═e2eddba9-3c11-43f6-bfd8-722c66bfbf9d
# ╠═afbeff8b-f0f1-4ed2-a09e-3848db5c47fb
# ╠═8e60dfbc-9a80-41e9-9091-67910940aeec
# ╠═0e5ddc64-3943-45d5-ac8d-e00e7e10c1cf
# ╠═0323a47c-3dcc-4162-b82b-d0504a7b029d
# ╠═5446ce8a-5143-4209-a5ea-bd2b5bd20e1a
# ╠═39595f13-7a31-401d-bc56-5a1106f567bc
# ╠═094184ad-325a-4b17-8e7a-89a29b5df473
# ╠═a0b12e5b-31e8-4af9-9abb-6e676d13a4cf
# ╠═f14e73f2-8d78-4d13-b7fd-3f7c2a729ebb
# ╠═9a793cd5-d7f2-496d-9a8e-2ec58716ec41
# ╠═948ceeda-6c71-41f7-a653-9ce961b586fd
# ╠═8de562f7-b519-40a6-a46b-250fd0c365d3
# ╠═2d24ddbe-58d1-4eb2-b176-08ac0842e083
# ╟─cf617485-1af4-4d0d-a9ef-1c3f986d413a
# ╠═0fbffd3c-e548-4344-b9fd-013911219a27
# ╠═2f6c7c75-0f53-42b1-8bcd-93fa489bf14b
# ╠═30050264-358d-46f6-b36e-ffc868034d06
# ╠═f9bbf516-c618-4ac8-a416-9716997d4358
# ╠═c4d7b9a5-9aaf-4350-bd31-61fe1d999916
# ╠═2b5f1501-430c-4ce3-bf89-827e906cc36e


text:
using Documenter
using MLStyle

makedocs(
    modules = [
        MLStyle.MatchCore,
        MLStyle.DataType,
        MLStyle.Err,
        MLStyle.Extension,
        MLStyle.Infras,
        MLStyle.Pervasives,

        # standard patterns
        MLStyle.TypeVarDecons,
        MLStyle.Active,
        MLStyle.Uncomprehensions,
        MLStyle.LambdaCases,
        MLStyle.WhenCases,

        MLStyle.Render
    ],
    clean = false,
    format = :html,
    sitename = "MLStyle.jl",
    linkcheck = !("skiplinks" in ARGS),
    analytics = "UA-89508993-1",
    pages = [
        "Home"   => "index.md",
        "Syntax" => Any[
            "syntax/adt.md",
            "syntax/pattern.md",
            "syntax/pattern-function.md",
            "syntax/when.md",
            "syntax/extension.md",
            "syntax/qualifier.md"
        ],
        "Tutorials" => Any[
            "tutorials/capture.md",
            "tutorials/query-lang.md",
        ],
        "Modules" => Any[
            "modules/ast.md",
            "modules/cond.md"
        ]
    ],
    html_prettyurls = !("local" in ARGS),
)

deploydocs(
    repo="github.com/thautwarm/MLStyle.jl",
    target="build",
    julia="1.0",
    deps=nothing,
    make=nothing)


text:
#=
Copyright Alec Hammond 2022

Simple fdtd in 2d excercise.

Assumes TMᶻ polarization. Relevant physics:
    −σₐHx − μ∂Hx/∂t = ∂Ez/∂y
     σₐHy + μ∂Hy/∂t = ∂Ez/∂x
     σEz + ε∂Ez/∂t = ∂Hy/∂x - ∂Hx/∂y
where
    σₐ is the magnetic conductivity
    σ is the electric conductivity
    ε is the permittivity
    μ is the permeability

Define Ez in lower left corner: Ez(m,n)
Define Hy just to the right:    Hy(m+1/2,n)
Define Hx just above:           Hx(m,n+1/2)

Nodes at top of grid lack Hx
Nodes at right of grid lack Hy

References
-------------------
* https://eecs.wsu.edu/~schneidj/ufdtd/chap8.pdf

=#

const USE_GPU = false
using ImplicitGlobalGrid
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
using Plots, Printf, Statistics

type = Float64
dims = 2
@static if USE_GPU
    @init_parallel_stencil(CUDA, type, dims);
else
    @init_parallel_stencil(Threads, type, dims);
end

@parallel function step_H!(Hx::Data.Array,Hy::Data.Array,Ez::Data.Array,dt::Data.Number,dx::Data.Number,dy::Data.Number)
    # assume σₐ=0
    # assume μ=1
    @all(Hx) = @all(Hx) - dt/dy*@d_ya(Ez)
    @all(Hy) = @all(Hy) + dt/dx*@d_xa(Ez)
    return
end

@parallel function step_E!(Hx::Data.Array,Hy::Data.Array,Ez::Data.Array,ε::Data.Array,dt::Data.Number, dx::Data.Number, dy::Data.Number)
    # assume σ=0
    @inn(Ez) = @inn(Ez) + (dt/dx./@inn(ε)).*(@d_xi(Hy)) - (dt/dy./@inn(ε)).*(@d_yi(Hx))
    return
end

@parallel_indices (ix,iy) function step_E_sources!(Ez::Data.Array,n::Int64,dt::Data.Number, dx::Data.Number, dy::Data.Number)
    # point source in center of cell
    if (ix==size(Ez,1)/2 && iy==size(Ez,2)/2)
        t = n*dt
        Ez[ix,iy] = cos(20*t)*exp((t-5)^2)
    end
    return
end

function fdtd2D_step!(Hx::Data.Array,Hy::Data.Array,Ez::Data.Array,ε::Data.Array,n::Int64,dt::Data.Number, dx::Data.Number, dy::Data.Number)
    @parallel step_H!(Hx,Hy,Ez,dt,dx,dy)
    @parallel step_E!(Hx,Hy,Ez,ε,dt,dx,dy)
    @parallel step_E_sources!(Ez,n,dt,dx,dy)
    return
end

function fdtd2D()

    # Physics
    lx, ly = 1.0, 1.0;

    # Numerics
    C      = 0.5;        # Courant factor
    nx, ny = 256, 256;
    nt     = 250;

    # Derived numerics
    dx, dy    = lx/(nx-1), ly/(ny-1) # cell sizes
    dt        = min(dx,dy)*C

    # Array initializations
    Hx       = @zeros(nx, ny-1);
    Hy       = @zeros(nx-1, ny);
    Ez       = @zeros(nx, ny);
    ε        = @ones(nx, ny);
    @printf("Starting...\n")
    # Time loop
    for it = 1:nt
        if (it==11)  global wtime0 = Base.time()  end
        fdtd2D_step!(Hx,Hy,Ez,ε,it,dt,dx,dy)
    end
    @printf("Stopping...\n")

    # Performance
    wtime    = Base.time() - wtime0
    A_eff    = (2*2+2)/1e9*nx*ny*sizeof(Data.Number)  # Effective main memory access per iteration [GB] (Lower bound of required memory access: H and dHdτ have to be read and written (dHdτ for damping): 4 whole-array memaccess; B has to be read: 1 whole-array memaccess)
    wtime_it = wtime/(nt-10)                          # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                         # Effective memory throughput [GB/s]

    @printf("Total steps=%d, time=%1.3e sec (@ T_eff = %1.2f GB/s) \n", nt, wtime, round(T_eff, sigdigits=2))
    gr()
    hm = heatmap(Ez, color = :bluesreds, aspect_ratio=:equal, show = true)
    png("Test")
    return
end

fdtd2D()

max_stars_repo_name: smartalecH/em-julia

max_issues_repo_name: smartalecH/em-julia

max_forks_repo_name: smartalecH/em-julia


text:
struct NaiveRigidBody <: AbstractRigidBody
    nn
    qdim
    udim
    thetamask
end
function NaiveRigidBody(qdim::Int, udim::Int, hidden_sizes, thetamask)

    return NaiveRigidBody(Chain(
        Dense(qdim*2+udim, hidden_sizes[1], tanh),
        (Dense(hidden_sizes[i], hidden_sizes[i+1], tanh) for i in 1:(length(hidden_sizes)-1))...,
        Dense(hidden_sizes[end], qdim, identity)
    ) |> f64, qdim, udim, thetamask)
end

function (f::NaiveRigidBody)(x, u)
    q = x[1:f.qdim, :]
    v = x[f.qdim+1:2*f.qdim, :]
    return vcat(v, f.nn(vcat(x, u)))
end
Flux.@functor NaiveRigidBody


text:
using FlexibleFunctors; const FF = FlexibleFunctors

# Standard struct with constant parameters
struct Foo{X,Y,PS<:Tuple}
    x::X
    y::Y
    ps::PS
end
FF.parameters(f::Foo) = f.ps

@testset "Parameters" begin
    @test FF.parameters([1, 2]) === ()
    A = [Ref(1), Ref(2)]
    @test FF.parameters(A) === A
    @test FF.parameters((1, 2)) == (1, 2)
    @test FF.parameters((a=1, b=2)) == (a=1, b=2)

    p0, re = ffunctor([1, 2])
    @test p0 === ()

    p0, re = ffunctor(A)
    @test p0 === A
end

@testset "FlexibleFunctors" begin
    # No parameters
    mf1 = Foo(1, 2, ())
    p0, re = ffunctor(mf1)
    @test re(()) == mf1

    # One parameters
    mf2 = Foo(1.0, 2, (:x, ))
    p0, re = ffunctor(mf2)
    @test re(p0) == mf2                         # Reconstruct self

    _, re2 = destructure(mf2)
    @test re2([1]) == Foo(1, 2, (:x,))          # Change type on Foo.x from Float64 to Int

    # Both fields as parameters
    mf3 = Foo(22.0, 3.14, (:x, :y))
    p0, re = destructure(mf3)
    @test re(p0) == mf3
    @test re(reverse(p0)) == Foo(3.14, 22.0, (:x, :y)) # Swap the fields

    # Nested parameters
    mf4 = Foo(Foo(1, 2, (:y,)), 3, (:x, :y))
    p0, re = destructure(mf4)
    @test re(p0) == mf4
    @test re(["hello", -1]) == Foo(Foo(1, "hello", (:y,)), -1, (:x, :y))
end

@testset "fieldmap" begin
    mf1 = Foo(3.0, 2, (:x,))
    @test fieldmap((x) -> x^2, mf1) == Foo(9.0, 2, (:x,))

    mf2 = Foo(-22.0, 3.14, (:x, :y))
    @test fieldmap(abs, mf2) == Foo(22.0, 3.14, (:x, :y))

    mf3 = Foo(Foo(1, 2, (:y,)), 3, (:x, :y))
    @test fieldmap(sqrt, mf3) == Foo(Foo(1, sqrt(2), (:y,)), sqrt(3), (:x, :y))
end

module MyMod
    import FlexibleFunctors
    struct Bar{A,B,PS<:Tuple}
        a::A
        b::B
        parameters::PS
    end
    FlexibleFunctors.parameters(b::Bar) = b.parameters
end
import .MyMod

@testset "Reconstruction of Module-Specific Type" begin
    mb = MyMod.Bar(1, 2, (:a,))
    p0, re = destructure(mb)
    @test re(p0) == mb
    @test re([3]) == MyMod.Bar(3, 2, (:a,))
end


text:
module ExprTools

export splitdef, combinedef

include("function.jl")

end


text:
g = basic_osm_graph_stub()

node_ids = keys(g.nodes)
n1_id, state = iterate(node_ids) # doesn't matter what one it is
n2_id, _ = iterate(node_ids, state) # doesn't matter what one it is
n1 = g.nodes[n1_id]
n2 = g.nodes[n2_id]

# Test GeoLocation methods, expect same node to be returned
idxs, dists = nearest_node(g, n1.location)
@test idxs[1][1] == n1.id
@test dists[1][1] == 0.0

# Test vector methods, expect same node to be returned
point1 = [n1.location.lat, n1.location.lon, n1.location.alt]
point2 = [n2.location.lat, n2.location.lon, n2.location.alt]
idxs, dists = nearest_node(g, point1)
@test idxs[1][1] == n1.id
@test dists[1][1] == 0.0
idxs, dists = nearest_node(g, [point1, point2])
@test all(x -> x == n1.id, idxs[1])
@test all(x -> x == n2.id, idxs[2])
@test all(x -> x == 0.0, dists[1])
@test all(x -> x == 0.0, dists[2])

# Test Node methods, expect different node to be returned
idxs, dists = nearest_node(g, n1)
@test idxs[1][1] !== n1.id
@test dists[1][1] !== 0.0
idxs, dists = nearest_node(g, [n1, n2], 2)
@test all(x -> x !== n1.id, idxs[1])
@test all(x -> x !== n2.id, idxs[2])
@test all(x -> x !== 0.0, dists[1])
@test all(x -> x !== 0.0, dists[2])
@test all(length.(idxs) .== 2) # Two points returned
@test all(x -> x[2] > x[1], dists) # 2nd point further away

# Test two nodes we know are closest in the stub graph
n1_id = 1005
idxs, dists = nearest_node(g, n1_id)
@test idxs == [[1004]]


text:
##############################
## Audio decomposition demo using NMF

using PyPlot, PyCall
using DataFrames
using Distributions

push!(LOAD_PATH, ".")
import NMF
@pyimport scipy.io.wavfile as wf

# load data
wavfile = "../data/organ.wav"
fs, data = wf.read(wavfile)

figure("data")
clf()
Pxx, freqs, t, pl = specgram(data[10000:318000,2], Fs=fs, NFFT=256, noverlap=0)
xlabel("time [sec]")
ylabel("frequency [Hz]")
ylim([0,22000])

# model
D, N = size(Pxx)
K = 2
a_t = 1.0
b_t = 1.0
a_v = 1.0
b_v = 100.0
prior = NMF.NMFModel(a_t*ones(D,K), b_t*ones(D, K), a_v, b_v)

# inference
max_iter = 100
posterior, S_est, T_est, V_est = NMF.VI(Int64.(round.(Pxx)), prior, max_iter)
X = T_est * V_est

# visualize
figure("T")
clf()
for k in 1 : K
    subplot(K,1,k)
    plot(T_est[:,k], linewidth=1.0)
    xlim([0, D])
    ylim([0, ylim()[2]])
end

figure("V")
clf()
for k in 1 : K
    subplot(K,1,k)
    plot(V_est[k,:], linewidth=1.0)
    xlim([0,N])
    ylim([0, ylim()[2]])
end
show()


text:
using Test

import ZMQ
import ZMQChannels

@testset "Releasing port when closing channels" begin
  @testset "ZMQPullChannel" begin
      try
        channel = ZMQChannels.ZMQPullChannel("tcp://*:6666")
        close(channel)

        sleep(0.1)

        channel = ZMQChannels.ZMQPullChannel("tcp://*:6666")
        close(channel)
      catch e
        @error "Should now throw any exceptions. Got: $e"
        @test false
      end
  end

  @testset "ZMQPushChannel" begin
    try
      channel = ZMQChannels.ZMQPushChannel("tcp://localhost:6666")
      close(channel)
      channel = ZMQChannels.ZMQPushChannel("tcp://localhost:6666")
      close(channel)
    catch e
      @error "Should now throw any exceptions. Got: $e"
      @test false
    end
  end
end


text:
"""
Helper infrastructure to compile and sample models using `cmdstan`.

[`StanModel`](@ref) wraps a model definition (source code), while [`stan_sample`](@ref) can
be used to sample from it.

[`stan_compile`](@ref) can be used to pre-compile a model without sampling. A
[`StanModelError`](@ref) is thrown if this fails, which contains the error messages from
`stanc`.
"""
module StanDiagnose

using Reexport

@reexport using StanBase

using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF

import StanRun: stan_sample, stan_cmd_and_paths, default_output_base
import StanBase: cmdline
import Base.show

include("stanmodel/diagnose_types.jl")
include("stanmodel/DiagnoseModel.jl")
include("stanrun/cmdline.jl")
include("stansamples/read_diagnose.jl")

stan_diagnose = stan_sample

export
  DiagnoseModel,
  stan_diagnose,
  read_diagnose

end # module


text:
@testset "Utility & Objective" begin
    D = 10
    N = 15
    lookback = 30

    financial_data = get_financial_data(D=D)

    lower_bound = 0.0
    upper_bound = 25.0
    wealth = 100.0
    max_risk = 500.0

    Δ_init = ones(D)
    Γ_init = 1.0

    μs = []
    sqrtΣs = []
    samples = []

    for n = 1:N
        d = Matrix(financial_data[n:lookback+n-1, 1:D])
        push!(μs, mean(Matrix(d), dims=1)[:])
        push!(sqrtΣs, Matrix(cholesky(cov(Matrix(d), dims=1)).U))
        push!(samples, Vector(financial_data[lookback+n, 1:D]))
    end

    POs = []
    for n = 1:N
        weights, constraints = DiffBasicWeightConvexSet(
            D,
            wealth,
            lower_bound,
            upper_bound,
            max_risk,
            μs[n],
            sqrtΣs[n]
        )
        push!(POs, DiffBertsimas(D, weights, constraints, μs[n], sqrtΣs[n]))
    end

    utility_function = utility(POs, samples, expected_return())

    function f(params::Vector{Float64})
        return utility_function.value(params)[1]
    end

    function g!(storage::Vector, params::Vector{Float64})
        storage[:] .= utility_function.grad(params)
    end

    result = Optim.maximize(
        f,
        g!,
        vcat(Δ_init, Γ_init),
        LBFGS(alphaguess = LineSearches.InitialStatic(scaled=true), linesearch = LineSearches.BackTracking()),
        Optim.Options(
            x_tol = 1.0e-12,
            f_tol = 1.0e-12,
            g_tol = 1.0e-2,
            f_calls_limit = 1000,
            g_calls_limit = 1000,
            iterations = 1000,
            show_trace = true,
        ),
    )

    @test norm(utility_function.grad(vcat(Δ_init, Γ_init))) > norm(utility_function.grad(result.res.minimizer))
end

max_stars_repo_name: molet/DiffDRO

max_issues_repo_name: molet/DiffDRO

max_forks_repo_name: molet/DiffDRO


text:
#
# camera.jl
#
# general camera status functions
# APIs
#

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

#==
    Setter functions
==#

"""
    Set exposure time
""" set_exposuretime
set_exposuretime(camera::Camera, exposuretime::Clong) = begin
        @checked_call(:FLISetExposureTIme,(CameraHandle, Clong), handle(camera), exposuretime)
end

"""
    Set_ROI
""" set_ROI
function set_ROI(camera::Camera, roi::Vector{Clong})
    if length(roi) != 4
        #overwrite the information
         roi = [1600,1200,0,0]
    end
	width = roi[1]
	height = roi[2]
	offsetx = roi[3]
	offsety = roi[4]

	max_width = 1600
	max_height = 1200

	upper_x = (max_width - width)/2 + offsetx;
	upper_y = (max_height - height)/2+ offsety;
	lower_x = upper_x + width + offsetx;
	lower_y = upper_y + height + offsety;

	lower_x <= max_width || error("X offset is too large")
	lower_y <= max_height || error("Y offset is too large")

	 #
     @checked_call(:FLISetImageArea,(CameraHandle, Clong, Clong, Clong, Clong), handle(camera),upper_x,upper_y, lower_x, lower_y)

     nothing
end

"""
    Set temperature
""" set_temperature
set_temperature(camera::Camera, temperature::Float64) = begin
    @checked_call(:FLISetTemperature,(CameraHandle, Float64), handle(camera), temperature)
end
#===
    Getter functions
===#

"""
    Get temperature
""" get_temperature
function get_temperature(camera::Camera)
	ref = Ref{Float64}(0)
	@checked_call(:FLIGetTemperature, (CameraHandle, Ptr{Float64}), handle(camera), ref)
	return  ref[]
end

"""
    Get timestamp to second using time() C library
    return: Float64
""" get_timestamp
@inline get_sys_timestamp() = time()

"""
    Flush CCD to clear the data before taking a new image

""" flushCCD

@inline flushCCD(camera::Camera) = @checked_call(:FLIFlushRow,(CameraHandle,Clong,Clong,),
                                    handle(camera), 1200 , 1)

"""
    Expose camera CCD to take image as long as specified by exposure time
"""
@inline exposeFrame(camera::Camera) = @checked_call(:FLIExposeFrame,(CameraHandle,),
                                                    handle(camera))

#===
Camera opearations
===#
"""
	initialize the camera by setting the temperature and wait until it cools down

"""
function initialize(camera::Camera; temperature::Float64 = -20.0)
	@info "Camera is initializing..."
	set_temperature(camera, temperature)
	temp = 30.0
	t0 = time()
	while temp > temperature
		temp = get_temperature(camera)
		t1 = time()
		if t1-t0 >= 5.0
			println("Camera temperature = ", temp, " ° C" )
			t0 = time()
		end
	end
	nothing

end

"""
    acquire and save images
"""
function acquire_n_save_image(camera::Camera, nSave::Int,imgConfig:: ImageConfigContext,
                                format::String)
    flushCCD(camera)
    width = imgConfig.width
    height = imgConfig.height
    widthInByte = imgConfig.width * imgConfig.bytePerPixel
    buff = Ref{UInt8}(0)
    img_cache = Matrix{UInt8}(udnef,height,width)
    for j in 1:nSave

        exposeFrame(camera)
        while !isexposing(camera) end

        for row in 1:height
            @checked_call(:FLIGrabRow, (CameraHandle, Ptr{UInt8}, Clong,),
                            handle(camera), buff, widthInByte)

            copyto!(img_cache[row,:], buff[])
        end

        flushCCD(camera)
        fname = @sprintf "FLI_image_%d" j
        fname = join([fname,format])
        save(fname, colorview(Gray,map(_convert,img_cache)))
        @printf "%s is saved\n" fname
    end

end

"""
    begin acquisition

""" beginAcquisition
function beginAcquisition(camera::Camera, sig::Channel{CameraStatus},
                            new_img::Condition, empty_buff::Condition)
    camera.status = Continue()
    flushCCD(camera)

    while true

        exposeFrame(camera)
        while !isexposing(camera) end
        notify(new_img)

        wait(empty_buff) # wait for buffer to finish grabing image

        flushCCD(camera)

        isready(sig) && (signal = take!(sig))
        if signal == Stop()
            camera.status = Stop()
            break
        end
    end
    nothing

end

"""
    Get the next image from camera buffer
    1.  img_cache provided by function caller to avoid allocating
    2.  imgConfig contains info about the image format
    3.  new_img conditional variable

""" getNextImage
function getNextImage(img_cache::Matrix{T}, timestamp::Vector{Float64},
                      imgConfig:: ImageConfigContext, new_img::Condition,
                      empty_buff::Condition) where T<:Number

    width = imgConfig.width
    height = imgConfig.height
    widthInByte = imgConfig.width * imgConfig.bytePerPixel
    wait(new_img)          # wait to grab image

    timestamp[1] = get_sys_timestamp()

    buff = Ref{UInt8}(0)
    for row in 1:height
        @checked_call(:FLIGrabRow, (CameraHandle, Ptr{UInt8}, Clong,),
                        handle(camera), buff, widthInByte)
        copyto!(img_cache[row,:], buff[])
    end
    notify(empty_buff)
    nothing

end

"""
end acquisition
""" endAcquisition
function endAcquisition(camera::Camera, sig::Channel{CameraStatus})
    camera.status = Continue() || error("Camera is not running")
    put!(sig, Stop());
    nothing
end

## Functions for remote workers

"""
    Work
    create image buffer, and start image acquisition.
""" work
function working(camNum::Int64)
      w = workers()
      try
         remote_do(FLICameras.work,w[1],camNum)

     catch e
         @info e
      end
     return w[1]
 end

"""
    Work function
"""
function work(camNum::Int64)
    Base.exit_on_sigint(false)
    interface = FLICameras.FLIDOMAIN_USB
    deviceType = FLICameras.FLIDEVICE_CAMERA

    camList = FLICameras.CameraList(interface, deviceType)
    camera = FLICameras.Camera(camList, camNum)

    # attach shared arrays
    img_array, imgTime_array = FLICameras.attach_remote_process()

    # operational variables
    new_img = Condition()
    empty_buff = Condition()
    sig = Channel{FLICameras.CameraStatus}(1)
    img_cache = Matrix{UInt8}(undef,imgConfig.width, imgConfig.height)
    timestamp_cache = [0.0];

    # start acquisition thread
    main_routine = @async beginAcquisition(camera, sig, new_img, empty_buff)

    # start grabbing image
    try
        while true

            getNextImage(img_cache, timestamp_cache, imgConfig, new_img, empty_buff)

            wrlock(img_array,1.0)do
                copyto!(img_array, img_cache)
            end
            wrlock(imgTime_array,1.0) do
                imgTime_array[1] = timestamp_cache
            end

        end

    catch e
        #
        FLICameras.endAcquisition(camera, sig)
        if e isa InterruptException
            @info "Acquisition loop is terminated"
            return nothing
        else
            rethrow(e)
            return nothing
        end

    end
end


text:
"""
	grdclip(cmd0::String="", arg1=[], kwargs...)

Clip the range of grid values. will set values < low to below and/or values > high to above.
You can also specify one or more intervals where all values should be set to ``between``,
or replace individual values.

Full option list at [`grdclip`](http://gmt.soest.hawaii.edu/doc/latest/grdclip.html)

Parameters
----------

- **cmd0** -- Str --

    Either the input file name or the full monolitic options string. Do not use this
    when the grid (a GMTgrid type) is passed via the ``arg1`` argument.
- **G** : **outgrid** : -- Str --

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdclip(....) form.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/grdclip.html#g)
- $(GMT.opt_R)
- **above** : **high** -- Array or Str --

    Two elements array with ``high`` and ``above`` or a string with "high/above".
    It sets all data[i] > ``high`` to ``above``.
- **below** : **low** -- Array or Str --

    Two elements array with ``low`` and ``below`` or a string with "low/below".
    It sets all data[i] < ``low`` to ``below``.
- **between** -- Array or Str --

    Three elements array with ``low, high`` and ``between`` or a string with "low/high/between".
    It sets all data[i] >= ``low`` and <= ``high`` to ``between``.
- **old** : **new** -- Array or Str --

    Two elements array with ``old`` and ``new`` or a string with "old/new".
    It sets all data[i] == ``old`` to ``new``.
- **S** -- Str --

    Condense all replacement options above in a single string.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/grdclip.html#t)
- $(GMT.opt_V)

-   Examples:

        G=gmt("grdmath", "-R0/10/0/10 -I1 X");
        G2=grdclip(G, above=[5 6], low=[2 2], between="3/4/3.5")

    or (note the use of -S for second on options because we can't repeat a kwarg name)

        G2=grdclip(G, S="a5/6 -Sb2/2 -Si3/4/3.5")
"""
function grdclip(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("grdclip", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd, = parse_R("", d)
	cmd = parse_V_params(cmd, d)

    cmd = add_opt(cmd, 'G', d, [:G :outgrid])
    cmd = add_opt(cmd, 'S', d, [:S])

	cmd = opt_S(d, cmd, [:high :above], 'a')
	cmd = opt_S(d, cmd, [:low :below], 'b')
	cmd = opt_S(d, cmd, [:old :new], 'r')
	cmd = opt_S(d, cmd, [:between], 'i')

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	if (isa(arg1, Array{<:Number}))		arg1 = mat2grid(arg1)	end
	return common_grd(d, cmd, got_fname, 1, "grdclip", arg1)	# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdclip(arg1=[], cmd0::String=""; kw...) = grdclip(cmd0, arg1; kw...)

# ---------------------------------------------------------------------------------------------------
function opt_S(d::Dict, cmd::String, symbs, flag::Char)
	# This is common to the 4 cases
	for sym in symbs
		if (haskey(d, sym))
			cmd = string(cmd, " -S", flag)
			if (isa(d[sym], String))
				cmd = cmd * d[sym]
			elseif (isa(d[sym], Array))
				if (sym == :between && length(d[sym]) == 3)
					cmd = @sprintf("%s%f/%f/%f", cmd, d[sym][1], d[sym][2], d[sym][3])
				elseif (length(d[sym]) == 2)
					cmd = @sprintf("%s%f/%f", cmd, d[sym][1], d[sym][2])
				else
					error("Wrong number of elements in S option")
				end
			else
				error(string("Argument of :", sym, " must be a string or a two elements array."))
			end
			break
		end
	end
	return cmd
end


text:
######################################################################
#
#     Atom Integration
#
#     - Redefines Media.render for type(s) passed to @redirect
#
######################################################################

import Media, Juno
import Media: render

# type RotoloDisplay <: Media.Display end
# const rd = RotoloDisplay()
# Media.@media RotoloThing
# Media.setdisplay(RotoloThing, rd)

@traitfn Media.render{X; ShouldRedirect{X}}(e::Juno.Editor, x::X) = begin
  Media.render(Juno.Inline(), Atom.icon("check"))
  showmsg(x)
end

@traitfn Media.render{X; !ShouldRedirect{X}}(e::Juno.Editor, x::X) = begin
  Media.render(Juno.Inline(), Juno.Copyable(x))
end

# redefinition of 'redirect'
# macro redirect(args...)
#   currentSession == nothing &&
#     error("A session needs to be created before calling @redirect")
#
#   for a in args
#     t = try
#           eval(current_module(), a)
#         catch e
#           error("can't evaluate $a, error $e")
#         end
#     isa(t, Type) || error("$a does not evaluate to a type")
#
#     push!(currentSession.redirected_types, t)
#     loadmsg(t)
#
#     @eval function Media.render(e::Atom.Editor, x::($t))
#         Media.render(e, nothing)
#         showmsg(x)
#       end
#
#   end
# end


text:
println("Dg) Test of the PhotoExcitationFluores module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
#
setDefaults("print summary: open", "zzz-PhotoExcitationFluores.sum")
wa = Atomic.Computation("xx",  Nuclear.Model(26.),
                        initialConfigs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                        intermediateConfigs=[Configuration("1s 2s^2 2p"), Configuration("1s 2s 2p^2") ],
                        finalConfigs  =[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p") ],
                        process = JAC.PhotoExcFluor,
                        processSettings=PhotoExcitationFluores.Settings()  )

wb = perform(wa)
setDefaults("print summary: close", "")

## [E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], false, false, false, true, true, Tuple{Int64,Int64,Int64}[(1,1,0), (1,2,0), (1,2,0)], 0., 0., 0., "Coulomb"


text:
include("tree.jl")

type KIFMM
    tree::Tree

    # surfaces variables
    nSurface
    stdUpEquivSurface
    stdUpCheckSurface
    stdDownEquivSurface
    stdDownCheckSurface
    stdNormal

    # kernel variables
    kernel

    # charge variables
    chargeTree

    KIFMM() = new(Tree(), 0, [], [], [], [], Nil, [])
end

function initialize!(k::KIFMM, nSurface, source, normal, target, charge, nSource, nTarget, rank, maxLevel, kernel::Function)
    k.tree = Tree()
    k.nSurface = nSurface
    # build tree
    t = k.tree
    populate!(t, source, normal, target, nSource, nTarget, rank, maxLevel)

    k.chargeTree = charge
    k.kernel = kernel

    k.stdUpEquivSurface = []
    k.stdUpCheckSurface = []
    k.stdDownEquivSurface = []
    k.stdDownCheckSurface = []
    k.stdNormal = []

    dr = 0.1

    getStdSurfaces!(√2 + dr, nSurface, k.stdUpEquivSurface)
    getStdSurfaces!(4.0 - √2 - 2 * dr, nSurface, k.stdUpCheckSurface)
    getStdSurfaces!(√2 + dr, nSurface, k.stdDownCheckSurface)
    getStdSurfaces!(4.0 - √2 - 2 * dr, nSurface, k.stdDownEquivSurface)
    getStdSurfaces!(1.0, nSurface, k.stdNormal)
end

function FMM!(k::KIFMM)
    upPass!(k, 1)
    potential = zeros(k.tree.nTarget)
    downPass!(k, 1, potential)
    return potential
end

function getStdSurfaces!(r::Float64, nSurface::Int, surface)
    θ = 2 * π/nSurface
    for l = 0:nSurface-1
        push!(surface, [r * cos(l*θ), r * sin(l*θ)])
    end
end

function getLocalSurface!(origin, center, radius, remote)
    nSurface = length(origin)
    for l = 1:nSurface
        push!(remote, [radius[1] * origin[l][1], radius[2] * origin[l][2]] + center)
    end
end

function upPass!(k::KIFMM, rootId)
    # not relevant to target, only related to source terms
    # S2M, M2M
    # bottom up, postorder traversal
    t = k.tree
    d = t.dict
    node = d[rootId]
    # at current node, initialize all surfaces

    # if there is no target points in this node then
    # this is no Down-* surfaces
    getLocalSurface!(k.stdDownEquivSurface, node.center, node.radius, node.shiftedDownEquivalentSurface)
    getLocalSurface!(k.stdDownCheckSurface, node.center, node.radius, node.shiftedDownCheckSurface)
    node.downwardEquivalent = zeros(k.nSurface)

    # if there is no source points in this node then
    # there is no Up-* surfaces
    getLocalSurface!(k.stdUpEquivSurface, node.center, node.radius, node.shiftedUpEquivalentSurface)
    getLocalSurface!(k.stdUpCheckSurface, node.center, node.radius, node.shiftedUpCheckSurface)
    node.upwardEquivalent = zeros(k.nSurface)

    if (node.isLeaf)
        # at leaf, S2M
        getCharge!(k, rootId)
        # upward transfer to equivalent surface (target)
        contribution = getAccumulation(k, node.source, node.normal, node.shiftedUpCheckSurface, node.charge)
        node.upwardEquivalent = getInversion(k, node.shiftedUpEquivalentSurface, k.stdNormal, node.shiftedUpCheckSurface, contribution)
    else
        # at non-leaf, M2M
        contribution = zeros(k.nSurface)
        for l = 1:4
            upPass!(k, node.child[l])
            if (!d[node.child[l]].isEmpty)
                contribution += getAccumulation(k, d[node.child[l]].shiftedUpEquivalentSurface, k.stdNormal, node.shiftedUpCheckSurface, d[node.child[l]].upwardEquivalent)
            end
        end
        node.upwardEquivalent = getInversion(k, node.shiftedUpEquivalentSurface, k.stdNormal, node.shiftedUpCheckSurface, contribution)
    end
end

function downPass!(k::KIFMM, rootId, potential)
    # preorder traversal, first root, then children
    # M2L, L2L, L2T
    t = k.tree
    d = t.dict
    node = d[rootId]

    if (node.parent != -1)
        # V List
        contribution = zeros(k.nSurface)
        for l = 1: node.nVList
            if (!d[node.vList[l]].isEmpty)
                contribution += getAccumulation(k, d[node.vList[l]].shiftedUpEquivalentSurface, k.stdNormal, node.shiftedDownCheckSurface, d[node.vList[l]].upwardEquivalent)
            end
        end

        node.downwardEquivalent += getInversion(k, node.shiftedDownEquivalentSurface, k.stdNormal, node.shiftedDownCheckSurface, contribution)

        # X list
        contribution = zeros(k.nSurface)
        for l = 1: node.nXList
            if (!d[node.xList[l]].isEmpty)
                contribution += getAccumulation(k, d[node.xList[l]].shiftedUpEquivalentSurface, k.stdNormal, node.shiftedDownCheckSurface, d[node.xList[l]].upwardEquivalent)
            end
        end

        node.downwardEquivalent += getInversion(k, node.shiftedDownEquivalentSurface, k.stdNormal, node.shiftedDownCheckSurface, contribution)

        # L2L
        # from parent
        parentNode = d[node.parent]
        contribution = getAccumulation(k, parentNode.shiftedDownEquivalentSurface, k.stdNormal, node.shiftedDownCheckSurface, parentNode.downwardEquivalent)
        node.downwardEquivalent += getInversion(k, node.shiftedDownEquivalentSurface, k.stdNormal, node.shiftedDownCheckSurface, contribution)
    end

    if (node.isLeaf && node.nTarget != 0)
        node.potential = zeros(node.nTarget)
        # get U , W, Parent
        for l = 1:node.nUList
            # direct interactions
            getCharge!(k, node.uList[l])
            if (!d[node.uList[l]].isEmpty)
                node.potential += getAccumulation(k, d[node.uList[l]].source, d[node.uList[l]].normal, node.target, d[node.uList[l]].charge)
            end
        end

        for l = 1:node.nWList
            if (!d[node.wList[l]].isEmpty)
                node.potential += getAccumulation(k, d[node.wList[l]].shiftedUpEquivalentSurface, k.stdNormal, node.target, d[node.wList[l]].upwardEquivalent)
            end
        end

        # from down-* surface
        node.potential += getAccumulation(k, node.shiftedDownEquivalentSurface, k.stdNormal, node.target, node.downwardEquivalent)

        for l = 1:node.nTarget
            potential[node.targetIndex[l]] += node.potential[l]
        end
    end

    if !node.isLeaf
        for l = 1:4
            downPass!(k, node.child[l], potential)
        end
    end

end

function getCharge!(k::KIFMM, id)
    # lazy evaluation
    node = k.tree.dict[id]
    if (node.chargeComputed == true)

    else
        node.chargeComputed = true
        node.charge = zeros(node.nSource)
        for l = 1:node.nSource
            node.charge[l] = k.chargeTree[node.sourceIndex[l]]
        end
    end
end

function kernelEval(source, normal, target, f)
    sourceSize = length(source)
    targetSize = length(target)
    K = zeros(targetSize, sourceSize)
    for s = 1:sourceSize
        for t = 1: targetSize
            K[t, s] = f(source[s], normal[s], target[t])
        end
    end
    return K
end

function getAccumulation(k::KIFMM, source, normal, target, charge)
    K = kernelEval(source, normal, target, k.kernel)
    # K(target, source) * charge
    # density of length(target)
    # charge of lenth(source)
    return K * charge
end

function getInversion(k::KIFMM, source, normal, target, density)
    K = kernelEval(source, normal, target, k.kernel)
    return pinv(K)*density
end


text:
incx = abs(incx)
dy = iy - ix
dz = iz - ix
@inbounds for i = ix:incx:ix+(n-1)*incx
    x[i] = $op( y[dy+i], z[dz+i] )
end


text:
# Canonical form of multivariate normal

### Generic types

"""
    MvNormalCanon

Multivariate normal distribution is an [exponential family distribution](http://en.wikipedia.org/wiki/Exponential_family),
with two *canonical parameters*: the *potential vector* ``\\mathbf{h}`` and the *precision matrix* ``\\mathbf{J}``.
The relation between these parameters and the conventional representation (*i.e.* the one using mean ``\\boldsymbol{\\mu}`` and
covariance ``\\boldsymbol{\\Sigma}``) is:

```math
\\mathbf{h} = \\boldsymbol{\\Sigma}^{-1} \\boldsymbol{\\mu}, \\quad \\text{ and } \\quad \\mathbf{J} = \\boldsymbol{\\Sigma}^{-1}
```

The canonical parameterization is widely used in Bayesian analysis. We provide a type `MvNormalCanon`,
which is also a subtype of `AbstractMvNormal` to represent a multivariate normal distribution using
canonical parameters. Particularly, `MvNormalCanon` is defined as:

```julia
struct MvNormalCanon{P<:AbstractPDMat,V<:Union{Vector,ZeroVector}} <: AbstractMvNormal
    μ::V    # the mean vector
    h::V    # potential vector, i.e. inv(Σ) * μ
    J::P    # precision matrix, i.e. inv(Σ)
end
```

We also define aliases for common specializations of this parametric type:

```julia
const FullNormalCanon = MvNormalCanon{PDMat,    Vector{Float64}}
const DiagNormalCanon = MvNormalCanon{PDiagMat, Vector{Float64}}
const IsoNormalCanon  = MvNormalCanon{ScalMat,  Vector{Float64}}

const ZeroMeanFullNormalCanon = MvNormalCanon{PDMat,    ZeroVector{Float64}}
const ZeroMeanDiagNormalCanon = MvNormalCanon{PDiagMat, ZeroVector{Float64}}
const ZeroMeanIsoNormalCanon  = MvNormalCanon{ScalMat,  ZeroVector{Float64}}
```

A multivariate distribution with canonical parameterization can be constructed using a common constructor `MvNormalCanon` as:

    MvNormalCanon(h, J)

Construct a multivariate normal distribution with potential vector `h` and precision matrix represented by `J`.

    MvNormalCanon(J)

Construct a multivariate normal distribution with zero mean (thus zero potential vector) and precision matrix represented by `J`.

    MvNormalCanon(d, J)

Construct a multivariate normal distribution of dimension `d`, with zero mean and a precision matrix as `J * eye(d)`.

# Arguments
- `d::Int`: dimension of distribution
- `h::Vector{T<:Real}`: the potential vector, of type `Vector{T}` with `T<:Real`.
- `J`: the representation of the precision matrix, which can be in either of the following forms (`T<:Real`):
    1. an instance of a subtype of `AbstractPDMat`
    2. a square matrix of type `Matrix{T}`
    3. a vector of type `Vector{T}`: indicating a diagonal precision matrix as `diagm(J)`.
    4. a real number: indicating an isotropic precision matrix as `J * eye(d)`.

**Note:** `MvNormalCanon` share the same set of methods as `MvNormal`.
"""
struct MvNormalCanon{T<:Real,P<:AbstractPDMat,V<:AbstractVector} <: AbstractMvNormal
    μ::V    # the mean vector
    h::V    # potential vector, i.e. inv(Σ) * μ
    J::P    # precision matrix, i.e. inv(Σ)
end

const FullNormalCanon = MvNormalCanon{Float64, PDMat{Float64,Matrix{Float64}},Vector{Float64}}
const DiagNormalCanon = MvNormalCanon{Float64,PDiagMat{Float64,Vector{Float64}},Vector{Float64}}
const IsoNormalCanon  = MvNormalCanon{Float64,ScalMat{Float64},Vector{Float64}}

const ZeroMeanFullNormalCanon = MvNormalCanon{Float64,PDMat{Float64,Matrix{Float64}},ZeroVector{Float64}}
const ZeroMeanDiagNormalCanon = MvNormalCanon{Float64,PDiagMat{Float64,Vector{Float64}},ZeroVector{Float64}}
const ZeroMeanIsoNormalCanon  = MvNormalCanon{Float64,ScalMat{Float64},ZeroVector{Float64}}

### Constructors

function MvNormalCanon(μ::AbstractVector{T}, h::AbstractVector{T}, J::AbstractPDMat{T}) where T<:Real
    length(μ) == length(h) == dim(J) || throw(DimensionMismatch("Inconsistent argument dimensions"))
    if typeof(μ) == typeof(h)
        return MvNormalCanon{T,typeof(J),typeof(μ)}(μ, h, J)
    else
        return MvNormalCanon{T,typeof(J),Vector{T}}(collect(μ), collect(h), J)
    end
end

function MvNormalCanon(μ::AbstractVector{T}, h::AbstractVector{T}, J::P) where {T<:Real, P<:AbstractPDMat}
    R = promote_type(T, eltype(J))
    MvNormalCanon(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, h), convert(AbstractArray{R}, J))
end

function MvNormalCanon(μ::AbstractVector{T}, h::AbstractVector{S}, J::P) where {T<:Real, S<:Real, P<:AbstractPDMat}
    R = Base.promote_eltype(μ, h, J)
    MvNormalCanon(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, h), convert(AbstractArray{R}, J))
end

function MvNormalCanon(J::P) where P<:AbstractPDMat
    z = ZeroVector(eltype(J), dim(J))
    MvNormalCanon{eltype(J),P,ZeroVector{eltype(J)}}(z, z, J)
end

function MvNormalCanon(h::AbstractVector{T}, J::P) where {T<:Real, P<:AbstractPDMat}
    length(h) == dim(J) || throw(DimensionMismatch("Inconsistent argument dimensions"))
    R = Base.promote_eltype(h, J)
    hh, JJ = collect(convert(AbstractArray{R}, h)), convert(AbstractArray{R}, J)
    MvNormalCanon{eltype(hh),typeof(JJ),typeof(hh)}(JJ \ hh, hh, JJ)
end

MvNormalCanon(h::AbstractVector{<:Real}, J::AbstractMatrix{<:Real}) = MvNormalCanon(h, PDMat(J))
MvNormalCanon(h::AbstractVector{<:Real}, prec::AbstractVector{<:Real}) = MvNormalCanon(h, PDiagMat(prec))
MvNormalCanon(h::AbstractVector{<:Real}, prec::Real) = MvNormalCanon(h, ScalMat(length(h), prec))

MvNormalCanon(J::AbstractMatrix) = MvNormalCanon(PDMat(J))
MvNormalCanon(prec::AbstractVector) = MvNormalCanon(PDiagMat(prec))
MvNormalCanon(d::Int, prec) = MvNormalCanon(ScalMat(d, prec))

### Show

distrname(d::IsoNormalCanon) = "IsoNormalCanon"
distrname(d::DiagNormalCanon) = "DiagNormalCanon"
distrname(d::FullNormalCanon) = "FullNormalCanon"

distrname(d::ZeroMeanIsoNormalCanon) = "ZeroMeanIsoNormalCanon"
distrname(d::ZeroMeanDiagNormalCanon) = "ZeroMeanDiagormalCanon"
distrname(d::ZeroMeanFullNormalCanon) = "ZeroMeanFullNormalCanon"

### Conversion
function convert(::Type{MvNormalCanon{T}}, d::MvNormalCanon) where {T<:Real}
    MvNormalCanon(convert(AbstractArray{T}, d.μ), convert(AbstractArray{T}, d.h), convert(AbstractArray{T}, d.J))
end
function convert(::Type{MvNormalCanon{T}}, μ::AbstractVector{<:Real}, h::AbstractVector{<:Real}, J::AbstractPDMat) where {T<:Real}
    MvNormalCanon(convert(AbstractArray{T}, μ), convert(AbstractArray{T}, h), convert(AbstractArray{T}, J))
end

### conversion between conventional form and canonical form

meanform(d::MvNormalCanon) = MvNormal(d.μ, inv(d.J))
# meanform{C, T<:Real}(d::MvNormalCanon{T,C,Vector{T}}) = MvNormal(d.μ, inv(d.J))
# meanform{C, T<:Real}(d::MvNormalCanon{T,C,ZeroVector{T}}) = MvNormal(inv(d.J))

function canonform(d::MvNormal{T,C,<:AbstractVector{T}}) where {C, T<:Real}
    J = inv(d.Σ)
    return MvNormalCanon(d.μ, J * collect(d.μ), J)
end
canonform(d::MvNormal{T,C,ZeroVector{T}}) where {C, T<:Real} = MvNormalCanon(inv(d.Σ))

### Basic statistics

length(d::MvNormalCanon) = length(d.μ)
mean(d::MvNormalCanon) = convert(Vector{eltype(d.μ)}, d.μ)
params(d::MvNormalCanon) = (d.μ, d.h, d.J)
@inline partype(d::MvNormalCanon{T}) where {T<:Real} = T
eltype(::MvNormalCanon{T}) where {T} = T

var(d::MvNormalCanon) = diag(inv(d.J))
cov(d::MvNormalCanon) = Matrix(inv(d.J))
invcov(d::MvNormalCanon) = Matrix(d.J)
logdetcov(d::MvNormalCanon) = -logdet(d.J)

### Evaluation

sqmahal(d::MvNormalCanon, x::AbstractVector) = quad(d.J, broadcast(-, x, d.μ))
sqmahal!(r::AbstractVector, d::MvNormalCanon, x::AbstractMatrix) = quad!(r, d.J, broadcast(-, x, d.μ))

# Sampling (for GenericMvNormal)

unwhiten_winv!(J::AbstractPDMat, x::AbstractVecOrMat) = unwhiten!(inv(J), x)
unwhiten_winv!(J::PDiagMat, x::AbstractVecOrMat) = whiten!(J, x)
unwhiten_winv!(J::ScalMat, x::AbstractVecOrMat) = whiten!(J, x)

_rand!(rng::AbstractRNG, d::MvNormalCanon, x::AbstractVector) =
    add!(unwhiten_winv!(d.J, randn!(rng,x)), d.μ)
_rand!(rng::AbstractRNG, d::MvNormalCanon, x::AbstractMatrix) =
    add!(unwhiten_winv!(d.J, randn!(rng,x)), d.μ)


text:
"""
Attitude objects contains Euler angles (rad) and quaternions.

# Constructors
    Attitude(psi, theta, phi)
    Attitude(q0, q1, q2, q3)
"""
struct Attitude
    # Euler angles
    psi::Number
    theta::Number
    phi::Number
    # Quaternions
    q0::Number
    q1::Number
    q2::Number
    q3::Number
end

# Constructor from Euler angles
Attitude(psi, theta, phi) = Attitude(psi, theta, phi,
                                     euler2quaternion(psi, theta, phi)...)

# Constructor from quaternions
Attitude(q0, q1, q2, q3) = Attitude(quaternion2euler(q0, q1, q2, q3)...,
                                    q0, q1, q2, q3)

# Getters
"""
    get_euler_angles(att::Attitude)

Return Euler angles ψ, θ, ϕ (yaw, pith, roll).
"""
get_euler_angles(att::Attitude) = [att.psi, att.theta, att.phi]

"""
    get_quaternions(att::Attitude)

Return quaternions.
"""
get_quaternions(att::Attitude) = [att.q0, att.q1, att.q2, att.q3]


text:
function deal!{T<:Integer}(d::Deck, hlen::T)
    if hlen < length(d)
        hand = Deck(d.cards[1:hlen], d.design)
        d.cards = d.cards[hlen+1:end]
    else
        hand = d
        empty!(d)
    end
    return hand
end

function deal!(d::Deck)
    deal!(d, d.design.hlen)
end

function pretty(d::Deck)
    s = ""
    llen = d.design.rlen
    dlen = length(d)
    for i in 1:llen:dlen
        j = min(i+llen-1, dlen)
        s *= print(d[i:j])*"\n"
    end
    chop(s)
end


text:
using Test
using Seis, SeisSplit
import Random, Dates

# Seed the random number generator 'randomly'
Random.seed!(Dates.second(Dates.now()))

read_test_data() = Seis.read_sac.(joinpath(@__DIR__, "data", "wave.BH").*("E", "N"))

"Test that all fields are approximately equal for two `SeisSplit.Result`s"
function splits_are_equal(s, s′, atol=0.01)
    for f in fieldnames(SeisSplit.Result)
        f in (:trace1, :trace2, :xcorr_map) && continue # Don't compare traces themselves
        @test getfield(s, f) ≈ getfield(s′, f) atol=atol
    end
end

@testset "All tests" begin
    include("minimum_eigenvalue.jl")
    include("rotation_correlation.jl")
    include("qc.jl")
end

max_stars_repo_name: anowacki/SeisSplit

max_issues_repo_name: anowacki/SeisSplit

max_forks_repo_name: anowacki/SeisSplit


text:
using ContGridWorld
using Test

@testset "ContGridWorld.jl" begin
    # Write your tests here.
end


text:
module ExampleEntity

using ..Ahorn, Maple

#=
Minimum valid declaration.

This defines a type `Example` with a constructor that takes in `x` and `y`,
which then creates a unique `Entity` with the name `"ExampleMod/ExampleEntity"`
and a Dict containing the keys x and y, set to the values passed to the `Example` constructor.
=#
@mapdef Entity "ExampleMod/ExampleEntity" Example(x::Number, y::Number)
#= equivalent to:
    @pardef Example(x::Number, y::Number) = Entity("ExampleMod/ExampleEntity", x=x, y=y)
=#
#= evaluates to:
    const Example = Entity{Symbol("ExampleMod/ExampleEntity")}
    Example(x::Number, y::Number) = Entity{Symbol("ExampleMod/ExampleEntity")}("ExampleMod/ExampleEntity", x = x, y = y)
=#

#=
Entity data has the following reserved attributes: x, y, width, height, nodes.
All attributes other than x and y must be given a default value.
Allowed attribute types (ones that won't break the map parser) include:
    String, Bool, Char, Int, Float
=#
@mapdef Entity "ExampleMod/ExampleEntity2" Example2(x::Number, y::Number, width::Number=16, height::Number=16, nodes::Vector{Tuple{Int, Int}}=[(0,0)], attrib::Bool=false)
#= equivalent to:
    @pardef Example2(x::Number, y::Number, width::Number=16, height::Number=16, nodes::Vector{Tuple{Int, Int}}=[(0,0)], attrib::Bool=false) =
        Entity("ExampleMod/ExampleEntity", x=x, y=y, width=width, height=height, nodes=nodes, attrib=attrib)
=#
#= evaluates to:
    const Example2 = Entity{Symbol("ExampleMod/ExampleEntity2")}
    Example2(x::Number, y::Number, width::Number = 16, height::Number = 16, nodes::Vector{Tuple{Int, Int}} = [(0, 0)], attrib::Bool = false) =
        Entity{Symbol("ExampleMod/ExampleEntity2")}("ExampleMod/ExampleEntity2", x = x, y = y, width = width, height = height, nodes = nodes, attrib = attrib)
=#

# On load, EntityPlacements are retrieved from the placements field of each module and combined into a list.
const placements = Ahorn.PlacementDict(
    "Example Entity" => Ahorn.EntityPlacement(
        Example,                    # Your Entity type
        "point",                    # One of: "point", "rectangle", "line"
        Dict{String, Any}(),        # Data specific to this placement
        function finalizer(args...) # Called when the entity is created
            # arguments can be any of:
            #    (target, map, room)
            #    (target, room)
            #    (target)
        end
    )
)

#= Render the entity in Ahorn.
render is relative to `Ahorn.position(entity)`, renderAbs is relative to 0,0

Drawing helper methods provided by Ahorn can be found here:
https://github.com/CelestialCartographers/Ahorn/blob/12f6ac29677dc1f3c05b468bd205b757cbd4ec1e/src/drawing.jl#L223
=#
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Example) = false
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Example, room::Room) = false
Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Example) = false
Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Example, room::Room) = false

# Render the entity while it is selected.
Ahorn.renderSelected(ctx::Ahorn.Cairo.CairoContext, entity::Example) = false
Ahorn.renderSelected(ctx::Ahorn.Cairo.CairoContext, entity::Example, room::Maple.Room) = false
Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Example) = false
Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Example, room::Maple.Room) = false

#=
Describe the selection box for your entity.
Supported return types:
    Ahorn.Rectangle, Ahorn.Rectangle[], Bool (does nothing)
=#
Ahorn.selection(entity::Example) = nothing
Ahorn.selection(entity::Example, room::Maple.Room) = nothing

# The minimum width and height allowed for the entity. Can be overridden with debug settings.
Ahorn.minimumSize(entity::Example) = 8, 8
# Whether to allow resizing the x and y axis. Can be overridden with debug settings.
Ahorn.resizable(entity::Example) = false, false

# The minimum and maximum number of allowed nodes. `-1` for unbounded.
Ahorn.nodeLimits(entity::Example) = 0, 0

#=
Add specific editing options for attributes.
Supported options include:
    Array, Dict
=#
Ahorn.editingOptions(entity::Example) = Dict{String, Any}(
    "attributeName" => ["opt1", "opt2", "opt3"]
)

# The order in which attributes will be displayed when editing.
# Attributes that aren't listed are appended in alphabetical order.
Ahorn.editingOrder(entity::Example) = String["x", "y", "width", "height"]

# Attributes that should not be displayed when editing.
# `multiple` is true when more than one entity is being edited simultaneously.
Ahorn.editingIgnored(entity::Example, multiple::Bool=false) = multiple ? String["x", "y", "width", "height", "nodes"] : String[]

# `deleted`, `moved`, and `resized` are callbacks provided by Ahorn for the corresponding actions.
Ahorn.deleted(entity::Example, node::Int) = nothing
Ahorn.moved(entity::Example) = nothing
Ahorn.moved(entity::Example, x::Int, y::Int) = nothing
Ahorn.resized(entity::Example) = nothing
Ahorn.resized(entity::Example, width::Int, height::Int) = nothing

# If an entity is returned, it replaces the current one
Ahorn.flipped(entity::Example, horizontal::Bool) = nothing

# If an entity is returned, it replaces the current one
Ahorn.rotated(entity::Example, steps::Int) = nothing

end


text:
include("imports.jl")

@testset ExtendedTestSet "Selectors" begin
    node = XNode(:body, Dict(:class => "content"), [XLeaf(1), XLeaf(2)])
    leaf = children(node)[1]
    @test matches(SelectTag(:body), node)
    @test !matches(SelectTag(:div), node)
    @test matches(!SelectTag(:div), node)
    @test matches(SelectTag(:div) | SelectTag(:body), node)
    @test matches(SelectLeaf(), leaf)
    @test matches(SelectNode(), node)
    @test matches(SelectHasAttr(:class), node)
    @test matches(SelectAttrEq(:class, "content"), node)
end


text:
@testset "XOR reversibility UInt" begin
    for T in (UInt8,UInt16,UInt32,UInt64)
        A = rand(T,1000)
        @test A == unxor_delta(xor_delta(A))
        @test A == xor_delta(unxor_delta(A))

        B = copy(A)
        xor_delta!(A)
        unxor_delta!(A)
        @test B == A

        unxor_delta!(A)
        xor_delta!(A)
        @test B == A
    end
end

@testset "XOR reversibility Float" begin
    for T in (Float32,Float64)
        A = rand(T,1000)
        @test A == unxor_delta(xor_delta(A))
        @test A == xor_delta(unxor_delta(A))
    end
end

@testset "UInt: Backtranspose of transpose" begin
    for T in (UInt8,UInt16,UInt32,UInt64)
        for s  in (10,999,2048,9999,123123)
            A = rand(T,s)
            @test A == bitbacktranspose(bittranspose(A))
        end
    end
end

@testset "Float: Backtranspose of transpose" begin
    for T in (Float32,Float64)
        for s in (10,999,2048,9999,123123)
            A = rand(T,s)
            @test A == bitbacktranspose(bittranspose(A))
        end
    end
end

@testset "N-dimensional arrays" begin
    A = rand(UInt32,123,234)
    @test A == bitbacktranspose(bittranspose(A))

    A = rand(UInt32,12,23,34)
    @test A == bitbacktranspose(bittranspose(A))

    A = rand(Float32,123,234)
    @test A == bitbacktranspose(bittranspose(A))

    A = rand(Float32,12,23,34)
    @test A == bitbacktranspose(bittranspose(A))
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct LogNormal <: AbstractDistDelayed
    mu::NestedFloat
    sigma::NestedFloat
end

lognormal(mu::Float64, sigma::Float64)::Float64 = exp(normal(mu, sigma))


text:
# Dantzig-Wolfe Reformulation and Column Generation
# Edward J. Xu
# 2019.5.28
########################################################################################################################

function SSCFLRand(m::Int64,n::Int64, maxXY::Int64, minF::Int64, maxF::Int64, minD::Int64, maxD::Int64, seed)
    # m:  number of facilities
    # n: number of customers
    Random.seed!(seed)
    # Cost of opening facitlities
    f = rand(minF:maxF, m)
    # Coordinates for facilities (chosen in the box [0,maxXY]x[0,maxXY])
    fac = rand(0:maxXY, m, 2)
    # Coordinates for customers (chosen in the box [0,maxXY]x[0,maxXY])
    cust = rand(0:maxXY, n, 2)
    d = rand(minD:maxD, n)
    totDemand = sum(d)
    capPerFacility = ceil(1.3 * totDemand / m)
    q = capPerFacility*ones(m)

    c = zeros(m, n)
    # compute distances
    for i = 1:m
        for j = 1:n
            dist = round(norm(fac[i,1:2] - cust[j,1:2]))  # norm function is in package LinearAlgebra
            c[i,j] = dist
        end
    end
    return (c, f, q, d)
end

function solveSSCFL(c, f, q, d)
    (m,n) = size(c)
    timeStart = time()
    mod = Model(solver = GurobiSolver())
    #mod = Model(solver = CplexSolver())
    @variable(mod, x[1:m, 1:n], Bin )
    @variable(mod, y[1:m], Bin )
    @objective(mod, Min, sum(c[i,j] * x[i,j] for i=1:m for j=1:n) + sum(f[i]*y[i] for i=1:m))
    @constraint(mod, [j = 1:n], sum(x[i,j] for i=1:m) == 1)
    #@constraint(mod, [i=1:m,j=1:n], x[i,j] <= y[i])
    @constraint(mod, [i = 1:m], sum(d[j]*x[i,j] for j=1:n) <= q[i]y[i])
    # writeLP(mod, "SSCFLP.lp", genericnames = false)
    status = solve(mod, relaxation = true)
    println("#### Objective value, relaxation: ", getobjectivevalue(mod))
    println("#### y = ", getvalue(y))
    println("#### Elapsed time is $(time() - timeStart) seconds.\n",
            "################################################################################") ########################
    status = solve(mod)
    println("#### Objective value: ", getobjectivevalue(mod))
    #println("x = ", getvalue(x))
    #println("y = ", getvalue(y))
    println("#### Number of facilities opened: ", sum(getvalue(y)))
    println("#### Elapsed time is $(time() - timeStart) seconds.\n",
            "#### End #######################################################################") ########################
end

function getData(whi::Float64)
    if whi == 1  # Convexify􏰁 constraints into one sub-problem
        vec_c = - [292  453  359  219  268  736  291  443  403  498  400  967]
        mat_a = [
            Matrix{Float64}(I, 5, 5) zeros(5, 1) Matrix{Float64}(I, 5, 5) zeros(5, 1);
            14 20 6 16 10 -43 zeros(1, 6);
            zeros(1, 6) 14 20 6 16 10 -43
            ]
        vec_b = hcat([ones(5,1); 0; 0])
        matIndexSub = hcat(collect(1:12))'
        numSub = 1
        numXInSub = 12
        m_mat_h = 5
        vecRowMatD = [2]
    elseif whi == 2  # Convexify􏰁 constraints into two sub-problems
        vec_c = - [292  453  359  219  268  736  291  443  403  498  400  967]
        mat_a = [
            Matrix{Float64}(I, 5, 5) zeros(5, 1) Matrix{Float64}(I, 5, 5) zeros(5, 1);
            14 20 6 16 10 -43 zeros(1, 6);
            zeros(1, 6) 14 20 6 16 10 -43
            ]
        vec_b = hcat([ones(5,1); 0; 0])
        matIndexSub = Array{Int64}(undef, 2, 6)
        matIndexSub[1, :] = collect(1:6)
        matIndexSub[2, :] = collect(7:12)
        numSub = 2
        numXInSub = 6
        m_mat_h = 5
        vecRowMatD = [1, 1]
    elseif whi == 3  # Convexify􏰁 constraints into one sub-problems, restrain the y
        vec_c = - [292  453  359  219  268  736  291  443  403  498  400  967]
        mat_a = [
            Matrix{Float64}(I, 6, 6) Matrix{Float64}(I, 6, 6);
            14 20 6 16 10 -43 zeros(1, 6);
            zeros(1, 6) 14 20 6 16 10 -43
            ]
        vec_b = hcat([ones(5,1); 2; 0; 0])
        matIndexSub = Array{Int64}(undef, 2, 6)
        matIndexSub[1, :] = collect(1:6)
        matIndexSub[2, :] = collect(7:12)
        numSub = 1
        numXInSub = 12
        m_mat_h = 6
        vecRowMatD = [2]
    elseif whi == 4  # Convexify􏰁 constraints into two sub-problems, restrain the y
        vec_c = - [292  453  359  219  268  736  291  443  403  498  400  967]
        mat_a = [
            Matrix{Float64}(I, 6, 6) Matrix{Float64}(I, 6, 6);
            14 20 6 16 10 -43 zeros(1, 6);
            zeros(1, 6) 14 20 6 16 10 -43
            ]
        vec_b = hcat([ones(5,1); 2; 0; 0])
        matIndexSub = Array{Int64}(undef, 2, 6)
        matIndexSub[1, :] = collect(1:6)
        matIndexSub[2, :] = collect(7:12)
        numSub = 2
        numXInSub = 6
        m_mat_h = 6
        vecRowMatD = [1, 1]
    elseif whi >= 5
        if round(whi) == 5
            (c, f, q, d) = SSCFLRand(2, 5, 500, 500, 1000, 1, 20, 3)
        elseif round(whi) == 6
            (c, f, q, d) = SSCFLRand(10, 25, 500, 500, 1000, 1, 20, 3)
        elseif round(whi) == 7
            (c, f, q, d) = SSCFLRand(20, 50, 500, 500, 1000, 1, 20, 3)
        end
        vec_c = - reshape(hcat(c, f)', 1, :)
        (m, n) = size(c)
        mat_a_atom = zeros(m, m + n, n + 1)
        for i = 1:m
            mat_a_atom[i, :, :] = [
                Matrix{Float64}(I, n, n) zeros(n, 1);
                zeros(m, n + 1)
                ]
            mat_a_atom[i, n + i, :] = vcat(d, - q[i])
        end
        mat_a = zeros(n + m, m * (n + 1))
        for i = 1:m
            mat_a[:, ((i - 1) * (n + 1) + 1):(i * (n + 1))] = mat_a_atom[i, :, :]
        end
        vec_b = hcat(zeros(n + m, 1))
        vec_b[1:n, 1] = ones(n, 1)
        if ((whi == 5.1) | (whi == 6.1) | (whi == 7.1))
            matIndexSub = hcat(collect(1:(m * (n + 1))))'
            numSub = 1
            numXInSub = m * (n + 1)
            m_mat_h = n
            vecRowMatD = Array{Int64}(undef, 1)
            vecRowMatD = [m]
        elseif ((whi == 5.2) | (whi == 6.2) | (whi == 7.2))
            matIndexSub = Array{Int64}(undef, m, n + 1)
            for i = 1:m
                matIndexSub[i, :] = collect(((i - 1) * (n + 1) + 1):(i * (n + 1)))
            end
            numSub = m
            numXInSub = n + 1
            m_mat_h = n
            vecRowMatD = Array{Int64}(undef, 1)
            vecRowMatD = ones(1, m)
            println("#### 0,  Solve SSCFL Model using Gurobi#########################################") ####################
            solveSSCFL(c, f, q, d)
            println("#### End #######################################################################") ####################
        else
            throw("Wrong \"whi\" input!")
        end
    else
        throw("Wrong \"whi\" input!")
    end
    return (vec_c, mat_a, vec_b, matIndexSub, numSub, numXInSub, m_mat_h, vecRowMatD)
end


text:
using LinearAlgebra, SparseArrays, Test

using MultivariatePolynomials
using SumOfSquares

# Taken from JuMP/test/solvers.jl
function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

if try_import(:DynamicPolynomials)
    import DynamicPolynomials.@polyvar
else
    if try_import(:TypedPolynomials)
        import TypedPolynomials.@polyvar
    else
        error("No polynomial implementation installed : Please install TypedPolynomials or DynamicPolynomials")
    end
end

include("gram_matrix.jl")

include("variable.jl")
include("constraint.jl")

include("mock_tests.jl")

include("solvers.jl")

include("certificate.jl")

include("motzkin.jl")

# SOSTools demos
include("sospoly.jl")
include("lyapunov.jl")
include("sosdemo3.jl")
include("sosdemo5.jl")
include("sosdemo6.jl")
include("domain.jl")
include("sosmatrix.jl")
include("equalitypolyconstr.jl")
include("dsos.jl")
include("extract.jl")


text:
using RRTMGP.OpticalProps
using RRTMGP.CloudOptics

"""
    load_cld_lutcoeff(cloud_spec, ds_cld_coeff)

read cloud optical property LUT coefficients from NetCDF file
"""
function load_cld_lutcoeff(::Type{FT}, ds_cld_coeff, icergh) where {FT}
    # Open cloud optical property coefficient file
    # Read LUT coefficient dimensions
    band_lims_wvn = Array{FT}(ds_cld_coeff["bnd_limits_wavenumber"][:])

    # Read LUT constants
    liq = LookUpTable(
        ds_cld_coeff["radliq_lwr"][:],
        ds_cld_coeff["radliq_upr"][:],
        ds_cld_coeff["radliq_fac"][:],
        Array{FT}(ds_cld_coeff["lut_extliq"][:]),
        Array{FT}(ds_cld_coeff["lut_ssaliq"][:]),
        Array{FT}(ds_cld_coeff["lut_asyliq"][:]),
    )

    ice = LookUpTable(
        ds_cld_coeff["radice_lwr"][:],
        ds_cld_coeff["radice_upr"][:],
        ds_cld_coeff["radice_fac"][:],
        Array{FT}(ds_cld_coeff["lut_extice"][:]),
        Array{FT}(ds_cld_coeff["lut_ssaice"][:]),
        Array{FT}(ds_cld_coeff["lut_asyice"][:]),
    )

    cloud_spec = CloudOpticsLUT{FT, Int}(OpticalPropsBase("RRTMGP cloud optics LUT", band_lims_wvn), icergh, liq, ice)
    return cloud_spec
end

"""
    load_cld_padecoeff

read cloud optical property Pade coefficients from NetCDF file
"""
function load_cld_padecoeff(::Type{FT}, ds_cld_coeff, icergh) where {FT}
    # # Spectral discretization
    # Read Pade coefficient dimensions
    band_lims_wvn = Array{FT}(ds_cld_coeff["bnd_limits_wavenumber"][:])

    # Allocate cloud property Pade coefficient input arrays

    liq = PadeMethod(
        Array{FT}(ds_cld_coeff["pade_extliq"][:]),
        Array{FT}(ds_cld_coeff["pade_ssaliq"][:]),
        Array{FT}(ds_cld_coeff["pade_asyliq"][:]),
        Array{FT}(ds_cld_coeff["pade_sizreg_extliq"][:]),
        Array{FT}(ds_cld_coeff["pade_sizreg_ssaliq"][:]),
        Array{FT}(ds_cld_coeff["pade_sizreg_asyliq"][:]),
    )

    ice = PadeMethod(
        Array{FT}(ds_cld_coeff["pade_extice"][:]),
        Array{FT}(ds_cld_coeff["pade_ssaice"][:]),
        Array{FT}(ds_cld_coeff["pade_asyice"][:]),
        Array{FT}(ds_cld_coeff["pade_sizreg_extice"][:]),
        Array{FT}(ds_cld_coeff["pade_sizreg_ssaice"][:]),
        Array{FT}(ds_cld_coeff["pade_sizreg_asyice"][:]),
    )

    cloud_spec = CloudOpticsPade{FT, Int}(OpticalPropsBase("RRTMGP cloud optics Pade", band_lims_wvn), icergh, liq, ice)
    return cloud_spec
end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
"""
escapeHTML(i::String)

Returns a string with special HTML characters escaped: &, <, >, ", '
"""
function escapeHTML(i::String)
    # Refer to http://stackoverflow.com/a/7382028/3822752 for spec. links
    o = replace(i, "&", "&")
    o = replace(o, "\"", """)
    o = replace(o, "'", "'")
    o = replace(o, "<", "<")
    o = replace(o, ">", ">")
    return o
end

macro retry(expr)
    :(@retry 2 $(esc(expr)))
end

macro retry(N, expr)
    :(@retryif Any $N $(esc(expr)))
end

macro retryif(cond, expr)
    :(@retryif $(esc(cond)) 2 $(esc(expr)))
end

macro retryif(cond, N, expr)
    quote
        local __r__
        for i = 1:$N
            try
                __r__ = $(esc(expr))
                break
            catch e
                typeof(e) <: $(esc(cond)) || rethrow(e)
                i == $N && rethrow(e)
                sleep(0.1)
            end
        end
        __r__
    end
end

"""
@timeout secs expr then pollint

Start executing `expr`; if it doesn't finish executing in `secs` seconds,
then execute `then`. `pollint` controls the amount of time to wait in between
checking if `expr` has finished executing (short for polling interval).
"""
macro timeout(t, expr, then, pollint=0.01)
    return quote
        tm = Float64($(esc(t)))
        start = time()
        tsk = @async $(esc(expr))
        while !istaskdone(tsk) && (time() - start < tm)
            sleep($pollint)
        end
        istaskdone(tsk) || $(esc(then))
        tsk.result
    end
end

macro src()
    @static if VERSION >= v"0.7-" && length(:(@test).args) == 2
        esc(quote
            (__module__,
             __source__.file == nothing ? "?" : String(__source__.file),
             __source__.line)
        end)
    else
        esc(quote
            (current_module(),
             (p = Base.source_path(); p == nothing ? "REPL" : p),
             Int(unsafe_load(cglobal(:jl_lineno, Cint))))
        end)
    end
end

"""
    @debug DEBUG expr
    @debug DEBUG "message"

A macro to aid when needing to turn on extremely verbose output for debugging.
Set `const DEBUG = true` in HTTP.jl and re-compile the package to see
debug-level output from the package. When `DEBUG = false`, all `@debug` statements
compile to `nothing`.
"""
macro debug(should, expr)
    m, f, l = @src()
    if typeof(expr) == String
        e = esc(:(println("[DEBUG - ", $m, '.', $f, ":", $(rpad(l, 5, ' ')), "]: ", $(escape_string(expr)))))
    else
        e = esc(:(println("[DEBUG - ", $m, '.', $f, ":", $(rpad(l, 5, ' ')), "]: ", $(sprint(Base.show_unquoted, expr)), " = ", escape_string(string($expr)))))
    end
    return quote
        @static if $should
            $e
        end
    end
end

macro log(stmt)
    # "[HTTP]: Connecting to remote host..."
    return esc(:(verbose && (write(logger, "[HTTP - $(rpad(now(), 23, ' '))]: $($stmt)\n"); flush(logger))))
end

# parsing utils
macro anyeq(var, vals...)
    ret = e = Expr(:||)
    for (i, v) in enumerate(vals)
        x = :($var == $v)
        push!(e.args, x)
        i >= length(vals) - 1 && continue
        ne = Expr(:||)
        push!(e.args, ne)
        e = ne
    end
    return esc(ret)
end

@inline islower(b::UInt8) = UInt8('a') <= b <= UInt8('z')
@inline isupper(b::UInt8) = UInt8('A') <= b <= UInt8('Z')
@inline lower(c::UInt8) = c | 0x20
@inline lower(c) = Char(UInt32(c) | 0x20)
@inline isurlchar(c) =  c > '\u80' ? true : normal_url_char[Int(c) + 1]
@inline ismark(c) = @anyeq(c, '-', '_', '.', '!', '~', '*', '\'', '(', ')')
@inline isalpha(c) = 'a' <= lower(c) <= 'z'
@inline isnum(c) = '0' <= c <= '9'
@inline isalphanum(c) = isalpha(c) || isnum(c)
@inline isuserinfochar(c) = isalphanum(c) || ismark(c) || @anyeq(c, '%', ';', ':', '&', '=', '+', '$', ',')
@inline ishex(c) =  isnum(c) || ('a' <= lower(c) <= 'f')
@inline ishostchar(c) = isalphanum(c) || @anyeq(c, '.', '-', '_', '~')
@inline isheaderchar(c) = c == CR || c == LF || c == Char(9) || (c > Char(31) && c != Char(127))

macro shifted(meth, i, char)
    return esc(:(Int($meth) << Int(16) | Int($i) << Int(8) | Int($char)))
end

macro errorif(cond, err)
    return esc(quote
        $cond && @err($err)
    end)
end

macro err(e)
    return esc(quote
        errno = $e
        @goto error
    end)
end

macro strictcheck(cond)
    return esc(:(strict && @errorif($cond, HPE_STRICT)))
end

# ensure the first character and subsequent characters that follow a '-' are uppercase
function canonicalize!(s::String)
    toUpper = UInt8('A') - UInt8('a')
    bytes = Vector{UInt8}(s)
    upper = true
    for i = 1:length(bytes)
        @inbounds b = bytes[i]
        if upper
            islower(b) && (bytes[i] = b + toUpper)
        else
            isupper(b) && (bytes[i] = lower(b))
        end
        upper = b == UInt8('-')
    end
    return s
end


text:
using LinearAlgebra
using PyCall
@pyimport ehtim as eh

function amp_debias(amp, σ; force_nonzero=false)
    amp = abs(amp)
    σ = abs(σ)
    if amp < σ
        if force_nonzero
            deb² = σ^2
        else
            deb² = 0
        end
    else
        deb² = amp^2 - σ^2
    end
    return √deb²
end

function logcamp_debias(log_camp, snr1, snr2, snr3, snr4)
    log_camp_debias = log_camp + 0.5*(1/snr1^2 + 1/snr2^2 - 1/snr3^2 -1/snr4^2)
    return log_camp_debias
end

"""
    Stealing this from eht-imaging, as long as
    we use the quadrangles from ehtim then
    the definitions of numerators and denominators
    will work out.
"""
function make_log_closure_amplitude(n1amp, n2amp, d1amp, d2amp, n1err, n2err, d1err, d2err; debias=true)
    p1 = abs(n1amp)
    p2 = abs(n2amp)
    p3 = abs(d1amp)
    p4 = abs(d2amp)

    if debias
        p1 = amp_debias(p1, n1err)
        p2 = amp_debias(p2, n2err)
        p3 = amp_debias(p3, d1err)
        p4 = amp_debias(p4, d2err)
    end
    snr1 = p1/n1err
    snr2 = p2/n2err
    snr3 = p3/d1err
    snr4 = p4/d2err

    logcamp = log(p1)+log(p2)-log(p3)-log(p4)
    logcamp_err = √(1/snr1^2 + 1/snr2^2 + 1/snr3^2 + 1/snr4^2)
    if debias
        logcamp = logcamp_debias(logcamp, snr1, snr2, snr3, snr4)
    end
    return logcamp, logcamp_err
end

function amp_add_syserr(amp, amp_error; fractional=0, additive=0)
    σ = √(amp_error^2+(fractional*amp)^2+additive^2)
    return amp, σ
end

function vis_add_syserr(vis, amp_error; fractional=0, additive=0)
    return amp_add_syserr(abs(vis), amp_error, fractional=fractional, additive=additive)
end

function logcamp_add_syserr(n1amp, n2amp, d1amp, d2amp, n1err, n2err, d1err, d2err; fractional=0, additive = 0, debias=true)
    n1amp, n1err = amp_add_syserr(n1amp, n1err, fractional=fractional, additive=additive)
    n2amp, n2err = amp_add_syserr(n2amp, n2err, fractional=fractional, additive=additive)
    d1amp, d1err = amp_add_syserr(d1amp, d1err, fractional=fractional, additive=additive)
    d2amp, d2err = amp_add_syserr(d2amp, d2err, fractional=fractional, additive=additive)
    return make_log_closure_amplitude(n1amp, n2amp, d1amp, d2amp, n1err, n2err, d1err, d2err, debias=debias)
end

function make_bispectrum(v1, v2, v3, v1err, v2err, v3err)
    bi = v1*v2*v3
    bisig = abs(bi)*√((v1err/abs(v1))^2+(v2err/abs(v2))^2+(v3err/abs(v3))^2)
    return bi, bisig
end

function closure_phase_from_bispectrum(bi, bisig)
    cphase = angle(bi)
    cphase_error = bisig / abs(bi)
    return cphase, cphase_error
end

function bispectrum_add_syserr(v1, v2, v3, v1err, v2err, v3err; fractional=0, additive=0)
    v1, v1err =vis_add_syserr(v1, v1err, fractional=fractional, additive=additive)
    v2, v2err =vis_add_syserr(v2, v2err, fractional=fractional, additive=additive)
    v3, v3err =vis_add_syserr(v3, v3err, fractional=fractional, additive=additive)
    return make_bispectrum(v1, v2, v3, v1err, v2err, v3err)
end

function cphase_add_syserr(v1, v2, v3, v1err, v2err, v3err; fractional=0, additive=0)
    bi, bisig = bispectrum_add_syserr(v1, v2, v3, v1err, v2err, v3err, fractional=fractional, additive=additive)
    return closure_phase_from_bispectrum(bi, bisig)
end


text: const enable_opengl_debugging = false