text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module SelfConvolutive

using Nemo

export ModuleSelfConvolutive, SelfConvRec
export L000698, L001710, L003319, L005411, L005412, L006012, L006318
export L047891, L062980, L082298, L082301, L082302, L105523, L107716
export L111529, L111530, L111531, L111532, L111533, L146559, L167872

"""
Some exactly solvable self-convolutive recurrences.

* SelfConvRec, L000698, L001710, L003319, L005411, L005412, L006012, L006318
L047891, L062980, L082298, L082301, L082302, L105523, L107716
L111529, L111530, L111531, L111532, L111533, L146559, L167872
"""
const ModuleSelfConvolutive = ""

# R. J. Martin and M. J. Kearney,
# An exactly solvable self-convolutive recurrence
# http://arXiv.org/abs/1103.4936
# Integral representation of certain combinatorial recurrences
# see also http://oeis.org/A172455

Convolution(A, n) = sum(A[j+1] * A[n-j+1] for j in 0:n)

"""
An exactly solvable self-convolutive recurrence studied by R. J. Martin and M. J. Kearney.
"""
function SelfConvRec(len, a, b, c)
    A = fill(ZZ(0), len)
    len == 0 && return A
    A[1] = 1

    for n in 1:len - 1
        A[n+1] = (a * (n + 1) + b) * A[n] + c * Convolution(A, n-1)
    end
    A
end

# """
# Return the powers of 2.
# """
# L000079(len::Int) = SelfConvRec(len, 1, 1, -1)
# """
# Return the Catalan numbers.
# """
# L000108(len::Int) = SelfConvRec(len, 0, 0, 1)
# """
# Return the factorial numbers.
# """
# L000142(len::Int) = SelfConvRec(len, 1, -1, 0)
# """
# Return the powers of 3.
# """
# L000244(len::Int) = SelfConvRec(len, 2, 1, -2)
# """
# Return the powers of 5.
# """
# L000351(len::Int) = SelfConvRec(len, 4, 1, -4)
# """
# Return the powers of 6.
# """
# L000400(len::Int) = SelfConvRec(len, 5, 1, -5)
# """
# Return the powers of 7.
# """
# L000420(len::Int) = SelfConvRec(len, 6, 1, -6)

"""
Return the number of indecomposable perfect matchings on ``[2n]``.
"""
L000698(len::Int) = SelfConvRec(len, 2, -3, 1)
"""
Return the order of alternating group ``A_n``, or number of even permutations of ``n`` letters.
"""
L001710(len::Int) = SelfConvRec(len, 1, 1, 0)

# """
# Return n!/3!.
# """
# L001715(len::Int) = SelfConvRec(len, 1, 2, 0)
# """
# Return n!/4!.
# """
# L001720(len::Int) = SelfConvRec(len, 1, 3, 0)
# """
# Return n!/5!.
# """
# L001725(len::Int) = SelfConvRec(len, 1, 4, 0)
# """
# Return n!/6!.
# """
# L001730(len::Int) = SelfConvRec(len, 1, 5, 0)

"""
Return the number of connected permutations of ``[1..n]``. Also called indecomposable permutations.
"""
L003319(len::Int) = SelfConvRec(len, 1, -2, 1)
"""
Return the number of non-vanishing Feynman diagrams of order ``2n`` for the electron or the photon propagators in quantum electrodynamics.
"""
L005411(len::Int) = SelfConvRec(len, 2, -4, 1)
"""
Return the number of non-vanishing Feynman diagrams of order ``2n`` for the vacuum polarization in quantum electrodynamics.
"""
L005412(len::Int) = SelfConvRec(len, 2, -2, 1)
"""
Counting some sets of permutations.
"""
L006012(len::Int) = SelfConvRec(len, -1, 2, 2)
"""
Return the large Schröder numbers.
"""
L006318(len::Int) = SelfConvRec(len, 0, 1, 1)
"""
Return the number of planar rooted trees with ``n`` nodes and tricolored end nodes.
"""
L047891(len::Int) = SelfConvRec(len, 0, 2, 1)
# """
# Return (n+7)!/7!.
# """
# L049388(len::Int) = SelfConvRec(len, 1, 6, 0)
# """
# Return (3*n+4)!!!/4!!!.
# """
# L051604(len::Int) = SelfConvRec(len, 3, 1, 0)
# """
# Return (3*n+5)!!!/5!!!.
# """
# L051605(len::Int) = SelfConvRec(len, 3, 2, 0)
# """
# Return (3*n+6)!!!/6!!!.
# """
# L051606(len::Int) = SelfConvRec(len, 3, 3, 0)
# """
# Return (3*n+7)!!!/7!!!.
# """
# L051607(len::Int) = SelfConvRec(len, 3, 4, 0)
# """
# Return (3*n+8)!!!/8!!!.
# """
# L051608(len::Int) = SelfConvRec(len, 3, 5, 0)
# """
# Return (3*n+9)!!!/9!!!.
# """
# L051609(len::Int) = SelfConvRec(len, 3, 6, 0)
# """
# Return (4*n+5)(!^4)/5(!^4).
# """
# L051617(len::Int) = SelfConvRec(len, 4, 1, 0)
# """
# Return  (4*n+6)(!^4)/6(!^4).
# """
# L051618(len::Int) = SelfConvRec(len, 4, 2, 0)
# """
# Return (4*n+7)(!^4)/7(!^4).
# """
# L051619(len::Int) = SelfConvRec(len, 4, 3, 0)
# """
# Return (4*n+8)(!^4)/8(!^4).
# """
# L051620(len::Int) = SelfConvRec(len, 4, 4, 0)
# """
# Return (4*n+9)(!^4)/9(!^4).
# """
# L051621(len::Int) = SelfConvRec(len, 4, 5, 0)
# """
# Return (4*n+10)(!^4)/10(!^4).
# """
# L051622(len::Int) = SelfConvRec(len, 4, 6, 0)
# """
# Return (5*n+6)(!^5)/6.
# """
# L051687(len::Int) = SelfConvRec(len, 5, 1, 0)
# """
# Return (5*n+7)(!^5)/7(!^5).
# """
# L051688(len::Int) = SelfConvRec(len, 5, 2, 0)
# """
# Return (5*n+8)(!^5)/8(!^5).
# """
# L051689(len::Int) = SelfConvRec(len, 5, 3, 0)
# """
# Return (5*n+9)(!^5)/9(!^5).
# """
# L051690(len::Int) = SelfConvRec(len, 5, 4, 0)
# """
# Return (5*n+10)(!^5)/10(!^5).
# """
# L051691(len::Int) = SelfConvRec(len, 5, 5, 0)
# """
# Return ((6*n+7)(!^6))/7.
# """
# L053100(len::Int) = SelfConvRec(len, 6, 1, 0)
# """
# Return ((6*n+8)(!^6))/8(!^6).
# """
# L053101(len::Int) = SelfConvRec(len, 6, 2, 0)
# """
# Return ((6*n+9)(!^6))/9(!^6).
# """
# L053102(len::Int) = SelfConvRec(len, 6, 3, 0)
# """
# Return ((6*n+10)(!^6))/10(!^6).
# """
# L053103(len::Int) = SelfConvRec(len, 6, 4, 0)
# """
# Return ((7*n+8)(!^7))/8.
# """
# L053104(len::Int) = SelfConvRec(len, 7, 1, 0)
# """
# Return ((7*n+9)(!^7))/9(!^7).
# """
# L053105(len::Int) = SelfConvRec(len, 7, 2, 0)
# """
# Return ((7*n+10)(!^7))/10(1^7).
# """
# L053106(len::Int) = SelfConvRec(len, 7, 3, 0)

"""
Return the number of rooted unlabeled connected triangular maps on a compact closed oriented surface with ``2n`` faces.
"""
L062980(len::Int) = SelfConvRec(len, 6, -8, 1)
"""
Return the number of lattice paths from ``(0,0)`` to ``(n+1,n+1)`` that consist of steps ``(i,0)`` and ``(0,j)`` with ``i,j≥1`` and stay strictly below the diagonal line ``y=x`` except at the endpoints.
"""
L082298(len::Int) = SelfConvRec(len, 0, 3, 1)
"""
Return the number of Schröder paths of semilength n in which the (2,0)-steps come in 4 colors.
"""
L082301(len::Int) = SelfConvRec(len, 0, 4, 1)
"""
Return the coefficients in the expansion of ``(1-5x-√(25x^2-14x+1))/(2x)``.
"""
L082302(len::Int) = SelfConvRec(len, 0, 5, 1)

# """
# Return the coefficients in the expansion of ``(1-6x-√(36x^2-16x+1))/(2x)``.
# """
# L082305(len::Int) = SelfConvRec(len, 0, 6, 1)
# """
# Return the coefficients in the expansion of ``(1-7x-√(49x^2-18x+1))/(2x)``.
# """
# L082366(len::Int) = SelfConvRec(len, 0, 7, 1)
# """
# Return the coefficients in the expansion of ``(1-8x-√(64x^2-20x+1))/(2*x)``.
# """
# L082367(len::Int) = SelfConvRec(len, 0, 8, 1)

"""
Return the coefficients in the expansion of ``(1+2x-√(1+4x^2))/(2x)``.
"""
L105523(len::Int) = SelfConvRec(len, 0, -2, 1)
"""
A transform the of triple factorial numbers.
"""
L107716(len::Int) = SelfConvRec(len, 3, -4, 1)
"""
A convolutory inverse of the factorial sequence.
"""
L111529(len::Int) = SelfConvRec(len, 1, -3, 2)
"""
A convolutory inverse of the factorial sequence.
"""
L111530(len::Int) = SelfConvRec(len, 1, -4, 3)
"""
A convolutory inverse of the factorial sequence.
"""
L111531(len::Int) = SelfConvRec(len, 1, -5, 4)
"""
A convolutory inverse of the factorial sequence.
"""
L111532(len::Int) = SelfConvRec(len, 1, -6, 5)
"""
A convolutory inverse of the factorial sequence.
"""
L111533(len::Int) = SelfConvRec(len, 1, -7, 6)

#"""
#"""
#L111546(len::Int) = SelfConvRec(len, 1, 0, 1)
#"""
#"""
#L111556(len::Int) = SelfConvRec(len, 1, 1, 1)

"""
Return the expansion of ``(1-x)/(1 - 2x + 2x^2)``.
"""
L146559(len::Int) = SelfConvRec(len, 1, 1, -2)
"""
Return half the number of Feynman diagrams of order ``2(n+1)``, for the electron self-energy in quantum electrodynamics.
"""
L167872(len::Int) = SelfConvRec(len, 2, -3, 2)

# """
# """
# L172450(len::Int) = SelfConvRec(len, 2, 0, -1)
# """
# """
# L172455(len::Int) = SelfConvRec(len, 6, -4, -1)
# """
# """
# L172485(len::Int) = SelfConvRec(len, -1, -2, 3)
# """
# """
# L177354(len::Int) = SelfConvRec(len, 1, 2, 1)

#START-TEST-########################################################
using Test, SeqTests

function test()
    @testset "SelfConvRec" begin
        if is_oeis_installed()

            L = [L005411, L005412, L006012, L006318, L047891, L062980, L082298,
                L082301, L082302, L105523, L107716, L111529, L111530, L111531,
                L111532, L111533, L146559, L167872]
            SeqTest(L, 'L')
        end
    end
end

function demo()

    for n in 0:6 println(n, " ↦ ", L005411(n)) end

    f(len) = SelfConvRec(len, 2, -1, -1)
    println(f(11))
end

"""
L005411(500) :: 0.238653 seconds (252.99 k allocations: 3.887 MB)
"""
function perf()
    GC.gc()
    @time L005411(500)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
L --> A000079
Nemo.fmpz[1, 2, 4, 8, 16, 32]
Nemo.fmpz[1, 2, 4, 8, 16, 32]
L --> A000108
Nemo.fmpz[1, 1, 2, 5, 14, 42]
Nemo.fmpz[1, 1, 2, 5, 14, 42]
L --> A000142
Nemo.fmpz[1, 1, 2, 6, 24, 120]
Nemo.fmpz[1, 1, 2, 6, 24, 120]
L --> A000244
Nemo.fmpz[1, 3, 9, 27, 81, 243]
Nemo.fmpz[1, 3, 9, 27, 81, 243]
L --> A000351
Nemo.fmpz[1, 5, 25, 125, 625, 3125]
Nemo.fmpz[1, 5, 25, 125, 625, 3125]
L --> A000400
Nemo.fmpz[1, 6, 36, 216, 1296, 7776]
Nemo.fmpz[1, 6, 36, 216, 1296, 7776]
L --> A000420
Nemo.fmpz[1, 7, 49, 343, 2401, 16807]
Nemo.fmpz[1, 7, 49, 343, 2401, 16807]
#############################
L --> A000698
Nemo.fmpz[1, 1, 2, 10, 74, 706]
Nemo.fmpz[1, 2, 10, 74, 706, 8162]
#####################
seq((n+2)!/2, n=0..9);
L --> A001710
Nemo.fmpz[1, 1, 1, 3, 12, 60]
Nemo.fmpz[1, 3, 12, 60, 360, 2520]
######################
L --> A001715
Nemo.fmpz[1, 4, 20, 120, 840, 6720]
Nemo.fmpz[1, 4, 20, 120, 840, 6720]
L --> A001720
Nemo.fmpz[1, 5, 30, 210, 1680, 15120]
Nemo.fmpz[1, 5, 30, 210, 1680, 15120]
L --> A001725
Nemo.fmpz[1, 6, 42, 336, 3024, 30240]
Nemo.fmpz[1, 6, 42, 336, 3024, 30240]
L --> A001730
Nemo.fmpz[1, 7, 56, 504, 5040, 55440]
Nemo.fmpz[1, 7, 56, 504, 5040, 55440]
##################
L --> A003319
Nemo.fmpz[0, 1, 1, 3, 13, 71]
Nemo.fmpz[1, 1, 3, 13, 71, 461]
##################
L --> A005411
Nemo.fmpz[1, 1, 4, 25, 208, 2146]
Nemo.fmpz[1, 1, 4, 25, 208, 2146]
L --> A005412
Nemo.fmpz[1, 3, 18, 153, 1638, 20898]
Nemo.fmpz[1, 3, 18, 153, 1638, 20898]
L --> A006012
Nemo.fmpz[1, 2, 6, 20, 68, 232]
Nemo.fmpz[1, 2, 6, 20, 68, 232]
L --> A006318
Nemo.fmpz[1, 2, 6, 22, 90, 394]
Nemo.fmpz[1, 2, 6, 22, 90, 394]
L --> A047891
Nemo.fmpz[1, 3, 12, 57, 300, 1686]
Nemo.fmpz[1, 3, 12, 57, 300, 1686]
L --> A049388
Nemo.fmpz[1, 8, 72, 720, 7920, 95040]
Nemo.fmpz[1, 8, 72, 720, 7920, 95040]
L --> A051604
Nemo.fmpz[1, 7, 70, 910, 14560, 276640]
Nemo.fmpz[1, 7, 70, 910, 14560, 276640]
L --> A051605
Nemo.fmpz[1, 8, 88, 1232, 20944, 418880]
Nemo.fmpz[1, 8, 88, 1232, 20944, 418880]
L --> A051606
Nemo.fmpz[1, 9, 108, 1620, 29160, 612360]
Nemo.fmpz[1, 9, 108, 1620, 29160, 612360]
L --> A051607
Nemo.fmpz[1, 10, 130, 2080, 39520, 869440]
Nemo.fmpz[1, 10, 130, 2080, 39520, 869440]
L --> A051608
Nemo.fmpz[1, 11, 154, 2618, 52360, 1204280]
Nemo.fmpz[1, 11, 154, 2618, 52360, 1204280]
L --> A051609
Nemo.fmpz[1, 12, 180, 3240, 68040, 1632960]
Nemo.fmpz[1, 12, 180, 3240, 68040, 1632960]
L --> A051617
Nemo.fmpz[1, 9, 117, 1989, 41769, 1044225]
Nemo.fmpz[1, 9, 117, 1989, 41769, 1044225]
L --> A051618
Nemo.fmpz[1, 10, 140, 2520, 55440, 1441440]
Nemo.fmpz[1, 10, 140, 2520, 55440, 1441440]
L --> A051619
Nemo.fmpz[1, 11, 165, 3135, 72105, 1946835]
Nemo.fmpz[1, 11, 165, 3135, 72105, 1946835]
L --> A051620
Nemo.fmpz[1, 12, 192, 3840, 92160, 2580480]
Nemo.fmpz[1, 12, 192, 3840, 92160, 2580480]
L --> A051621
Nemo.fmpz[1, 13, 221, 4641, 116025, 3364725]
Nemo.fmpz[1, 13, 221, 4641, 116025, 3364725]
L --> A051622
Nemo.fmpz[1, 14, 252, 5544, 144144, 4324320]
Nemo.fmpz[1, 14, 252, 5544, 144144, 4324320]
L --> A051687
Nemo.fmpz[1, 11, 176, 3696, 96096, 2978976]
Nemo.fmpz[1, 11, 176, 3696, 96096, 2978976]
L --> A051688
Nemo.fmpz[1, 12, 204, 4488, 121176, 3877632]
Nemo.fmpz[1, 12, 204, 4488, 121176, 3877632]
L --> A051689
Nemo.fmpz[1, 13, 234, 5382, 150696, 4972968]
Nemo.fmpz[1, 13, 234, 5382, 150696, 4972968]
L --> A051690
Nemo.fmpz[1, 14, 266, 6384, 185136, 6294624]
Nemo.fmpz[1, 14, 266, 6384, 185136, 6294624]
L --> A051691
Nemo.fmpz[1, 15, 300, 7500, 225000, 7875000]
Nemo.fmpz[1, 15, 300, 7500, 225000, 7875000]
L --> A053100
Nemo.fmpz[1, 13, 247, 6175, 191425, 7082725]
Nemo.fmpz[1, 13, 247, 6175, 191425, 7082725]
L --> A053101
Nemo.fmpz[1, 14, 280, 7280, 232960, 8852480]
Nemo.fmpz[1, 14, 280, 7280, 232960, 8852480]
L --> A053102
Nemo.fmpz[1, 15, 315, 8505, 280665, 10945935]
Nemo.fmpz[1, 15, 315, 8505, 280665, 10945935]
L --> A053103
Nemo.fmpz[1, 16, 352, 9856, 335104, 13404160]
Nemo.fmpz[1, 16, 352, 9856, 335104, 13404160]
L --> A053104
Nemo.fmpz[1, 15, 330, 9570, 344520, 14814360]
Nemo.fmpz[1, 15, 330, 9570, 344520, 14814360]
L --> A053105
Nemo.fmpz[1, 16, 368, 11040, 408480, 17973120]
Nemo.fmpz[1, 16, 368, 11040, 408480, 17973120]
L --> A053106
Nemo.fmpz[1, 17, 408, 12648, 480624, 21628080]
Nemo.fmpz[1, 17, 408, 12648, 480624, 21628080]
L --> A062980
Nemo.fmpz[1, 5, 60, 1105, 27120, 828250]
Nemo.fmpz[1, 5, 60, 1105, 27120, 828250]
L --> A082298
Nemo.fmpz[1, 4, 20, 116, 740, 5028]
Nemo.fmpz[1, 4, 20, 116, 740, 5028]
L --> A082301
Nemo.fmpz[1, 5, 30, 205, 1530, 12130]
Nemo.fmpz[1, 5, 30, 205, 1530, 12130]
L --> A082302
Nemo.fmpz[1, 6, 42, 330, 2814, 25422]
Nemo.fmpz[1, 6, 42, 330, 2814, 25422]
L --> A082305
Nemo.fmpz[1, 7, 56, 497, 4760, 48174]
Nemo.fmpz[1, 7, 56, 497, 4760, 48174]
L --> A082366
Nemo.fmpz[1, 8, 72, 712, 7560, 84616]
Nemo.fmpz[1, 8, 72, 712, 7560, 84616]
L --> A082367
Nemo.fmpz[1, 9, 90, 981, 11430, 140058]
Nemo.fmpz[1, 9, 90, 981, 11430, 140058]
L --> A105523
Nemo.fmpz[1, -1, 0, 1, 0, -2]
Nemo.fmpz[1, -1, 0, 1, 0, -2]
L --> A107716
Nemo.fmpz[1, 3, 21, 219, 2973, 49323]
Nemo.fmpz[1, 3, 21, 219, 2973, 49323]
L --> A111529
Nemo.fmpz[1, 1, 4, 22, 148, 1156]
Nemo.fmpz[1, 1, 4, 22, 148, 1156]
L --> A111530
Nemo.fmpz[1, 1, 5, 33, 261, 2361]
Nemo.fmpz[1, 1, 5, 33, 261, 2361]
L --> A111531
Nemo.fmpz[1, 1, 6, 46, 416, 4256]
Nemo.fmpz[1, 1, 6, 46, 416, 4256]
L --> A111532
Nemo.fmpz[1, 1, 7, 61, 619, 7045]
Nemo.fmpz[1, 1, 7, 61, 619, 7045]
L --> A111533
Nemo.fmpz[1, 1, 8, 78, 876, 10956]
Nemo.fmpz[1, 1, 8, 78, 876, 10956]
L --> A111546
Nemo.fmpz[1, 3, 15, 99, 783, 7083]
Nemo.fmpz[1, 3, 15, 99, 783, 7083]
L --> A111556
Nemo.fmpz[1, 4, 24, 184, 1664, 17024]
Nemo.fmpz[1, 4, 24, 184, 1664, 17024]
L --> A146559
Nemo.fmpz[1, 1, 0, -2, -4, -4]
Nemo.fmpz[1, 1, 0, -2, -4, -4]
L --> A167872
Nemo.fmpz[1, 3, 21, 207, 2529, 36243]
Nemo.fmpz[1, 3, 21, 207, 2529, 36243]
L --> A172450
Nemo.fmpz[1, 3, 12, 63, 432, 3798]
Nemo.fmpz[1, 3, 12, 63, 432, 3798]
L --> A172485
Nemo.fmpz[1, -1, -1, 3, 3, -21]
Nemo.fmpz[1, -1, -1, 3, 3, -21]
#################################
L --> A177354
Nemo.fmpz[5, 35, 305, 3095, 35225, 439775]
Nemo.fmpz[1, 5, 35, 305, 3095, 35225]
###########################

L = [ L000079, L000108, L000142, L000244, L000351, L000400,
 L000420, L001715, L001720, L001725, L001730,
 L005411, L005412, L006012, L006318, L047891, L049388, L051604,
 L051605, L051606, L051607, L051608, L051609, L051617, L051618,
 L051619, L051620, L051621, L051622, L051687, L051688, L051689,
 L051690, L051691, L053100, L053101, L053102, L053103, L053104,
 L053105, L053106, L062980, L082298, L082301, L082302, L082305,
 L082366, L082367, L105523, L107716, L111529, L111530, L111531,
 L111532, L111533, L111546, L111556, L146559, L167872,
 L172450, L172485 ]
# L000698, L001710, L003319, L177354
SeqTest(L, 'L')

=#


text:
# NOTE: these are older implementations/functions which were used in the past
# but have now been optimised heavily and ultimately replaced by other
# functions they're only kept here for the tests.

function sum_var_sij(Z::AbstractMatrix, CovZ::AbstractMatrix, n::Int,
                     corrected=false; with_diag=true)
    scale = ifelse(corrected, n-1, n)
    Z²    = Z.^2
    π̂mat  = ((Z²'*Z²)/n - CovZ.^2) / scale
    with_diag && return sum(π̂mat)
    return sum(π̂mat) - sum(diag(π̂mat))
end

function centercols(X::AbstractMatrix, μ=nothing; dims=1)
    Xc = (dims == 1) ? copy(X) : transpose(X)
    if μ === nothing
        Xc .-= mean(X, dims=1)
    elseif iszero(μ)
        # assume X is centered already
    elseif μ isa AbstractVector
        @assert length(μ) == size(Xc, 2) "Dimension mismatch, given mean has length $(length(μ)) while p=$(size(Xc, 2))"
        Xc .-= μ
    else
        throw(ArgumentError("Incorrect assignment for `mean`, expected either 0, nothing or a vector, got $μ."))
    end
    Xc
end

ext: jl


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCellDG
include("../../test/useful_routines.jl")
include("circular_bc_transformation_elasticity_solver.jl")

function compute_interface_maxnorm_stress_error(
    polyoder,
    nelmts,
    penaltyfactor,
)
    width = 1.0
    K1, K2 = 247.0, 192.0
    mu1, mu2 = 126.0, 87.0
    lambda1 = lame_lambda(K1, mu1)
    lambda2 = lame_lambda(K2, mu2)
    stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)

    theta0 = -0.067
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    meshwidth = [width, width]
    numqp = required_quadrature_order(polyorder) + 2

    dx = width / nelmts
    penalty = penaltyfactor / dx * 0.5 * (lambda1 + mu1 + lambda2 + mu2)

    basis = TensorProductBasis(2, polyorder)
    interfacecenter = [0.5, 0.5]
    interfaceradius = width / 3.0
    outerradius = 2.0
    analyticalsolution = AnalyticalSolution(
        interfaceradius,
        outerradius,
        interfacecenter,
        lambda1,
        mu1,
        lambda2,
        mu2,
        theta0,
    )

    mesh, cellquads, facequads, interfacequads = construct_mesh_and_quadratures(
        meshwidth,
        nelmts,
        basis,
        interfacecenter,
        interfaceradius,
        numqp,
    )
    nodaldisplacement = nodal_displacement(
        mesh,
        basis,
        cellquads,
        facequads,
        interfacequads,
        stiffness,
        theta0,
        analyticalsolution,
        penalty,
    )

    parenterror, producterror = interface_stress_maxnorm_error(
        nodaldisplacement,
        basis,
        interfacequads,
        stiffness,
        transfstress,
        theta0,
        mesh,
        x -> core_stress(analyticalsolution),
        x -> shell_stress(analyticalsolution, x),
    )

    parentmaxnorm = maxnorm_on_interface(
        x -> core_stress(analyticalsolution),
        interfacequads,
        -1,
        mesh,
        4,
    )
    productmaxnorm = maxnorm_on_interface(
        x -> shell_stress(analyticalsolution, x),
        interfacequads,
        +1,
        mesh,
        4,
    )

    return parenterror ./ parentmaxnorm, producterror ./ productmaxnorm
end

polyorder = 2
penaltyfactor = 1e2
powers = [3,4,5]
nelmts = [2^p for p in powers]

stresserror2 = [compute_interface_maxnorm_stress_error(polyorder, ne, penaltyfactor) for ne in nelmts]

parenterror2 = [st[1] for st in stresserror2]
producterror2 = [st[2] for st in stresserror2]

parenterror2 = Array(transpose(hcat(parenterror2...)))
producterror2 = Array(transpose(hcat(producterror2...)))

dx = 1.0 ./ nelmts

parentrates2 = mapslices(x->convergence_rate(dx,x),parenterror2,dims=1)
parentrates2 = parentrates2[:,[1,2,4]]
productrates2 = mapslices(x->convergence_rate(dx,x),producterror2,dims=1)


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

module ShootoutPerf
import Perftests: @perf, meta

# Find relative paths easily
rpath(filename) = joinpath(dirname(@__FILE__), filename)
function download_if_needed(url)
    path = joinpath(dirname(@__FILE__),basename(url))
    if !isfile(path)
        println("Downloading $(basename(url))")
        download(url, path)
    end
end

# If we don't already have the input files we need, download them!
cache_url = "https://cache.e.ip.saba.us/http://benchmarksgame.alioth.debian.org/download"
in_and_outs = ["knucleotide", "regexdna", "revcomp"]
just_outs = ["spectralnorm",  "fasta", "mandelbrot", "nbody", "fannkuchredux", "pidigits", "meteor"]
for file in in_and_outs
    download_if_needed("$cache_url/$file-input.txt")
    download_if_needed("$cache_url/$file-output.txt")
end
for file in just_outs
    download_if_needed("$cache_url/$file-output.txt")
end

include("binary_trees.jl")
@perf binary_trees(10) meta("binary_trees", "Allocate and deallocate many many binary trees")

include("fannkuch.jl")
@perf fannkuch(7) meta("fannkuch", "Indexed-access to tiny integer-sequence")

include("fasta.jl")
@perf fasta(100) meta("fasta", "Generate and write random DNA sequences")

include("k_nucleotide.jl")
infile = rpath("knucleotide-input.txt")
@perf k_nucleotide(infile) meta("k_nucleotide", "Hashtable update and k-nucleotide strings")

include("mandelbrot.jl")
outfile = rpath("mandelbrot-output-julia.txt")
@perf mandelbrot(200, outfile) meta("mandelbrot", "Generate Mandelbrot set portable bitmap file")

include("meteor_contest.jl")
@perf meteor_contest() meta("meteor_contest", "Search for solutions to shape packing puzzle")

include("nbody.jl")
@perf NBody.nbody() meta("nbody", "Double-precision N-body simulation")

include("nbody_vec.jl")
@perf NBodyVec.nbody_vec() meta("nbody_vec", "A vectorized double-precision N-body simulation")

include("pidigits.jl")
@assert pidigits(1000) == 9216420198
@perf pidigits(1000) meta("pidigits", "Streaming arbitrary-precision arithmetic")

include("regex_dna.jl")
infile = rpath("regexdna-input.txt")
@perf regex_dna(infile) meta("regex_dna", "Match DNA 8-mers and substitute nucleotides for IUB codes")

include("revcomp.jl")
infile = rpath("revcomp-input.txt")
@perf revcomp(infile) meta("revcomp", "Read DNA sequences - write their reverse-complement")

include("spectralnorm.jl")
@perf spectralnorm() meta("spectralnorm", "Eigenvalue using the power method")

end # module


text:
using QCMaterial
using PyCall
using LinearAlgebra

scipy_opt = pyimport("scipy.optimize")
struct Hamiltonian
    ham
    n_sites
end

hfstate(n_qubit, n_electron) = parse(Int, repeat("0", n_qubit-n_electron) * repeat("1", n_electron), base=2)

function generate_ham(nsite)
    # ## 2-site Hubbard modelのハミルトニアンを定義する。
    #
    # $H=-t \sum_{\langle i, j>\sigma=(\uparrow, \downarrow)} \sum_{i, \sigma}\left(a_{i, \sigma}^{\dagger} a_{j, \sigma}+a_{j, \sigma}^{\dagger} a_{i, \sigma}\right)+U \sum_{i} n_{i, \uparrow} n_{i, \downarrow}-\mu \sum_{i} \sum_{\sigma=(\uparrow, \downarrow)} n_{i, \sigma}$
    #
    # $1_\uparrow$->0,$1_\downarrow$->1, $2_\uparrow$->2, $2_\downarrow$->3
    #
    # クーロン斥力(U=2): 2 * $n_{1,\uparrow}$ $n_{1,\downarrow}$ + 2 * $n_{2,\uparrow}$ $n_{2,\downarrow}$
    #                = 2 * $c^\dagger_0$ $c_0$ $c^\dagger_1$ $c_1$ + 2 * $c^\dagger_2$ $c_2$ $c^\dagger_3$ $c_3$
    #
    # ホッピング項(t=-0.1): -0.1 * $c^\dagger_{1,\uparrow}$ $c_{2,\uparrow}$ -0.1 * $c^\dagger_{2,\uparrow}$ $c_{1,\uparrow}$
    #                    -0.1 * $c^\dagger_{1,\downarrow}$ $c_{2,\downarrow}$ -0.1 * $c^\dagger_{2,\downarrow}$ $c_{1,\downarrow}$
    #   = -0.1 * $c^\dagger_0$ $c_2$ + -0.1 * $c^\dagger_2$ $c_0$ + -0.1 * $c^\dagger_1$ $c_3$ + -0.1 * $c_3^\dagger$ $c_1$
    #
    # ケミカルポテンシャル項(μ=U/2=1):-1 * $n_{1,\uparrow}$ - 1 * $n_{1,\downarrow}$ - 1 * $n_{2,\uparrow}$ - 1 * $n_{2,\downarrow}$
    #                = - 1 * $c^\dagger_0$ $c_0$ - 1 * $c^\dagger_1$ $c_1$ - 1 * $c^\dagger_2$ $c_2$ - 1 * $c^\dagger_3$ $c_3$

    #ハバードモデルハミルトニアンを定義し、対角化まで行う
    n_qubit = 2*nsite #量子ビットの数
    U = 10.0
    t = -1.0

    #ハーフフィリングを仮定(電子数 = サイトの数)するとケミカルポテンシャルはμ=U/2.マイナスがつくので以下。
    μ = U/2

    ham = FermionOperator()
    #斥力項
    for i in 1:nsite
        #up_index,down_indexの定義は、QC_materialを参照。
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^ $(down)^ $(up) $(down)", -U) #左側に生成演算子。右際に消滅演算子をもっていく過程で半交換関係が1回でマイナスをつける。
    end
    #ホッピング項
    for i in 1:nsite-1
        ham += FermionOperator("$(up_index(i+1))^ $(up_index(i))", t)
        ham += FermionOperator("$(up_index(i))^ $(up_index(i+1))", t)
        ham += FermionOperator("$(down_index(i+1))^ $(down_index(i))", t)
        ham += FermionOperator("$(down_index(i))^ $(down_index(i+1))", t)
    end

    #ケミカルポテンシャルの項
    for i in 1:nsite
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^  $(up) ", -μ)
        ham += FermionOperator("$(down)^ $(down)", -μ)
    end

    Hamiltonian(ham, nsite)
end

"""
Evaluate energy
"""
function eval_energy(circuit, ham_jw, theta_list, n_qubit)
    state = create_hf_state(n_qubit, n_electron) #|0000> を準備
    update_circuit_param!(circuit, theta_list) #量子回路にパラメータをセット
    update_quantum_state!(circuit, state) #量子回路を状態に作用
    get_expectation_value(ham_jw, state) #ハミルトニアンの期待値
end

function construct_circuit(hamiltonian)
    ham = hamiltonian.ham
    n_qubit = 2*hamiltonian.n_sites

    @assert mod(n_electron, 2) == 0
    @assert n_electron <= hamiltonian.n_sites

    #JW変換
    ham_jw = jordan_wigner(ham)

    # Prepare a circuit
    circuit = uccgsd(n_qubit, orbital_rot=true)
    circuit, ham_jw
end

function solve(hamiltonian, n_electron;theta_init=nothing, comm=nothing)
    if comm === nothing
        rank = 0
    else
        rank = MPI.Comm_rank(comm)
    end

    ham = hamiltonian.ham
    n_qubit = 2*hamiltonian.n_sites

    println(ham)　#ハバードモデルハミルトニアンの表示。
    @assert mod(n_electron, 2) == 0
    @assert n_electron <= hamiltonian.n_sites

    # Compute exact ground-state energy
    sparse_mat = get_number_preserving_sparse_operator(ham, n_qubit, n_electron)　#行列の取得
    enes_ed = eigvals(sparse_mat.toarray())　#対角化を行う
    EigVal_min = minimum(enes_ed)
    if rank == 0
        println("Exact ground-state energy = $(EigVal_min)")
    end

    # Construct circuit
    circuit, ham_jw = construct_circuit(hamiltonian)
    if rank == 0
        println("Number of Qubits:", n_qubit)
        println("Number of Electrons:", n_electron)
    end

    # Define a cost function
    function cost(theta_list)
        eval_energy(circuit, ham_jw, theta_list, n_qubit)
    end

    # Define the gradient of the cost function
    function grad_cost(theta_list)
        t1 = time_ns()
        if comm === nothing
            first_idx, size = 1, length(theta_list)
        else
            first_idx, size = distribute(length(theta_list), MPI.Comm_size(comm), MPI.Comm_rank(comm))
        end
        last_idx = first_idx + size - 1
        res = numerical_grad(cost, theta_list, first_idx=first_idx, last_idx=last_idx)
        if comm !== nothing
            res = MPI.Allreduce(res, MPI.SUM, comm)
        end
        t2 = time_ns()
        if rank == 0
            println("g: ", (t2-t1)*1e-9)
        end
        res
    end

    if theta_init === nothing
        theta_init = rand(size(circuit.theta_offsets)[1])
        if comm !== nothing
            # Make sure all processes use the same initial values
            MPI.Bcast!(theta_init, 0, comm)
        end
    end
    cost_history = Float64[] #コスト関数の箱
    init_theta_list = theta_init
    push!(cost_history, cost(init_theta_list))

    method = "BFGS"
    options = Dict("disp" => true, "maxiter" => 200, "gtol" => 1e-5)
    function callback(x)
        push!(cost_history, cost(x))
        if rank == 0
            println("iter ", length(cost_history), " ", cost_history[end])
        end
    end
    opt = scipy_opt.minimize(cost, init_theta_list, method=method, callback=callback, jac=grad_cost,
        options=options)

    return cost_history, circuit, EigVal_min, opt
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
module samplet3
function samplet3mesh()
    xyz = [0.0 0.0
0.5 0.0
1.0 0.0
0.0 0.3333333333333333
0.5 0.3333333333333333
1.0 0.3333333333333333
0.0 0.6666666666666666
0.5 0.6666666666666666
1.0 0.6666666666666666
0.0 1.0
0.5 1.0
1.0 1.0]
    c = [1 2 5
1 5 4
4 5 8
4 8 7
7 8 11
7 11 10
2 3 6
2 6 5
5 6 9
5 9 8
8 9 12
8 12 11
]
    return xyz, c
end
end


text:
#=
using RoguePkg

julia> using VariantVisualization
INFO: Loading HttpServer methods...
WARNING: Method definition ==(Base.Nullable{S}, Base.Nullable{T}) in module Base at nullable.jl:238 overwritten in module NullableArrays at /Users/George/.julia/test/v0.6/NullableArrays/src/operators.jl:99.

julia> Pkg.test(pkg_for"VariantVisualization")

=#

@testset "VCFUtils" begin

vcf_filename = "test_files/test_4X_191.vcf"
#vcf_filename_with_chr = "test_files/test_with_chr.vcf"

reader = VCF.Reader(open(vcf_filename, "r"))
#reader_with_chr = VCF.Reader(open(vcf_filename_with_chr, "r"))
sample_names = get_sample_names(reader)

dp_limit = 50

@testset "clean_column1!" begin
    df = Matrix(["X" 1 2; "Y" 2 3; 2 4 1])
    clean_column1!(df)
    @test df[1,1] == "23"
    @test df[2,1] == "24"
    @test df[3,1] == 2
end

#functions for variant filters

@testset "io_genomic_range_vcf_filter" begin
sub = io_genomic_range_vcf_filter("chr4:0-400000000",vcf_filename)
@test typeof(sub) == Array{Any,1}
@test size(sub,1) == 1012
#println("io_genomic_range_vcf_filter type is $(typeof(sub))")
#println("io_genomic_range_vcf_filter size is $(size(sub,1))")
end

@testset "filters_with_siglist" begin

    @testset "load_siglist" begin
    sig_list=load_siglist("test_files/positions_list.csv")

    @testset "pass_genomic_range_siglist_filter" begin
    sig_list=load_siglist("test_files/positions_list.csv")
    sub = pass_genomic_range_siglist_filter(vcf_filename,sig_list,"chr4:0-5000000000")

    @test (typeof(sub[1])) == GeneticVariation.VCF.Record
    @test (length(sub)) ==  5

    end

            @testset "io_sig_list_vcf_filter" begin

            sub = VariantVisualization.io_sig_list_vcf_filter(sig_list,vcf_filename)
            @test (typeof(sub[1])) == GeneticVariation.VCF.Record
            @test (length(sub)) ==  11
            end

            @testset "pass_siglist_filter" begin
            sub = pass_siglist_filter(vcf_filename, sig_list)
            @test (typeof(sub[1])) == GeneticVariation.VCF.Record
            @test (length(sub)) ==  10
            end

            @testset "genomic_range_siglist_filter" begin
            sub = genomic_range_siglist_filter(vcf_filename,sig_list,"chr4:0-400000000")
            @test (typeof(sub[1])) == GeneticVariation.VCF.Record
            @test (length(sub)) ==  5
            end

    end
end

@testset "io_pass_filter" begin
    sub = io_pass_filter(vcf_filename)
    @test (typeof(sub[1])) == GeneticVariation.VCF.Record
    @test (length(sub)) ==  1164
end

@testset "pass_genomic_range_filter" begin
    reader = VCF.Reader(open(vcf_filename, "r"))
    sub = pass_genomic_range_filter(reader,"chr4:0-400000000",vcf_filename)
    @test (typeof(sub[1])) == GeneticVariation.VCF.Record
    @test (length(sub)) ==  856
end

#functions for converting vcf record array to numerical array
@testset "combined_all_genotype_array_functions" begin

sub = io_pass_filter(vcf_filename)
number_rows = size(sub,1)

gt_num_array,gt_chromosome_labels=combined_all_genotype_array_functions(sub)
#println("combined_all_genotype_array_functions gt array is type: $(typeof(gt_num_array))")
#println("combined_all_genotype_array_functions gt array is length: $(length(gt_num_array))")
#println("combined_all_genotype_array_functions gt_chromosome_labels is typeof: $(typeof(gt_chromosome_labels))")
#println("combined_all_genotype_array_functions gt_chromosome_labels is length: $(length(gt_chromosome_labels))")
@test typeof(gt_num_array) == Array{Int64,2}
@test length(gt_num_array) == 222324
@test typeof(gt_chromosome_labels) == Array{Any,2}
@test length(gt_chromosome_labels) == 2328

    @testset "chromosome_label_generator" begin
    chrom_label_info = VariantVisualization.chromosome_label_generator(gt_chromosome_labels[:,1])
    #println("chromosome_label_generator chrom_label_info is type $(typeof(chrom_label_info))")
    #println("chromosome_label_generator chrom_label_info is length $(length(chrom_label_info))")
    @test typeof(chrom_label_info) == Tuple{Array{String,1},Array{Int64,1},String}
    @test size(chrom_label_info,1) == 3
    end

    @testset "generate_chromosome_positions_for_hover_labels" begin
    chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(gt_chromosome_labels)
    #println("generate_chromosome_positions_for_hover_labels chr_pos_tuple_list is type $(typeof(chr_pos_tuple_list))")
    #println("generate_chromosome_positions_for_hover_labels chr_pos_tuple_list is length $(length(chr_pos_tuple_list))")
    @test typeof(chr_pos_tuple_list) == Array{String,1}
    @test size(chr_pos_tuple_list,1) == 1164
    end

    @testset "generate_genotype_array" begin
    sub = io_pass_filter(vcf_filename)
    genotype_array=generate_genotype_array(sub,"GT")

    #println("generate_genotype_array is $(typeof(genotype_array))")
    #println("generate_genotype_array is $(size(genotype_array,1))")
    @test typeof(genotype_array) == Array{String,2}
    @test size(genotype_array,1) == 1164

    @testset "define_geno_dict" begin
    geno_dict = define_geno_dict()
    #println("define_geno_dict is type is $(typeof(geno_dict))")
    #println("define_geno_dict is length is $(length(geno_dict))")
    @test typeof(geno_dict) == Dict{Any,Any}
    @test length(geno_dict) == 101

    @testset "translate_genotype_to_num_array" begin
    gt_num_array,gt_chromosome_labels=translate_genotype_to_num_array(genotype_array,geno_dict)
    #println("translate_genotype_to_num_array gt_num_array type is $(typeof(gt_num_array))")
    #println("translate_genotype_to_num_array gt_num_array length is $(length(gt_num_array))")
    #println("translate_genotype_to_num_array gt_chromosome_labels typeof is $(typeof(gt_chromosome_labels))")
    #println("translate_genotype_to_num_array gt_chromosome_labels length is $(length(gt_chromosome_labels))")
    @test typeof(gt_num_array) == Array{Int64,2}
    @test length(gt_num_array) == 222324
    @test typeof(gt_chromosome_labels) == Array{String,2}
    @test length(gt_chromosome_labels) == 2328
    end
    end
    end

end

@testset "combined_all_read_depth_array_functions" begin

sub = io_pass_filter(vcf_filename)
dp_num_array,dp_chromosome_labels=combined_all_read_depth_array_functions(sub)

#println("combined_all_read_depth_array_functions dp_num_array type is $(typeof(dp_num_array))")
#println("combined_all_read_depth_array_functions dp_num_array length is $(length(dp_num_array))")
#println("combined_all_read_depth_array_functions gt_chromosome_labels typeof is $(typeof(dp_chromosome_labels))")
#println("combined_all_read_depth_array_functions gt_chromosome_labels length is $(length(dp_chromosome_labels))")
@test typeof(dp_num_array) == Array{Int64,2}
@test length(dp_num_array) == 222324
@test typeof(dp_chromosome_labels) == Array{Any,2}
@test length(dp_chromosome_labels) == 2328

    @testset "get_sample_names" begin
    reader = VCF.Reader(open(vcf_filename, "r"))
    sample_names=get_sample_names(reader)
    #println("get_sample_names")
    #println("get_sample_names sample_names type is $(typeof(sample_names))")
    #println("get_sample_names sample_names length is $(length(sample_names))")
    @test typeof(sample_names) == Array{Symbol,2}
    @test length(sample_names) == 191

            @testset "read_depth_threshhold" begin
            dp_num_array=read_depth_threshhold(dp_num_array)
            #println("read_depth_threshhold dp_num_array is type $(typeof(dp_num_array))")
            #println("read_depth_threshhold dp_num_array is size $(size(dp_num_array,1))")
            @test typeof(dp_num_array) == Array{Int64,2}
            @test size(dp_num_array,1) == 1164
            end

        @testset "avg_dp_samples" begin
        avg_sample_list=avg_dp_samples(dp_num_array)
        #println("avg_dp_samples avg_sample_list type is $(typeof(avg_sample_list))")
        #println("avg_dp_samples avg_sample_list length is $(length(avg_sample_list))")
        @test typeof(avg_sample_list) == Array{Float64,1}
        @test size(avg_sample_list,1) == 191

        @testset "list_sample_names_low_dp" begin
        list=list_sample_names_low_dp(avg_sample_list,sample_names)
        #println("list_sample_names_low_dp list type is $(typeof(list))")
        #println("list_sample_names_low_dp list length is $(length(list))")
        @test typeof(list) == Array{String,1}
        @test size(list,1) == 6
        end

    end

    @testset "avg_dp_variant" begin
    avg_variant_list=avg_dp_variant(dp_num_array)
    #println("avg_dp_variant avg_variant_list type is $(typeof(avg_variant_list))")
    #println("avg_dp_variant avg_variant_list length is $(length(avg_variant_list))")
    @test typeof(avg_variant_list) == Array{Float64,1}
    @test size(avg_variant_list,1) == 1164

        @testset  "list_variant_positions_low_dp" begin
        list=list_variant_positions_low_dp(avg_variant_list,dp_chromosome_labels)
        #println("list_variant_positions_low_dp list type is $(typeof(list))")
        #println("list_variant_positions_low_dp list length is $(length(list))")
        @test typeof(list) == Array{Tuple{Any,Int64},1}
        @test size(list,1) == 33

        end

    end

    @testset "sortcols_by_phenotype_matrix" begin
    vcf,group_label_pack=sortcols_by_phenotype_matrix("test_files/sample_metadata_matrix.csv","case,control", dp_num_array, sample_names)
    #println("sortcols_by_phenotype_matrix vcf type is $(typeof(vcf))")
    #println("sortcols_by_phenotype_matrix vcf size is $(size(vcf,1))")
    #println("sortcols_by_phenotype_matrix group_label_pack type is $(typeof(group_label_pack))")
    #println("sortcols_by_phenotype_matrix group_label_pack size is $(size(group_label_pack,1))")
    @test typeof(vcf) == Array{Int64,2}
    @test size(vcf,1) == 1164
    @test typeof(group_label_pack) == Array{Any,1}
    @test size(group_label_pack,1) == 5

#=
        @testset "find_group_label_indices" begin
        pheno = DelimitedFiles.readdlm("test_files/fixed_pheno_matrix_test.csv", ',')
        row_to_sort_by = find(x -> x == "control,case", pheno)
        row_to_sort_by = row_to_sort_by[1]
        group_label_pack=find_group_label_indices(pheno,"control,case",row_to_sort_by)
        #println("find_group_label_indices group_label_pack type is $(typeof(group_label_pack))")
        #println("find_group_label_indices group_label_pack length is $(length(group_label_pack))")
        @test typeof(group_label_pack) == Array{Any,1}
        @test size(group_label_pack,1) == 5
        end

        =#

        @testset "select_columns" begin
        dp_num_array=select_columns("test_files/select_samples_list.txt", dp_num_array, sample_names)
        #println("select_columns dp_num_array type is $(typeof(dp_num_array))")
        #println("select_columns dp_num_array size is $(size(dp_num_array,1))")
        @test typeof(dp_num_array) == Tuple{Array{Int64,2},Array{Any,2}}
        @test size(dp_num_array,1) == 2
        end

    end

    end
    end

end


text:
module IsOdd

using PyCall
@pyimport pyiseven

function odd(x)
  if x == 1
    true
  elseif x == 0
    false
  else
    pyiseven.even(x-1)
  end
end

end

max_stars_repo_name: astrieanna/madpug

max_issues_repo_name: astrieanna/madpug

max_forks_repo_name: astrieanna/madpug


text:
# Put in a separate page so it can be used by SciMLDocs.jl

pages=[
    "Home" => "index.md",
    "matrix_exponentials.md",
    "expv.md",
    "arnoldi.md"
]

max_stars_repo_name: JuliaDiffEq/ExponentialUtilities

max_issues_repo_name: JuliaDiffEq/ExponentialUtilities

max_forks_repo_name: JuliaDiffEq/ExponentialUtilities


text:
using SimpleSolvers
using Test
using GeometricIntegrators
using GeometricIntegrators.Integrators.VPRK
using GeometricIntegrators.Utils
using GeometricProblems.LotkaVolterra4d
using GeometricProblems.LotkaVolterra4d: Δt, nt, reference_solution

SimpleSolvers.set_config(:nls_atol, 8eps())
SimpleSolvers.set_config(:nls_rtol, 2eps())

@testset "$(rpad("Lotka-Volterra 4D",80))" begin
    ode  = lotka_volterra_4d_ode()
    iode = lotka_volterra_4d_iode()
    idae = lotka_volterra_4d_idae()

    int = Integrator(ode, TableauGauss(2), Δt)
    sol = integrate(ode, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 8E-4

    int = IntegratorVPRKpMidpoint(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 2E-3

    int = IntegratorVPRKpSymmetric(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 8E-4

    int = Integrator(idae, TableauVSPARKGLRKpMidpoint(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 2E-3

    int = Integrator(idae, TableauVSPARKGLRKpSymmetric(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 8E-4
end


text:
using DifferentialEquations
using SimpleDiffEq
using DiffEqCallbacks
using Random
using Distributions
using Plots

function sir_ode!(du,u,p,t)
    (S,I,R,C) = u
    (β,c,γ) = p
    N = S+I+R
    infection = β*c*I/N*S
    recovery = γ*I
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
        du[4] = infection
    end
    nothing
end;

tmax = 40.0
δt = 1.0
tspan = (0.0,tmax)
obstimes = 1.0:δt:tmax;
u0 = [990.0,10.0,0.0,0.0]; # S,I.R,C
p = [0.05,10.0,0.25]; # β,c,γ

prob_ode = ODEProblem(sir_ode!,u0,tspan,p)
sol_ode_cumulative = solve(prob_ode,Tsit5(),saveat=δt);

out = Array(sol_ode_cumulative)
C = out[4,:];

X = C[2:end] .- C[1:(end-1)];

Random.seed!(1234);

Y = rand.(Poisson.(X));

bar(obstimes,Y)
plot!(obstimes,X)

S = out[1,:]
Cpred = 990.0 .- S
Cdiff = Cpred .- C
plot(obstimes,Cdiff[2:end])

affect!(integrator) = integrator.u[4] = 0.0
cb_zero = PresetTimeCallback(obstimes,affect!);

sol_ode_cb = solve(prob_ode,Tsit5(),saveat=δt,callback=cb_zero);

X_cb = sol_ode_cb(obstimes)[4,:];

Random.seed!(1234);

Y_cb = rand.(Poisson.(X_cb));

X_diff_cb = X_cb .- X
plot(obstimes,X_diff_cb)

Y_diff_cb = Y_cb .- Y
plot(obstimes,Y_diff_cb)

function sir_dde!(du,u,h,p,t)
    (S,I,R,C) = u
    (β,c,γ) = p
    N = S+I+R
    infection = β*c*I/N*S
    recovery = γ*I
    e = oneunit(t)
    history = h(p, t-e)*inv(e)
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
        du[4] = infection - history[4]
    end
    nothing
end;

function sir_history(p, t; idxs = 5)
    zero(t)
end;

prob_dde = DDEProblem(DDEFunction(sir_dde!),
        u0,
        sir_history,
        tspan,
        p;
        constant_lags = [1.0]);

sol_dde = solve(prob_dde,MethodOfSteps(Tsit5()));

X_dde = sol_dde(obstimes)[4,:];

Random.seed!(1234)
Y_dde = rand.(Poisson.(X_dde));

X_diff_dde = X_dde .- X
plot(X_diff_dde)

Y_diff_dde = Y_dde .- Y
plot(obstimes, Y_diff_dde)


text:
module TensorFlow

# N.B. This is only a stub used to run tests in the generated module,
#      and does not get installed to TensorFlow/src.

include("CoreTypes.jl")
using .CoreTypes

include("API.jl")
using .API

end


text:
#----------------------------------------------------
function generatelatentZ(; S = S, D = D, seed=1)
#----------------------------------------------------

    rg = MersenneTwister(seed)

    Z  = [randn(rg, D) for s in 1:S]

    return Z

end

#----------------------------------------------------
function eigendecomposition(Σ)
#----------------------------------------------------

    # decompose covariance
    EigenDecomposition = eigen(Σ)

    # get root of eigenvalues
    Esqrt = sqrt.(max.(EigenDecomposition.values, 1e-8))

    # get eigenvectors that stay fixed throughout
    V     = EigenDecomposition.vectors

    return V, Esqrt

end


text:
using RemoteFiles

stats_base = "https://julialang-logs.s3.amazonaws.com/public_outputs/current/"
ext = ".csv.gz"

all_rollups = [
    "client_types",
    "client_types_by_region",
    "client_types_by_date",
    "client_types_by_region_by_date",
    "julia_systems",
    "julia_systems_by_region",
    "julia_systems_by_date",
    "julia_systems_by_region_by_date",
    "julia_versions",
    "julia_versions_by_region",
    "julia_versions_by_date",
    "julia_versions_by_region_by_date",
    "resource_types",
    "resource_types_by_region",
    "resource_types_by_date",
    "resource_types_by_region_by_date",
    "package_requests",
    "package_requests_by_region",
    "package_requests_by_date",
    "package_requests_by_region_by_date",
]


text:
#=
  Combograms in Julia.

  List all the possible words that can created given a source word,
  with some options:
  - as many characters as in the word, or as many as wanted
  - plain subwords
  - anagrams
  - minimum length
  - must contain certain characters

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

=#

#
# Return a Dict of the elements and their occurences
# in the collection a.
#
function make_hash(a)
    d = Dict{Char,Int}()
    for e in collect(a)
        if !haskey(d,e)
            d[e] = 0
        end
        d[e] += 1
    end
    return d
end

#
# combogram(word,
#           lang,         # :eng, :swe
#           type,         # :as_many_as_in_source,:as_many_as_wanted, :plain_subwords, :anagram
#           must_contain, # character that must be in the word
#           min_length,   # minimum length of the word
#           )
#
function combogram(word,
                   lang=:eng,
                   type=:as_many_as_in_source,
                   must_contain="",
                   min_length=0)

    word = lowercase(word)
    println("word:$word lang:$lang type:$type must_contain:$must_contain min_length:$min_length")

    wordlist = "eng_dict.txt" # :eng
    if lang == :swe
        # Swedish words
        wordlist = "sv_spelling_org_utf8.txt"
    end

    if must_contain != ""
        must_contain = collect(must_contain)
    end

    words = readlines(wordlist)
    println("Number of words: $(length(words))")

    d = make_hash(word) # Characters in source word

    res = []
    cc = 0
    for w in words

        w_len = length(w)
        if w_len == 0
            continue
        end

        if min_length > 0 && w_len < min_length
            continue
        end

        # Hash for this word
        wd = make_hash(w)

        # Anagrams
        if type == :anagram
            if wd != d
                continue
            end
        end

        # Plain subword
        if type == :plain_subwords
            if !occursin(w,word)
                continue
            end
        end

        good = true
        # Check all characters
        for c in keys(wd)

            # Char is not in source word
            if !haskey(d,c)
                good = false
                break
            end
            if !good continue end

            # Has not a mandatory char
            if must_contain != ""
                for mc in must_contain
                    if !haskey(wd,mc)
                        good = false
                        break
                    end
                end
            end
            if !good continue end

            if type == :as_many_as_in_source
                # Too many occurrences of this character
                if wd[c] > d[c]
                    good = false
                    break
                end
            end
            if !good continue end

        end

        if good
            push!(res,w)
        end

    end
    return res
end

# 747 Swedish words
# res = combogram("programming",:swe,:as_many_as_wanted,"",0)

# 615 English words
# res = combogram("programming",:eng,:as_many_as_wanted,"",0)

# Plain subwords
# res = combogram("kjellerstrand",:swe,:plain_subwords,"",0)

# Anagrams
# res = combogram("glare",:eng,:anagram,"",0)

# min length
# Any["darkeners", "slanderer", "terrellas", "treadlers"]
# res = combogram("kjellerstrand",:eng,:as_many_as_in_source,"",9)

len = length(res)
println(res)
println("found $len words")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# TODO: is there a way to have build_executable import the cli.jl, but also
# have the `main()` at the bottom? it errors now when it tries to run it....
# This feels kinds dumb now
include("./cli.jl")

main()


text:
# only need to be run once to install packages
#Pkg.clone("https://github.com/JuliaData/SplitApplyCombine.jl.git")
#Pkg.clone("https://github.com/xiaodaigh/FastGroupBy.jl.git")

using Revise, DataBench
using FastGroupBy

const M=100_000_000; const K=100
srand(1)
svec1 = rand(["i"*dec(k,7) for k in 1:M÷K], M)
@time radixsort!(svec1)
issorted(svec1)

using PooledArrays

const N = 100_000_000
# const N = Int(2^31-1) # 368 seconds to run
const K = 100

using Base.Threads
nthreads()

srand(1)
# generate string ids
function randstrarray1(pool, N)
    K = length(pool)
    PooledArray(PooledArrays.RefArray(rand(1:K, N)), pool)
end
const pool1 = [@sprintf "id%010d" k for k in 1:(N/K)]
const id3 = randstrarray1(pool1, N)
v1 =  rand(Int32(1):Int32(5), N)

# treat it as Pooledarray
@time sumby(id3, v1)
@time fastby(id3, v1, sum)

# treat by as strings and use dictionary method; REALLY SLOW
const id3_str = rand(pool1, N)
@time sumby(id3_str, v1)

@time Int.(getindex.(id3_str,1 ))

@time all(isascii.(id3_str))

@time sort(id3_str)

srand(1)
pool_str = rand([@sprintf "id%010d" k for k in 1:(N/K)], N)
id = rand(1:K, N)
valvec = rand(N)
@time all(isascii.(pool_str))
@time svec = sizeof.(pool_str)
by_vec = pool_str

using FastGroupBy
@time sumby_dict(zip(pool_str, id), valvec)

# fast group by for unicode strings
function fastby!(fn::Function, byvec::Vector{T}, valvec::Vector{S}, skip_sizeof_grouping = false, ascii_only = false) where {T<:AbstractString,S}
    l = length(byvec)

    # firstly sort the string by size
    if skip_sizeof_grouping
        ##
    else
        sizevec = sizeof.(svec)
        # typically the range of sizes for
        minsize, maxsize = extrema(svec)
        if  minsize != maxsize
            # if there is only one size then ignore
            indices = collect(1:l)
            # grouptwo!(svec, indices)
        else
        end
    end
end

function isgrouped(x, hashx)
    x1 = x[1]
    hashx1 = hashx[1]
    for i = 2:length(x)
        @inbounds if x1 != x[i]
            if hashx1 == hashx[i]
                return false
            else
                x1 = x[i]
                hashx1 = hashx[i]
            end
        end
    end
    return true
end

function roughhash(s::AbstractString)
    pp = pointer(s)
    sz = sizeof(s)
    hh = zero(UInt64) | Base.pointerref(pp, 1, 1)
    for j = 2:min(sz, 8)
        hh = (hh << 8) | Base.pointerref(pp, j, 1)
    end
    hh
end
@benchmark roughhash(x)

function radixgroup(fn::Function, svec::Vector{T}, valvec::Vector{S}) where {T <: AbstractString, S}
    a = zeros(UInt128, length(svec))
    @time for (i, s) in enumerate(svec)
        pp = pointer(s)
        hh = zero(UInt128) ⊻ Base.pointerref(pp, 1, 1)
        sz = sizeof(s)
        for j = 2:min(16,sz)
            hh = (hh << 8) ⊻ Base.pointerref(pp, j, 1)
            # hh += (hh << 10);
            # hh ^= (hh >> 6);
        end
        @inbounds a[i] = hh
    end
    a
    # @time grouptwo!(a, valvec)
    #
    # res = Dict{T, S}()
    # l = length(svec)
    #
    # j = 1
    # lastby = svec[1]
    # @time for i = 2:l
    #     @inbounds byval = svec[i]
    #     if byval != lastby
    #         viewvalvec = @view valvec[j:i-1]
    #         @inbounds res[lastby] = fn(viewvalvec)
    #         j = i
    #         @inbounds lastby = svec[i]
    #     end
    # end
    #
    # viewvalvec = @view valvec[j:l]
    # @inbounds res[svec[l]] = fn(viewvalvec)
    # return res
end

@time permi = radixgroup(sum, svec, valvec)
using StatsBase
countmap(permi)

hh = zero(UInt)
j=1
pp = pointer(svec[1])
bp = Base.pointerref(pp, j, 1)
hh = xor(hh << 8, bp)
j = j + 1

@code_warntype radixgroup!(svec)

function fastby(fn::Function, byvec::Vector{T}, valvec::Vector) where T <: AbstractString
    permi = radixgroup!(byvec)
    fastby_contiguous(fn, byvec[permi], valvec[permi])
end

function fastby_contiguous(fn::Function, byvec::Vector{T}, valvec::Vector{S}) where {T <: AbstractString, S}
    res = Dict{T, S}()
    l = length(byvec)

    j = 1
    lastby = byvec[1]
    for i = 2:l
        @inbounds byval = byvec[i]
        if byval != lastby
            viewvalvec = @view valvec[j:i-1]
            @inbounds res[lastby] = fn(viewvalvec)
            j = i
            @inbounds lastby = byvec[i]
        end
    end

    viewvalvec = @view valvec[j:l]
    @inbounds res[byvec[l]] = fn(viewvalvec)
    return res
end

@time fastby(sum, idstr, valvec)

using StatsBase
@time countmap(idstr)


text:
# Partitioning (localization of global objects) and assembly (global assembly of local objects)

partition(c::AbstractVector, Nc, Nr, r) = c[1 + (r-1) * Nc : 1 + Nc * r]
partition(c::Colon,          Nc, Nr, r) = Colon()
partition(c::Tuple,          Nc, Nr, r) = (c[1] + (r-1) * (c[2] - c[1]) / Nr,    c[1] + r * (c[2] - c[1]) / Nr)

function partition(c::UnitRange, Nc, Nr, r)
    g = (first(c), last(r))
    ℓ = partition(g, Nc, Nr, r)
    return UnitRange(ℓ[1], ℓ[2])
end

assemble(c::Tuple, Nc, Nr, r, arch) = (c[2] - r * (c[2] - c[1]), c[2] - (r - Nr) * (c[2] - c[1]))

"""
    assemble(c::AbstractVector, Nc, Nr, r, arch)

Build a linear global coordinate vector given a local coordinate vector `c_local`
a local number of elements `Nc`, number of ranks `Nr`, rank `r`,
and `arch`itecture.
"""
function assemble(c_local::AbstractVector, Nc, Nr, r, arch)
    c_global = zeros(eltype(c_local), Nc*Nr+1)
    c_global[1 + (r-1) * Nc : Nc * r] .= c[1:end-1]
    r == Nr && (c_global[end] = c[end])

    MPI.Allreduce!(c_global, +, arch.communicator)

    return c_global
end


text:
"""
    break_stick_ibp(ν)

Break stick (IBP).
"""
function break_stick_ibp(ν)
    Kmax = size(ν, 1)
    p = []
    push!(p, ν[1,:])
    for k = 2:Kmax
        push!(p, p[k-1] .* ν[k,:])
    end
    p = hcat(p...)'
    return p
end

"""
    break_logstick_ibp(logν)

Break stick (IBP) in logarithmic space.
"""
function break_logstick_ibp(logν)
    Kmax = size(logν, 1)
    logp = []
    push!(logp, logν[1,:])
    for k = 2:Kmax
        push!(logp, logp[k-1] .+ logν[k,:])
    end
    logp = hcat(logp...)'
    return logp
end


text:
# TODO: Connect with GitLink config
const _MIN_LOOP_FREC = 3.0
const _MAX_LOOP_FREC = 60.0

const _LOOP_FREC_KEY = :loop_frec_key
loop_frec(gl::GitLink) = clamp(
    get!(gl, _LOOP_FREC_KEY, _MIN_LOOP_FREC),
    _MIN_LOOP_FREC, _MAX_LOOP_FREC
)

loop_frec!(gl::GitLink, val) = set!(gl, _LOOP_FREC_KEY, clamp(float(val), _MIN_LOOP_FREC, _MAX_LOOP_FREC))
add_loop_frec!(gl::GitLink, val) = loop_frec!(gl, loop_frec(gl) + val)


text:
fs()

mkpath(F.path(:literate))

@testset "Literate-0" begin
    @test_throws ErrorException literate_folder("foo/")
    litpath = literate_folder("_literate/")
    @test litpath == literate_folder(F.path(:literate))
end

# @testset "Literate-a" begin
#     # Post processing: numbering of julia blocks
#     s = raw"""
#         A
#
#         ```julia
#         B
#         ```
#
#         C
#
#         ```julia
#         D
#         ```
#         """
#     @test F.literate_post_process(s) == """
#
#         A
#
#         ```julia:ex1
#         B
#         ```
#
#         C
#
#         ```julia:ex2
#         D
#         ```
#         """
# end

@testset "Literate-b" begin
    # Literate to Franklin
    s = raw"""
        # # Rational numbers
        #
        # In julia rational numbers can be constructed with the `//` operator.
        # Lets define two rational numbers, `x` and `y`:

        ## Define variable x and y
        x = 1//3
        y = 2//5

        # When adding `x` and `y` together we obtain a new rational number:

        z = x + y
        """
    path = joinpath(F.path(:literate), "tutorial.jl")
    write(path, s)
    opath, = F.literate_to_franklin("/_literate/tutorial")
    @test endswith(opath, joinpath(F.PATHS[:site], "assets", "literate", "tutorial.md"))
    out = read(opath, String)
    @test out == """

        # Rational numbers

        In julia rational numbers can be constructed with the `//` operator.
        Lets define two rational numbers, `x` and `y`:

        ```julia:ex1
        # Define variable x and y
        x = 1//3
        y = 2//5
        ```

        When adding `x` and `y` together we obtain a new rational number:

        ```julia:ex2
        z = x + y
        ```

        """

    # Use of `\literate` command
    h = raw"""
        @def hascode = true
        @def showall = true
        @def reeval = true

        \literate{/_literate/tutorial.jl}
        """ |> fd2html_td
    @test isapproxstr(h, """

Rational numbers

In julia rational numbers can be constructed with the // operator. Lets define two rational numbers, x and y:

$(F.htmlesc(raw"""
            # Define variable x and y
            x = 1//3
            y = 2//5
            """))

2//5

When adding x and y together we obtain a new rational number:

$(F.htmlesc(raw"""
            z = x + y
            """))

11//15

        """)

    # issue 592
    # Literate to Franklin
    s = raw"""
        # # Rational numbers
        # ```julia
        # const a = 1
        # ```
        a = 5
        """
    path = joinpath(F.path(:literate), "tutorial.jl")
    write(path, s)
    opath, = F.literate_to_franklin("/_literate/tutorial")
    @test endswith(opath, joinpath(F.PATHS[:site], "assets", "literate", "tutorial.md"))
    out = read(opath, String)
    @test out // """

        # Rational numbers
        ```julia
        const a = 1
        ```

        ```julia:ex1
        a = 5
        ```
        """

    # Use of `\literate` command
    h = raw"""
        @def hascode = true
        @def showall = true
        @def reeval = true

        \literate{/_literate/tutorial.jl}
        """ |> fd2html_td
    @test h // """

Rational numbers

$(F.htmlesc(raw"""const a = 1"""))

$(F.htmlesc(raw"""a = 5"""))

5

        """
end

@testset "Literate-c" begin
    # s = raw"""
    #     \literate{/foo/bar/foo}
    #     """
    # @test_throws F.LiterateRelativePathError (s |> fd2html_td)
    s = raw"""
        \literate{/foo}
        """

    global r = ""; s = @capture_out begin
        global r
        r = s |> fd2html_td
    end
    @test r == """
// Literate file matching '/foo' not found. //
\n"""
    @test occursin("File not found when", s)
end


text:
# IIF #485 --

# using Revise

using Test
using Logging
using Statistics
using DistributedFactorGraphs
using IncrementalInference
using TensorCast

@testset "test basic three variable graph with prior" begin

VAR1 = :a
VAR2 = :b
VAR3 = :c

# logger = SimpleLogger(stdout, Logging.Debug)
# global_logger(logger)
dfg = initfg() #LocalDFG{SolverParams}(solverParams=SolverParams())
# Add some nodes.
v1 = addVariable!(dfg, VAR1, ContinuousScalar, labels = [:POSE])
v2 = addVariable!(dfg, VAR2, ContinuousScalar, labels = [:POSE])
v3 = addVariable!(dfg, VAR3, ContinuousScalar, labels = [:LANDMARK])
f1 = addFactor!(dfg, [VAR1; VAR2], LinearRelative(Normal(50.0,2.0)) )
f2 = addFactor!(dfg, [VAR2; VAR3], LinearRelative(Normal(50.0,2.0)) )

addFactor!(dfg, [VAR1], Prior(Normal()))

# drawGraph(dfg, show=true)

# tree = buildTreeReset!(dfg)
# # drawTree(tree, show=true)
#
# getCliqFactors(tree, VAR3)
# getCliqFactors(tree, VAR1)

initAll!(dfg)

# cliq= getClique(tree, VAR3)
# getCliqueData(cliq)
#
# cliq= getClique(tree, VAR1)
# getCliqueData(cliq)

getSolverParams(dfg).limititers = 50
# getSolverParams(dfg).drawtree = true
# getSolverParams(dfg).showtree = true
# getSolverParams(dfg).dbg = true
## getSolverParams(dfg).async = true

tree = solveTree!(dfg) #, recordcliqs=ls(dfg))

pts_ = getBelief(dfg, :c) |> getPoints
TensorCast.@cast pts[i,j] := pts_[j][i]

@test 70 < Statistics.mean(pts) < 130

# #
# using Gadfly, Cairo, Fontconfig
# drawTree(tree, show=true, imgs=true)

end

#


text:
export rotate, rotate!

"""
    rotate(arr, θ, rotation_plane=(1,2))

Rotate an `arr` in the plane `rotation_plane` with an angle `θ` in degree
around the center pixel.

`rotate!` is also available.
"""
function rotate(arr, θ, rotation_plane=(1, 2))
    return rotate!(copy(arr), θ, rotation_plane)
end

"""
    rotate!(arr, θ, rotation_plane=(1,2))

Rotate an `arr` in the plane `rotation_plane` with an angle `θ` in degree
around the center pixel.
"""
function rotate!(arr, θ, rotation_plane=(1, 2))
    # convert to radiants
    θ = deg2rad(θ)
    a,b = rotation_plane

    # parameters for shearing
    α = -tan(θ/2) * size(arr, a)
	β = sin(θ) * size(arr, b)

    # do the three step shearing
    shear!(arr, α, a, b)
    shear!(arr, β, b, a)
    shear!(arr, α, a, b)

    return arr
end


text:
"""
Render an inline R script, substituting invalid "\$" signs for Julia symbols
"""
function render(script::String)
    symdict = OrderedDict{String,Any}()
    local k = 0
    local lastex = RParseError()
    local line
    local col
    local c

    if !isascii(script)
        if !rcopy(reval("isTRUE(l10n_info()\$`UTF-8`)"))
            throw(RParseError("unicode script is not supported"))
        end
    end
    while true
        parse_error = false
        sf = protect(rcall_p(:srcfile, "xx"))
        try
            rparse_p(script, sf)
        catch ex
            lastex = ex
            if isa(ex, RParseError)
                parse_error = true
            else
                unprotect(1)
                throw(ex)
            end
        end
        # break if parse complete
        if !parse_error
            unprotect(1)
            break
        end

        try
            parsedata_p = rcall_p(:getParseData, sf)
            try
                parsedata = protect(parsedata_p)
                n = length(parsedata[1])
                line = parsedata[1][n]
                col = parsedata[2][n]
                c = rcopy(String, parsedata[:text][n])[1]
            finally
                unprotect(1)
            end
        catch
            throw(lastex)
        finally
            unprotect(1)
        end

        # break if the parse error is not caused by $
        if c != '\$'
            throw(lastex)
        end

        index = 0
        for i in 1:(line-1)
            index = something(findnext(isequal('\n'), script, index+1), 0)
        end
        for j in 1:col
            index = nextind(script, index)
        end

        try
            c = script[index]
        catch e
            c = ' '
        end
        if c != '\$'
            throw(lastex)
        end

        ast, i = Meta.parse(script, index+1, greedy=false)

        if isa(ast,Symbol)
            sym = "$ast"
        elseif isa(ast, Expr) && !(ast.head == :error || ast.head == :continue || ast.head == :incomplete)
            sym = "($ast)"
            # if an expression has already appeared, we generate a new symbol so it will be evaluated twice (e.g. `R"$(rand(10)) == $(rand(10))"`)
            if haskey(symdict, sym)
                sym *= "##$k"
                k = k + 1
            end
        elseif isa(ast, Expr) && (ast.head == :incomplete || ast.head == :continue)
            throw(RParseIncomplete("incomplete julia expression"))
        else
            throw(lastex)
        end

        symdict[sym] = ast
        script = string(script[1:index-1], "`#JL`\$`",sym,'`', script[i:end])
    end

    return script, symdict
end

"""
Prepare code for evaluating the julia expressions. When the code is execulated,
the results are stored in the R environment `#JL`.
"""
function prepare_inline_julia_code(symdict, escape::Bool=false)
    new_env = Expr(:(=), :env, Expr(:call, reval_p, Expr(:call, rparse_p, "`#JL` <- new.env()")))
    blk = Expr(:block)
    for (rsym, expr) in symdict
        push!(blk.args, Expr(:(=), Expr(:ref, :env, rsym), escape ? esc(expr) : expr))
    end
    Expr(:let, new_env, blk)
end


text:
using Comrade
using Distributions
using Pathfinder
using AdvancedHMC
using Plots
# load eht-imaging we use this to load eht data
load_ehtim()
# To download the data visit https://doi.org/10.25739/g85n-f134
obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "SR1_M87_2017_096_lo_hops_netcal_StokesI.uvfits"))
obs.add_scans()
# kill 0-baselines since we don't care about
# large scale flux and make scan-average data
obs = obs.flag_uvdist(uv_min=0.1e9).avg_coherent(0.0, scan_avg=true)
# extract log closure amplitudes and closure phases
dlcamp = extract_lcamp(obs; count="min")
dcphase = extract_cphase(obs, count="min")
# form the likelihood
lklhd = RadioLikelihood(dlcamp, dcphase)
# build the model here we fit a ring with a azimuthal
#brightness variation and a Gaussian
function model(θ)
  @unpack rad, wid, a, b, f, sig, asy, pa, x, y = θ
  ring = f*smoothed(stretched(MRing((a,), (b,)), rad, rad), wid)
  g = (1-f)*shifted(rotated(stretched(Gaussian(), sig*asy, sig), pa), x, y)
  return ring + g
end
# define the priors
prior = (
          rad = Uniform(μas2rad(10.0), μas2rad(30.0)),
          wid = Uniform(μas2rad(1.0), μas2rad(10.0)),
          a = Uniform(-0.5, 0.5), b = Uniform(-0.5, 0.5),
          f = Uniform(0.0, 1.0),
          sig = Uniform(μas2rad(1.0), μas2rad(40.0)),
          asy = Uniform(0.0, 0.75),
          pa = Uniform(0.0, 1π),
          x = Uniform(-μas2rad(80.0), μas2rad(80.0)),
          y = Uniform(-μas2rad(80.0), μas2rad(80.0))
        )
# Now form the posterior
post = Posterior(lklhd, prior, model)
# We will use HMC to sample the posterior.
# First to reduce burn in we use pathfinder
q, phi, _ = multipathfinder(post, 100)
# now we sample using hmc
metric = DiagEuclideanMetric(dimension(post))
chain, stats = sample(post, HMC(;metric), 2000; nadapts=1000, init_params=phi[1])
# plot a draw from the posterior
plot(model(chain[end-1]), xlims=(-80.0, 80.0), ylims=(-80.0,80.0), colorbar=nothing)


text:
__precompile__(true)

module Optim
    using Calculus
    using PositiveFactorizations
    using Compat
    using ForwardDiff
    using LineSearches

    import Compat.String
    import Compat.view

    import Base.length,
           Base.push!,
           Base.show,
           Base.getindex,
           Base.setindex!

    export optimize,
           DifferentiableFunction,
           TwiceDifferentiableFunction,
           OptimizationOptions,
           OptimizationState,
           OptimizationTrace,

           AcceleratedGradientDescent,
           BFGS,
           Brent,
           ConjugateGradient,
           Fminbox,
           GoldenSection,
           GradientDescent,
           LBFGS,
           MomentumGradientDescent,
           NelderMead,
           Newton,
           NewtonTrustRegion,
           SimulatedAnnealing,
           ParticleSwarm

    # Types
    include("types.jl")

    # API
    include("api.jl")

    # Generic stuff
    include("utilities/generic.jl")

    # Maxdiff
    include("utilities/maxdiff.jl")

    # Tracing
    include("utilities/update.jl")

    # Grid Search
    include("grid_search.jl")

    # preconditioning functionality
    include("precon.jl")

    # Gradient Descent
    include("gradient_descent.jl")
    include("accelerated_gradient_descent.jl")
    include("momentum_gradient_descent.jl")

    # Conjugate gradient
    include("cg.jl")

    # Newton and Quasi-Newton Methods
    include("newton.jl")
    include("newton_trust_region.jl")
    include("bfgs.jl")
    include("l_bfgs.jl")

    # Constrained optimization
    include("fminbox.jl")

    # trust region methods
    include("levenberg_marquardt.jl")

    # Heuristic Optimization Methods
    include("nelder_mead.jl")
    include("simulated_annealing.jl")
    include("particle_swarm.jl")

    # Univariate methods
    include("golden_section.jl")
    include("brent.jl")

    # Backward compatibility
    include("deprecate.jl")

    # End-User Facing Wrapper Functions
    include("optimize.jl")

    # Convergence
    include("utilities/assess_convergence.jl")

    # Traces
    include("utilities/trace.jl")

    # Examples for testing
    include(joinpath("problems", "unconstrained.jl"))
    include(joinpath("problems", "univariate.jl"))
end


text:
using AutomatonEdification

@testset "mse" begin
    @test AutomatonEdification.LossFunctions.mse([0], [0]) == 0
    @test AutomatonEdification.LossFunctions.mse([1], [-1]) == 4

    @test AutomatonEdification.LossFunctions.mse([-1, -1, -1], [1, 2, 3]) ==
        ((1 - (-1))^2 + (2 - (-1))^2 + (3 - (-1))^2) / 3
end

@testset "cross_entropy" begin
    # TODO: add tests for cross_entropy
end

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
"""
    parseqps([t::Type{T}, ]filename::AbstractString) -> mps

Parse QPS formatted file `filename`. `T` must be an `AbstractFloat` type. When
only `filename` is provided, the first argument defaults to `Float64`.

`mps` is of `MPSDescription` type. It contains the description of the the
following optimization problem formuliation:

```
minimize    0.5*x'*Q*x + q₁'*x + q₂
subject to       A*x = b
            c₁ ≤ C*x ≤ c₂
            lb ≤  x  ≤ ub
```

where all the matrices and vectors are of type `T`.

See also: `name`, `variables`, `objective`, `equalities`, `inequalities`,
  `bounds` and `canonical`.
"""
function parseqps(::Type{T}, filename::AbstractString) where {T<:AbstractFloat}
  qpname, ineqrowidx, eqrowidx, objgiven, objsymbol, vars, rows, cols, rhs, ranges, bounds, quadobj, varvec = _parseqps(T,filename)

  # Create the matrices
  n, m₁, m₂ = length(vars), ineqrowidx, eqrowidx
  Q         = zeros(T, n, n )
  q₁        = zeros(T, n)
  q₂        = zero(T)
  c₁        = fill(convert(T, -Inf), m₁)
  C         = zeros(T, m₁, n)
  c₂        = fill(convert(T, Inf), m₁)
  A         = zeros(T, m₂, n)
  b         = zeros(T, m₂)
  lb        = zeros(T, n)
  ub        = fill(convert(T, Inf), n)

  # Fill (in)equalities with ranges
  for (key, val) in cols
    colidx = vars[key]
    for (rowsymbol, rowval) in val
      rowidx, rowtype = rows[rowsymbol]
      if rowidx == 0
        q₁[colidx] = rowval
      elseif rowtype == :eq
        A[rowidx, colidx] = rowval
        b[rowidx]         = get(rhs, rowsymbol, zero(T))
      elseif rowtype == :leq
        C[rowidx, colidx] = rowval
        c₂[rowidx]        = get(rhs, rowsymbol, zero(T))
        if haskey(ranges, rowsymbol)
          c₁[rowidx]      = c₂[rowidx] - abs(ranges[rowsymbol])
        end
      else
        C[rowidx, colidx] = rowval
        c₁[rowidx]        = get(rhs, rowsymbol, zero(T))
        if haskey(ranges, rowsymbol)
          c₂[rowidx]      = c₁[rowidx] + abs(ranges[rowsymbol])
        end
      end
    end
  end

  if objgiven
    q₂ = -get(rhs, objsymbol, zero(T))
  end

  # Fill bounds
  for (key, val) in bounds
    colidx = vars[key]
    for (bndtype, bndval) in val
      if bndtype == :lower
        lb[colidx] = bndval
      elseif bndtype == :upper
        ub[colidx] = bndval
      elseif bndtype == :fixed
        lb[colidx] = ub[colidx] = bndval
      else
        lb[colidx] = convert(T, -Inf)
      end
    end
  end

  # Fill Q
  for (colkey, val) in quadobj
    colidx = vars[colkey]
    for (rowkey, rowval) in val
      rowidx = vars[rowkey]
      Q[rowidx, colidx] = Q[colidx, rowidx] = rowval
    end
  end

  # Return the problem
  return MPSDescription(T, Q, q₁, q₂, c₁, C, c₂, A, b, lb, ub, varvec, qpname)
end

parseqps(filename::AbstractString) = parseqps(Float64, filename)

"""
    parse_sparseqps([t::Type{T}, ]filename::AbstractString) -> mps

Parse QPS formatted file `filename`, storing matrices in a sparse format. `T` must be an `AbstractFloat` type. When
only `filename` is provided, the first argument defaults to `Float64`.

`mps` is of `MPSDescription` type. It contains the description of the the
following optimization problem formuliation:

```
minimize    0.5*x'*Q*x + q₁'*x + q₂
subject to       A*x = b
            c₁ ≤ C*x ≤ c₂
            lb ≤  x  ≤ ub
```

where all the matrices and vectors are of type `T`.

See also: `name`, `variables`, `objective`, `equalities`, `inequalities`,
  `bounds` and `canonical`.
"""
function parse_sparseqps(::Type{T}, filename::AbstractString) where {T<:AbstractFloat}
  qpname, ineqrowidx, eqrowidx, objgiven, objsymbol, vars, rows, cols, rhs, ranges, bounds, quadobj, varvec = _parseqps(T,filename)

  # Create the vectors
  n, m₁, m₂ = length(vars), ineqrowidx, eqrowidx
  q₁        = zeros(T, n)
  q₂        = zero(T)
  c₁        = fill(convert(T, -Inf), m₁)
  c₂        = fill(convert(T, Inf), m₁)
  b         = zeros(T, m₂)
  lb        = zeros(T, n)
  ub        = fill(convert(T, Inf), n)

  # Prepare matrices
  Qᵢ        = Vector{Int}()
  Qⱼ        = Vector{Int}()
  Qᵥ        = Vector{T}()
  Cᵢ        = Vector{Int}()
  Cⱼ        = Vector{Int}()
  Cᵥ        = Vector{T}()
  Aᵢ        = Vector{Int}()
  Aⱼ        = Vector{Int}()
  Aᵥ        = Vector{T}()

  # Fill (in)equalities with ranges
  for (key, val) in cols
    colidx = vars[key]
    for (rowsymbol, rowval) in val
      rowidx, rowtype = rows[rowsymbol]
      if rowidx == 0
        q₁[colidx] = rowval
      elseif rowtype == :eq
        push!(Aᵢ,rowidx)
        push!(Aⱼ,colidx)
        push!(Aᵥ,rowval)
        b[rowidx]         = get(rhs, rowsymbol, zero(T))
      elseif rowtype == :leq
        push!(Cᵢ,rowidx)
        push!(Cⱼ,colidx)
        push!(Cᵥ,rowval)
        c₂[rowidx]        = get(rhs, rowsymbol, zero(T))
        if haskey(ranges, rowsymbol)
          c₁[rowidx]      = c₂[rowidx] - abs(ranges[rowsymbol])
        end
      else
        push!(Cᵢ,rowidx)
        push!(Cⱼ,colidx)
        push!(Cᵥ,rowval)
        c₁[rowidx]        = get(rhs, rowsymbol, zero(T))
        if haskey(ranges, rowsymbol)
          c₂[rowidx]      = c₁[rowidx] + abs(ranges[rowsymbol])
        end
      end
    end
  end
  C = sparse(Cᵢ,Cⱼ,Cᵥ,m₁,n)
  A = sparse(Aᵢ,Aⱼ,Aᵥ,m₂,n)

  if objgiven
    q₂ = -get(rhs, objsymbol, zero(T))
  end

  # Fill bounds
  for (key, val) in bounds
    colidx = vars[key]
    for (bndtype, bndval) in val
      if bndtype == :lower
        lb[colidx] = bndval
      elseif bndtype == :upper
        ub[colidx] = bndval
      elseif bndtype == :fixed
        lb[colidx] = ub[colidx] = bndval
      else
        lb[colidx] = convert(T, -Inf)
      end
    end
  end

  # create Q
  for (colkey, val) in quadobj
    colidx = vars[colkey]
    for (rowkey, rowval) in val
      rowidx = vars[rowkey]
      push!(Qᵢ,rowidx)
      push!(Qⱼ,colidx)
      push!(Qᵥ,rowval)
      if rowidx != colidx
        push!(Qᵢ,colidx)
        push!(Qⱼ,rowidx)
        push!(Qᵥ,rowval)
      end
    end
  end
  Q = sparse(Qᵢ,Qⱼ,Qᵥ,n,n)

  # Return the problem
  return MPSDescription(T, Q, q₁, q₂, c₁, C, c₂, A, b, lb, ub, varvec, qpname)
end

parse_sparseqps(filename::AbstractString) = parse_sparseqps(Float64, filename)

function _parseqps(::Type{T}, filename::AbstractString) where {T<:AbstractFloat}
  # Initialize auxiliary variables
  qpname                        = "QP"
  mode                          = :NAME
  ineqrowidx, eqrowidx, varidx  = 0, 0, 1
  objgiven                      = false
  objsymbol                     = :obj

  # Define sections
  vars    = Dict{Symbol,Int}()
  rows    = Dict{Symbol,Pair{Int,Symbol}}()
  cols    = Dict{Symbol,Vector{Pair{Symbol,T}}}()
  rhs     = Dict{Symbol,T}()
  ranges  = Dict{Symbol,T}()
  bounds  = Dict{Symbol,Vector{Pair{Symbol,T}}}()
  quadobj = Dict{Symbol,Vector{Pair{Symbol,T}}}()

  # Parse the file
  open(filename) do io
    firstline = split(readline(io))
    if Symbol(firstline[1]) == :NAME
      qpname = join(firstline[2:end], " ")
    else
      throw(ArgumentError("`NAME` field is expected on the first line"))
    end
    for line in eachline(io)
      words = split(line)
      length(words) == 1 && (mode = Symbol(words[1]); continue)
      if mode == :ROWS
        rowsym = Symbol(words[2])
        rows[rowsym] =
          words[1] == "N" ? (objgiven = true; objsymbol = rowsym; Pair(0, :obj)) :
        words[1] == "L" ? (ineqrowidx += 1; Pair(ineqrowidx, :leq))             :
          words[1] == "G" ? (ineqrowidx += 1; Pair(ineqrowidx, :geq))           :
          (eqrowidx   += 1; Pair(eqrowidx, :eq))
      elseif mode == :COLUMNS
        var = Symbol(words[1])
        if get!(vars, var, 0) == 0
          vars[var] = varidx
          varidx    += 1
        end
        for idx = 2:2:length(words)
          push!(get!(cols, var, Pair{Symbol,T}[]),
                Pair(Symbol(words[idx]), convert(T, parse(Float64, words[idx+1]))))
        end
      elseif mode == :RHS
        # NOTE: Discarding rhs name, *i.e.*, assuming only one rhs.
        for idx = 2:2:length(words)
          rhs[Symbol(words[idx])] = convert(T, parse(Float64, words[idx+1]))
        end
      elseif mode == :RANGES
        # NOTE: Discarding range name, *i.e.*, assuming only one range.
        for idx = 2:2:length(words)
          ranges[Symbol(words[idx])] = convert(T, parse(Float64, words[idx+1]))
        end
      elseif mode == :BOUNDS
        var = Symbol(words[3])
        bnd = words[1] == "LO" ? :lower :
          words[1] == "UP" ? :upper :
          words[1] == "FR" ? :free  : :fixed
        push!(get!(bounds, var, Pair{Symbol,T}[]),
              Pair(bnd, convert(T, parse(Float64, bnd == :free ? "0" : words[4]))))
      elseif mode == :QUADOBJ
        var = Symbol(words[1])
        for idx = 2:2:length(words)
          push!(get!(quadobj, var, Pair{Symbol,T}[]),
                Pair(Symbol(words[idx]), convert(T, parse(Float64, words[idx+1]))))
        end
      elseif mode == :ENDATA
        break
      else
        throw(ArgumentError("$(mode) is not a valid word"))
      end
    end

    varvec    = Vector{Symbol}(undef,length(vars))
    for (varname, varidx) in vars
      varvec[varidx] = varname
    end

    return qpname, ineqrowidx, eqrowidx, objgiven, objsymbol, vars, rows, cols, rhs, ranges, bounds, quadobj, varvec
  end
end


text:
using Plots, LinearAlgebra
pyplot()

function draw_covariance_ellipse!(pose, cov, n)
  pose_xy = pose[1:2]
  cov_xy = cov[1:2, 1:2]

  eig_vals, eig_vecs = eigen(cov_xy)
  if eig_vals[1] >= eig_vals[2]
    big_idx = 1
    small_idx = 2
  else
    big_idx = 2
    small_idx = 1
  end

  # 3 sigma covariance, x-y
  t = 0:0.1:(pi * 2 + 0.1)
  a = sqrt(eig_vals[big_idx])
  b = sqrt(eig_vals[small_idx])
  x = [a * n * cos(it) for it in t]
  y = [b * n * sin(it) for it in t]
  angle = atan(eig_vecs[big_idx, 2], eig_vecs[big_idx, 1])
  rot_mat = [cos(angle) sin(angle); -sin(angle) cos(angle)]
  rot_xy = rot_mat * [x y]'
  px = rot_xy[1, begin:end] .+ pose_xy[1]
  py = rot_xy[2, begin:end] .+ pose_xy[2]
  plot!(px, py, color="blue", legend=false,
        aspect_ratio=true, alpha=0.5)

  # 3 sigma covariance, theta
  if size(pose) == (1, 3) && size(cov) == (3, 3)
    pose_theta = pose[3]
    theta_3_sigma = sqrt(cov[3, 3]) * n
    xs = [pose_xy[1] + cos(pose_theta - theta_3_sigma),
          pose_xy[1],
          pose_xy[1] + cos(pose_theta + theta_3_sigma)]
    ys = [pose_xy[2] + sin(pose_theta - theta_3_sigma),
          pose_xy[2],
          pose_xy[2] + sin(pose_theta + theta_3_sigma)]
    plot!(xs, ys, color="blue", legend=false,
          aspect_ratio=true, alpha=0.5)
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
module Optimisers

using Functors: functor, fmap, isleaf

include("interface.jl")
include("rules.jl")

export Descent, ADAM, Momentum, Nesterov, RMSProp,
       ADAGrad, AdaMax, ADADelta, AMSGrad, NADAM, ADAMW, RADAM, OADAM, AdaBelief,
       WeightDecay, OptimiserChain

end # module


text:
#= src/linop.jl
=#

abstract type AbstractLinearOperator{T} end
const AbstractLinOp = AbstractLinearOperator
const AbstractMatOrLinOp{T} = Union{AbstractMatrix{T}, AbstractLinOp{T}}

mutable struct LinearOperator{T, Tmul<:Function, Tmulc<:Function} <: AbstractLinOp{T}
  m::Int
  n::Int
  mul!::Tmul
  mulc!::Tmulc
  _tmp::Nullable{Array{T}}
  LinearOperator{T}(m, n, mul!, mulc!, _tmp) where T = new{T, typeof(mul!), typeof(mulc!)}(m, n, mul!, mulc!, _tmp)
end
const LinOp = LinearOperator

mutable struct HermitianLinearOperator{T, Tmul<:Function} <: AbstractLinOp{T}
  n::Int
  mul!::Tmul
  _tmp::Nullable{Array{T}}
  HermitianLinearOperator{T}(m, n, mul!, _tmp) where T = new{T, typeof(mul!)}(m, n, mul!, _tmp)
end
const HermLinOp = HermitianLinearOperator

function LinOp(A)
  ishermitian(A) && return HermLinOp(A)
  T = eltype(A)
  m, n = size(A)
  ml!  = (y, _, x) ->  mul!(y, A, x)
  mulc! = (y, _, x) -> mul!(y, A', x)
  LinOp{T}(m, n, ml!, mulc!, nothing)
end

function HermLinOp(A)
  T = eltype(A)
  m, n = size(A)
  m == n || throw(DimensionMismatch)
  ml! = (y, _, x) ->  mul!(y, A, x)
  HermLinOp{T}(n, ml!, nothing)
end

convert(::Type{Array}, A::AbstractLinOp) = Matrix(A)
convert(::Type{Array{T}}, A::AbstractLinOp) where {T} = convert(Array{T}, Matrix(A))
convert(::Type{LinOp}, A::HermLinOp) = convert(LinOp{eltype(A)}, A)
convert(::Type{LinOp{T}}, A::HermLinOp{T}) where T = LinOp{T}(A.n, A.n, A.mul!, A.mul!, A._tmp)

adjoint(A::LinOp{T}) where {T} = LinOp{T}(A.n, A.m, A.mulc!, A.mul!, nothing)
adjoint(A::HermLinOp) = A

eltype(::AbstractLinOp{T}) where {T} = T

Matrix(A::AbstractLinOp{T}) where {T} = A*Matrix{T}(I, size(A)...)

getindex(A::AbstractLinOp, ::Colon, ::Colon) = Matrix(A)
function getindex(A::AbstractLinOp{T}, ::Colon, cols) where T
  k = length(cols)
  S = zeros(T, size(A,2), k)
  for i = 1:k
    S[cols[i],i] = 1
  end
  A*S
end
function getindex(A::AbstractLinOp{T}, rows, ::Colon) where T
  k = length(rows)
  S = zeros(T, size(A,1), k)
  for i = 1:k
    S[rows[i],i] = 1
  end
  (A'*S)'
end
function getindex(A::AbstractLinOp{T}, rows, cols) where T
  if length(rows) >= length(cols)  return A[:,cols][rows,:]
  else                             return A[rows,:][:,cols]
  end
end

ishermitian(::LinOp) = false
ishermitian(::HermLinOp) = true
issymmetric(::LinOp) = false
issymmetric(A::HermLinOp) = isreal(A)

isreal(::AbstractLinOp{T}) where {T} = T <: Real

size(A::LinOp) = (A.m, A.n)
size(A::LinOp, dim::Integer) = dim == 1 ? A.m : (dim == 2 ? A.n : 1)
size(A::HermLinOp) = (A.n, A.n)
size(A::HermLinOp, dim::Integer) = (dim == 1 || dim == 2) ? A.n : 1

function transpose(A::LinOp{T}) where T
  n, m = size(A)
  mul!  = (y, L, x) -> (A.mulc!(y, L, conj(x)); conj!(y))
  mulc! = (y, L, x) -> ( A.mul!(y, L, conj(x)); conj!(y))
  LinOp{T}(m, n, mul!, mulc!, nothing)
end
function transpose(A::HermLinOp{T}) where T
  n = size(A, 1)
  mul! = (y, L, x) -> (A.mul!(y, L, conj(x)); conj!(y))
  HermLinOp{T}(n, mul!, nothing)
end

# matrix multiplication

^(A::AbstractLinOp, p::Integer) = Base.power_by_squaring(A, p)

mul!(C, A::AbstractLinOp, B::AbstractVecOrMat) = A.mul!(C, A, B)

  mul!(C, A::Adjoint{<:Any,<:LinOp}, B::AbstractVecOrMat) = parent(A).mulc!(C, parent(A), B)
  mul!(C, A::Adjoint{<:Any,<:HermLinOp}, B::AbstractVecOrMat) = mul!(C, parent(A), B)
  mul!(C, A::AbstractMatrix, B::Adjoint{<:Any,<:AbstractLinOp}) = adjoint!(C, parent(B)*A')

mul!(C, A::AbstractMatrix, B::AbstractLinOp) = adjoint!(C, B'*A')

*(A::AbstractLinOp{T}, x::AbstractVector) where {T} =
  (y = Array{T}(undef, size(A,1)); mul!(y, A, x))
*(A::AbstractLinOp{T}, B::AbstractMatrix) where {T} =
  (C = Array{T}(undef, size(A,1), size(B,2)); mul!(C, A, B))
*(A::AbstractMatrix, B::AbstractLinOp{T}) where {T} =
  (C = Array{T}(undef, size(A,1), size(B,2)); mul!(C, A, B))

# scalar multiplication/division

for (f, g) in ((:(A::LinOp), :(c::Number)), (:(c::Number), :(A::LinOp)))
  @eval begin
    function *($f, $g)
      T = eltype(A)
      m, n = size(A)
      sc_mul!  = (y, _, x) -> ( mul!(y, A, x); lmul!(c, y))
      sc_mulc! = (y, _, x) -> (mul!(y, A', x); lmul!(c, y))
      LinOp{T}(m, n, sc_mul!, sc_mulc!, nothing)
    end
  end
end

for (f, g) in ((:(A::HermLinOp), :(c::Number)), (:(c::Number), :(A::HermLinOp)))
  @eval begin
    function *($f, $g)
      T = eltype(A)
      n = size(A, 1)
      sc_mul! = (y, _, x) -> (mul!(y, A, x); lmul!(c, y))
      HermLinOp{T}(n, sc_mul!, nothing)
    end
  end
end

-(A::AbstractLinOp) = -1*A

/(A::AbstractLinOp, c::Number) = A*(1/c)
\(c::Number, A::AbstractLinOp) = (1/c)*A

# operator addition/subtraction

for (f, a) in ((:+, 1), (:-, -1))
  @eval begin
    function $f(A::AbstractLinOp{T}, B::AbstractLinOp{T}) where T
      size(A) == size(B) || throw(DimensionMismatch)
      m, n = size(A)
      alpha = T($a)
      mul!  =  gen_linop_axpy(A, B, alpha)
      mulc! = gen_linop_axpyc(A, B, alpha)
      LinOp{T}(m, n, mul!, mulc!, nothing)
    end

    function $f(A::HermLinOp{T}, B::HermLinOp{T}) where T
      size(A) == size(B) || throw(DimensionMismatch)
      n = size(A, 1)
      alpha = T($a)
      mul! = gen_linop_axpy(A, B, alpha)
      HermLinOp{T}(n, mul!, nothing)
    end
  end
end

function gen_linop_axpy(A::AbstractLinOp{T}, B::AbstractLinOp{T}, alpha::T) where T
  function linop_axpy!(
      y::AbstractVecOrMat{T}, L::AbstractLinOp{T}, x::AbstractVecOrMat{T}) where T
    if isnull(L._tmp) || size(get(L._tmp)) != size(y)
      L._tmp = similar(y)
    end
    tmp = get(L._tmp)
    mul!( y , A, x)
    mul!(tmp, B, x)
    BLAS.axpy!(alpha*one(T), tmp, y)
  end
end

  function gen_linop_axpyc(A::AbstractLinOp{T}, B::AbstractLinOp{T}, alpha::T) where T
    function linop_axpyc!(
        y::AbstractVecOrMat{T}, L::AbstractLinOp{T}, x::AbstractVecOrMat{T}) where T
      if isnull(L._tmp) || size(get(L._tmp)) != size(y)
        L._tmp = similar(y)
      end
      tmp = get(L._tmp)
      mul!( y , A', x)
      mul!(tmp, B', x)
      BLAS.axpy!(alpha*one(T), tmp, y)
    end
  end

# operator composition

function *(A::AbstractLinOp{T}, B::AbstractLinOp{T}) where T
  mA, nA = size(A)
  mB, nB = size(B)
  nA == mB || throw(DimensionMismatch)
  mul!  =  gen_linop_comp(A, B)
  mulc! = gen_linop_compc(A, B)
  LinOp{T}(mA, nB, mul!, mulc!, nothing)
end

function gen_linop_comp(A::AbstractLinOp{T}, B::AbstractLinOp{T}) where T
  function linop_comp!(
      y::AbstractVector{T}, L::AbstractLinOp{T}, x::AbstractVector{T}) where T
    n = size(B, 1)
    if isnull(L._tmp) || length(get(L._tmp)) != n
      L._tmp = Array{T}(undef, n)
    end
    tmp = get(L._tmp)
    mul!(tmp, B,  x )
    mul!( y , A, tmp)
  end
  function linop_comp!(
      Y::AbstractMatrix{T}, L::AbstractLinOp{T}, X::AbstractMatrix{T}) where T
    m = size(B, 1)
    n = size(X, 2)
    if isnull(L._tmp) || size(get(L._tmp)) != (m, n)
      L._tmp = Array{T}(undef, m, n)
    end
    tmp = get(L._tmp)
    mul!(tmp, B,  X )
    mul!( Y , A, tmp)
  end
end

  function gen_linop_compc(A::AbstractLinOp{T}, B::AbstractLinOp{T}) where T
    function linop_compc!(
        y::AbstractVector{T}, L::AbstractLinOp{T}, x::AbstractVector{T}) where T
      n = size(B, 1)
      if isnull(L._tmp) || length(get(L._tmp)) != n
        L._tmp = Array{T}(undef, n)
      end
      tmp = get(L._tmp)
      mul!(tmp, B',  x )
      mul!( y , A', tmp)
    end
    function linop_compc!(
        Y::AbstractMatrix{T}, L::AbstractLinOp{T}, X::AbstractMatrix{T}) where T
      m = size(B, 1)
      n = size(X, 2)
      if isnull(L._tmp) || size(get(L._tmp)) != (m, n)
        L._tmp = Array{T}(undef, m, n)
      end
      tmp = get(L._tmp)
      mul!(tmp, B',  X )
      mul!( Y , A', tmp)
    end
  end


text:
export @init

function initm(ex)
  quote
    if !@isdefined __inits__
      const __inits__ = []
    end
    if !@isdefined __init__
      __init__() = @init
    end
    push!(__inits__, () -> $ex)
    nothing
  end |> esc
end

function initm()
  :(for f in __inits__
      f()
    end) |> esc
end

macro init(args...)
  initm(args...)
end