text:
@testset "Stack Operators" begin
    @testset "encode/decode num" begin
        tests = [([0x75, 0x39, 0xd7, 0xf9], -2044148085),
        ([0x1b, 0x44, 0x8d, 0x17, 0x61, 0x36, 0x45, 0x51, 0x44, 0x38, 0x15, 0x20, 0xeb, 0x3f, 0xd1, 0xbc], -80840166210881125725200074314630382619),
        ([0xd3, 0x2c, 0xd5, 0xb9, 0x8b, 0xfa, 0x37, 0xd5, 0x24, 0xd0, 0xd2, 0xe8, 0x81, 0x8c, 0xcb, 0x15, 0xe5, 0xfe, 0xdd, 0x95, 0x55, 0xe4, 0x3d, 0x01, 0x82, 0x67, 0xc2, 0x23, 0xc5, 0xd6, 0x0c, 0x90, 0xe9, 0x72, 0xff, 0x9a, 0xd1, 0xcf, 0x16, 0xfa, 0x83], -8494629442587523471751141882831378743107401446683360409874466988135963066279425110931949604908243),
        ([0x7e, 0x5e, 0x14, 0x7b, 0x73, 0x6b, 0x57, 0x64, 0x67, 0x6d, 0x0b, 0x70, 0xa1, 0x6e, 0x71], 588973400481137931526979715425656446),
        ([0x7b], 123)]
        for i in 1:length(tests)
            @test Bitcoin.decode_num(tests[i][1]) == tests[i][2]
            @test Bitcoin.encode_num(tests[i][2]) == tests[i][1]
        end
    end
    @testset "Push value onto stack" begin
        test_stack = [[0x75, 0x39, 0xd7, 0xf9],
        [0x86, 0x8f, 0xf9, 0x73, 0x0a, 0x3c, 0x0c, 0xb9],
        [0xe2, 0xd6],
        [0x7b],
        [0x92, 0x8b, 0xae, 0xed, 0x1e, 0xfb, 0x04, 0x58],
        [0xdb, 0xaf, 0x44],
        [0x84, 0xd7, 0x35, 0xb7],
        [0xaf, 0x80, 0x91]]
        empty_stack = UInt8[]
        @testset "op_0" begin
            stack = test_stack
            want = stack
            push!(want, UInt8[])
            @test Bitcoin.op_0(stack) == true
            @test stack == want
        end
        @testset "op_1negate" begin
            stack = test_stack
            want = stack
            push!(want, [0x81])
            @test Bitcoin.op_1negate(stack) == true
            @test stack == want
        end
        @testset "op_2" begin
            stack = test_stack
            want = stack
            push!(want, [0x02])
            @test Bitcoin.op_2(stack) == true
            @test stack == want
        end
        @testset "op_3" begin
            stack = test_stack
            want = stack
            push!(want, [0x03])
            @test Bitcoin.op_3(stack) == true
            @test stack == want
        end
        @testset "op_4" begin
            stack = test_stack
            want = stack
            push!(want, [0x01])
            @test Bitcoin.op_4(stack) == true
            @test stack == want
        end
        @testset "op_5" begin
            stack = test_stack
            want = stack
            push!(want, [0x05])
            @test Bitcoin.op_5(stack) == true
            @test stack == want
        end
        @testset "op_6" begin
            stack = test_stack
            want = stack
            push!(want, [0x06])
            @test Bitcoin.op_6(stack) == true
            @test stack == want
        end
        @testset "op_7" begin
            stack = test_stack
            want = stack
            push!(want, [0x07])
            @test Bitcoin.op_7(stack) == true
            @test stack == want
        end
        @testset "op_8" begin
            stack = test_stack
            want = stack
            push!(want, [0x08])
            @test Bitcoin.op_8(stack) == true
            @test stack == want
        end
        @testset "op_9" begin
            stack = test_stack
            want = stack
            push!(want, [0x09])
            @test Bitcoin.op_9(stack) == true
            @test stack == want
        end
        @testset "op_10" begin
            stack = test_stack
            want = stack
            push!(want, [0x0a])
            @test Bitcoin.op_10(stack) == true
            @test stack == want
        end
        @testset "op_11" begin
            stack = test_stack
            want = stack
            push!(want, [0x0b])
            @test Bitcoin.op_11(stack) == true
            @test stack == want
        end
        @testset "op_12" begin
            stack = test_stack
            want = stack
            push!(want, [0x0c])
            @test Bitcoin.op_12(stack) == true
            @test stack == want
        end
        @testset "op_13" begin
            stack = test_stack
            want = stack
            push!(want, [0x0d])
            @test Bitcoin.op_13(stack) == true
            @test stack == want
        end
        @testset "op_14" begin
            stack = test_stack
            want = stack
            push!(want, [0x0e])
            @test Bitcoin.op_14(stack) == true
            @test stack == want
        end
        @testset "op_15" begin
            stack = test_stack
            want = stack
            push!(want, [0x0f])
            @test Bitcoin.op_15(stack) == true
            @test stack == want
        end
        @testset "op_16" begin
            stack = test_stack
            want = stack
            push!(want, [0x10])
            @test Bitcoin.op_16(stack) == true
            @test stack == want
        end
    end
    @testset "Conditional control flow" begin
        @testset "op_nop" begin
            @test Bitcoin.op_nop([UInt8[]]) == true
            @test Bitcoin.op_nop([[0x64],
                [0x86, 0x8f, 0xf9, 0x73, 0x0a, 0x3c, 0x0c, 0xb9],
                [0xe2, 0xd6]]) == true
        end
        # TODO review those tests
        @testset "op_if" begin
            want = [[0x64],
                [0x86, 0x8f, 0xf9, 0x73, 0x0a, 0x3c, 0x0c, 0xb9],
                [0xe2, 0xd6]]
            stack = Array{UInt8,1}[]
            items = [0x63,0x0c,0xe2]
            @test Bitcoin.op_if(stack,items) == false
            @test stack == []
            @test items == [0x63,0x0c,0xe2]

            stack = want
            @test Bitcoin.op_if(stack,items) == false
            @test stack == want
            @test items == UInt8[]

            items == [0x67,0x0d,0x69]
            @test Bitcoin.op_if(stack,items) == false
            @test stack == want
            @test items == UInt8[]

            items == [0x63,0x68,0x6b]
            @test Bitcoin.op_if(stack,items) == false
            @test stack == want
            @test items == UInt8[]

            items == [0x0d,0x0f,0x7a]
            @test Bitcoin.op_if(stack,items) == false
            @test stack == want
            @test items == UInt8[]

            push!(stack,UInt8[])
            items = [0x68,0xc4,0x07]
            @test Bitcoin.op_if(stack,items) == true
            @test stack == want
            @test items == [0xc4,0x07]

            items = [0x68,0xc4,0x07]
            @test Bitcoin.op_if(stack,items) == true
            pop!(want)
            @test stack == want
            @test items == [0xc4,0x07]
        end
        # TODO create those tests
        @testset "op_notif" begin
        end
        @testset "op_verify" begin
            stack = Array{UInt8,1}[]
            @test Bitcoin.op_verify(stack) == false
            @test stack ==  Array{UInt8,1}[]
            stack = [[0x64], [0x86, 0x8f, 0xf9], [0xe2, 0xd6]]
            @test Bitcoin.op_verify(stack) == true
            @test stack ==  [[0x64], [0x86, 0x8f, 0xf9]]
            stack = [[0x64], [0x86, 0x8f, 0xf9], [0xe2, 0xd6], UInt8[]]
            @test Bitcoin.op_verify(stack) == false
            @test stack ==  [[0x64], [0x86, 0x8f, 0xf9], [0xe2, 0xd6]]
        end
        @testset "op_return" begin
            @test Bitcoin.op_return([UInt8[]]) == false
            @test Bitcoin.op_return([[0x64],
                [0x86, 0x8f, 0xf9, 0x73, 0x0a, 0x3c, 0x0c, 0xb9],
                [0xe2, 0xd6]]) == false
        end
    end
    @testset "Timelock operations" begin
        empty = Array{UInt8,1}[]
        stack = Array{UInt8,1}[[0x03], [0x04]]
        @testset "op_checklocktimeverify" begin
            @test Bitcoin.op_checklocktimeverify(stack, 1, 0xffffffff) == false
            @test Bitcoin.op_checklocktimeverify(empty, 100, 1) == false
            push!(stack, [0x81])
            @test Bitcoin.op_checklocktimeverify(stack, 100, 1) == false
            pop!(stack)
            @test Bitcoin.op_checklocktimeverify(stack, 500000001, 1) == false
            @test Bitcoin.op_checklocktimeverify(stack, 1, 1) == false
            @test Bitcoin.op_checklocktimeverify(stack, 100, 1) == true
        end
        # TODO
        @testset "op_checksequenceverify" begin
        end
    end
    @testset "Stack operations" begin
        @testset "op_toaltstack" begin
            stack = Array{UInt8,1}[]
            altstack = stack
            @test Bitcoin.op_toaltstack(stack, altstack) == false
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            altstack = [[0xf4], [0xf6, 0xff, 0xf9], [0xf2, 0xf6]]
            @test Bitcoin.op_toaltstack(stack, altstack) == true
            @test stack == [[0x04], [0x06, 0x0f, 0x09]] && altstack == [[0xf4], [0xf6, 0xff, 0xf9], [0xf2, 0xf6], [0x02, 0x06]]
        end
        @testset "op_fromaltstack" begin
            stack = Array{UInt8,1}[]
            altstack = stack
            @test Bitcoin.op_toaltstack(stack, altstack) == false
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            altstack = [[0xf4], [0xf6, 0xff, 0xf9], [0xf2, 0xf6]]
            @test Bitcoin.op_fromaltstack(stack, altstack) == true
            @test stack == [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06],[0xf2, 0xf6]] && altstack == [[0xf4], [0xf6, 0xff, 0xf9]]
        end
        @testset "op_2drop" begin
            stack = [[0x04]]
            @test Bitcoin.op_2drop(stack) == false
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            @test Bitcoin.op_2drop(stack) == true
            @test stack == [[0x04]]
        end
        @testset "op_2dup" begin
            stack = [[0x04]]
            @test Bitcoin.op_2dup(stack) == false
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            @test Bitcoin.op_2dup(stack) == true
            @test stack == [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06], [0x06, 0x0f, 0x09], [0x02, 0x06]]
        end
        @testset "op_3dup" begin
            stack = [[0x04], [0x06, 0x0f, 0x09]]
            @test Bitcoin.op_3dup(stack) == false
            stack = [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
            @test Bitcoin.op_3dup(stack) == true
            @test stack == [[0x04], [0x06, 0x0f, 0x09], [0x02, 0x06], [0x04], [0x06, 0x0f, 0x09], [0x02, 0x06]]
        end
        @testset "op_2over" begin
            stack = [[0x01], [0x02], [0x03]]
            @test Bitcoin.op_2over(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04], [0x05]]
            @test Bitcoin.op_2over(stack) == true
            @test stack == [[0x01], [0x02], [0x03], [0x04], [0x05], [0x02], [0x03]]
        end
        @testset "op_2rot" begin
            stack = [[0x01], [0x02], [0x03], [0x04], [0x05]]
            @test Bitcoin.op_2rot(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04], [0x05], [0x06], [0x07]]
            @test Bitcoin.op_2rot(stack) == true
            @test stack == [[0x01], [0x02], [0x03], [0x04], [0x05], [0x06], [0x07], [0x02], [0x03]]
        end
        @testset "op_2swap" begin
            stack = [[0x01], [0x02], [0x03]]
            @test Bitcoin.op_2swap(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04], [0x05]]
            @test Bitcoin.op_2swap(stack) == true
            @test stack == [[0x01], [0x04], [0x05], [0x02], [0x03]]
        end
        @testset "op_ifdup" begin
            stack = Array{UInt8,1}[]
            @test Bitcoin.op_ifdup(stack) == false
            stack = [[0x01], [0x00]]
            @test Bitcoin.op_ifdup(stack) == true
            @test stack == [[0x01], [0x00]]
            stack = [[0x01], [0x02]]
            @test Bitcoin.op_ifdup(stack) == true
            @test stack == [[0x01], [0x02], [0x02]]
        end
        @testset "op_depth" begin
            stack = Array{UInt8,1}[]
            @test Bitcoin.op_depth(stack) == true
            @test stack == [[]]
            stack = [[0x01], [0x02]]
            @test Bitcoin.op_depth(stack) == true
            @test stack == [[0x01], [0x02], [0x02]]
        end
        @testset "op_dup" begin
            stack = Array{UInt8,1}[]
            @test Bitcoin.op_dup(stack) == false
            stack = [[0x01], [0x02], [0x03]]
            @test Bitcoin.op_dup(stack) == true
            @test stack == [[0x01], [0x02], [0x03], [0x03]]
        end
        @testset "op_nip" begin
            stack = Array{UInt8,1}[[0x01]]
            @test Bitcoin.op_nip(stack) == false
            stack = [[0x01], [0x02], [0x03]]
            @test Bitcoin.op_nip(stack) == true
            @test stack == [[0x01], [0x03]]
        end
        @testset "op_over" begin
            stack = [[0x01]]
            @test Bitcoin.op_over(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test Bitcoin.op_over(stack) == true
            @test stack == [[0x01], [0x02], [0x03], [0x04], [0x03]]
        end
        @testset "op_pick" begin
            stack = Array{UInt8,1}[[0x01]]
            @test Bitcoin.op_pick(stack) == false
            stack = [[0x01], [0x0b], [0x03], [0x05]]
            @test Bitcoin.op_pick(stack) == false
            stack = [[0x01], [0x0b], [0x03], [0x02]]
            @test Bitcoin.op_pick(stack) == true
            @test stack == [[0x01], [0x0b], [0x03], [0x0b]]
        end
        @testset "op_roll" begin
            stack = Array{UInt8,1}[]
            @test Bitcoin.op_roll(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test Bitcoin.op_roll(stack) == false
            @test stack == [[0x01], [0x02], [0x03]]
            stack = [[0x01], [0x0b], [0x03], [0x00]]
            @test Bitcoin.op_roll(stack) == true
            @test stack == [[0x01], [0x0b], [0x03]]
            stack = [[0x01], [0x0b], [0x03], [0x02]]
            @test Bitcoin.op_roll(stack) == true
            @test stack == [[0x01], [0x03], [0x0b]]
        end
        @testset "op_rot" begin
            stack = [[0x01], [0x02]]
            @test Bitcoin.op_rot(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test Bitcoin.op_rot(stack) == true
            @test stack == [[0x01], [0x03], [0x04], [0x02]]
        end
        @testset "op_swap" begin
            stack = [[0x01]]
            @test Bitcoin.op_swap(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test Bitcoin.op_swap(stack) == true
            @test stack == [[0x01], [0x02], [0x04], [0x03]]
        end
        @testset "op_tuck" begin
            stack = [[0x01]]
            @test Bitcoin.op_tuck(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test Bitcoin.op_tuck(stack) == true
            @test stack == [[0x01], [0x02], [0x04], [0x03], [0x04]]
        end
    end
    @testset "String splice operations" begin
        @testset "op_size" begin
            stack = Array{UInt8,1}[]
            @test Bitcoin.op_size(stack) == false
            stack = [[0x01], [0x02], [0x01, 0x08]]
            @test Bitcoin.op_size(stack) == true
            @test stack == [[0x01], [0x02], [0x01, 0x08], [0x02]]
        end
    end
    @testset "Binary arithmetic and conditionals" begin
        @testset "op_equal" begin
            stack = Array{UInt8,1}[]
            @test Bitcoin.op_equal(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test Bitcoin.op_equal(stack) == true
            @test stack == [[0x01], [0x02], []]
            stack = [[0x01], [0x02], [0x03], [0x03]]
            @test Bitcoin.op_equal(stack) == true
            @test stack == [[0x01], [0x02], [0x01]]
        end
        @testset "op_equalverify" begin
            stack = Array{UInt8,1}[]
            @test Bitcoin.op_equalverify(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test Bitcoin.op_equalverify(stack) == false
            stack = [[0x01], [0x02], [0x03], [0x03]]
            @test Bitcoin.op_equalverify(stack) == true
        end
    end
    @testset "Numeric operators" begin
        empty = Array{UInt8,1}[]
        stack = Array{UInt8,1}[[0x03]]
        @testset "op_1add" begin
            @test Bitcoin.op_1add(empty) == false
            @test Bitcoin.op_1add(stack) == true
            @test stack == [[0x04]]
        end
        @testset "op_1sub" begin
            @test Bitcoin.op_1sub(empty) == false
            @test Bitcoin.op_1sub(stack) == true
            @test stack == [[0x03]]
        end
        @testset "op_negate" begin
            @test Bitcoin.op_negate(empty) == false
            @test Bitcoin.op_negate(stack) == true
            @test stack == [[0x83]]
        end
        @testset "op_abs" begin
            @test Bitcoin.op_abs(empty) == false
            @test Bitcoin.op_abs(stack) == true
            @test stack == [[0x03]]
            @test Bitcoin.op_abs(stack) == true
            @test stack == [[0x03]]
        end
        @testset "op_not" begin
            @test Bitcoin.op_not(empty) == false
            @test Bitcoin.op_not(stack) == true
            @test stack == [UInt8[]]
            @test Bitcoin.op_not(stack) == true
            @test stack == [[0x01]]
        end
        @testset "op_0notequal" begin
            @test Bitcoin.op_0notequal(empty) == false
            stack = [[0x1a]]
            @test Bitcoin.op_0notequal(stack) == true
            @test stack == [[0x01]]
            stack = [[0x00]]
            @test Bitcoin.op_0notequal(stack) == true
            @test stack == [UInt8[]]
        end
        @testset "op_add" begin
            @test Bitcoin.op_add(empty) == false
            stack = [[0x09], [0x02], [0x03]]
            @test Bitcoin.op_add(stack) == true
            @test stack == [[0x09], [0x05]]
        end
        @testset "op_sub" begin
            @test Bitcoin.op_sub(empty) == false
            @test Bitcoin.op_sub(stack) == true
            @test stack == [[0x04]]
        end
        @testset "op_booland" begin
            @test Bitcoin.op_booland(empty) == false
            stack = [[0x00], [0x00]]
            @test Bitcoin.op_booland(stack) == true
            @test stack == [UInt8[]]
            stack = [[0x01], [0x00]]
            @test Bitcoin.op_booland(stack) == true
            @test stack == [UInt8[]]
            stack = [[0x00], [0x01]]
            @test Bitcoin.op_booland(stack) == true
            @test stack == [UInt8[]]
            stack = [[0x01], [0x01]]
            @test Bitcoin.op_booland(stack) == true
            @test stack == [UInt8[0x01]]
        end
        @testset "op_boolor" begin
            @test Bitcoin.op_boolor(empty) == false
            stack = [[0x00], [0x00]]
            @test Bitcoin.op_boolor(stack) == true
            @test stack == [UInt8[]]
            stack = [[0x00], [0x01]]
            @test Bitcoin.op_boolor(stack) == true
            @test stack == [UInt8[0x01]]
            stack = [[0x01], [0x00]]
            @test Bitcoin.op_boolor(stack) == true
            @test stack == [UInt8[0x01]]
            stack = [[0x01], [0x01]]
            @test Bitcoin.op_boolor(stack) == true
            @test stack == [UInt8[0x01]]
        end
        @testset "op_numequal" begin
            @test Bitcoin.op_numequal(empty) == false
            stack = [[0x01], [0x02], [0x00], [0x00]]
            @test Bitcoin.op_numequal(stack) == true
            @test stack == [[0x01], [0x02], [0x01]]
            @test Bitcoin.op_numequal(stack) == true
            @test stack == [[0x01], UInt8[]]
        end
        @testset "op_numequalverify" begin
            @test Bitcoin.op_numequalverify(empty) == false
            stack = [[0x01], [0x02], [0x00], [0x00]]
            @test Bitcoin.op_numequalverify(stack) == true
            @test stack == [[0x01], [0x02]]
            @test Bitcoin.op_numequalverify(stack) == false
            @test stack == empty
        end
        @testset "op_numnotequal" begin
            @test Bitcoin.op_numnotequal(empty) == false
            stack = [[0x01], [0x02], [0x03], [0x03]]
            @test Bitcoin.op_numnotequal(stack) == true
            @test stack == [[0x01], [0x02], UInt8[]]
            stack = [[0x01], [0x02], [0x03], [0x04]]
            @test Bitcoin.op_numnotequal(stack) == true
            @test stack == [[0x01], [0x02], [0x01]]
        end
        @testset "op_lessthan" begin
            @test Bitcoin.op_lessthan(empty) == false
            stack = [[0x01], [0x01], [0x03], [0x04]]
            @test Bitcoin.op_lessthan(stack) == true
            @test stack == [[0x01], [0x01], [0x01]]
            @test Bitcoin.op_lessthan(stack) == true
            @test stack == [[0x01], UInt8[]]
        end
        @testset "op_greaterthan" begin
            @test Bitcoin.op_greaterthan(empty) == false
            stack = [[0x04], [0x01], [0x02], [0x01]]
            @test Bitcoin.op_greaterthan(stack) == true
            @test stack == [[0x04], [0x01], [0x01]]
            @test Bitcoin.op_greaterthan(stack) == true
            @test stack == [[0x04], UInt8[]]
        end
        @testset "op_lessthanorequal" begin
            @test Bitcoin.op_lessthanorequal(empty) == false
            stack = [[0x01], [0x02], [0x04], [0x04]]
            @test Bitcoin.op_lessthanorequal(stack) == true
            @test stack == [[0x01], [0x02], [0x01]]
            @test Bitcoin.op_lessthanorequal(stack) == true
            @test stack == [[0x01], UInt8[]]
        end
        @testset "op_greaterthanorequal" begin
            @test Bitcoin.op_greaterthanorequal(empty) == false
            stack = [[0x04], [0x00], [0x02], [0x02]]
            @test Bitcoin.op_greaterthanorequal(stack) == true
            @test stack == [[0x04], [0x00], [0x01]]
            @test Bitcoin.op_greaterthanorequal(stack) == true
            @test stack == [[0x04], UInt8[]]
        end
        @testset "op_min" begin
            @test Bitcoin.op_min(empty) == false
            stack = [[0x01], [0x08], [0x03], [0x04]]
            @test Bitcoin.op_min(stack) == true
            @test stack == [[0x01], [0x08], [0x03]]
            @test Bitcoin.op_min(stack) == true
            @test stack == [[0x01], [0x03]]
        end
        @testset "op_max" begin
            @test Bitcoin.op_max(empty) == false
            stack = [[0x01], [0x08], [0x03], [0x04]]
            @test Bitcoin.op_max(stack) == true
            @test stack == [[0x01], [0x08], [0x04]]
            @test Bitcoin.op_max(stack) == true
            @test stack == [[0x01], [0x08]]
        end
        @testset "op_within" begin
            @test Bitcoin.op_within(empty) == false
            stack = [[0xa1], [0x02], [0x02], [0x04]]
            @test Bitcoin.op_within(stack) == true
            @test stack == [[0xa1], [0x01]]
            stack = [[0xa1], [0x08], [0x02], [0x04]]
            @test Bitcoin.op_within(stack) == true
            @test stack == [[0xa1], UInt8[]]
        end
    end
    @testset "Cryptographic and hashing operations" begin
        empty = Array{UInt8,1}[]
        stack = Array{UInt8,1}[[0x03], [0x04]]
        @testset "op_ripemd160" begin
            @test Bitcoin.op_ripemd160(empty) == false
            @test Bitcoin.op_ripemd160(stack) == true
            @test stack == Array{UInt8,1}[[0x03], [0x44, 0x9b, 0x34, 0xb6, 0xa3, 0x41, 0x19, 0x43, 0xe3, 0x3a, 0x25, 0x87, 0xeb, 0xf2, 0x81, 0xca, 0xff, 0x16, 0x74, 0x98]]
        end
        stack = Array{UInt8,1}[[0x03], [0x04]]
        @testset "op_sha1" begin
            @test Bitcoin.op_sha1(empty) == false
            @test Bitcoin.op_sha1(stack) == true
            @test stack == Array{UInt8,1}[[0x03], [0xa4, 0x2c, 0x6c, 0xf1, 0xde, 0x3a, 0xbf, 0xde, 0xa9, 0xb9, 0x5f, 0x34, 0x68, 0x7c, 0xbb, 0xe9, 0x2b, 0x9a, 0x73, 0x83]]
        end
        stack = Array{UInt8,1}[[0x03], [0x04]]
        @testset "op_sha256" begin
            @test Bitcoin.op_sha256(empty) == false
            @test Bitcoin.op_sha256(stack) == true
            @test stack == Array{UInt8,1}[[0x03], [0xe5, 0x2d, 0x9c, 0x50, 0x8c, 0x50, 0x23, 0x47, 0x34, 0x4d, 0x8c, 0x07, 0xad, 0x91, 0xcb, 0xd6, 0x06, 0x8a, 0xfc, 0x75, 0xff, 0x62, 0x92, 0xf0, 0x62, 0xa0, 0x9c, 0xa3, 0x81, 0xc8, 0x9e, 0x71]]
        end
        stack = Array{UInt8,1}[[0x03], [0x04]]
        @testset "op_hash160" begin
            @test Bitcoin.op_hash160(empty) == false
            @test Bitcoin.op_hash160(stack) == true
            @test stack == Array{UInt8,1}[[0x03], [0x6d, 0x4c, 0x0a, 0xa9, 0x72, 0xc3, 0x14, 0x84, 0x0a, 0xc0, 0x7b, 0xe9, 0x6c, 0x5d, 0xde, 0x9c, 0x71, 0x4c, 0x9c, 0xa4]]
        end
        stack = Array{UInt8,1}[[0x03], [0x04]]
        @testset "op_hash256" begin
            @test Bitcoin.op_hash256(empty) == false
            @test Bitcoin.op_hash256(stack) == true
            @test stack == Array{UInt8,1}[[0x03], [0x21, 0x4e, 0x63, 0xbf, 0x41, 0x49, 0x0e, 0x67, 0xd3, 0x44, 0x76, 0x77, 0x8f, 0x67, 0x07, 0xaa, 0x6c, 0x8d, 0x2c, 0x8d, 0xcc, 0xdf, 0x78, 0xae, 0x11, 0xe4, 0x0e, 0xe9, 0xf9, 0x1e, 0x89, 0xa7]]
        end
        # TODO
        @testset "op_codeseparator" begin
        end
        # TODO
        @testset "op_checksig" begin
            z = parse(BigInt, "7c076ff316692a3d7eb3c3bb0f8b1488cf72e1afcd929e29307032997a838a3d", base=16)
            sec = hex2bytes("04887387e452b8eacc4acfde10d9aaf7f6d9a0f975aabb10d006e4da568744d06c61de6d95231cd89026e286df3b6ae4a894a3378e393e93a0f45b666329a0ae34")
            sig = hex2bytes("3045022000eff69ef2b1bd93a66ed5219add4fb51e11a840f404876325a1e8ffe0529a2c022100c7207fee197d27c618aea621406f6bf5ef6fca38681d82b2f06fddbdce6feab601")
            stack = [sig, sec]
            @test Bitcoin.op_checksig(stack, z) == true
            @test Bitcoin.decode_num(stack[1]) == 1
        end
        # TODO
        @testset "op_checksigverify" begin
        end
        # TODO
        @testset "op_checkmultisig" begin
        end
        # TODO
        @testset "op_checkmultisigverify" begin
        end
    end
    @testset "Nonoperators" begin
    end
    @testset "Reserved OP codes for internal use by the parser operators" begin
    end
end


text:
using ..Galerkin:
    galerkin_matrix,
    galerkin_projection!,
    galerkin_projection  # for Documenter only...

using LinearAlgebra:
    cholesky!,
    Cholesky,
    ldiv!

@doc raw"""
    MinimiseL2Error <: AbstractApproxMethod

Approximate a given function ``f(x)`` by minimisation of the ``L^2`` distance
between ``f`` and its spline approximation ``g(x)``.

# Extended help

Minimises the ``L^2`` distance between the two functions:

```math
{\left\lVert f - g \right\rVert}^2 = \left< f - g, f - g \right>,
```

where

```math
\left< u, v \right> = ∫_a^b u(x) \, v(x) \, \mathrm{d}x
```

is the inner product between two functions, and ``a`` and ``b`` are the
boundaries of the prescribed B-spline basis.
Here, ``g`` is the spline ``g(x) = ∑_{i = 1}^N c_i \, b_i(x)``, and
``\{ b_i \}_{i = 1}^N`` is a prescribed B-spline basis.

One can show that the optimal coefficients ``c_i`` minimising the ``L^2`` error
are the solution to the linear system ``\bm{M} \bm{c} = \bm{φ}``,
where ``M_{ij} = \left< b_i, b_j \right>`` and ``φ_i = \left< b_i, f \right>``.
These two terms are respectively computed by [`galerkin_matrix`](@ref) and
[`galerkin_projection`](@ref).

The integrals associated to ``\bm{M}`` and ``\bm{φ}`` are computed via
Gauss--Legendre quadrature.
The number of quadrature nodes is chosen as a function of the order ``k`` of the
prescribed B-spline basis, ensuring that ``\bm{M}`` is computed exactly (see
also [`galerkin_matrix`](@ref)).
In the particular case where ``f`` is a polynomial of degree ``k - 1``, this
also results in an exact computation of ``\bm{φ}``.
In more general cases, as long as ``f`` is smooth enough, this is still expected
to yield a very good approximation of the integral, and thus of the optimal coefficients ``c_i``.

"""
struct MinimiseL2Error <: AbstractApproxMethod end

function approximate(f, B::AbstractBSplineBasis, m::MinimiseL2Error)
    T = typeof(f(first(knots(B))))
    S = Spline{T}(undef, B)
    M = galerkin_matrix(B)  # by default it's a BandedMatrix

    # We annotate the return type to avoid inference issue in ArrayLayouts...
    # https://github.com/JuliaMatrices/ArrayLayouts.jl/issues/66
    Mfact = cholesky!(M) :: Cholesky{eltype(M), typeof(parent(M))}

    data = (; M = Mfact)
    A = SplineApproximation(m, S, data)
    approximate!(f, A)
end

function _approximate!(f, A, m::MinimiseL2Error)
    @assert method(A) === m
    S = spline(A)
    cs = coefficients(S)
    galerkin_projection!(f, cs, basis(S))  # computes rhs onto cs
    ldiv!(data(A).M, cs)  # now cs = M \ rhs
    A
end

ext: jl

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
import .ArchGDAL: RasterDataset, AbstractRasterBand,
  getgeotransform, width, height, getname, getcolorinterp,
  getband, nraster, getdataset

function dimname(a::RasterDataset, i)
    if i == 1
        return :Y
    elseif i == 2
        return :X
    elseif i == 3
        return :Band
    else
        error("RasterDataset only has 3 dimensions")
    end
end
function dimvals(a::RasterDataset, i)
    if i == 1
        geo=getgeotransform(a)
        latr = range(geo[1],length=width(a), step=geo[2])
    elseif i == 2
        geo=getgeotransform(a)
        range(geo[4],length=height(a), step=geo[6])
    elseif i == 3
        colnames = map(ib -> getname(getcolorinterp(getband(a,ib))),1:nraster(a))
        if !allunique(colnames)
            colnames = string.("Band_",1:nraster(a))
        end
        colnames
    else
        error("RasterDataset only has 3 dimensions")
    end
end
iscontdim(a::RasterDataset, i) = i < 3 ? true : nraster(a)<8
function getattributes(a::RasterDataset)
    globatts = Dict{String,Any}(
        "projection_PROJ4"=>ArchGDAL.toPROJ4(ArchGDAL.newspatialref(ArchGDAL.getproj(a))),
        "projection_WKT"=>ArchGDAL.toWKT(ArchGDAL.newspatialref(ArchGDAL.getproj(a))),
    )
    bands = (getbandattributes(ArchGDAL.getband(a, i)) for i in 1:size(a, 3))
    allbands = mergewith(bands...) do a1,a2
        isequal(a1,a2) ? a1 : missing
    end
    merge(globatts, allbands)
end

function dimname(::AbstractRasterBand, i)
    if i == 1
        return :Y
    elseif i == 2
        return :X
    else
        error("RasterDataset only has 3 dimensions")
    end
end
function dimvals(b::AbstractRasterBand, i)
    geo = getgeotransform(getdataset(b))
    if i == 1
        range(geo[1],length=width(b), step=geo[2])
    elseif i == 2
        range(geo[4],length=height(b), step=geo[6])
    else
        error("RasterDataset only has 3 dimensions")
    end
end
iscontdim(a::AbstractRasterBand, i) = true
function getattributes(a::AbstractRasterBand)
  atts = getattributes(ArchGDAL.RasterDataset(ArchGDAL.getdataset(a)))
  bandatts = getbandattributes(a)
  merge(atts, bandatts)
end

function insertattifnot!(attrs, val, name, condition)
    if !condition(val)
        attrs[name] = val
    end
end
function getbandattributes(a::AbstractRasterBand)
  atts = Dict{String,Any}()
  catdict = Dict((i-1)=>v for (i,v) in enumerate(AG.getcategorynames(a)))
  insertattifnot!(atts, AG.getnodatavalue(a), "missing_value", isnothing)
  insertattifnot!(atts, catdict, "labels", isempty)
  insertattifnot!(atts, AG.getunittype(a), "units", isempty)
  insertattifnot!(atts, AG.getoffset(a), "add_offset", iszero)
  insertattifnot!(atts, AG.getscale(a), "scale_factor", x->isequal(x, one(x)))
  atts
end


text:
# Types #
# ==========================
include("attributes.jl")
include("optimizer.jl")
include("solution_cache.jl")
include("scenario.jl")
include("sampler.jl")
include("stage.jl")
include("decisions/decisions.jl")
include("stochasticmodel.jl")
include("stochasticstructure.jl")
include("deterministic_equivalent/deterministicequivalent.jl")
include("block_structured/block_structured.jl")
include("stochasticprogram.jl")
include("decision_variable.jl")
include("objective.jl")
include("constraints.jl")
include("confidenceinterval.jl")
include("helpers.jl")


text:
using Debug
import Plots

@inline function mul3!(T,a,b)
    for j = 1:3, i = 1:3
        T[i,j] = a[i,1]*b[1,j] + a[i,2]*b[2,j] + a[i,3]*b[3,j]
    end
    return T
end
@inline function mul4!(T,a,b)
    for j = 1:4, i = 1:3
        T[i,j] = a[i,1]*b[1,j] + a[i,2]*b[2,j] + a[i,3]*b[3,j] + a[i,4]*b[4,j]
    end
    return T
end

@inline function mul3(a,b)
    T = Matrix{Float64}(3,3)
    for j = 1:3, i = 1:3
        T[i,j] = a[i,1]*b[1,j] + a[i,2]*b[2,j] + a[i,3]*b[3,j]
    end
    return T
end
@inline function mul4(a,b)
    T = Matrix{Float64}(4,4)
    for j = 1:4, i = 1:3
        T[i,j] = a[i,1]*b[1,j] + a[i,2]*b[2,j] + a[i,3]*b[3,j] + a[i,4]*b[4,j]
    end
    T[4,1] = T[4,2] = T[4,3] = 0
    T[4,4] = 1
    return T
end

# """
# Calculates the weight update for a bunch of particles
# x_mu: difference between expected seam location and actual seam location,
# SigmaV: measurement noise covariance matrix
# """

function measUpdate(x_μ, theta, x,traj_nom, σV, index, γ, T_TF_M)
    N = size(x_μ,1)
    xtheta = zeros(N)
    for i = 1:N
        if index[i] == size(traj_nom,3)
            xtheta[i] = 0
        else
            frameP = γ*traj_nom[:,:,index[i]+1] + (1-γ)*traj_nom[:,:,index[i]]
            xtheta[i] = xyθ(frameP*T_TF_M,x[:,:,i]*T_TF_M)
        end
    end
    σ = 2σV.^2
    theta_mu2 = (theta-xtheta).^2
    r = -(x_μ[:,1].^2)./σ[1] - (x_μ[:,2].^2)./σ[2] - theta_mu2./σ[3]

    return r
end

# """
# Calculates the distance between where particle `i` thinks the seam should
# be and where it actually is, accoring to sensor.
# `m`: measurement in ℜ², (x,y)-coordinates from sensor
# `traj_nom`: nominal trajectory, in ℜ4x4xT
# `xp`: particles, in ℜ4x4xN
# `T_TCP_M`: tool to sensor transformation matrix, in ℜ4x4
# """

function gamma(i1::Int, i2, n, m̂i, search_traj)
    p₂ = search_traj[i2,:][:]
    p₁ = search_traj[i1,:][:]
    v = p₂-p₁
    γ = (n⋅(m̂i-p₁))/(n[1]v[1]+n[2]v[2]+n[3]v[3])
    p = p₁ + γ*v
    γ, v, p, i1
end

macro OKγ(γ)
    return :(γ <= 1 && γ >= 0)
end

@debug function distToSeam(m, search_traj, xp, T_TF_M; tries = 2)

    Npart = size(xp,3)
    Ntn = size(search_traj,1)

    # m̂ is where the seam should be according to the hypothesis
    m̂ = zeros(Npart,3)
    xpM = similar(xp)
    for n = 1:Npart
        xpM[:,:,n] = xp[:,:,n]*T_TF_M
        m̂[n,:] = (xpM[:,:,n]*[m[1], m[2], 0, 1])[1:3]'# The z-measurement is effectively 0 since we know the point is in the laser plane

    end

    # TODO: change distance from euclidian 2-norm to Mahalanobis distance
    ni, nd = knnsearch(m̂,xpM,search_traj)
    # @bp
    md = Array(Float64,Npart,3)
    # assert(all(all(isfinite(ni))))
    # assert(all(all(isfinite(nd))))
    index = zeros(Int,size(nd))
    # @bp
    for i = 1:Npart
        n = xpM[1:3,3,i] # Normal is the z-axis
        m̂i = (m̂[i,:]')[:]
        if ni[i] == Ntn # End of seam case
            γ, v, p, i1 = gamma(Ntn-1, Ntn, n, m̂i, search_traj)
            if γ > 1 # No measurement aquired
                nd[i] = 0
                md[i,:] = 0
                index[i] = Ntn
                continue
            end
        else
            γ, v, p, i1 = gamma(ni[i], min(ni[i]+1,Ntn), n, m̂i, search_traj)
        end
        if !@OKγ(γ) && ni[i] != 1# Go forward
            for j = 1:tries-1 # Forward had one try for free
                γ, v, p, i1 = gamma(min(ni[i]+j,Ntn), min(ni[i]+j+1,Ntn), n, m̂i, search_traj)
                if @OKγ(γ)
                    break
                end
            end
        end
        if !@OKγ(γ) && ni[i] != 1# # Go Backwards
            for j = 1:tries
                γ, v, p, i1 = gamma(max(ni[i]-j,1), max(ni[i]-j+1,1), n, m̂i, search_traj)
                if @OKγ(γ)
                    break
                end
            end
        end

        if (!@OKγ(γ) && ni[i] != 1) #|| (n⋅v)/norm(v) < 0.1 # Failed to find a measurement OR Angle is too small
            # @show (n⋅v)/norm(v)
            # @show γ
            # @bp norm(v) != 0
            nd[i] = 0
            md[i,:] = 0
            index[i] = i1
            continue
        end

        # Proceed to calculate measurement
        e = p-m̂i
        e = (xpM[1:3,1:3,i])'e
        nd[i] = e[1]^2 + e[2]^2
        md[i,:] = e
        index[i] = i1
        if abs(e[3]) > 1e-6 # Make sure that the error lies in the laser plane
            @show e[3]
            @bp
        end
    end

    x_μ2 = nd
    # @show sum(x_μ2 .== 0)
    return x_μ2, md, index, γ
end

# """
# `knnsearch(m̂,search_traj, range=0, center = 0)`
# returns the index and the distance²
# """

function knnsearch(m̂,xpM,search_traj)

    Npart,M = size(m̂)
    L = size(search_traj,1)
    idx = zeros(Int64,Npart)
    D = zeros(Npart)
    d² = Array(Float64,L)
    # Loop for each query point
    for n=1:Npart
        @inbounds @fastmath for l = 1:L # Transform back to the sensor frame so we are sure only to compare distance in the xy-plane, this hard coded loop is the most essential place for optimization in the entire filter, it reduces runtime to 1/7 over generic matrix multiplication
            a1 = search_traj[l,1] - m̂[n,1]
            a2 = search_traj[l,2] - m̂[n,2]
            a3 = search_traj[l,3] - m̂[n,3]
            b1 = a1*xpM[1,1,n] + a2*xpM[2,1,n] + a3*xpM[3,1,n]
            b2 = a1*xpM[1,2,n] + a2*xpM[2,2,n] + a3*xpM[3,2,n]
            b3 = a1*xpM[1,3,n] + a2*xpM[2,3,n] + a3*xpM[3,3,n]
            d²[l] = b1*b1 + b2*b2 + b3*b3
        end
        D[n],idx[n] = findmin(d²) # Find the smallest distance
    end
    return idx,D
end

function measFKUpdate(x,T_RB_TCP,f,sigmaVFK)
    N = size(x,3)
    d = sigmaVFK[3] + sigmaVFK[4]*norm(f)
    dR = sigmaVFK[5]
    x_mu2 = zeros(N)
    theta = zeros(N)
    # TODO: update this according to paper
    for i = 1:N
        x_mu2[i] = sum((T_RB_TCP[1:3,4] - x[1:3,4,i]).^2)  # this has been optimized to remove norm calculation which was squared in the next step
        costheta = (trace(T_RB_TCP[1:3,1:3]' * x[1:3,1:3,i])-1)/2
        if costheta > 1 # will never be < -1 because theta is small
            costheta = 1;
        end
        theta[i] = acos(costheta)
    end
    w = zeros(theta)
    outsidePlateauR = theta .> dR^2
    w[outsidePlateauR] -= ((theta[outsidePlateauR]-dR).^2)./(2sigmaVFK[2]^2) # Angle update
    outsidePlateau = x_mu2 .> d^2
    w[outsidePlateau] -= ((sqrt(x_mu2[outsidePlateau])-d).^2)./(2sigmaVFK[1]^2); # xyz-update
    return w

end

# """
# This function propagates all particles forward in time.
# xp: particles to be propagated
# increment: 4x4 transformation matrix describing the move done by the
# robot during the last time step.
# T\_RB\_TF(t-1)\\T\_RB\_TF(t) = T\_TF(t-1)\_TF(t)
# """

function timeUpdate!(xp, increment, sigmaW,sigmaWr, h)
    N = size(xp,3)
    I = eye(3)
    noise = randn(3,N).*(√(h)*sigmaW[1:3])
    R = eye(3)
    for n = 1:N
        xp[:,:,n] = increment*xp[:,:,n]
        xp[1:3,4,n] += noise[:,n] # Add state noise

        #         w = 0.01*pi/180*randn(3,1);
        #         W = [0, -w(3), w(2); w(3), 0, -w(1); -w(2), w(1), 0];
        #         R = expm(W); # this takes a lot of time!
        k = randn(3)
        k /= norm(k) # Axis of rotation
        theta = √(h)*sigmaWr*randn()
        kx = skew!(R,k)
        R = I + sin(theta)*kx + (1-cos(theta))*(k*k'-I)
        @assert !any(isnan(R[:]))
        xp[1:3,1:3,n] = R*xp[1:3,1:3,n]
    end

    return xp
end

@inline function skew!(R,s)
    R[1,1] = 0
    R[2,1] = s[3]
    R[3,1] = -s[2]
    R[1,2] = -s[3]
    R[2,2] = 0
    R[3,2] = s[1]
    R[1,3] = s[2]
    R[2,3] = -s[1]
    R[3,3] = 0
    return R
end

function resample_systematic_exp(w,M)
    N = size(w,1)
    bins = cumsum(w) # devec
    s = (rand()/M+0):1/M:bins[end]
    j = Array(Int64,M)
    bo = 1
    for i = 1:M
        for b = bo:N
            if s[i] <= bins[b]
                j[i] = b
                bo = b
                break
            end
        end
    end
    return j
end


text:
using GapFindFill, Test
# using JuMP, Gurobi, CPLEX, GLPK

println("Starting tests....................")
@time begin include("testGapFind.jl") end
@time begin include("testGapFill.jl") end


text:
# Problem Set I
#
# Student: Luan Mugarte
# Group with: Gabriel Petrini

# Loading packages

using Pkg

import JLD2
import Random
import LinearAlgebra
import Statistics
import CSV
import DataFrames
import FreqTables
import Distributions

using JLD2
using Random
using LinearAlgebra
using Statistics
using CSV
using DataFrames
using FreqTables
using Distributions

#### Item 1 ####

function q1()

    # a)
    Random.seed!(1234)

    # i.
    A = rand(Uniform(-5,10),10,7)
    # ii.
    B = rand(Normal(-2,15),10,7)
    # iii.
    C = [A[1:5,1:5] B[1:5,6:7]]
    # iv,
    D = A
    D[D .> 0] .= 0
    D

    # b)
    length(A)

    # c)
    length(unique(A))

    # d)
    E = reshape(B,(70,1))

    ## I did not know an easier way to accomplish the vectorization,
    ## without having to use a loop.

    # e)
    F = cat(A,B, dims = 3)

    # f)
    F = permutedims(F,(3,1,2))

    # g)
    G = kron(B,C)
    #G2 = kron(C,F)

    ## Trying to use a Kronecker product between C and F turns a error due
    ## to different dimensions.

    # h)
    @save "matrixpractice.jld" A,B,C,D,E,F,G

    # i)
    @save "firstmatrix.jld" A,B,C,D

    # j)
    C = DataFrame(C)
    CSV.write("Cmatrix.csv",C)

    # k)
    D = DataFrame(D)
    CSV.write("Dmatrix.dat",D;delim="   ")

    return A,B,C,D

end;

A,B,C,D = q1()

# l)
q1()

#### Item 2 ####

function q2(A,B,C)

    # a)
    AB = [x*y for (x,y) in zip(A, B)]

    AB2 = A .* B

    # b)
    Cprime = Vector{Float64}()
    for i in C
        if  -5 <= i  <= 5.0
            append!(Cprime, i)
        end
    end

    Cprime2 = C[-5 .<= C .<= 5]

    # c)
    N = 15169
    K = 6
    T = 5
    X = zeros(N,K,T)
    for t in 1:T
        X[:,1,t] .= ones(N) # stationary over time.
        X[:,2,t] .= rand(Bernoulli(.75*(6-t)/5), N)
        X[:,3,t] .= rand(Normal(15+t-1, 5*(t-1)), N)
        X[:,4,t] .= rand(Normal(π*(6-t)/3, 1/ℯ), N)
        X[:,5,t] .= rand(Binomial(20, 0.6), N)
        X[:,6,t] .= rand(Binomial(20,0.5),N) # Discrete normal stationary over time.
    end

    # d)

    β = zeros(K,T)
    β[1,:] .= [1 + t*(.25) for t in 1:T] # d1
    β[2,:] .= [log(t) for t in 1:T] # d2
    β[3,:] .= [-sqrt(t) for t in 1:T] # d3
    β[4,:] .= [(ℯ^(t)  - ℯ^(t+1)) for t in 1:T] # d4
    β[5,:] .= [t for t in 1:T] # d5
    β[6,:] .= [t/3 for t in 1:T] # d6

    # e)
    Y = zeros(N,T)
    for t in 1:T # How to use comprehensions in this case?
       Y[:,t] .= X[:,:,t]*β[:,t] + rand(Normal(0, .36), N)
    end
end

q2()

# 3
function q3()
    # a)
    df = CSV.read("./PS1-julia-intro/nlsw88.csv")# |> dropmissing
    save("./PS1-julia-intro/nlsw88.jld.jld", "df", df)
    # b)
    println((df.never_married |> sum)*100/nrow(df) |> round,"% sample has never been married")
    println((df.collgrad |> sum)*100/nrow(df) |> round,"% are college graduates")

    # c)
    freqtable(df.race) |> prop

    # d)
    summarystats = describe(df)
    println(df.grade .|> ismissing |> sum, " grades observations are missing")

    # e)
    freqtable(df.occupation, df.industry) |> prop

    # f)
    gdf = groupby(df, ["occupation","industry"], sort=true)
    combine(gdf, :wage => mean => :avg_ind_occ_wage) |> dropmissing
    return nothing
end

q3()

#### Item 4 ####

# a)

function q4()
    load("./PS1-julia-intro/firstmatrix.jld")

    function matrixops(m1,m2)
        if size(m1) != size(m2)
            error("inputs must have the same size.")
        end
        #=
        Takes as inputs the matrices A and B from question (a) of problem 1 and has three outputs:
            (i) the element-by-element product of the inputs,
            (ii) the product A'B, and
            (iii) the sum of all the elements of A + B.
        =#
        return m1 .* m2, m1' * m2, m1 + m2
    end

    a,b, c = matrixops(A, B);

    df = CSV.read("./PS1-julia-intro/nlsw88.csv") |> dropmissing
    convert(Array,df.ttl_exp); convert(Array,df.wage);
    matrixops(
        df.ttl_exp,
        df.wage);
    return nothing
end

q4()


text:
module SolarizedLogging
using Logging

function __init__()
    # Only modify console loggers
    cur_logger = global_logger()
    typeof(cur_logger) == Logging.ConsoleLogger || return

    new_logger = ConsoleLogger(cur_logger.stream, cur_logger.min_level;
        meta_formatter=solarized_metafmt, show_limited=cur_logger.show_limited,
        right_justify=cur_logger.right_justify)

    global_logger(new_logger)
end

function solarized_logcolor(level)
    level < Logging.Info  ? 4 :
    level < Logging.Warn  ? 6 :
    level < Logging.Error ? 3 :
                            1
end

function solarized_metafmt(level, _module, group, id, file, line)
    prefix_color = solarized_logcolor(level)
    prefix = (level == Logging.Warn ? "Warning" : string(level))*':'
    suffix_color = 10
    suffix = ""
    Logging.Info <= level < Logging.Warn && return prefix_color, prefix, suffix_color, suffix
    _module !== nothing && (suffix *= "$(_module)")
    if file !== nothing
        _module !== nothing && (suffix *= " ")
        suffix *= Base.contractuser(file)
        if line !== nothing
            suffix *= ":$(isa(line, UnitRange) ? "$(first(line))-$(last(line))" : line)"
        end
    end
    !isempty(suffix) && (suffix = "@ " * suffix)
    return prefix_color, prefix, suffix_color, suffix
end

end # module


text:
using Random
using POMDPSimulators
using PMDPs.LP

@testset "HistoryReplayer.jl" begin
    mg, me = simple_mdps()
    pp = PMDPs.pp(mg)

    # History replayer takes either full SimulationHistory structure
    simhistory = PMDPs.simulate_trace(mg, MersenneTwister(1))
    # Or more limited structure containing Abstract Arrays of NamedTuples that contain state and info fields.
    trace = simple_short_trace(pp)

    @test isa(PMDPs.HistoryReplayer(mg, simhistory), PMDPs.PMDP)
    @test isa(PMDPs.HistoryReplayer(mg, trace), PMDPs.PMDP)

    hrpl = PMDPs.HistoryReplayer(mg, trace)

    # Test some basic properties
    @test POMDPs.actions(hrpl) == POMDPs.actions(mg)

    # Test gen method
    s = trace[1].s
    (sp, r, info) = POMDPs.gen(hrpl, s, 1000., MersenneTwister(99))
    @test sp == trace[2].s
    @test r == 0.
    @test info == trace[1].info

    # s = PMDPs.State(SA[2,3], 5, SA[true, true])

    # test that the new trace from replayer matches the input trace
    hrpl_sim = PMDPs.HistoryReplayer(mg, simhistory)
    t = PMDPs.simulate_trace(hrpl_sim, MersenneTwister(321))
    @test t == simhistory

    # TODO: The following should likely be tested in tests of the specific policies

    hrec = HistoryRecorder(max_steps = PMDPs.selling_period_end(mg), rng = MersenneTwister(4321))

    # test VI with HistoryReplayer
    policy = PMDPs.get_VI_policy(me)
    hᵥ = simulate(hrec, hrpl, policy)
    @test length(hᵥ) == length(trace)
    @test sum(collect(hᵥ[:r])) > 0.

    # test MCTS with historyReplayer
    planner = PMDPs.get_MCTS_planner(mg; params_mcts = Dict(:rng=>MersenneTwister(1)))
    hₘ = simulate(hrec, hrpl, planner)
    @test length(hₘ) == length(trace)
    @test sum(collect(hₘ[:r])) > 0.

    # test hindsight with historyReplayer
    hindsight = PMDPs.LP.get_MILP_hindsight_policy(mg, trace)
    hₕ = simulate(hrec, hrpl, hindsight)
    @test length(hₕ) == length(trace)
    @test sum(collect(hₕ[:r])) > 0.

    # test flatrate with historyReplayer
    # R, U = PMDPs.optimize_flatrate_policy(mg, [trace, trace])
    flatrate = PMDPs.get_flatrate_policy(mg, [trace, trace])
    hᵣ = simulate(hrec, hrpl, flatrate)
    @test length(hᵣ) == length(trace)
    @test sum(collect(hᵣ[:r])) > 0.

end


text:
using ChainRulesTestUtils: rand_tangent

# Test struct for `rand_tangent`.
struct Foo
   a::Float64
   b::Int
   c::Any
end

@testset "generate_tangent" begin
    rng = MersenneTwister(123456)

    foreach([
        ("hi", DoesNotExist),
        ('a', DoesNotExist),
        (:a, DoesNotExist),
        (true, DoesNotExist),
        (4, DoesNotExist),
        (5.0, Float64),
        (5.0 + 0.4im, Complex{Float64}),
        (randn(Float32, 3), Vector{Float32}),
        (randn(Complex{Float64}, 2), Vector{Complex{Float64}}),
        (randn(5, 4), Matrix{Float64}),
        (randn(Complex{Float32}, 5, 4), Matrix{Complex{Float32}}),
        ([randn(5, 4), 4.0], Vector{Any}),
        ((4.0, ), Composite{Tuple{Float64}}),
        ((5.0, randn(3)), Composite{Tuple{Float64, Vector{Float64}}}),
        ((a=4.0, ), Composite{NamedTuple{(:a,), Tuple{Float64}}}),
        ((a=5.0, b=1), Composite{NamedTuple{(:a, :b), Tuple{Float64, Int}}}),
        (sin, typeof(NO_FIELDS)),
        (Foo(5.0, 4, rand(rng, 3)), Composite{Foo}),
        (Foo(4.0, 3, Foo(5.0, 2, 4)), Composite{Foo}),
    ]) do (x, T_tangent)
        @test rand_tangent(rng, x) isa T_tangent
        @test rand_tangent(x) isa T_tangent
        @test x + rand_tangent(rng, x) isa typeof(x)
    end

    # Ensure struct fallback errors for non-struct types.
    @test_throws ArgumentError invoke(rand_tangent, Tuple{AbstractRNG, Any}, rng, 5.0)
end


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

index(A, i) = A[i]
emit(index, Array{Int64, 1}, Int64)
#

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.index), Vector{Int64}, Int64}, svec(), index(A, i) in Main.Main at /home/mccoy/Dev/brutus/test/Codegen/lower/index_linear_1D.jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.index), Array{Int64, 1}, Int64}"(%arg0: !jlir<"typeof(Main.index)">, %arg1: !jlir<"Array{Int64, 1}">, %arg2: !jlir.Int64) -> !jlir.Int64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir.Core.arrayref} : () -> !jlir<"typeof(Core.arrayref)">
# CHECK-NEXT:     %1 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %2 = "jlir.call"(%0, %1, %arg1, %arg2) : (!jlir<"typeof(Core.arrayref)">, !jlir.Bool, !jlir<"Array{Int64, 1}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     "jlir.return"(%2) : (!jlir.Int64) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.index), Array{Int64, 1}, Int64}"(%arg0: !jlir<"typeof(Main.index)">, %arg1: memref<?xi64>, %arg2: i64) -> i64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %c1 = constant 1 : index
# CHECK-NEXT:     %0 = "jlir.convertstd"(%arg2) : (i64) -> index
# CHECK-NEXT:     %1 = subi %0, %c1 : index
# CHECK-NEXT:     %2 = load %arg1[%1] : memref<?xi64>
# CHECK-NEXT:     return %2 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.index), Array{Int64, 1}, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: !llvm.ptr<i64>, %arg2: !llvm.ptr<i64>, %arg3: i64, %arg4: i64, %arg5: i64, %arg6: i64) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.mlir.undef : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %1 = llvm.insertvalue %arg1, %0[0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %2 = llvm.insertvalue %arg2, %1[1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %3 = llvm.insertvalue %arg3, %2[2] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %4 = llvm.insertvalue %arg4, %3[3, 0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %5 = llvm.insertvalue %arg5, %4[4, 0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %6 = llvm.mlir.constant(1 : index) : i64
# CHECK-NEXT:     %7 = llvm.sub %arg6, %6  : i64
# CHECK-NEXT:     %8 = llvm.extractvalue %5[1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %9 = llvm.getelementptr %8[%7] : (!llvm.ptr<i64>, i64) -> !llvm.ptr<i64>
# CHECK-NEXT:     %10 = llvm.load %9 : !llvm.ptr<i64>
# CHECK-NEXT:     llvm.return %10 : i64
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.index), Array{Int64, 1}, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: !llvm.ptr<struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>>, %arg2: i64) -> i64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.load %arg1 : !llvm.ptr<struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>>
# CHECK-NEXT:     %1 = llvm.extractvalue %0[0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %2 = llvm.extractvalue %0[1] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %3 = llvm.extractvalue %0[2] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %4 = llvm.extractvalue %0[3, 0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %5 = llvm.extractvalue %0[4, 0] : !llvm.struct<(ptr<i64>, ptr<i64>, i64, array<1 x i64>, array<1 x i64>)>
# CHECK-NEXT:     %6 = llvm.call @"Tuple{typeof(Main.index), Array{Int64, 1}, Int64}"(%arg0, %1, %2, %3, %4, %5, %arg2) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, !llvm.ptr<i64>, !llvm.ptr<i64>, i64, i64, i64, i64) -> i64
# CHECK-NEXT:     llvm.return %6 : i64
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: JuliaLabs/brutus

max_issues_repo_name: JuliaLabs/brutus

max_forks_repo_name: JuliaLabs/brutus


text:
@testset "SDP - scalar parameter" begin
    cached = MOI.Bridges.full_bridge_optimizer(
        MOIU.CachingOptimizer(
            MOIU.UniversalFallback(MOIU.Model{Float64}()),
            SCS.Optimizer(),
        ),
        Float64,
    )
    optimizer = POI.Optimizer(cached)

    m = direct_model(optimizer)
    set_silent(m)

    @variable(m, p in POI.Parameter(0))
    @variable(m, x[1:2, 1:2], Symmetric)
    @objective(m, Min, tr(x))
    @constraint(m, con, Symmetric(x .- [1+p 0; 0 1+p]) in PSDCone())

    optimize!(m)
    @test all(isapprox.(value.(x), [1 0; 0 1], atol = ATOL))

    MOI.set(m, POI.ParameterValue(), p, 1)
    optimize!(m)

    @test all(isapprox.(value.(x), [2 0; 0 2], atol = ATOL))
end

@testset "SDP - Matrix parameter" begin
    cached = MOI.Bridges.full_bridge_optimizer(
        MOIU.CachingOptimizer(
            MOIU.UniversalFallback(MOIU.Model{Float64}()),
            SCS.Optimizer(),
        ),
        Float64,
    )
    optimizer = POI.Optimizer(cached)

    m = direct_model(optimizer)
    set_silent(m)

    P1 = [1 2; 2 3]
    @variable(m, p[1:2, 1:2] in POI.Parameter.(P1))
    @variable(m, x[1:2, 1:2], Symmetric)
    @objective(m, Min, tr(x))
    @constraint(m, con, Symmetric(x - p) in PSDCone())

    optimize!(m)

    @test all(isapprox.(value.(x), P1, atol = ATOL))

    P2 = [1 2; 2 1]
    MOI.set.(m, POI.ParameterValue(), p, P2)
    optimize!(m)

    @test all(isapprox.(value.(x), P2, atol = ATOL))
end


text:
# Partition mesh using a static domain decomposition algorithm based on leaf cell count alone
function partition!(mesh)
  # Determine number of leaf cells per rank
  leaves = leaf_cells(mesh.tree)
  @assert length(leaves) > mpi_nranks()
  n_leaves_per_rank = OffsetArray(fill(div(length(leaves), mpi_nranks()), mpi_nranks()),
                                  0:(mpi_nranks() - 1))
  for d in 0:(rem(length(leaves), mpi_nranks()) - 1)
    n_leaves_per_rank[d] += 1
  end
  @assert sum(n_leaves_per_rank) == length(leaves)

  # Assign MPI ranks to all cells such that all ancestors of each cell - if not yet assigned to a
  # rank - belong to the same rank
  mesh.first_cell_by_rank = similar(n_leaves_per_rank)
  mesh.n_cells_by_rank = similar(n_leaves_per_rank)

  leaf_count = 0
  last_id = leaves[n_leaves_per_rank[0]]
  mesh.first_cell_by_rank[0] = 1
  mesh.n_cells_by_rank[0] = last_id
  mesh.tree.mpi_ranks[1:last_id] .= 0
  for d in 1:(length(n_leaves_per_rank)-1)
    leaf_count += n_leaves_per_rank[d-1]
    last_id = leaves[leaf_count + n_leaves_per_rank[d]]
    mesh.first_cell_by_rank[d] = mesh.first_cell_by_rank[d-1] + mesh.n_cells_by_rank[d-1]
    mesh.n_cells_by_rank[d] = last_id - mesh.first_cell_by_rank[d] + 1
    mesh.tree.mpi_ranks[mesh.first_cell_by_rank[d]:last_id] .= d
  end

  @assert all(x->x >= 0, mesh.tree.mpi_ranks[1:length(mesh.tree)])
  @assert sum(mesh.n_cells_by_rank) == length(mesh.tree)

  return nothing
end

function load_mesh(restart_filename, mpi_parallel::Val{true})
  # Get number of spatial dimensions
  ndims_ = parameter("ndims")

  # Get maximum number of cells that should be supported
  n_cells_max = parameter("n_cells_max")

  # Create mesh
  @timeit timer() "creation" mesh = TreeMesh(ParallelTree{ndims_}, n_cells_max)

  # Determine mesh filename
  filename = get_restart_mesh_filename(restart_filename, Val(true))
  mesh.current_filename = filename
  mesh.unsaved_changes = false

  # Read mesh file
  if mpi_isroot()
    h5open(filename, "r") do file
      # Set domain information
      mesh.tree.center_level_0 = read(attrs(file)["center_level_0"])
      mesh.tree.length_level_0 = read(attrs(file)["length_level_0"])
      mesh.tree.periodicity    = Tuple(read(attrs(file)["periodicity"]))
      MPI.Bcast!(collect(mesh.tree.center_level_0), mpi_root(), mpi_comm())
      MPI.Bcast!(collect(mesh.tree.length_level_0), mpi_root(), mpi_comm())
      MPI.Bcast!(collect(mesh.tree.periodicity),    mpi_root(), mpi_comm())

      # Set length
      n_cells = read(attrs(file)["n_cells"])
      MPI.Bcast!(Ref(n_cells), mpi_root(), mpi_comm())
      resize!(mesh.tree, n_cells)

      # Read in data
      mesh.tree.parent_ids[1:n_cells] = read(file["parent_ids"])
      mesh.tree.child_ids[:, 1:n_cells] = read(file["child_ids"])
      mesh.tree.neighbor_ids[:, 1:n_cells] = read(file["neighbor_ids"])
      mesh.tree.levels[1:n_cells] = read(file["levels"])
      mesh.tree.coordinates[:, 1:n_cells] = read(file["coordinates"])
      @views MPI.Bcast!(mesh.tree.parent_ids[1:n_cells],      mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.child_ids[:, 1:n_cells],    mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.neighbor_ids[:, 1:n_cells], mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.levels[1:n_cells],          mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.coordinates[:, 1:n_cells],  mpi_root(), mpi_comm())
    end
  else # non-root ranks
    # Set domain information
    mesh.tree.center_level_0 = MPI.Bcast!(collect(mesh.tree.center_level_0), mpi_root(), mpi_comm())
    mesh.tree.length_level_0 = MPI.Bcast!(collect(mesh.tree.length_level_0), mpi_root(), mpi_comm())[1]
    mesh.tree.periodicity    = Tuple(MPI.Bcast!(collect(mesh.tree.periodicity),    mpi_root(), mpi_comm()))

    # Set length
    n_cells = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
    resize!(mesh.tree, n_cells)

    # Read in data
    @views MPI.Bcast!(mesh.tree.parent_ids[1:n_cells],      mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.child_ids[:, 1:n_cells],    mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.neighbor_ids[:, 1:n_cells], mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.levels[1:n_cells],          mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.coordinates[:, 1:n_cells],  mpi_root(), mpi_comm())
  end

  # Partition mesh
  partition!(mesh)

  return mesh
end

function get_restart_mesh_filename(restart_filename, mpi_parallel::Val{true})
  # Get directory name
  dirname, _ = splitdir(restart_filename)

  if mpi_isroot()
    # Read mesh filename from restart file
    mesh_file = ""
    h5open(restart_filename, "r") do file
      mesh_file = read(attrs(file)["mesh_file"])
    end

    buffer = Vector{UInt8}(mesh_file)
    MPI.Bcast!(Ref(length(buffer)), mpi_root(), mpi_comm())
    MPI.Bcast!(buffer, mpi_root(), mpi_comm())
  else # non-root ranks
    count = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())
    buffer = Vector{UInt8}(undef, count[])
    MPI.Bcast!(buffer, mpi_root(), mpi_comm())
    mesh_file = String(buffer)
  end

  # Construct and return filename
  return joinpath(dirname, mesh_file)
end


text:
using BinaryBuilder

name = "Jansson"
version = v"2.14"

# Collection of sources required to build ZMQ
sources = [
    GitSource("https://github.com/akheron/jansson.git", "684e18c927e89615c2d501737e90018f4930d6c5"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/jansson

autoreconf -fi
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libjansson", :libjansson),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
export run_unittfopf, run_ac_unittfopf, run_dc_unittfopf

""
function run_ac_unittfopf(data::Dict{String,Any}, solver; kwargs...)
    return run_unittfopf(data, ACPPowerModel, solver; kwargs...)
end

""
function run_dc_unittfopf(data::Dict{String,Any}, solver; kwargs...)
    return run_unittfopf(data, DCPPowerModel, solver; kwargs...)
end

""
function run_unittfopf(data::Dict{String,Any}, model_constructor, solver; kwargs...)
    process_additional_data!(data)
    pm = PowerModels.build_generic_model(data, model_constructor, post_unittfopf; kwargs...)
    return PowerModels.solve_generic_model(pm, solver; solution_builder = get_solution_tf)
end

""
function post_unittfopf(pm::GenericPowerModel)
    add_load_model!(pm) # To add load data
    PowerModels.variable_voltage(pm)
    PowerModels.variable_generation(pm)
    PowerModels.variable_branch_flow(pm)
    PowerModels.variable_dcline_flow(pm)
    variable_transformation(pm)
    variable_node_aggregation(pm)
    variable_load(pm)
    variable_action_indicator(pm)
    variable_auxiliary_power(pm)

    objective_min_redispatch_cost(pm)

    PowerModels.constraint_voltage(pm)
    for i in PowerModels.ids(pm, :ref_buses)
        PowerModels.constraint_theta_ref(pm, i)
    end

    for i in PowerModels.ids(pm, :bus)
        constraint_kcl_shunt_aggregated(pm, i)
        constraint_load_gen_aggregation_sheddable(pm, i)
    end

    for i in PowerModels.ids(pm, :load)
        constraint_flexible_load(pm, i)
        constraint_redispatch_power_load(pm, i)
    end

    for i in PowerModels.ids(pm, :gen)
        constraint_flexible_gen(pm, i)
        constraint_redispatch_power_gen(pm, i)
    end

    for i in PowerModels.ids(pm, :branch)
        branch = PowerModels.ref(pm, :branch, i)

        if branch["shiftable"] == false && branch["tappable"] == false
            PowerModels.constraint_ohms_yt_from(pm, i)
            PowerModels.constraint_ohms_yt_to(pm, i)
            constraint_link_voltage_magnitudes(pm, i)
        else
            constraint_variable_transformer_y_from(pm, i)
            constraint_variable_transformer_y_to(pm, i)
        end
        PowerModels.constraint_voltage_angle_difference(pm, i)

        PowerModels.constraint_thermal_limit_from(pm, i)
        PowerModels.constraint_thermal_limit_to(pm, i)
    end
    for i in PowerModels.ids(pm, :dcline)
        PowerModels.constraint_dcline(pm, i)
    end
end

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
include(joinpath("gammaReg", "chosenVariables_inverse_test.jl"))
include(joinpath("gammaReg", "chosenVariables_log_test.jl"))
include(joinpath("gammaReg", "research_inverse_test.jl"))
include(joinpath("gammaReg", "research_inverse_test.jl"))

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
"""
CORENUMS
---------
    compute the core number for each vertex in the graph and returns the core
    numbers for each vertex of the graph A along with the removal order of the vertex in the
    tuple (d,rt). This function works on directed graphs but gives the in-degree core number.
    To get the out-degree core numbers call corenums(A')

Functions
---------
- (d,rt) = corenums(A::MatrixNetwork)
- (d,rt) = corenums{T}(A::SparseMatrixCSC{T,Int64})
- (d,rt) = corenums(ei::Vector{Int64},ej::Vector{Int64})

Example
-------
~~~
A = load_matrix_network("cores_example")
(d,rt) = corenums(A)
~~~
"""

function corenums(A::MatrixNetwork)
    (rp,ci) = (A.rp,A.ci)
    n=length(rp)-1
    nz=length(ci)

    d=zeros(Int64,n)
    maxd=0
    rt=zeros(Int64,n)
    for k=1:nz
        newd=d[ci[k]]+1
        d[ci[k]]=newd
        if newd>maxd
            maxd=newd
        end
    end

    # compute the bucket sort
    dp=zeros(Int64,maxd+2)
    vs=zeros(Int64,n)
    vi=zeros(Int64,n)# degree position, vertices
    for i=1:n
        dp[d[i]+2]=dp[d[i]+2]+1
    end # plus 2 because degrees start at 0

    dp=cumsum(dp)
    dp=dp+1
    for i=1:n
        vs[dp[d[i]+1]]=i
        vi[i]=dp[d[i]+1]
        dp[d[i]+1]=dp[d[i]+1]+1
    end

    for i=maxd:-1:1
        dp[i+1]=dp[i]
    end

    # start the algorithm
    t=1
    for i=1:n
        v = vs[i]
        dv = d[v]
        rt[v]=t
        t=t+1

        for rpi=rp[v]:rp[v+1]-1
            w=ci[rpi]
            dw=d[w]

            if dw<=dv # we already removed w
            else #need to remove edge (v,w), which decreases d(w)
                # swap w with the vertex at the head of its degree
                pw=vi[w] # get the position of w
                px=dp[dw+1] #get the pos of the vertex at the head of dw list
                x=vs[px]
                #swap w, x
                vs[pw]=x
                vs[px]=w
                vi[w]=px
                vi[x]=pw

                # decrement the degree of w and increment the start of dw
                d[w]=dw-1
                dp[dw+1]=px+1
            end
        end
    end
    return (d,rt)
end

############################
### Additional functions ###
############################

## CSC:
function corenums{T}(A::SparseMatrixCSC{T,Int64})
    return corenums(MatrixNetwork(A))
end

## Triplet Format:
function corenums(ei::Vector{Int64},ej::Vector{Int64})
    return corenums(MatrixNetwork(ei,ej))
end

## CSR sparse matrices:
function corenums{T}(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64)
    return corenums(MatrixNetwork(n,rp,ci,vals))
end


text:
"""
	ImplicitTrajectory{T}
This structure holds the trajectory of evaluations and Jacobians of the implicit dynamics.
These evaluations and Jacobians are computed using a linearizedimation computed around `lin`.
"""
mutable struct ImplicitTrajectory{T,R,RZ,Rθ,NQ}
	H::Int
	lin::Vector{LinearizedStep{T}}
	d::Vector{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int}},true}} # dynamics violation
	dq2::Vector{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int}},true}}
	dγ1::Vector{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int}},true}}
	db1::Vector{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int}},true}}
	δq0::Vector{SubArray{T,2,Array{T,2},Tuple{UnitRange{Int},UnitRange{Int}},false}} # q0 solution gradient length=H
	δq1::Vector{SubArray{T,2,Array{T,2},Tuple{UnitRange{Int},UnitRange{Int}},false}}  # q1 solution gradient length=H
	δu1::Vector{SubArray{T,2,Array{T,2},Tuple{UnitRange{Int},UnitRange{Int}},false}}  # u1 solution gradient length=H
	ip::Vector{InteriorPoint{T,R,RZ,Rθ}}
	mode::Symbol
	iq2::SVector{NQ,Int}
end

function ImplicitTrajectory(ref_traj::ContactTraj, s::Simulation;
	κ = ref_traj.κ[1],
	max_time = 1e5,
	mode = :configurationforce,
	opts = InteriorPointOptions(
			undercut = 5.0,
			γ_reg = 0.1,
			κ_tol = κ[1],
			r_tol = 1.0e-8,
			diff_sol = true,
			solver = :empty_solver,
			max_time = max_time))

	model = s.model
	env = s.env

	iq2 = index_q2(model, env, quat = false)
	nq2 = length(iq2)

	H = ref_traj.H

	nq = model.nq
	nu = model.nu
	nc = model.nc
	nb = nc * friction_dim(env)
	if mode == :configurationforce
		nd = nq + nc + nb
	elseif mode == :configuration
		nd = nq
	else
		@error "invalid mode"
	end
	nz = num_var(model, env)
	nθ = num_data(model)

	lin = [LinearizedStep(s, ref_traj.z[t], ref_traj.θ[t], κ) for t = 1:H]

	ip =  [interior_point(
			 deepcopy(ref_traj.z[t]),
			 deepcopy(ref_traj.θ[t]),
			 idx = IndicesOptimization(model, env),
			 r! = rlin!,
			 rz! = rzlin!,
			 rθ! = rθ!,
			 r  = RLin(s, lin[t].z, lin[t].θ, lin[t].r, lin[t].rz, lin[t].rθ),
			 rz = RZLin(s, lin[t].rz),
			 rθ = RθLin(s, lin[t].rθ),
			 opts = opts) for t = 1:H]

	# views
	d = [view(ip[t].z, 1:nd) for t = 1:H]
	dq2 = [view(ip[t].z, 1:nq) for t = 1:H]
	if mode == :configurationforce
		dγ1 = [view(ip[t].z, nq .+ (1:nc)) for t = 1:H]
		db1 = [view(ip[t].z, nq + nc .+ (1:nb)) for t = 1:H]
	elseif mode == :configuration
		dγ1 = [view(ip[t].z, 1:0) for t = 1:H]
		db1 = [view(ip[t].z, 1:0) for t = 1:H]
	else
		@error "invalid mode"
	end

	off = 0
	δq0 = [view(ip[t].δz, 1:nd, off .+ (1:nq)) for t = 1:H]; off += nq
	δq1 = [view(ip[t].δz, 1:nd, off .+ (1:nq)) for t = 1:H]; off += nq
	δu1 = [view(ip[t].δz, 1:nd, off .+ (1:nu)) for t = 1:H]; off += nu

	return ImplicitTrajectory{typeof.([ip[1].z[1], ip[1].r, ip[1].rz, ip[1].rθ])...,nq2}(H, lin, d, dq2, dγ1, db1, δq0, δq1, δu1, ip, mode,
		SVector{nq2,Int}(iq2))
end

function update!(im_traj::ImplicitTrajectory{T,R,RZ,Rθ,nq}, ref_traj::ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ},
	s::Simulation{T,W,FC}, alt::Vector{T}, κ::T, H::Int) where {T,R,RZ,Rθ,nq,nu,nw,nc,nb,nz,nθ,W,FC}

	for t = 1:H
		# central-path parameter
		im_traj.ip[t].κ[1] = κ

		t > H-1 && continue

		# # residual
		# im_traj.lin[t] = im_traj.lin[t+1]
		# im_traj.ip[t].r = im_traj.ip[t+1].r
		# im_traj.ip[t].rz = im_traj.ip[t+1].rz
		# im_traj.ip[t].rθ = im_traj.ip[t+1].rθ

		# # altitude
		# im_traj.ip[t].r.alt = alt
	end

	# update!(im_traj.lin[H], s, ref_traj.z[H], ref_traj.θ[H])

	# z0  = im_traj.lin[H].z
	# θ0  = im_traj.lin[H].θ
	# r0  = im_traj.lin[H].r
	# rz0 = im_traj.lin[H].rz
	# rθ0 = im_traj.lin[H].rθ

	# update!(im_traj.ip[H].r, z0, θ0, r0, rz0, rθ0)
	# update!(im_traj.ip[H].rz, rz0)
	# update!(im_traj.ip[H].rθ, rθ0)

	# # altitude
	# im_traj.ip[H].r.alt = alt
	return nothing
end

function set_implicit_trajectory!(im_traj::ImplicitTrajectory, im_traj_cache::ImplicitTrajectory)
	H = im_traj.H

	for t = 1:H
		im_traj.lin[t] = im_traj_cache.lin[t]
		im_traj.ip[t].r = im_traj_cache.ip[t].r
		im_traj.ip[t].rz = im_traj_cache.ip[t].rz
		im_traj.ip[t].rθ = im_traj_cache.ip[t].rθ
	end

	return nothing
end

function set_altitude!(im_traj::ImplicitTrajectory, alt::Vector)
	H = im_traj.H
	for t = 1:H
		# altitude
		set_altitude!(im_traj.ip[t], alt)
	end
	return nothing
end

function set_altitude!(ip::InteriorPoint, alt::Vector)
	# altitude
	ip.r.alt = alt
	return nothing
end

function implicit_dynamics!(im_traj::ImplicitTrajectory, traj::ContactTraj;
	threads=false,
	window=collect(1:traj.H + 2))

	for (i, t) in enumerate(window[1:end-2])
		# initialized solver
		z_initialize!(im_traj.ip[t].z, im_traj.iq2, traj.q[i+2]) #TODO: try alt. schemes
		im_traj.ip[t].θ .= traj.θ[i]
	end

	if threads
		Threads.@threads for t in window[1:end-2]
			# solve
			status = interior_point_solve!(im_traj.ip[t])
			!status && (@warn "implicit dynamics failure (t = $t)")
		end
	else
		for t in window[1:end-2]
			# solve
			status = interior_point_solve!(im_traj.ip[t])
			!status && (@warn "implicit dynamics failure (t = $t)")
		end
	end

	for (i, t) in enumerate(window[1:end-2])
		# compute dynamics violation
		im_traj.dq2[t] .-= traj.q[i+2]

		if im_traj.mode == :configurationforce
			im_traj.dγ1[t] .-= traj.γ[i]
			im_traj.db1[t] .-= traj.b[i]
		elseif im_traj.mode == :configuration
			nothing
		end
	end
	return nothing
end


text:
module Suri
using Hecke

#= follows Sebastian Posur's idea
  Start with a pseudo_matrix with ideals A_i and rows alpha_i
  plus an additional pseudo-element      A            alpha
  assume, at least for now, the alpha_i are independent.
  assume sum gamma_i alpha_i + gamma alpha = 0
  (gamma are the dependency)
  Now P :=
  (inv(A_i) gamma_i)
  ( inv(A)  gamma  )
  is a (n+1) x 1 pseudo matrix. Using the (pseudo) HNF we find a matrix
  T in Gl(n+1, K) and ideals C_i s.th.
  T*P = H is a pseudo-HNF. By construction, H has 1. row 1, all other rows 0
  and the coeff. ideals C_i
  Then T is an isomorphism between
  sum inv(A_i) + inv(A) and sum C_i
  as zk-modules
  The inverse-transpose of T is an isomorphism between
  sum inv(C_i) and sum A_i + A
  the dual modules to the ones above

  The idea:
  T is an iso between all modules given by the ideals
  The construction forces inv(T) to have the relation as 1. row/col
  The dual of an abstract module is giving by the inverse ideals
  The corresponding map is transpose(inv(T)) = inv(T)'
  Thus it works...
=#

function extend(M::Hecke.PMat, b::Generic.MatSpaceElem{nf_elem}, gamma::Generic.MatSpaceElem{nf_elem})

  zk = base_ring(M)
  nc = ncols(gamma)
  @assert nc == ncols(b)
  p = pseudo_matrix(gamma, vcat(map(inv, coefficient_ideals(M)), [1*zk for i=1:nc]))
  h, T = Hecke.pseudo_hnf_with_transform(p)
  #for n x 1 matrices, the trasform, especially the inverse can more easily be computed
#  @assert prod(coefficient_ideals(p)) ==  (det(T))*(prod(coefficient_ideals(h)))
#  @assert all(T[j,i] in (coefficient_ideals(p)[i])*inv(coefficient_ideals(h)[j]) for i=1:nrows(M) for j=1:nrows(M))
  Ti = inv(T)
#  @assert all(Ti[i,j] in inv(coefficient_ideals(p)[i])*(coefficient_ideals(h)[j]) for i=1:nrows(M) for j=1:nrows(M))
  e = pseudo_matrix((Ti'*hcat(M.matrix', b)')[nc+1:nrows(M)+nc, :], map(inv, coefficient_ideals(h)[nc+1:nrows(M)+nc]))
  return e
end

function Hecke.denominator(P::Hecke.PMat, M::NfOrd)
  l = fmpz(1)
  p = matrix(P)
  for i=1:nrows(P)
    I = coefficient_ideals(P)[i]
    Hecke.assure_2_normal(I.num)
    for j=1:ncols(P)
      l = lcm(l, denominator(simplify(p[i,j]*I)))
    end
  end
  return l
end

function Hecke.simplify(P::Hecke.PMat)
  c = copy(coefficient_ideals(P))
  m = deepcopy(matrix(P))
  for i=1:length(c)
    a, b = Hecke.reduce_ideal2(c[i])
    m[i, :] *= inv(b)
    c[i] = fractional_ideal(order(a), a)
  end
  return pseudo_matrix(m, c)
end

function Hecke.dual(P::Hecke.PMat)
  return pseudo_matrix(inv(P.matrix)', map(inv, coefficient_ideals(P)))
end

end

#= example

Hecke.example("Suri.jl")
Hecke.revise("Suri.jl")

k, a = wildanger_field(3, 13)
m = rand(MatrixSpace(k, 10, 10), 1:10);
m = cat(m,m, dims=(1,2));
b = rand(MatrixSpace(k, 20, 1), 1:10);
S = kernel(hcat(m, b));

m1 = Suri.extend(pseudo_matrix(m'), b, S[2]);

norm(det(m))
norm(det(m1))

b = rand(MatrixSpace(k, 20, 1), 1:10);
S = kernel(hcat(m1.matrix', b));
m2 = Suri.extend(m1, b, S[2]);

norm(det(m2))

=#


text:
#=
The code below was derived from the scipy implementation of the LOBPCG algorithm in https://github.com/scipy/scipy/blob/v1.1.0/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py#L109-L568.

Since the link above mentions the license is BSD license, the notice for the BSD license 2.0 is hereby given below giving credit to the authors of the Python implementation.

Copyright (c) 2018, Robert Cimrman, Andrew Knyazev
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * The names of the contributors may not be used to endorse or promote products
      derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

struct Blocks{Generalized, T, TA<:AbstractArray{T}}
    block::TA # X, R or P
    A_block::TA # AX, AR or AP
    B_block::TA # BX, BR or BP
end
Blocks(X, AX) = Blocks{false, eltype(X), typeof(X)}(X, AX, X)
Blocks(X, AX, BX) = Blocks{true, eltype(X), typeof(X)}(X, AX, BX)
function A_mul_X!(b::Blocks, A)
    A_mul_B!(b.A_block, A, b.block)
    return
end
function A_mul_X!(b::Blocks, A, n)
    A_mul_B!(view(b.A_block, :, 1:n), A, view(b.block, :, 1:n))
    return
end
function B_mul_X!(b::Blocks{true}, B)
    A_mul_B!(b.B_block, B, b.block)
    return
end
function B_mul_X!(b::Blocks{true}, B, n)
    A_mul_B!(view(b.B_block, :, 1:n), B, view(b.block, :, 1:n))
    return
end
function B_mul_X!(b::Blocks{false}, B, n = 0)
    return
end

struct Constraint{T, TA<:AbstractArray{T}, TC}
    Y::TA
    BY::TA
    gram_chol::TC
    gramYBV::TA # to be used in view
    tmp::TA # to be used in view
end
function Constraint(::Void, B, X)
    return Constraint{Void, Matrix{Void}, Void}(Matrix{Void}(0,0), Matrix{Void}(0,0), nothing, Matrix{Void}(0,0), Matrix{Void}(0,0))
end
function Constraint(Y, B, X)
    T = eltype(X)
    if B isa Void
        B = Y
    else
        BY = similar(B)
        A_mul_B!(BY, B, Y)
    end
    gramYBY = At_mul_B(Y, BY)
    gramYBY_chol = cholfact!(Hermitian(gramYBY))
    gramYBV = zeros(T, size(Y, 2), size(X, 2))
    tmp = similar(gramYBV)

    return Constraint(Y, BY, gramYBY_chol, gramYBV, tmp)
end

function (constr!::Constraint{Void})(X)
    nothing
end

function (constr!::Constraint)(X)
    sizeX = size(X, 2)
    sizeY = size(constr!.Y, 2)
    gramYBV_view = view(constr!.gramYBV, 1:sizeY, 1:sizeX)
    At_mul_B!(gramYBV_view, constr!.BY, X)
    tmp_view = view(constr!.tmp, 1:sizeY, 1:sizeX)
    A_ldiv_B!(tmp_view, gram_chol, gramYBV_view)
    A_mul_B!(X, constr!.Y, tmp_view)

    nothing
end

struct RPreconditioner{TM, T, TA<:AbstractArray{T}}
    M::TM
    buffer::TA
    RPreconditioner{TM, T, TA}(M, X) where {TM, T, TA<:AbstractArray{T}} = new(M, similar(X))
end
RPreconditioner(M, X) = RPreconditioner{typeof(M), eltype(X), typeof(X)}(M, X)

function (precond!::RPreconditioner{Void})(X)
    nothing
end
function (precond!::RPreconditioner)(X)
    bs = size(X, 2)
    A_ldiv_B!(view(precond!.buffer, :, 1:bs), precond!.M, X)
    # Just returning buffer would be cheaper but struct at call site must be mutable
    X .= view(precond!.buffer, :, 1:bs)
    nothing
end

struct BlockGram{Generalized, TA}
    XAX::TA
    XAR::TA
    XAP::TA
    RAR::TA
    RAP::TA
    PAP::TA
end
function BlockGram(XBlocks::Blocks{Generalized, T}) where {Generalized, T}
    sizeX = size(XBlocks.block, 2)
    XAX = zeros(T, sizeX, sizeX)
    XAP = zeros(T, sizeX, sizeX)
    XAR = zeros(T, sizeX, sizeX)
    RAR = zeros(T, sizeX, sizeX)
    RAP = zeros(T, sizeX, sizeX)
    PAP = zeros(T, sizeX, sizeX)
    return BlockGram{Generalized, Matrix{T}}(XAX, XAR, XAP, RAR, RAP, PAP)
end
XAX!(BlockGram, XBlocks) = At_mul_B!(BlockGram.XAX, XBlocks.block, XBlocks.A_block)
XAP!(BlockGram, XBlocks, PBlocks, n) = At_mul_B!(view(BlockGram.XAP, :, 1:n), XBlocks.block, view(PBlocks.A_block, :, 1:n))
XAR!(BlockGram, XBlocks, RBlocks, n) = At_mul_B!(view(BlockGram.XAR, :, 1:n), XBlocks.block, view(RBlocks.A_block, :, 1:n))
RAR!(BlockGram, RBlocks, n) = At_mul_B!(view(BlockGram.RAR, 1:n, 1:n), view(RBlocks.block, :, 1:n), view(RBlocks.A_block, :, 1:n))
RAP!(BlockGram, RBlocks, PBlocks, n) = At_mul_B!(view(BlockGram.RAP, 1:n, 1:n), view(RBlocks.A_block, :, 1:n), view(PBlocks.block, :, 1:n))
PAP!(BlockGram, PBlocks, n) = At_mul_B!(view(BlockGram.PAP, 1:n, 1:n), view(PBlocks.block, :, 1:n), view(PBlocks.A_block, :, 1:n))
XBP!(BlockGram, XBlocks, PBlocks, n) = At_mul_B!(view(BlockGram.XAP, :, 1:n), XBlocks.block, view(PBlocks.B_block, :, 1:n))
XBR!(BlockGram, XBlocks, RBlocks, n) = At_mul_B!(view(BlockGram.XAR, :, 1:n), XBlocks.block, view(RBlocks.B_block, :, 1:n))
RBP!(BlockGram, RBlocks, PBlocks, n) = At_mul_B!(view(BlockGram.RAP, 1:n, 1:n), view(RBlocks.B_block, :, 1:n), view(PBlocks.block, :, 1:n))
XBX!(BlockGram, XBlocks) = At_mul_B!(BlockGram.XAX, XBlocks.block, XBlocks.B_block)
RBR!(BlockGram, RBlocks, n) = At_mul_B!(view(BlockGram.RAR, 1:n, 1:n), view(RBlocks.block, :, 1:n), view(RBlocks.B_block, :, 1:n))
PBP!(BlockGram, PBlocks, n) = At_mul_B!(view(BlockGram.PAP, 1:n, 1:n), view(PBlocks.block, :, 1:n), view(PBlocks.B_block, :, 1:n))

function I!(G, xr)
    for j in xr, i in xr
        G[i, j] = ifelse(i==j, 1, 0)
    end
    return
end

function (g::BlockGram)(gram, lambda, n1::Int, n2::Int, n3::Int)
    xr = 1:n1
    rr = n1+1:n1+n2
    pr = n1+n2+1:n1+n2+n3
    if n1 > 0
        #gram[xr, xr] .= view(g.XAX, 1:n1, 1:n1)
        gram[xr, xr] .= Diagonal(view(lambda, 1:n1))
    end
    if n2 > 0
        gram[rr, rr] .= view(g.RAR, 1:n2, 1:n2)
        gram[xr, rr] .= view(g.XAR, 1:n1, 1:n2)
        transpose!(view(gram, rr, xr), view(g.XAR, 1:n1, 1:n2))
    end
    if n3 > 0
        gram[pr, pr] .= view(g.PAP, 1:n3, 1:n3)
        gram[rr, pr] .= view(g.RAP, 1:n2, 1:n3)
        gram[xr, pr] .= view(g.XAP, 1:n1, 1:n3)
        transpose!(view(gram, pr, rr), view(g.RAP, 1:n2, 1:n3))
        transpose!(view(gram, pr, xr), view(g.XAP, 1:n1, 1:n3))
    end
    return
end
function (g::BlockGram)(gram, n1::Int, n2::Int, n3::Int, normalized::Bool=true)
    xr = 1:n1
    rr = n1+1:n1+n2
    pr = n1+n2+1:n1+n2+n3
    if n1 > 0
        if normalized
            I!(gram, xr)
        else
            gram[xr, xr] .= view(g.XAX, 1:n1, 1:n1)
        end
    end
    if n2 > 0
        if normalized
            I!(gram, rr)
        else
            gram[rr, rr] .= view(g.RAR, 1:n2, 1:n2)
        end
        gram[xr, rr] .= view(g.XAR, 1:n1, 1:n2)
        transpose!(view(gram, rr, xr), view(g.XAR, 1:n1, 1:n2))
    end
    if n3 > 0
        if normalized
            I!(gram, pr)
        else
            gram[pr, pr] .= view(g.PAP, 1:n3, 1:n3)
        end
        gram[rr, pr] .= view(g.RAP, 1:n2, 1:n3)
        gram[xr, pr] .= view(g.XAP, 1:n1, 1:n3)
        transpose!(view(gram, pr, rr), view(g.RAP, 1:n2, 1:n3))
        transpose!(view(gram, pr, xr), view(g.XAP, 1:n1, 1:n3))
    end
    return
end

abstract type AbstractOrtho end
struct CholQR{TA} <: AbstractOrtho
    gramVBV::TA # to be used in view
end

function A_rdiv_B!(A, B::UpperTriangular)
    s = size(A, 2)
    A[:,1] .= view(A, :, 1) ./ B[1,1]
    for i in 2:s
        for j in 1:i-1
            A[:,i] .= view(A, :, i) - view(A, :, j) .* B[j,i]
        end
        A[:,i] .= view(A, :, i) ./ B[i,i]
    end
    return
end

function (ortho!::CholQR)(XBlocks::Blocks{Generalized}, sizeX = -1; update_AX=false, update_BX=false) where Generalized
    useview = sizeX != -1
    if sizeX == -1
        sizeX = size(XBlocks.block, 2)
    end
    X = XBlocks.block
    BX = XBlocks.B_block # Assumes it is premultiplied
    AX = XBlocks.A_block
    gram_view = view(ortho!.gramVBV, 1:sizeX, 1:sizeX)
    if useview
        At_mul_B!(gram_view, view(X, :, 1:sizeX), view(BX, :, 1:sizeX))
    else
        At_mul_B!(gram_view, X, BX)
    end
    cholf = cholfact!(Hermitian(gram_view))
    R = cholf.factors
    if useview
        A_rdiv_B!(view(X, :, 1:sizeX), UpperTriangular(R))
        update_AX && A_rdiv_B!(view(AX, :, 1:sizeX), UpperTriangular(R))
        Generalized && update_BX && A_rdiv_B!(view(BX, :, 1:sizeX), UpperTriangular(R))
    else
        A_rdiv_B!(X, UpperTriangular(R))
        update_AX && A_rdiv_B!(AX, UpperTriangular(R))
        Generalized && update_BX && A_rdiv_B!(BX, UpperTriangular(R))
    end

    return
end

struct RayleighRitz{Generalized, T, TA, TB, TL<:AbstractVector{T}, TVec<:AbstractVector{Int}, TV<:AbstractArray{T}, TBlocks<:Blocks{Generalized, T}, TO<:AbstractOrtho, TP, TC, TG, TH, TM}
    A::TA
    B::TB
    λ::TL
    λperm::TVec
    V::TV
    residuals::TL
    largest::Bool
    XBlocks::TBlocks
    tempXBlocks::TBlocks
    PBlocks::TBlocks
    activePBlocks::TBlocks # to be used in view
    RBlocks::TBlocks
    activeRBlocks::TBlocks # to be used in view
    iteration::Base.RefValue{Int}
    currentBlockSize::Base.RefValue{Int}
    ortho!::TO
    precond!::TP
    constr!::TC
    gramABlock::TG
    gramBBlock::TG
    gramA::TV
    gramB::TV
    residualNormsHistory::TH
    activeMask::TM
end
function RayleighRitz(A, B, M, Y, X, largest)
    T = eltype(X)
    constr! = Constraint(Y, B, X)
    precond! = RPreconditioner(M, X)

    nev = size(X, 2)
    if B isa Void
        XBlocks = Blocks(copy(X), similar(X))
        tempXBlocks = Blocks(copy(X), similar(X))
        RBlocks = Blocks(similar(X), similar(X))
        activeRBlocks = Blocks(similar(X), similar(X))
        PBlocks = Blocks(similar(X), similar(X))
        activePBlocks = Blocks(similar(X), similar(X))
    else
        XBlocks = Blocks(copy(X), similar(X), similar(X))
        tempXBlocks = Blocks(copy(X), similar(X), similar(X))
        RBlocks = Blocks(similar(X), similar(X), similar(X))
        activeRBlocks = Blocks(similar(X), similar(X), similar(X))
        PBlocks = Blocks(similar(X), similar(X), similar(X))
        activePBlocks = Blocks(similar(X), similar(X), similar(X))
    end
    λ = zeros(T, nev*3)
    λperm = zeros(Int, nev*3)
    V = zeros(T, nev*3, nev*3)
    residuals = zeros(T, nev)
    iteration = Ref(1)
    currentBlockSize = Ref(nev)
    generalized = !(B isa Void)
    ortho! = CholQR(zeros(T, nev, nev))

    gramABlock = BlockGram(XBlocks)
    gramBBlock = BlockGram(XBlocks)

    residualNormsHistory = Vector{Float64}[]
    activeMask = ones(Bool, nev)

    gramA = zeros(T, 3*nev, 3*nev)
    gramB = zeros(T, 3*nev, 3*nev)

    return RayleighRitz{generalized, T, typeof(A), typeof(B), typeof(λ), typeof(λperm), typeof(V), typeof(XBlocks), typeof(ortho!), typeof(precond!), typeof(constr!), typeof(gramABlock), typeof(residualNormsHistory), typeof(activeMask)}(A, B, λ, λperm, V, residuals, largest, XBlocks, tempXBlocks, PBlocks, activePBlocks, RBlocks, activeRBlocks, iteration, currentBlockSize, ortho!, precond!, constr!, gramABlock, gramBBlock, gramA, gramB, residualNormsHistory, activeMask)
end

function ortho_AB_mul_X!(blocks::Blocks, ortho!, A, B, bs=-1)
    # Finds BX
    bs == -1 ? B_mul_X!(blocks, B) : B_mul_X!(blocks, B, bs)
    # Orthonormalizes X and updates BX
    bs == -1 ? ortho!(blocks, update_BX=true) : ortho!(blocks, bs, update_BX=true)
    # Updates AX
    bs == -1 ? A_mul_X!(blocks, A) : A_mul_X!(blocks, A, bs)
    return
end
function residuals!(rr)
    sizeX = size(rr.XBlocks.block, 2)
    A_mul_B!(rr.RBlocks.block, rr.XBlocks.B_block, Diagonal(view(rr.λ, 1:sizeX)))
    rr.RBlocks.block .= rr.XBlocks.A_block .- rr.RBlocks.block
    # Finds residual norms
    for j in 1:size(rr.RBlocks.block, 2)
        rr.residuals[j] = 0
        for i in 1:size(rr.RBlocks.block, 1)
            x = rr.RBlocks.block[i,j]
            rr.residuals[j] += conj(x)*x
        end
        rr.residuals[j] = sqrt(rr.residuals[j])
    end
    return
end

function update_mask!(rr, residualTolerance)
    sizeX = size(rr.XBlocks.block, 2)
    # Update active vectors mask
    rr.activeMask .*= view(rr.residuals, 1:sizeX) .> residualTolerance
    rr.currentBlockSize[] = sum(rr.activeMask)
    return
end

function update_active!(mask, bs::Int, blockPairs...)
    for (activeblock, block) in blockPairs
        activeblock[:, 1:bs] .= view(block, :, mask)
    end
    return
end

function precond_constr!(block, bs, precond!, constr!)
    precond!(view(block, :, 1:bs))
    # Constrain the active residual vectors to be B-orthogonal to Y
    constr!(view(block, :, 1:bs))
    return
end
function block_grams_1x1!(rr)
    # Finds gram matrix X'AX
    XAX!(rr.gramABlock, rr.XBlocks)
    return
end
function block_grams_2x2!(rr, bs)
    sizeX = size(rr.XBlocks.block, 2)
    #XAX!(rr.gramABlock, rr.XBlocks)
    XAR!(rr.gramABlock, rr.XBlocks, rr.activeRBlocks, bs)
    RAR!(rr.gramABlock, rr.activeRBlocks, bs)
    XBR!(rr.gramBBlock, rr.XBlocks, rr.activeRBlocks, bs)
    rr.gramABlock(rr.gramA, view(rr.λ, 1:sizeX), sizeX, bs, 0)
    rr.gramBBlock(rr.gramB, sizeX, bs, 0, true)

    return
end
function block_grams_3x3!(rr, bs)
    # Find R'AR, P'AP, X'AR, X'AP and R'AP
    sizeX = size(rr.XBlocks.block, 2)
    #XAX!(rr.gramABlock, rr.XBlocks)
    XAR!(rr.gramABlock, rr.XBlocks, rr.activeRBlocks, bs)
    XAP!(rr.gramABlock, rr.XBlocks, rr.activePBlocks, bs)
    RAR!(rr.gramABlock, rr.activeRBlocks, bs)
    RAP!(rr.gramABlock, rr.activeRBlocks, rr.activePBlocks, bs)
    PAP!(rr.gramABlock, rr.activePBlocks, bs)
    # Find X'BR, X'BP and P'BR
    XBR!(rr.gramBBlock, rr.XBlocks, rr.activeRBlocks, bs)
    XBP!(rr.gramBBlock, rr.XBlocks, rr.activePBlocks, bs)
    RBP!(rr.gramBBlock, rr.activeRBlocks, rr.activePBlocks, bs)
    # Update the gram matrix [X R P]' A [X R P]
    rr.gramABlock(rr.gramA, view(rr.λ, 1:sizeX), sizeX, bs, bs)
    # Update the gram matrix [X R P]' B [X R P]
    rr.gramBBlock(rr.gramB, sizeX, bs, bs, true)

    return
end

function sub_problem!(rr, sizeX, bs1, bs2)
    subdim = sizeX+bs1+bs2
    if bs1 == 0
        gramAview = view(rr.gramABlock.XAX, 1:subdim, 1:subdim)
        # Source of type instability
        eigf = eigfact!(Hermitian(gramAview))
    else
        gramAview = view(rr.gramA, 1:subdim, 1:subdim)
        gramBview = view(rr.gramB, 1:subdim, 1:subdim)
        # Source of type instability
        eigf = eigfact!(Hermitian(gramAview), Hermitian(gramBview))
    end
    # Selects extremal eigenvalues and corresponding vectors
    selectperm!(view(rr.λperm, 1:subdim), eigf.values, 1:subdim, rev=rr.largest)
    rr.λ[1:sizeX] .= view(eigf.values, view(rr.λperm, 1:sizeX))
    rr.V[1:subdim, 1:sizeX] .= view(eigf.vectors, :, view(rr.λperm, 1:sizeX))

    return
end

function update_X_P!(rr::RayleighRitz{Generalized}, bs1, bs2) where Generalized
    sizeX = size(rr.XBlocks.block, 2)
    x_eigview = view(rr.V, 1:sizeX, 1:sizeX)
    r_eigview = view(rr.V, sizeX+1:sizeX+bs1, 1:sizeX)
    p_eigview = view(rr.V, sizeX+bs1+1:sizeX+bs1+bs2, 1:sizeX)
    r_blockview = view(rr.activeRBlocks.block, :, 1:bs1)
    ra_blockview = view(rr.activeRBlocks.A_block, :, 1:bs1)
    p_blockview = view(rr.activePBlocks.block, :, 1:bs2)
    pa_blockview = view(rr.activePBlocks.A_block, :, 1:bs2)
    if Generalized
        rb_blockview = view(rr.activeRBlocks.B_block, :, 1:bs1)
        pb_blockview = view(rr.activePBlocks.B_block, :, 1:bs2)
    end
    if bs1 > 0
        A_mul_B!(rr.PBlocks.block, r_blockview, r_eigview)
        A_mul_B!(rr.PBlocks.A_block, ra_blockview, r_eigview)
        if Generalized
            A_mul_B!(rr.PBlocks.B_block, rb_blockview, r_eigview)
        end
    end
    if bs2 > 0
        A_mul_B!(rr.tempXBlocks.block, p_blockview, p_eigview)
        A_mul_B!(rr.tempXBlocks.A_block, pa_blockview, p_eigview)
        if Generalized
            A_mul_B!(rr.tempXBlocks.B_block, pb_blockview, p_eigview)
        end
        rr.PBlocks.block .= rr.PBlocks.block .+ rr.tempXBlocks.block
        rr.PBlocks.A_block .= rr.PBlocks.A_block .+ rr.tempXBlocks.A_block
        if Generalized
            rr.PBlocks.B_block .= rr.PBlocks.B_block .+ rr.tempXBlocks.B_block
        end
    end
    block = rr.XBlocks.block
    tempblock = rr.tempXBlocks.block
    A_mul_B!(tempblock, block, x_eigview)
    block = rr.XBlocks.A_block
    tempblock = rr.tempXBlocks.A_block
    A_mul_B!(tempblock, block, x_eigview)
    if Generalized
        block = rr.XBlocks.B_block
        tempblock = rr.tempXBlocks.B_block
        A_mul_B!(tempblock, block, x_eigview)
    end
    if bs1 > 0
        rr.XBlocks.block .= rr.tempXBlocks.block .+ rr.PBlocks.block
        rr.XBlocks.A_block .= rr.tempXBlocks.A_block .+ rr.PBlocks.A_block
        if Generalized
            rr.XBlocks.B_block .= rr.tempXBlocks.B_block .+ rr.PBlocks.B_block
        end
    else
        rr.XBlocks.block .= rr.tempXBlocks.block
        rr.XBlocks.A_block .= rr.tempXBlocks.A_block
        if Generalized
            rr.XBlocks.B_block .= rr.tempXBlocks.B_block
        end
    end
    return
end

function (rr::RayleighRitz{Generalized})(residualTolerance) where Generalized
    sizeX = size(rr.XBlocks.block, 2)
    iteration = rr.iteration[]
    if iteration == 1
        ortho_AB_mul_X!(rr.XBlocks, rr.ortho!, rr.A, rr.B)
        # Finds gram matrix X'AX
        block_grams_1x1!(rr)
        sub_problem!(rr, sizeX, 0, 0)
        # Updates Ritz vectors X and updates AX and BX accordingly
        update_X_P!(rr, 0, 0)
    elseif iteration == 2
        residuals!(rr)
        # Store history of norms
        push!(rr.residualNormsHistory, rr.residuals[1:sizeX])
        update_mask!(rr, residualTolerance)
        rr.currentBlockSize[] == 0 && return
        bs = rr.currentBlockSize[]
        # Update active R blocks
        update_active!(rr.activeMask, bs, (rr.activeRBlocks.block, rr.RBlocks.block))
        # Precondition and constrain the active residual vectors
        precond_constr!(rr.activeRBlocks.block, bs, rr.precond!, rr.constr!)
        # Orthonormalizes R[:,1:bs] and finds AR[:,1:bs] and BR[:,1:bs]
        ortho_AB_mul_X!(rr.activeRBlocks, rr.ortho!, rr.A, rr.B, bs)
        # Find [X R] A [X R] and [X R]' B [X R]
        block_grams_2x2!(rr, bs)
        # Solve the Rayleigh-Ritz sub-problem
        sub_problem!(rr, sizeX, bs, 0)
        update_X_P!(rr, bs, 0)
    else
        residuals!(rr)
        # Store history of norms
        push!(rr.residualNormsHistory, rr.residuals[1:sizeX])
        update_mask!(rr, residualTolerance)
        rr.currentBlockSize[] == 0 && return
        # Update active blocks
        bs = rr.currentBlockSize[]
        # Update active R and P blocks
        update_active!(rr.activeMask, bs, (rr.activeRBlocks.block, rr.RBlocks.block),
                                          (rr.activePBlocks.block, rr.PBlocks.block),
                                          (rr.activePBlocks.A_block, rr.PBlocks.A_block),
                                          (rr.activePBlocks.B_block, rr.PBlocks.B_block))
        # Precondition and constrain the active residual vectors
        precond_constr!(rr.activeRBlocks.block, bs, rr.precond!, rr.constr!)
        # Orthonormalizes R[:,1:bs] and finds AR[:,1:bs] and BR[:,1:bs]
        ortho_AB_mul_X!(rr.activeRBlocks, rr.ortho!, rr.A, rr.B, bs)
        # Orthonormalizes P and updates AP
        rr.ortho!(rr.activePBlocks, bs, update_AX=true, update_BX=true)
        # Update the gram matrix [X R P]' A [X R P] and [X R P]' B [X R P]
        block_grams_3x3!(rr, bs)
        # Solve the Rayleigh-Ritz sub-problem
        sub_problem!(rr, sizeX, bs, bs)
        # Updates Ritz vectors X and updates AX and BX accordingly
        # And updates P, AP and BP
        update_X_P!(rr, bs, bs)
    end

    return
end

function dense_solver(A, B, X, largest)
    warn("The problem size is small compared to the block size. Using dense eigensolver instead of LOBPCG.")
    # Define the closed range of indices of eigenvalues to return.
    n, sizeX = size(X)
    eigvals = largest ? (n - sizeX + 1, n) : (1, sizeX)
    A_dense = A*eye(n)
    if B isa Void
        return eig(Hermitian(A_dense))
    else
        B_dense = B*eye(n)
        return eig(Hermitian(A_dense), Hermitian(B_dense))
    end
end

function lobpcg(A, X::AbstractMatrix, largest::Bool=true, ::Type{Val{residualhistory}} = Val{false}; preconditioner=nothing, constraint=nothing, maxiter=100, tol::Number=nothing) where {residualhistory}
    lobpcg(A, nothing, X, largest, Val{residualhistory}; tol=tol, maxiter=maxiter, preconditioner=preconditioner, constraint=constraint)
end

"""
    The Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG)

    function lobpcg(A, B, X, largest=true, ::Type{Val{residualhistory}}=Val{false};
        preconditioner=nothing, constraint=nothing,
        tol=nothing, maxiter=100) where {residualhistory}

    Finds the k extremal eigenvalues and their corresponding eigenvectors satisfying `AX = λBX` returning `λ` and `X`, using `X` as an initial guess, where `k = size(X,2)`. The input `X` is overwritten with the solution.

    `A` and `B` may be generic types but `Base.A_mul_B!(C, AorB, X)` and `AorB * X` must be defined for vectors and strided matrices `X` and `C`.

    Example:

    T = Float64;
    n = 1000; # Problem dimension
    k = 3; # Number of eigenvalues

    A = 10*sprand(T, n, n, 10/n);
    A = A' + A + 50I;
    B = 10*sprand(T, n, n, 10/n);
    B = B' + B + 500I;
    X0 = zeros(T, n, 3);
    tol = T(1e-6);
    maxiter = 200;

    # Smallest generalized eigenvalues
    λ, X = lobpcg(A, B, X0, false, tol=tol, maxiter=maxiter);
    # Smallest simple eigenvalues
    λ, X = lobpcg(A, X0, false, tol=tol, maxiter=maxiter);
    # Largest generalized eigenvalues
    λ, X = lobpcg(A, B, X0, true, tol=tol, maxiter=maxiter);
    # Largest simple eigenvalues
    λ, X = lobpcg(A, X0, true, tol=tol, maxiter=maxiter);

"""
function lobpcg(A, B, X0, largest=true, ::Type{Val{residualhistory}}=Val{false};
                preconditioner=nothing, constraint=nothing,
                tol=nothing, maxiter=100) where {residualhistory}

    X = copy(X0)
    T = eltype(X)
    M = preconditioner
    Y = constraint
    n, sizeX = size(X)
    if Y isa Void
        n < 5 * sizeX && return dense_solver(A, B, X, largest)
    else
        sizeY = size(Y, 2)
        (n - sizeY) < 5 * sizeX && throw("The dense eigensolver does not support constraints.")
    end
    sizeX > n && throw("X column dimension exceeds the row dimension")

    for j in 1:size(X,2)
        if all(x -> x==0, view(X, :, j))
            X[:,j] .= rand.()
        end
    end
    rr = RayleighRitz(A, B, M, Y, X, largest)
    residualTolerance = (tol isa Void) ? sqrt(1e-15)*n : tol
    maxiter = min(n, maxiter)
    for iteration in 1:maxiter
        rr.iteration[] = iteration
        rr(residualTolerance)
        rr.currentBlockSize[] == 0 && break
    end
    if residualhistory
        return rr.λ[1:sizeX], rr.XBlocks.block, rr.residualNormsHistory
    else
        return rr.λ[1:sizeX], rr.XBlocks.block
    end
end