text:
using OpenCL

struct OCLContext
    ocl_device::OpenCL.cl.Device
    ocl_ctx::OpenCL.cl.Context
    ocl_queue::OpenCL.cl.CmdQueue
end

function OCLContext()
    device, ctx, queue = OpenCL.cl.create_compute_context()
    return OCLContext(device, ctx, queue)
end

function OCLContext(device::OpenCL.cl.Device)
    ctx   = cl.Context(device);
    queue = cl.CmdQueue(ctx)
    return OCLContext(device,ctx,queue)
end

"""
FWSim_OpenCL contains the created OpenCL code for the FW-Sim

Structure of the code is a single kernel with 3 parameters:
__kernel void fwsim(__global const float *semus,
                  __global const float *fluxes,
                  __global float *mds_out)

semus[ i*length(semu_mass_ratios) ] is the first element of the semus data
for sample i.
fluxes[ i*length(fluxes) ] is the first element of the flux rate data
for sample i.
"""
mutable struct FWSim_OpenCL <: FWSim
    code::String
    fluxes::Vector{String} # defines the order of fluxes which have to be passed to the function, is equal to the flux order in the supplied network
    semu_mass_ratios::Vector{ASTEMUMassRatio} # defines the order of defined semus
    target_mass_ratios::Vector{ASTEMUMassRatio} # defines the order of result semus , is equal to the supplied vector

    required_matrix_size_x::Int
    required_matrix_size_y::Int

    ocl_context::Nullable{OCLContext}
    k::Nullable{OpenCL.cl.Kernel}
    p::Nullable{OpenCL.cl.Program}
end
function FWSim_OpenCL( code::String, fluxes::Vector{String}, semu_mass_ratios::Vector{ASTEMUMassRatio} , target_mass_ratios::Vector{ASTEMUMassRatio} , required_matrix_size_x::Int , required_matrix_size_y::Int  )
    return FWSim_OpenCL(code,fluxes,semu_mass_ratios,target_mass_ratios, required_matrix_size_x,required_matrix_size_y,Nullable{OCLContext}() , Nullable{OpenCL.cl.Kernel}() , Nullable{OpenCL.cl.Program}() )
end

function ast_to_opencl_01( net::Network , ast::ASTEMUSim , target_smrs::Vector{ASTEMUMassRatio} )

    flag_include_debug::Bool = false

    # 1. create function "header"
    code_header::IOBuffer = IOBuffer()
    print(code_header ,"     __kernel void fwsim(__global const float *semus,
                          __global const float *fluxes,
                          __global float *mds_out)
        {
          const int gid = get_group_id(0) + get_group_id(1);
          const int idx = get_local_id(0);
          const int idy = get_local_id(1);

          //printf(\"\\n group_id=%d idx=%d idy=%d glob_x=%d glob_y=%d \\n\",gid,idx,idy,get_global_id(0),get_global_id(1));

          // Here we store the sizes of the given matrices in a step
          int mAm = 0;
          int mAn = 0;
          int mBm = 0;
          int mBn = 0;
          int mXm = 0;
          int mXn = 0;
          int mYm = 0;
          int mYn = 0;

        \n\n")

    # 1.b Init the local matrices needed for the single steps..
    req_size_A = maximum( map( x -> maximum( collect(map( xi -> prod(size(xi.A.m)) , x))) , values(ast.sim_steps) ) )
    req_size_B = maximum( map( x -> maximum( collect(map( xi -> prod(size(xi.B.m)) , x))) , values(ast.sim_steps) ) )
    req_size_X = maximum( map( x -> maximum( collect(map( xi -> prod(size(xi.X.m)) , x))) , values(ast.sim_steps) ) )
    req_size_Y = maximum( map( x -> maximum( collect(map( xi -> prod(size(xi.Y.m)) , x))) , values(ast.sim_steps) ) )
    req_size_BY = maximum( map( x -> maximum( collect(map( xi -> size(xi.B.m,1) * size(xi.Y.m,2) , x))) , values(ast.sim_steps) ) )
    # 1.c Create local matrices:
    print(code_header , "__local float mA[$(req_size_A)]; \n")
    print(code_header , "__local float mA_L[$(req_size_A)]; \n") # lu decomp, l part
    print(code_header , "__local float mB[$(req_size_B)]; \n")
    print(code_header , "__local float mX[$(req_size_X)]; \n")
    print(code_header , "__local float mY[$(req_size_Y)]; \n")
    print(code_header , "__local float mBY[$(req_size_BY)]; \n")

    # here we add all that we find in the X matrix, in any step
    semu_mass_ratios_computed::Set{ASTEMUMassRatio} = Set{ASTEMUMassRatio}()

    # 3. step in, start with processing the sim steps..
    code_computation::IOBuffer = IOBuffer()

    req_constants::Dict{ASTNumNode,String} = Dict{ASTNumNode,String}()
    for siz in sort( collect( keys(ast.sim_steps) ) )
        ss_cnt = 0
        for ss in ast.sim_steps[siz]
            ss_cnt+=1
            print(code_computation, "//----Init A---------------------------------------------------------------------------------\n")
            print(code_computation, "if(idx==0 && idy==0){\n")
            print( code_computation , to_opencl_code_ast_matrix( "mA" , ss.A , req_constants ) )
            print(code_computation, "}\n")
            print(code_computation, "//----Init B---------------------------------------------------------------------------\n")
            print(code_computation, "if(idx==1 && idy==0){\n")
            print( code_computation , to_opencl_code_ast_matrix( "mB" , ss.B , req_constants ) )
            print(code_computation, "}\n")
            print(code_computation, "//----Init Y---------------------------------------------------------------------------\n")
            print(code_computation, "if(idx==2 && idy==0){\n")
            print( code_computation , to_opencl_code_ast_matrix( "mY" , ss.Y , req_constants ) )
            print(code_computation, "}\n")

            # !! We have to set the matrix size indicators on work item level!!
            print(code_computation,"mAm = $(size(ss.A.m,1)); \n")
            print(code_computation,"mAn = $(size(ss.A.m,2)); \n")
            print(code_computation,"mBm = $(size(ss.B.m,1)); \n")
            print(code_computation,"mBn = $(size(ss.B.m,2)); \n")
            print(code_computation,"mYm = $(size(ss.Y.m,1)); \n")
            print(code_computation,"mYn = $(size(ss.Y.m,2)); \n")

            # might not be needed..
            print(code_computation, "barrier(CLK_LOCAL_MEM_FENCE);\n")

            if(flag_include_debug) # some debug code.. :)
                print( code_computation ,  "\nif(idx==0&&idy==0){printf(\"\\nBefore Computation:\");}\n")
                print( code_computation ,  ocl_print_matrix( "mA" , size(ss.A.m,1) , size(ss.A.m,2)) )
                print( code_computation ,  ocl_print_matrix( "mB" , size(ss.B.m,1) , size(ss.B.m,2)) )
                print( code_computation ,  ocl_print_matrix( "mY" , size(ss.Y.m,1) , size(ss.Y.m,2)) )
            end

            print(code_computation, "\n//------COMPUTE!-----------------------------------------------------------------------\n")

            #print(code_computation, to_opencl_code_ludecomp_2("mA",size(ss.A.m,1)) )
            print(code_computation, to_opencl_code_ludecomp_2_1d("mA",size(ss.A.m,1)) )
            print(code_computation, "//----------------------------------\n")
            #print(code_computation, to_opencl_code_matrixmul("mBY","mB","mY", ss.B , ss.Y ) )
            print(code_computation, to_opencl_code_matrixmul_1d("mBY","mB","mY", ss.B , ss.Y ) )
            print(code_computation,"barrier(CLK_LOCAL_MEM_FENCE);\n")

            if(flag_include_debug) # some debug code.. :)
                print( code_computation ,  "\nif(idx==0&&idy==0){printf(\"\\nResult LU Decomp:\");}\n")
                print( code_computation ,  ocl_print_matrix( "mA" , size(ss.A.m,1) , size(ss.A.m,2)) )
                print( code_computation ,  ocl_print_matrix( "mA_L" , size(ss.A.m,1) , size(ss.A.m,2)) )
            end

            if(flag_include_debug) # some debug code.. :)
                print( code_computation ,  "\nif(idx==0&&idy==0){printf(\"\\nafter Multiply B and Y:\");}\n")
                print( code_computation ,  ocl_print_matrix( "mBY" , size(ss.B.m,1) , size(ss.Y.m,2)) )
            end

            print(code_computation, to_opencl_code_solve_L_system( "mA_L" , "mBY" , size(ss.A.m,1) , size(ss.A.m,2) , size(ss.Y.m,2) ))
            print(code_computation,"barrier(CLK_LOCAL_MEM_FENCE);\n")

            if(flag_include_debug) # some debug code.. :)
                print( code_computation ,  "\nif(idx==0&&idy==0){printf(\"\\nafter L Solve:\");}\n")
                print( code_computation ,  ocl_print_matrix( "mBY" , size(ss.B.m,1) , size(ss.Y.m,2)) )
            end

            print(code_computation, to_opencl_code_solve_U_system( "mA" , "mBY" , size(ss.A.m,1) , size(ss.A.m,2) , size(ss.Y.m,2) ))
            print(code_computation,"barrier(CLK_LOCAL_MEM_FENCE);\n")

            if(flag_include_debug) # some debug code.. :)
                print( code_computation ,  "\nif(idx==0&&idy==0){printf(\"\\nafter U Solve:\");}\n")
                print( code_computation ,  ocl_print_matrix( "mBY" , size(ss.B.m,1) , size(ss.Y.m,2)) )
            end

            #    print( code_computation , "( sA_L_$(siz)_$(ss_cnt) , sA_U_$(siz)_$(ss_cnt) ) = lu(sA_$(siz)_$(ss_cnt),Val{false})\n" )
            #    print( code_computation , "sX_pre_$(siz)_$(ss_cnt) = sA_L_$(siz)_$(ss_cnt) \\ ( sB_$(siz)_$(ss_cnt) * sY_$(siz)_$(ss_cnt) )\n" )
            #    print( code_computation , "sX_$(siz)_$(ss_cnt) = sA_U_$(siz)_$(ss_cnt) \\ sX_pre_$(siz)_$(ss_cnt) \n\n" )

            print(code_computation, "//-------------------------------------------------------------------------------------\n")
            #print(code_computation, "if(idx==0&&idy==0){ \n") # !! ALL WORK ITEMS HAVE TO DO THIS!!
            # now read out the newly computed values from X:
            for zib = 1:size(ss.X.m,2)
                for zia = 1:size(ss.X.m,1)
                    smri = ss.X.m[zia,zib]
                    @assert typeof(smri) == ASTEMUMassRatio
                    var_name_si = create_var_name(smri)
                    req_constants[smri] = var_name_si
                    push!( semu_mass_ratios_computed , smri )
                    #print(code_computation, "$(var_name_si) = sX_$(siz)_$(ss_cnt)[$(zia),$(zib)]\n")
                    print(code_computation, "$(var_name_si) = mBY[ $( (zia-1) + (zib-1) * size(ss.X.m,1)) ];\n") # !! -1
                end
            end
            #print(code_computation, "} \n")
            print(code_computation, "barrier(CLK_LOCAL_MEM_FENCE); \n")
            print(code_computation, "//-------------------------------------------------------------------------------------\n")
        end
    end

    # 4. now we add the required variable initializations
    # 4.1 collect the required mass distribution ratios:
    # code_varinit also sets the offsets for flux/semu/mratio data
    code_varinit::IOBuffer = IOBuffer()

    rc_flux = filter( x1 -> typeof(x1)==ASTFluxRate , collect(keys(req_constants)) )
    rc_emr  = filter( x1 -> typeof(x1)==ASTEMUMassRatio , collect(keys(req_constants)) )
    rc_flux = convert(Vector{ASTFluxRate},collect(rc_flux))
    rc_emr  = convert(Vector{ASTEMUMassRatio},collect(rc_emr))

    rc_emr_needed::Vector{ASTEMUMassRatio} = collect( setdiff( Set(rc_emr) , Set(semu_mass_ratios_computed) ) )
    rc_emr_needed                          = sort(rc_emr_needed)
    print("\n\nRequired EMRs: $(length(rc_emr_needed))  Required Fluxes: $(length(rc_flux))  Target MDs: $(length(target_smrs)) \n");
    print(code_varinit,"\n\n//Required EMRs: $(length(rc_emr_needed))  Required Fluxes: $(length(rc_flux))  Target MDs: $(length(target_smrs)) \n");
    print(code_varinit,"int offset_emrs   = 0; \n"); # CHANGE THIS IF EVERY WORK ITEM HAS DIFF. Substrate Config..
    print(code_varinit,"int offset_fluxes = gid * $(length(rc_flux)); \n");
    print(code_varinit,"int offset_target = gid * $(length(target_smrs)); \n");

    #print(code_varinit,"\n\nRequired EMRs: \n");
    print(code_varinit, "//-----Init EMRs-------------------------------------------------------------------------------------\n")
    emri_cnt = 0
    for emri in rc_emr_needed
        emri_cnt += 1
        println("$(emri)")
        print(code_varinit,"float $(req_constants[emri]) = semus[offset_emrs + $(emri_cnt-1)];\n") #!!!! -1
    end

    rc_flux = sort(rc_flux) # NOTE: this does not change anything, just makes the resulting code look nicer..
    #print("\n\nRequired Fluxes: \n");
    print(code_varinit, "//-----Init Fluxes-----------------------------------------------------------------------------------\n")
    for fluxii in rc_flux
        println("$(fluxii)")
        #fluxpos = find( xfi -> xfi==fluxii.flux , data.net.fluxes )
        fluxpos = find( xfi -> xfi==fluxii.flux , net.fluxes )
        @assert length(fluxpos)==1
        print(code_varinit,"float $(req_constants[fluxii]) = fluxes[ offset_fluxes + $(fluxpos[1]-1)];\n") # !!! -1
    end
    print(code_varinit, "//---------------------------------------------------------------------------------------------------\n")
    print(code_varinit, "//------Declare All Temporary EMR Float values--------------------------------------------------------------\n")
    all_emrs_without_initial_ones = setdiff( collect(keys(req_constants)) , union( rc_emr_needed , rc_flux ) )
    for emi in all_emrs_without_initial_ones
        print(code_varinit,"float $(req_constants[emi]);\n")
    end
    print(code_varinit, "//---------------------------------------------------------------------------------------------------\n")
    # 5. collect the results
    code_results::IOBuffer = IOBuffer()
    print(code_results, "//---RESULTS-----------------------------------------------------------------------------------------\n")
    print(code_results, "if(idx==0&&idy==0){ \n")
    cnt_tsr = 0
    #print(code_results , "target_smrs = Array{Float64,1}($(length(target_smrs)))\n")
    for t_semu in target_smrs
        cnt_tsr +=1
        print(code_results , "mds_out[ offset_target + $(cnt_tsr-1) ] = $(req_constants[t_semu]);\n") #!!!! -1
    end
    print(code_results, "}\n")
    print(code_results, "//---------------------------------------------------------------------------------------------------\n")
    #print(code_results,"return target_smrs\n")
    #print(code_results, "#---------------------------------------------------------------------------------------------------\n")

    #str_total = string( takebuf_string(code_header) ,  takebuf_string(code_header) ,  takebuf_string(code_computation) )
    code_total = string( String(take!(code_header)) ,  String(take!(code_varinit)) ,  String(take!(code_computation) ) , String(take!(code_results) ) , "}" )

    #return FWSim_Julia(code_total,data.net.fluxes,rc_emr_needed,[])
    #return code_total

    # find required matrix size..
    req_matrix_size = maximum( map( x -> maximum( collect(map( xi -> maximum( [ size(xi.A.m,1) ; size(xi.A.m,2) ;size(xi.B.m,1) ;size(xi.B.m,2) ;size(xi.Y.m,1) ;size(xi.Y.m,2) ] ) , x))) , values(ast.sim_steps) ) )

    return FWSim_OpenCL(code_total,net.fluxes,rc_emr_needed,target_smrs , req_matrix_size , req_matrix_size )
end

function ocl_mpos(m::ASTMatrix,x::Int,y::Int)
    return x*size(m.m,1)+y
end
function ocl_mpos(s::Int,x::Int,y::Int)
    return x*size(m.m,1)+y
end

function to_opencl_code_ast_matrix(matrix_name::String , m::ASTMatrix ,  req_constants::Dict{ASTNumNode,String} )
    sb::IOBuffer = IOBuffer()
    #print(sb,"$(matrix_name)::Array{Float64,2} = Array{Float64,2}($(size(m.m,1)),$(size(m.m,2))) \n")

    # this has to happen on work-item level!!
    #print(sb,"$(matrix_name)m = $(size(m.m,1)); \n")
    #print(sb,"$(matrix_name)n = $(size(m.m,2)); \n")

    # fill the matrix..
    for zj=0:size(m.m,2)-1
        for zi=0:size(m.m,1)-1
            # evaluate cell..
            if(isassigned(m.m,zi+1,zj+1))
                str_cell = to_opencl_code_ast_num(m.m[zi+1,zj+1],req_constants)
                print(sb,"$(matrix_name)[$(ocl_mpos(m,zj,zi))] = $(str_cell); \n")
            else
                print(sb,"$(matrix_name)[$(ocl_mpos(m,zj,zi))] = 0.0f; \n")
            end
        end
    end
    return String( take!( sb ) )
end

function to_opencl_code_ast_num( n::ASTSymmOp , req_constants::Dict{ASTNumNode,String})
    str = "( "
    op_str = "$(n.op)"
    for zi=1:length(n.args)
        str = string( str , to_opencl_code_ast_num(n.args[zi] , req_constants) , " " , op_str )
    end
    str = str[1:length(str)-1]
    str = string(str," ) ")
    return str
end

function to_opencl_code_ast_num( n::ASTRealValue , req_constants::Dict{ASTNumNode,String} )
    return " $(n.val)f "
end

function to_opencl_code_ast_num( n::ASTFluxRate , req_constants::Dict{ASTNumNode,String} )
    if( ! haskey( req_constants , n ) )
        req_constants[n] = create_var_name(n)
    end
    return " $(req_constants[n]) "
end
function to_opencl_code_ast_num( n::ASTEMUMassRatio , req_constants::Dict{ASTNumNode,String} )
    if( ! haskey( req_constants , n ) )
        req_constants[n] = create_var_name(n)
    end
    return " $(req_constants[n]) "
end

# the resulting code will compute lu decomposition of the given matrix.
# U will be in place of the original matrix, and L
#
# function assumes that we have a 2d group, at least as large as the matrix
#
function to_opencl_code_ludecomp( matrix_name::String , size::Int)
    sb::IOBuffer = IOBuffer()
    mn  = matrix_name;
    mnl = string(matrix_name,"_L")

    # 1. Init eye matrix for _L
    print(sb,"//LU-DECOMPOSITION (in place): $(matrix_name)_L * $(matrix_name) = $(matrix_name)\n")
    print(sb,"if(idx==0){
            for(int i=0; i<$(size*size); i++){ $(mnl)[i]=0f; }
            for(int i=0; i<$(size); i++){ $(mnl)[i*$(size)+i]=1f; }
        }\n")
    print(sb,"barrier(CLK_LOCAL_MEM_FENCE);\n")
    # 2. start with lu decomp:
    print(sb, "if(idx==0){\n")
    print(sb, "  for(int k=0;k<$(size);k++){\n")
    print(sb, "    for(int ii=k+1;ii<$(size);ii++){ $(mnl)[ k*$(size)+ii ] = $(mn)[ $(mnl)[ k*$(size)+ii ] / $(mn)[ k*$(size)+k ] ; } \n")

    print(sb, "    for(int ii=k+1;ii<$(size);ii++){ for(int ax=0;ax<$(size);ax++) { $(mn)[ ax*$(size)+ii ] -= $(mnl)[ k*$(size)+ii ] * $(mn)[ ax*$(size)+k ] ; } } \n")

    print(sb, "} }\n")

    return String( take!( sb ) )
end

# the resulting code will compute lu decomposition of the given matrix.
# U will be in place of the original matrix, and L
#
# function assumes that we have a 2d group, at least as large as the matrix
#
function to_opencl_code_ludecomp_2( matrix_name::String , size::Int)
    sb::IOBuffer = IOBuffer()
    mn  = matrix_name;
    mnl = string(matrix_name,"_L")
    s   = size

    # 1. Init eye matrix for _L
    print(sb,"//LU-DECOMPOSITION (in place): $(matrix_name)_L * $(matrix_name) = $(matrix_name)\n")
    print(sb,"if( idx<$(s) && idy<$(s) ){ $(mnl)[ idx+$(s)*idy ] = (idx==idy)?1.0f:0.0f; }\nbarrier(CLK_LOCAL_MEM_FENCE);\n")

    # 2. Loop over steps of lu decomp..
    for k in 0:s-1
        print(sb,"//LU-DECOMPOSITION STEP $(k) / $(s-1)\n")
        # process the L submatrix  (k+1:s) x (k)
        # process the U submatrix  (k+1:s) x (0:s)   //(k:s-1)

        # first handle L:
        # work item (0,idy) processes element (idy,k)
        print(sb, "if( idx==0 && idy > $(k) && idy < $(s) ){ \n");
        print(sb, "  $(mnl)[ idy + $(k*s) ] = $(mn)[ idy + $(k*s) ] / $(mn)[ $( k+k*s ) ]; \n");    # k*s is ok here
        print(sb, "} \n");
        print(sb, "barrier(CLK_LOCAL_MEM_FENCE);\n")

        # then handle U:
        # work item (idx,idy) processes element (idx,idy)
        print(sb, "if( idx>$(k) && idx < $(s) && idy < $(s) ){ \n");
        print(sb, "  $(mn)[ idx + $(s)*idy ] -= $(mnl)[ idx + $(s*k) ] * $(mn)[$(k)+ idy*$(s) ]; \n");
        print(sb, "} \n");
        print(sb, "barrier(CLK_LOCAL_MEM_FENCE);\n")

        if(false) # add some debug output..
            print(sb , "\nif(idx==0&&idy==0){ printf(\"After LU It:\\n\"); }\n")
            print(sb , ocl_print_matrix( mn , size , size))
            print(sb , ocl_print_matrix( mnl , size , size))
        end
    end
    print(sb,"//LU-DECOMPOSITION DONE!\n")
    return String( take!( sb ) )
end

# the resulting code will compute lu decomposition of the given matrix.
# U will be in place of the original matrix, and L
#
# function assumes that we have a 1d group, at least as large as the height of the matrix
#
function to_opencl_code_ludecomp_2_1d( matrix_name::String , size::Int)
    sb::IOBuffer = IOBuffer()
    mn  = matrix_name;
    mnl = string(matrix_name,"_L")
    s   = size

    # 1. Init eye matrix for _L
    print(sb,"//LU-DECOMPOSITION (in place): $(matrix_name)_L * $(matrix_name) = $(matrix_name)\n")
    #print(sb,"if( idx<$(s) && idy<$(s) ){ $(mnl)[ idx+$(s)*idy ] = (idx==idy)?1.0f:0.0f; }\nbarrier(CLK_LOCAL_MEM_FENCE);\n")
    print(sb,"if( idx<$(s)){\n")
    for idy_v in 0:s-1
         print(sb,"$(mnl)[ idx+$(s)*$(idy_v) ] = (idx==$(idy_v))?1.0f:0.0f; \n")
     end
     print(sb, "}\nbarrier(CLK_LOCAL_MEM_FENCE);\n" )

    # 2. Loop over steps of lu decomp..
    for k in 0:s-1
        print(sb,"//LU-DECOMPOSITION STEP $(k) / $(s-1)\n")
        # process the L submatrix  (k+1:s) x (k)
        # process the U submatrix  (k+1:s) x (0:s)   //(k:s-1)

        # first handle L:
        # work item (idx) processes element (idx,k)
        #print(sb, "if( idx==0 && idy > $(k) && idy < $(s) ){ \n");
        print(sb, "if( idx > $(k) && idx < $(s) && idy==0 ){ \n");
        print(sb, "  $(mnl)[ idx + $(k*s) ] = $(mn)[ idx + $(k*s) ] / $(mn)[ $( k+k*s ) ]; \n");    # k*s is ok here
        print(sb, "} \n");
        print(sb, "barrier(CLK_LOCAL_MEM_FENCE);\n")

        # then handle U:
        # work item (idx,idy) processes element (idx,idy)

        #print(sb, "if( idx>$(k) && idx < $(s) && idy < $(s) ){ \n");
        print(sb, "if( idx>$(k) && idx < $(s) && idy==0 ){ \n");
        for idy_v in 0:s-1
            #print(sb, "  $(mn)[ idx + $(s)*idy ] -= $(mnl)[ idx + $(s*k) ] * $(mn)[$(k)+ idy*$(s) ]; \n");
            print(sb, "  $(mn)[ idx + $(s)*$(idy_v) ] -= $(mnl)[ idx + $(s*k) ] * $(mn)[$(k)+ $(idy_v)*$(s) ]; \n");
        end
        print(sb, "} \n");
        print(sb, "barrier(CLK_LOCAL_MEM_FENCE);\n")

        if(false) # add some debug output..
            print(sb , "\nif(idx==0&&idy==0){ printf(\"After LU It:\\n\"); }\n")
            print(sb , ocl_print_matrix( mn , size , size))
            print(sb , ocl_print_matrix( mnl , size , size))
        end
    end
    print(sb,"//LU-DECOMPOSITION DONE!\n")
    return String( take!( sb ) )
end

# function assumes that we have a 2d group, at least as large as (mbm x myn)
#
function to_opencl_code_matrixmul( mby::String , mb::String , my::String , b::ASTMatrix , y::ASTMatrix )
    sb::IOBuffer = IOBuffer()
    mbm = size(b.m,1)
    myn = size(y.m,2)
    mbn = size(b.m,2)
    mym = size(y.m,1)

    # 1.
    print(sb,"//MATRIX MUL: $(mby) = $(mb) * $(my)   \n")
    print(sb,"//Work item (idx,idy) compute $(mby)[idx,idy]   \n")

    print(sb,"if( idx < $(mbm) && idy < $(myn) ) {  \n")
    # now this work item just computes this matrix element..
    # unroll.. ;)
    print(sb,"float mulsum = 0.0f;\n")
    for ni in 0:mbn-1
        print(sb,"mulsum += $(mb)[ $(mbm)*$(ni) + idx ] * $(my)[ idy * $(mym) + $(ni) ];\n")
    end
    print(sb,"$(mby)[idx+$(mbm)*idy] = mulsum;\n}\n")
    print(sb,"//MATRIX MUL COMPLETE! \n")

    return String( take!( sb ) )
end

# function assumes that we have a 1d group, at least as large as the height of b
#
function to_opencl_code_matrixmul_1d( mby::String , mb::String , my::String , b::ASTMatrix , y::ASTMatrix )
    sb::IOBuffer = IOBuffer()
    mbm = size(b.m,1)
    myn = size(y.m,2)
    mbn = size(b.m,2)
    mym = size(y.m,1)

    # 1.
    print(sb,"//MATRIX MUL: $(mby) = $(mb) * $(my)   \n")
    print(sb,"//Work item (idx) compute $(mby)[idx,:]   \n")

    #print(sb,"if( idx < $(mbm) && idy < $(myn) ) {  \n")
    print(sb,"if( idx < $(mbm) && idy < $(myn) ) {  \n")
    print(sb,"float mulsum = 0.0f;\n")
    #print(sb,"for(int idy_v=0; idy_v < $(myn);idy_v++ ) {  \n")
    for idy_v in 0:myn-1

        # now this iteration just computes this matrix element..
        # unroll.. ;)
        print(sb,"mulsum = 0.0f;\n")
        for ni in 0:mbn-1
            print(sb,"mulsum += $(mb)[ $(mbm)*$(ni) + idx ] * $(my)[ $(idy_v) * $(mym) + $(ni) ];\n")
        end
        print(sb,"$(mby)[idx+$(mbm)*$(idy_v)] = mulsum;\n")
    end
    print(sb,"}\n")

    print(sb,"//MATRIX MUL COMPLETE! \n")

    return String( take!( sb ) )
end

#
# Solves mL * x = mB
# solution will be in mB
#
# Not very parallel, work item (j) processes column j of mB
function to_opencl_code_solve_L_system( ml::String , mb::String ,  lm::Int  , ln::Int ,  bn::Int )
    sb::IOBuffer = IOBuffer()

    print(sb,"// Solve L System\n")
    #unroll..
    for i in 0:lm-1
        print(sb,"// Solve L System - Row $(i)\n")
        #print(sb,"if(idx==0 && idy<$(bn) ){\n")
        print(sb,"if(idx<$(bn) && idy==0 ){\n")
        print(sb,"  float sum_i = 0.0f;\n")
        for j in 0:(i-1)
            #print(sb,"  sum_i += $(mb)[ $(j) + $(lm)*idy ];\n")       # bm = lm
            #print(sb,"  sum_i += $(mb)[ $(j) + $(lm)*idy ] * $(ml)[ $(i) + $(j*lm) ] ;\n")       # bm = lm
            print(sb,"  sum_i += $(mb)[ $(j) + $(lm)*idx ] * $(ml)[ $(i) + $(j*lm) ] ;\n")       # bm = lm
        end
        #print(sb,"  $(mb)[ $(i) + $(lm)*idy ] = ($(mb)[ $(i) + $(lm)*idy ] - sum_i ) / $(ml)[ $(i+lm*i) ] ;\n")         # bm = lm
        print(sb,"  $(mb)[ $(i) + $(lm)*idx ] = ($(mb)[ $(i) + $(lm)*idx ] - sum_i ) / $(ml)[ $(i+lm*i) ] ;\n")         # bm = lm
        print(sb , "}\n")
        print(sb, "barrier(CLK_LOCAL_MEM_FENCE);\n")
    end
    print(sb,"// Solve L System DONE!\n")
    return String( take!( sb ) )
end

#
# Solves mU * x = mB
# In place, solution will be in mB
#
# Work item (j) processes column j of mB
function to_opencl_code_solve_U_system( mu::String , mb::String ,  um::Int  , un::Int ,  bn::Int )
    sb::IOBuffer = IOBuffer()

    print(sb,"// Solve U System\n")
    #unroll..
    for i_t in 0:(um-1) #(um-1):-1:0
        i = (um-1)-i_t
        print(sb,"// Solve U System - Row $(i)\n")
        #print(sb,"if(idx==0 && idy<$(bn) ){\n")
        print(sb,"if(idx<$(bn) && idy==0 ){\n")
        print(sb,"  float sum_i = 0.0f;\n")
        for j in ((un-1)-i_t+1):(un-1)  # un:-1:( (un-1)-i+1)  # 0:(i-1)
            #print(sb,"  sum_i += $(mb)[ $(j) + $(um)*idy ] * $(mu)[ $(i) + $(j*um) ] ;\n")       # bm = lm
            print(sb,"  sum_i += $(mb)[ $(j) + $(um)*idx ] * $(mu)[ $(i) + $(j*um) ] ;\n")       # bm = lm
        end
        #print(sb,"  $(mb)[ $(i) + $(um)*idy ] = ($(mb)[ $(i) + $(um)*idy ] - sum_i ) / $(mu)[ $(i+um*i) ] ;\n")         # bm = lm
        print(sb,"  $(mb)[ $(i) + $(um)*idx ] = ($(mb)[ $(i) + $(um)*idx ] - sum_i ) / $(mu)[ $(i+um*i) ] ;\n")         # bm = lm
        print(sb , "}\n")
        print(sb, "barrier(CLK_LOCAL_MEM_FENCE);\n")
    end
    print(sb,"// Solve U System DONE!\n")
    return String( take!( sb ) )
end

function ocl_print_matrix( name::String , sm::Int , sn::Int)
    return "if(idx==0&&idy==0){
              printf(\"\\n$(name) : $(sm)x$(sn) \\n\");
              for(int cm=0;cm<$(sm);cm++){ for( int cn=0;cn<$(sn);cn++){ printf( \" %6.3f \" , $(name)[cm+$(sm)*cn] ); }
                  printf(\"\\n\");
              }
              printf(\"\\n\");
             }
            "
end

function test_0()

    a = rand(Float32, 50_000)
    b = rand(Float32, 50_000)
    device, ctx, queue = OpenCL.cl.create_compute_context()

    a_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=a)
    b_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=b)
    c_buff = cl.Buffer(Float32, ctx, :w, length(a))

    p = cl.Program(ctx, source=test_source) |> cl.build!
    #cl.CL_device_info(1)

    k = cl.Kernel(p, "fwsim")
    queue(k, [2000], nothing, a_buff, b_buff, c_buff)
    r = cl.read(queue, c_buff)
end


text:
module PooksoftTextMagicAPI

    # include my codes -
    include("Include.jl")

    # export types -
    export PSResult
    export PSTextMagicAPIUserObject

    # export methods -
    export load_config_file
    export send_text_message
    export build_user_model_object
    export load_csv_data_file

    # filestore watcher methods -
    export register_file_store_watcher
    export unregister_file_store_watcher

end # module

ext: jl


text:
# TODO stop-gap string storage of Distrubtion types, should be upgraded to more efficient storage
function normalfromstring(str::AS) where {AS <: AbstractString}
  meanstr = match(r"μ=[+-]?([0-9]*[.])?[0-9]+", str).match
  mean = split(meanstr, '=')[2]
  sigmastr = match(r"σ=[+-]?([0-9]*[.])?[0-9]+", str).match
  sigma = split(sigmastr, '=')[2]
  Normal{Float64}(parse(Float64,mean), parse(Float64,sigma))
end

function mvnormalfromstring(str::AS) where {AS <: AbstractString}
  means = split(split(split(str, 'μ')[2],']')[1],'[')[end]
  mean = Float64[]
  for ms in split(means, ',')
    push!(mean, parse(Float64, ms))
  end
  sigs = split(split(split(str, 'Σ')[2],']')[1],'[')[end]
  sig = Float64[]
  for ms in split(sigs, ';')
    for m in split(ms, ' ')
      length(m) > 0 ? push!(sig, parse(Float64, m)) : nothing
    end
  end
  len = length(mean)
  sigm = reshape(sig, len,len)
  MvNormal(mean, sigm)
end

function categoricalfromstring(str::AS)::Distributions.Categorical where {AS <: AbstractString}
  # pstr = match(r"p=\[", str).match
  psubs = split(str, '=')[end]
  psubs = split(psubs, '[')[end]
  psubsub = split(psubs, ']')[1]
  pw = split(psubsub, ',')
  p = parse.(Float64, pw)
  return Categorical(p ./ sum(p))
end

function extractdistribution(str::AS)::Union{Nothing, SamplableBelief} where {AS <: AbstractString}
  # TODO improve use of multidispatch and packing of Distribution types
  if str == ""
    return nothing
  elseif (occursin(r"Normal", str) && !occursin(r"FullNormal", str))
    return normalfromstring(str)
  elseif occursin(r"FullNormal", str)
    return mvnormalfromstring(str)
  elseif occursin(r"Categorical", str)
    return categoricalfromstring(str)
  elseif occursin(r"DiscreteNonParametric", str)
    return categoricalfromstring(str)
  elseif occursin(r"KDE:", str)
    return convert(BallTreeDensity, str)
  elseif occursin(r"AliasingScalarSampler", str)
    return convert(AliasingScalarSampler, str)
  else
    error("Don't know how to extract distribution from str=$(str)")
  end
end

function extractdistributionJson(packedFactor::Dict{String, Any})::Union{Void, Distributions.Distribution}
  # TODO improve use of multidispatch and packing of Distribution types
  error("Not implemented yet")
end


text:
for op in (:(+),:(-),:(*),:(/))
  @eval begin
     function ($op){G1<:Grasp,G2<:Grasp,R1<:Real,R2<:Real}(a::Rvl{G1,R1}, b::Rvl{G2,R2})
         R = promote_type(R1,R2)
         alo,blo = promote(a.lo,b.lo)
         ahi,bhi = promote(a.hi,b.hi)
         ($op)(Rvl{G1,R}(alo,ahi), Rvl{G2,R}(blo,bhi))
     end
  end
end


text:
# # @__NAME__

# PREAMBLE

# PKG_SETUP

# ## Setup

# In addition to `RigidBodyDynamics`, we'll be using the `StaticArrays` package, used throughout `RigidBodyDynamics`, which provides stack-allocated, fixed-size arrays:

using RigidBodyDynamics
using LinearAlgebra
using StaticArrays

# ## Creating a double pendulum `Mechanism`

# We're going to create a simple `Mechanism` that represents a [double pendulum](https://en.wikipedia.org/wiki/Double_pendulum). The `Mechanism` type can be thought of as an interconnection of rigid bodies and joints.
#
# We'll start by creating a 'root' rigid body, representing the fixed world, and using it to create a new `Mechanism`:

g = -9.81 # gravitational acceleration in z-direction
world = RigidBody{Float64}("world")
doublependulum = Mechanism(world; gravity = SVector(0, 0, g))

# Note that the `RigidBody` type is parameterized on the 'scalar type', here `Float64`.
#
# We'll now add a second body, called 'upper link', to the `Mechanism`. We'll attach it to the world with a revolute joint, with the $y$-axis as the axis of rotation. We'll start by creating a `SpatialInertia`, which stores the inertial properties of the new body:

axis = SVector(0., 1., 0.) # joint axis
I_1 = 0.333 # moment of inertia about joint axis
c_1 = -0.5 # center of mass location with respect to joint axis
m_1 = 1. # mass
frame1 = CartesianFrame3D("upper_link") # the reference frame in which the spatial inertia will be expressed
inertia1 = SpatialInertia(frame1,
    moment=I_1 * axis * axis',
    com=SVector(0, 0, c_1),
    mass=m_1)

# Note that the created `SpatialInertia` is annotated with the frame in which it is expressed (in the form of a `CartesianFrame3D`). This is a common theme among `RigidBodyDynamics` objects. Storing frame information with the data obviates the need for the complicated variable naming conventions that are used in some other libraries to disambiguate the frame in which quantities are expressed. It also enables automated reference frame checks.

# We'll now create the second body:

upperlink = RigidBody(inertia1)

# and a new revolute joint called 'shoulder':

shoulder = Joint("shoulder", Revolute(axis))

# Creating a `Joint` automatically constructs two new `CartesianFrame3D` objects: a frame directly before the joint, and one directly after. To attach the new body to the world by this joint, we'll have to specify where the frame before the joint is located on the parent body (here, the world):

before_shoulder_to_world = one(Transform3D,
    frame_before(shoulder), default_frame(world))

# Now we can attach the upper link to the world:

attach!(doublependulum, world, upperlink, shoulder,
    joint_pose = before_shoulder_to_world)

# which changes the tree representation of the `Mechanism`.

# We can attach the lower link in similar fashion:

l_1 = -1. # length of the upper link
I_2 = 0.333 # moment of inertia about joint axis
c_2 = -0.5 # center of mass location with respect to joint axis
m_2 = 1. # mass
inertia2 = SpatialInertia(CartesianFrame3D("lower_link"),
    moment=I_2 * axis * axis',
    com=SVector(0, 0, c_2),
    mass=m_2)
lowerlink = RigidBody(inertia2)
elbow = Joint("elbow", Revolute(axis))
before_elbow_to_after_shoulder = Transform3D(
    frame_before(elbow), frame_after(shoulder), SVector(0, 0, l_1))
attach!(doublependulum, upperlink, lowerlink, elbow,
    joint_pose = before_elbow_to_after_shoulder)

# Now our double pendulum `Mechanism` is complete.

# **Note**: instead of defining the `Mechanism` in this way, it is also possible to load in a [URDF](http://wiki.ros.org/urdf) file (an XML file format used in ROS), using the `parse_urdf` function, e.g.:

srcdir = dirname(pathof(RigidBodyDynamics))
urdf = joinpath(srcdir, "..", "test", "urdf", "Acrobot.urdf")
parse_urdf(urdf)

# ## The state of a `Mechanism`

# A `Mechanism` stores the joint/rigid body layout, but no state information. State information is separated out into a `MechanismState` object:

state = MechanismState(doublependulum)

# Let's first set the configurations and velocities of the joints:

set_configuration!(state, shoulder, 0.3)
set_configuration!(state, elbow, 0.4)
set_velocity!(state, shoulder, 1.)
set_velocity!(state, elbow, 2.);

# **Important**: a `MechanismState` contains cache variables that depend on the configurations and velocities of the joints. These need to be invalidated when the configurations and velocities are changed. To do this, call

setdirty!(state)

# The joint configurations and velocities are stored as `Vector`s (denoted $q$ and $v$ respectively in this package) inside the `MechanismState` object:

q = configuration(state)
v = velocity(state)

# ## Kinematics

# We are now ready to do kinematics. Here's how you transform a point at the origin of the frame after the elbow joint to world frame:

transform(state, Point3D(frame_after(elbow), zero(SVector{3})),
    default_frame(world))

# Other objects like `Wrench`es, `Twist`s, and `SpatialInertia`s can be transformed in similar fashion.

# You can also ask for the homogeneous transform to world:

transform_to_root(state, frame_after(elbow))

# Or a relative transform:

relative_transform(state, frame_after(elbow), frame_after(shoulder))

# and here's the center of mass of the double pendulum:

center_of_mass(state)

# ## Dynamics

# A `MechanismState` can also be used to compute quantities related to the dynamics of the `Mechanism`. Here we compute the mass matrix:

mass_matrix(state)

# Note that there is also a zero-allocation version, `mass_matrix!` (the `!` at the end of a method is a Julia convention signifying that the function is 'in-place', i.e. modifies its input data).

# We can do inverse dynamics as follows (note again that there is a non-allocating version of this method as well):

v̇ = similar(velocity(state)) # the joint acceleration vector, i.e., the time derivative of the joint velocity vector v
v̇[shoulder][1] = 1
v̇[elbow][1] = 2
inverse_dynamics(state, v̇)

# ## Simulation

# Let's simulate the double pendulum for 5 seconds, starting from the state we defined earlier. For this, we can use the basic `simulate` function:

ts, qs, vs = simulate(state, 5., Δt = 1e-3);

# `simulate` returns a vector of times (`ts`) and associated joint configurations (`qs`) and velocities (`vs`). You can of course plot the trajectories using your favorite plotting package (see e.g. [Plots.jl](https://github.com/JuliaPlots/Plots.jl)). The [MeshCatMechanisms](https://github.com/JuliaRobotics/MeshCatMechanisms.jl) or [RigidBodyTreeInspector](https://github.com/rdeits/RigidBodyTreeInspector.jl) packages can also be used for 3D animation of the double pendulum in action. See also [RigidBodySim.jl](https://github.com/JuliaRobotics/RigidBodySim.jl) for a more full-fledge simulation environment.

# A lower level interface for simulation/ODE integration with more options is also available.
# Consult the documentation for more information.
# In addition, [RigidBodySim.jl](https://github.com/JuliaRobotics/RigidBodySim.jl) offers a more full-featured simulation environment.

max_stars_repo_path: examples/1. Quickstart - double pendulum/1. Quickstart - double pendulum.jl


text:
using MadMax
using Test

# BisectionEval
# RegulaFalsiEval
# NewtonEval
# GaussJacobiEval
# GaussSeidelEval
# SOREval
# SOREval2
# EulerODE
# ModifiedEulerODE
# RK2ODE
# RK4ODE

X = [5 5 15 -40
      5 15 35 -104
      15 35 99 -290]

X2 = [3  2  0 -5
      2  3 -1 -4
      0 -1  2 -1]

x, a = GaussSeidelEval(X2, [0.9;-3.1;0.9], N=3)
x, a = GaussJacobiEval(X, [0;0;0], N=7)
x, a = SOREval(X2, [0;0;0], 1.25, N=3)
x, a = SOREval2(X2, [0;0;0], 0.8, N=7)

f(x) = x^24-2.6x+1.6
RegulaFalsiEval(f, 0, 1, tol=0.5*10^-3)
NewtonEval(f, 0.05, N=4, tol=0.5*10^-4)

f(t, x) = -x*sin(t) + sin(t)
ModifiedEulerODE(f, 0, 2+ℯ, 2π, 2π/10)
RK2ODE(f, 0, 2+ℯ, 2π, 2π/10)


text:
module Multirate

include( "enum.jl" )

import Base: filt, filt!, reset
import Polynomials: Poly, polyval

using DSP.Windows
export  hanning,
        hammming,
        kaiser,
        blackman

include( "support.jl" )
include( "FIRDesign.jl" )
export  firdes,
        kaiserlength,
        FIRResponse,
        LOWPASS,
        HIGPASS,
        BANDPASS,
        BANDSTOP

include( "Filters.jl" )
export  FIRFilter,
        FIRInterpolator,
        FIRArbitrary,
        FIRDecimator,
        FIRFarrow,
        FIRRational,
        FIRStandard,
        filt!,
        filt,
        setphase,
        tapsforphase!,
        tapsforphase,
        taps2pfb,
        reset,
        outputlength,
        inputlength

include( "NaiveResamplers.jl" )

end # module


text:
@static if VERSION < v"0.7-"
    const IOError = Base.UVError
else
    using Base: IOError
end

Base.open(vis::Visualizer, args...; kw...) = open(vis.core, args...; kw...)

function Base.open(core::CoreVisualizer, port::Integer)
    @async WebIO.webio_serve(Mux.page("/", req -> core.scope), port)
    url = "http://127.0.0.1:$port"
    Compat.@info("Serving MeshCat visualizer at $url")
    open_url(url)
end

function open_url(url)
    try
        if Compat.Sys.iswindows()
            run(`start $url`)
        elseif Compat.Sys.isapple()
            run(`open $url`)
        elseif Compat.Sys.islinux()
            run(`xdg-open $url`)
        end
    catch e
        println("Could not open browser automatically: $e")
        println("Please open the following URL in your browser:")
        println(url)
    end
end

function Base.open(core::CoreVisualizer; default_port=8700, max_retries=500)
    for port in default_port:(default_port + max_retries)
        server = try
            listen(port)
        catch e
            if e isa IOError
                continue
            end
        end
        close(server)
        # It is *possible* that a race condition could occur here, in which
        # some other process grabs the given port in between the close() above
        # and the open() below. But it's unlikely and would not be terribly
        # damaging (the user would just have to call open() again).
        return open(core, port)
    end
end


text:
__precompile__()

module Struve

include("fast.jl")
include("integral.jl")
include("large_arg.jl")
include("small_arg.jl")

# define defaults, should be exact
const H0 = _H0_integral
const K0 = _K0_integral
const L0 = _L0_integral
const M0 = _M0_integral

const H = _H_integral
const K = _K_integral
const L = _L_integral
const M = _M_integral

# define fast, less exact alternatives
const H0_fast = _H0_fast
const H1_fast = _H1_fast

end # module


text:
function landau_mignotte(p::AbstractPolynomial)
    l = ceil(Int, hypot(coefficients(p)...))
    d = deg(p) ÷ 2
    maximum([binomial(d-1,j)*l + binomial(d-1,j-1)*leading(p) for j = 1:d])
end

function factor_distinct_degree(p::Polynomial{true, ℤₚ{n}}) where n
    f = FactoredPoly()
    p₀ = p
    x = modular(n, var(p))
    d = 1
    while deg(p) > 1 # && d <= deg(p) ÷ 2
        # q = modular(n, x^(n^d) - x)
        q = modpow(x, n^d, p₀) - x
        g = gcd(p, q)

        if !isone(g)
            add_factor!(f, g, d)
            p ÷= g
        end
        d += 1
    end
    # add_factor!(f, p, deg(p))
    f
end

@polyvar 𝑢

# function factor_equal_degree(p::Polynomial{true, ℤₚ{n}}, d, i₀=0) where n
#     f = FactoredPoly()
#     p₀ = p
#     # x = var(p)
#     x = var(p)
#     e = (n^d - 1) ÷ 2
#
#     for i = i₀:n
#         a = modular(n, x+i)
#         # a = modular(n, ((i-1)÷n)*x^2 + x + ((i-1)%n))
#         q = modpow(a, e, p₀) - one(p₀)
#         g = gcd(p, q)
#         # qᵤ = modular(n, 𝑢^e - 1)
#         # pᵤ = p(x => 𝑢 - i)
#         # gᵤ = gcd(pᵤ, qᵤ)
#         # g = gᵤ(𝑢 => x + i)
#         g ÷= leading(g)
#
#         if !isone(g)
#             if deg(g) == d
#                 add_factor!(f, g, 1)
#             else
#                 f₂ = factor_equal_degree(g, d, i+1)
#                 combine!(f, f₂)
#             end
#             p ÷= g
#         end
#
#         isone(p) && return f
#     end
#     return f
# end
#
# function factor_equal_degree(f::FactoredPoly)
#     h = FactoredPoly()
#
#     for w in factors(f)
#         p, d = first(w), last(w)
#         if deg(p) == d
#             add_factor!(h, p, 1)
#         else
#             combine!(h, factor_equal_degree(p, d))
#         end
#     end
#     h
# end

function is_square_free(p::AbstractPolynomial)
    p′ = derivative(p)
    g = gcd(p, p′)
    return isone(g)
end

# function find_integer_factorization(p::AbstractPolynomial, f)
#     m = length(f)
#     S = Set{Int}()
#     fs = []
#
#     d = 1
#     while length(S) < m && d <= deg(p)÷2
#         for a in Iterators.product([1:m for i=1:d]...)
#             if length(unique(a)) == d
#                 if !any(i ∈ S for i in a)
#                     ρ = prod(f[i] for i in a)
#                     if iszero(p % ρ)
#                         push!(fs, ρ)
#                         p = remove_factor(p, ρ)
#                         for i in a
#                             push!(S, i)
#                         end
#                     end
#                 end
#             end
#         end
#         d += 1
#     end
#
#     push!(fs, p)
#     fs
# end

function factor_modular(p::AbstractPolynomial, n)
    q = modular(n, p)
    # f = factor_distinct_degree(q)
    # f = factor_equal_degree(f)
    f = factor_modular(q)
    [demodular(v) for v in first.(factors(f))]
end

# p should be integer and monic
function factor_roundabout(p::AbstractPolynomial, n)
    !isprime(n) && error("$n is not a prime!")
    lc = leading(p)
    # p, undo = standard_form(p)
    undo = x->x
    λ = landau_mignotte(p)
    x = var(p)

    f = FactoredPoly()
    p₁ = one(p)

    for w in factors(decompose(p))
        v, k = first(w), last(w)

        if deg(v, x) > 0
            f₁ = factor_modular(v, n)
            f₂ = lift(v, f₁, n, 2λ)
            f₃ = find_integer_factorization(v, f₂)

            for u in f₃
                if deg(u) > 0
                    w₁ = undo(u)
                    add_factor!(f, w₁, k)
                    p₁ = p₁ * w₁^k
                end
            end
        end
    end

    ρ = lc // leading(p₁)
    if denominator(ρ) == 1
        ρ = numerator(ρ)
    end

    if !isone(ρ)
        add_factor!(f, ρ, 1)
    end

    f
end

function factor_roundabout(eq, n)
    p, v = prewrap(eq)
    unwrap(factor_roundabout(p, n), v)
end

function factor_roundabout(p::AbstractPolynomial; n=3, max_prime=50)
    n >= max_prime && return nothing
    n = nextprime(n)

    f₁ = factor_roundabout(p, n)

    i = 1
    while length(f₁) == 1 && i < 4
        n = nextprime(n+1)
        f₁ = factor_roundabout(p, n)
        i += 1
    end

    length(f₁) == 1 && return f₁

    f₂ = FactoredPoly()
    k = 1

    for w₁ in factors(f₁)
        v₁, k₁ = first(w₁), last(w₁)
        if deg(v₁) == 0
            k *= v₁^k₁
        elseif deg(v₁) == 1
            k *= cont(v₁)^k₁
            add_factor!(f₂, prim(v₁), k₁)
        else
            f₃ = factor_roundabout(v₁; n=n+1, max_prime=max_prime)
            if f₃ != nothing
                for w₃ in factors(f₃)
                    v₃, k₃ = first(w₃), last(w₃)
                    if deg(v₃) == 0
                        k *= v₃^(k₁*k₃)
                    else
                        add_factor!(f₂, v₃, k₁*k₃)
                    end
                end
            end
        end
    end

    if !isone(k)
        add_factor!(f₂, k, 1)
    end

    return f₂
end

factor_roundabout(eq) = wrap(factor_roundabout, eq)

##############################################################################

function remove_factor(p::AbstractPolynomialLike, f)
    q, r = divrem(rationalize(p), f)
    !iszero(r) && error("$f is not a proper factor of $p")
    integer_poly(q)
end

function lift(p::AbstractPolynomial, s₁, t₁, n, λ)
    sₙ = modular(n, s₁)
    tₙ = modular(n, t₁)
    _, σₙ, τₙ = gcdx(sₙ, tₙ)

    s, t = s₁, t₁
    d = p - s*t
    cₙ = modular(n, d ÷ n)

    i = 2
    while n^i < 2λ && !iszero(d)
        s̄ₙ = rem(τₙ*cₙ, modular(n, s))
        t̄ₙ = rem(σₙ*cₙ, modular(n, t))
        s += n^(i-1) * demodular(s̄ₙ)
        t += n^(i-1) * demodular(t̄ₙ)
        d = p - s*t
        cₘ = modular(n, d ÷ n^i)
        i += 1
    end

    f = FactoredPoly()
    add_factor!(f, s)
    add_factor!(f, t)
    f
end

function lift(p::AbstractPolynomial, f, m, λ)
    nf = length(f)
    if nf <= 1
        return f
    elseif nf == 2
        s, t = f[1], f[2]
        return lift(p, s, t, m, λ)
    else
        s, t = f[1], integer_poly(prod(f[i] for i=2:nf))
        h₁ = lift(p, s, t, m, λ)
        s = h₁[1]
        # h₂ = lift(remove_factor(p,s), f[2:nf], m, λ)
        h₂ = lift(t, f[2:nf], m, λ)
        add_factor!(h₂, s)
        return h₂
    end
end

function modpow(p::AbstractPolynomialLike, k::Integer, q::AbstractPolynomialLike)
    x = rem(p, q)
    r = one(p)
    while k > 0
        if isodd(k)
            r = rem(r*x, q)
        end
        x = rem(x*x, q)
        k >>= 1
    end
    r
end

##############################################################################

# p is a modular polynomial
function factor_modular(p::Polynomial{true, ℤₚ{n}}) where n
    x = modular(n, var(p))
    h = x
    v = p ÷ leading(p)
    f = FactoredPoly()
    i = 0

    while deg(v) >= 2*(i+1) && i < n
        i += 1
        h = modpow(h, n, p)
        g = gcd(h - x, v)

        if !isone(g)
            f₁ = factor_equal_degree(g, i)
            if f₁ == nothing return nothing end
            for gₛ in f₁
                if deg(gₛ) > 0
                    k = 0
                    while iszero(v % gₛ)
                        k += 1
                        v ÷= gₛ
                    end
                    if k > 0
                        add_factor!(f, gₛ, k)
                    end
                end
            end
        end
    end
    f
end

function factor_equal_degree(p::Polynomial{true, ℤₚ{n}}, d) where n
    m = deg(p)
    m == d && return [p]

    g = nothing
    i = 0
    while g == nothing || !(1 < deg(g) < m)
        g = split_equal_degree(p, d)
        i += 1
        if i == 10
            return nothing
        end
    end

    if deg(g) < m
        f₁ = factor_equal_degree(g, d)
        f₂ = factor_equal_degree(p÷g, d)
        if f₁ != nothing && f₂ != nothing
            return [f₁; f₂]
        elseif f₁ == nothing && f₂ != nothing
            return f₂
        elseif f₁ != nothing && f₂ == nothing
            return f₁
        else
            return nothing
        end
    end
    return [p]
end

function split_equal_degree(p::Polynomial{true, ℤₚ{n}}, d) where n
    x = modular(n, var(p))
    a = x + rand(0:n-1)
    g = gcd(p, a)
    !isone(g) && return g ÷ leading(g)
    b = modpow(a, (n^d-1)÷2, p)
    g = gcd(p, b - 1)
    !isone(g) && !iszero(g) && return g ÷ leading(g)
    nothing
end

function traverse_patterns(pat, mask, n, k, fun)
    bit = 1

    for i = 1:min(trailing_zeros(pat), n)
        pat₁ = pat | bit
        if (pat != pat₁) && (mask & pat₁ == 0)
            nb = count_ones(pat₁)
            if nb == k
                mask |= fun(pat₁)
            elseif nb < k
                mask = traverse_patterns(pat₁, mask, n, k, fun)
            end
        end
        bit <<= 1
    end
    return mask
end

function find_integer_factorization(p::AbstractPolynomial, f::FactoredPoly)
    p₀ = copy(p)
    x = var(p)
    Φ = first.(f.factors)
    η = last.(f.factors)
    n = length(η)
    f₁ = FactoredPoly()

    fun = function(pat)
        l = [i for i=1:n if testbit(pat, i)]
        q = demodular(prod(Φ[l]))
        k = 0
        while iszero(p % q)
            k += 1
            η[l] .-= 1
            p = integer_poly(p ÷ q)
        end
        if k > 0
            add_factor!(f₁, prim(q), k)
            return sum(1<<(i-1) for i=1:n if testbit(pat,i) && η[i]<=0; init=0)
        end
        return 0
    end

    mask = 0
    println("n = ", n)
    for k = 1:n
        # if binomial(n - count_ones(mask), k) > 100000
        #     break
        # end
        mask = traverse_patterns(0, mask, n, k, fun)
    end

    integer_poly(p₀ ÷ poly(f₁)), f₁
end

##############################################################################

function factor_combined(p::AbstractPolynomialLike; N=20, first_prime=3)
    c₀, p = integer_poly_coef(p)
    p₀ = copy(p)
    x = var(p)
    f₀ = FactoredPoly()
    fc = 0  # failure count
    a = first_prime

    while deg(p, x) > N
        println(p)
        a = nextprime(a+1)
        println(a)
        q = modular(a, p)
        # if !isone(gcd(q, derivative(q))) continue end
        println(q)
        f = factor_modular(q)
        println(f)

        if f != nothing
            p, f = find_integer_factorization(p, f)

            if length(f) == 0
                fc += 1
                if fc == 5
                    break
                end
            else
                fc = 0
            end

            combine!(f₀, f)
        end
    end

    if 1 < deg(p) <= N
        f = factor_roots_comb(p)
        combine!(f₀, f)
    end

    c, ρ = integer_poly_coef(p₀ ÷ poly(f₀))
    !isone(prim(ρ)) && add_factor!(f₀, prim(ρ), 1)
    c *= cont(ρ) * c₀

    if denominator(c) == 1
         c = numerator(c)
    end

    !isone(c) && add_factor!(f₀, c, 1)

    f₀
end


text:
"""
    coarseguess!(cache::PararealCache, problem::AbstractInitialValueProblem, parareal::Parareal)

computes the coarse solution of an [`AbstractInitialValueProblem`](@ref), and saves the outcome into a [`PararealCache`](@ref).
"""
function coarseguess!(cache::PararealCache, problem::AbstractInitialValueProblem, parareal::Parareal)
    @↓ U, T, G = cache
    @↓ u0, (t0, tN) ← tspan = problem
    @↓ coarsolver, P = parareal
    T[1] = t0
    for n = 1:P
        # stable sum
        T[n + 1] = (1 - n) * t0 / P + n * tN / P
    end
    G[1] = u0
    for n = 1:P
        chunkproblem = subproblemof(problem, G[n], T[n], T[n+1])
        chunksolution = coarsolver(chunkproblem)
        G[n + 1] = chunksolution.u[end]
    end
    U .= G
    return cache
end


text:
module TestStrategy

using Printf
import SearchBall, GeometryTypes, LinearAlgebra
const sb = SearchBall
const gt = GeometryTypes
const la = LinearAlgebra

function init_strategies(config; kwargs...)
    sb.init_text!(config, "dist:opp-off")
    sb.init_text!(config, "dist:opp-ball")
    sb.init_text!(config, "dist:off-ball")
    sb.add_slider!(config, "c1", 0:0.05:1, 0.7)
end

v = 0.1 # m/s
function calc_myPlayers(config::sb.Config, dt::Number, set_type::Symbol, moved_object::Union{sb.Object{<:Number}, Nothing}=nothing)
    ws = config.ws

    # push!(config.printables["dist:opp-off"], @sprintf("%.2fm", sb.distance(ws.opponents[1], ws.myPlayers[1]))) # debug output
    sb.show_text!(config, "dist:opp-off", @sprintf("%.2fm", sb.distance(ws.opponents[1], ws.myPlayers[1]))) # debug output
    sb.show_text!(config, "dist:opp-ball", @sprintf("%.2fm", sb.distance(ws.opponents[1], ws.obstacles[1]))) # debug output
    sb.show_text!(config, "dist:off-ball", @sprintf("%.2fm", sb.distance(ws.obstacles[1], ws.myPlayers[1]))) # debug output
    if sb.is_in_line_of_sight(ws, ws.myPlayers[1], ws.obstacles[1])
        sb.send_order!(ws.myPlayers[1], sb.GoToAbsPosition(ws.myPlayers[1] .+ la.normalize(sb.vector(ws.myPlayers[1], ws.obstacles[1])) .* dt .* v))
    else
        t = sb.outer_tangent(gt.Circle(gt.Point2(ws.myPlayers[1]), dt*v), gt.Circle(gt.Point2(ws.opponents[1]), sb.get_size(ws.opponents[1])))
        sb.send_order!(ws.myPlayers[1], sb.GoToAbsPosition([x for x in t[1]]))
    end

    return
end

# c1, c2 = 0.5, 0.5
function calc_opponend(config::sb.Config, dt::Number, set_type::Symbol, moved_object::Union{sb.Object{<:Number}, Nothing}=nothing)
    ws = config.ws

    # parameter decision
    if sb.is_in_line_of_sight(ws, ws.myPlayers[1], ws.obstacles[1])
        c1, c2 = 0.0, 1.0
    else
        # c1, c2 = 0.7, 0.3
        c1 = sb.get_manip_value(config, "c1")
        c2 = 1.0 - c1
    end
    # helping lines
    a = la.normalize(sb.vector(ws.opponents[1], ws.myPlayers[1]))
    offender_ball = sb.vector(ws.myPlayers[1], ws.obstacles[1])
    n_off_ball = sb.normal2d(offender_ball)
    crossing, isec = sb.intersection(gt.LineSegment(gt.Point2(ws.myPlayers[1]), gt.Point2(ws.obstacles[1])), gt.LineSegment(gt.Point2(ws.opponents[1]), gt.Point2(ws.opponents[1].+n_off_ball)))

    sb.show_graphic!(config, "opp1", [ws.opponents[1], isec], sb.GraphicArrow(arrowloc=1.0, color=[0.0, 0.8, 0.2, 0.6]))

    b = la.normalize(sb.vector(ws.opponents[1], isec))
    if any(isnan.(b))
        b = gt.Point2(0.0) # prevents error if `sb.vector(ws.opponents[1], isec)` is a zero sb.vector, becaus normalize evaluates zero sb.vectors to NaN
    end

    g = c1.* a .+ c2.* b
    sb.send_order!(ws.opponents[1], sb.GoToAbsPosition(ws.opponents[1] .+ la.normalize(g) .* dt .* v))

    return
end

function run_strategies(config::sb.Config, dt::Number=0, set_type::Symbol=:default, moved_object::Union{sb.Object{<:Number}, Nothing}=nothing)
    calc_myPlayers(config, dt, set_type, moved_object)
    calc_opponend(config, dt, set_type, moved_object)
    return
end

end


text:
###
### parsing
###

const OPTIONS = Parsers.Options(
    sentinel=missing,
    quoted=true,
    openquotechar='\'',
    closequotechar='\'',
    delim=',',
    stripwhitespace=true,
)

getbytes(source::Vector{UInt8}) = source, 1, length(source)
getbytes(source::IOBuffer) = source.data, source.ptr, source.size
getbytes(source) = getbytes(read(source))

"""
    parse_network(source) -> Network

Read a PSS/E-format `.raw` Power Flow Data file and return a [`Network`](@ref) object.

The version of the PSS/E format can be specified with the `v` keyword, like `v=33`,
or else it will be automatically detected when parsing the file.
"""
function parse_network(source; v::Union{Integer,Nothing}=nothing)
    @debug 1 "source = $source, v = $v"
    bytes, pos, len = getbytes(source)

    caseid, pos = parse_idrow(CaseID, bytes, pos, len, OPTIONS)
    @debug 1 "Parsed CaseID: rev = $(caseid.rev), pos = $pos"
    # when `v` not given, if `caseid.rev` missing we assume it is because data is v30 format
    version = something(v, coalesce(caseid.rev, 30))
    is_v33 = version == 33
    @debug 1 "Set version = $version"

    # Skip the 2 lines of comments
    pos = next_line(bytes, pos, len)
    pos = next_line(bytes, pos, len)
    @debug 1 "Parsed comments: pos = $pos"
    return if is_v33
        parse_network33(source, version, caseid, bytes, pos, len, OPTIONS)
    else
        parse_network30(source, version, caseid, bytes, pos, len, OPTIONS)
    end
end

function parse_network33(source, version, caseid, bytes, pos, len, options)
    buses, pos = parse_records!(Buses33(len÷1000), bytes, pos, len, options)
    nbuses = length(buses)
    loads, pos = parse_records!(Loads(nbuses), bytes, pos, len, options)
    fixed_shunts, pos = parse_records!(FixedShunts(), bytes, pos, len, options)
    gens, pos = parse_records!(Generators(nbuses÷10), bytes, pos, len, options)
    ngens = length(gens)
    branches, pos = parse_records!(Branches33(nbuses), bytes, pos, len, options)
    transformers, pos = parse_records!(Transformers(ngens*2), bytes, pos, len, options)
    interchanges, pos = parse_records!(AreaInterchanges(), bytes, pos, len, options)
    two_terminal_dc, pos = parse_records!(TwoTerminalDCLines33(), bytes, pos, len, options)
    vsc_dc, pos = parse_records!(VSCDCLines(), bytes, pos, len, options)
    impedance_corrections, pos = parse_records!(ImpedanceCorrections(), bytes, pos, len, options)
    multi_terminal_dc, pos = parse_records!(MultiTerminalDCLines{DCLineID33}(), bytes, pos, len, options)
    multi_section_lines, pos = parse_records!(MultiSectionLineGroups33(), bytes, pos, len, options)
    zones, pos = parse_records!(Zones(), bytes, pos, len, options)
    area_transfers, pos = parse_records!(InterAreaTransfers(), bytes, pos, len, options)
    owners, pos = parse_records!(Owners(), bytes, pos, len, options)
    facts, pos = parse_records!(FACTSDevices33(), bytes, pos, len, options)
    switched_shunts, pos = parse_records!(SwitchedShunts33(nbuses÷11), bytes, pos, len, options)
    return Network(
        version,
        caseid,
        buses,
        loads,
        fixed_shunts,
        gens,
        branches,
        transformers,
        interchanges,
        two_terminal_dc,
        vsc_dc,
        switched_shunts,
        impedance_corrections,
        multi_terminal_dc,
        multi_section_lines,
        zones,
        area_transfers,
        owners,
        facts,
    )
end

function parse_network30(source, version, caseid, bytes, pos, len, options)
    buses, pos = parse_records!(Buses30(len÷1000), bytes, pos, len, options)
    nbuses = length(buses)
    loads, pos = parse_records!(Loads(nbuses), bytes, pos, len, options)
    fixed_shunts = nothing
    gens, pos = parse_records!(Generators(nbuses÷10), bytes, pos, len, options)
    ngens = length(gens)
    branches, pos = parse_records!(Branches30(nbuses), bytes, pos, len, options)
    transformers, pos = parse_records!(Transformers(ngens*2), bytes, pos, len, options)
    interchanges, pos = parse_records!(AreaInterchanges(), bytes, pos, len, options)
    two_terminal_dc, pos = parse_records!(TwoTerminalDCLines30(), bytes, pos, len, options)
    vsc_dc, pos = parse_records!(VSCDCLines(), bytes, pos, len, options)
    switched_shunts, pos = parse_records!(SwitchedShunts30(nbuses÷11), bytes, pos, len, options)
    impedance_corrections, pos = parse_records!(ImpedanceCorrections(), bytes, pos, len, options)
    multi_terminal_dc, pos = parse_records!(MultiTerminalDCLines{DCLineID30}(), bytes, pos, len, options)
    multi_section_lines, pos = parse_records!(MultiSectionLineGroups30(), bytes, pos, len, options)
    zones, pos = parse_records!(Zones(), bytes, pos, len, options)
    area_transfers, pos = parse_records!(InterAreaTransfers(), bytes, pos, len, options)
    owners, pos = parse_records!(Owners(), bytes, pos, len, options)
    facts, pos = parse_records!(FACTSDevices30(), bytes, pos, len, options)
    return Network(
        version,
        caseid,
        buses,
        loads,
        fixed_shunts,
        gens,
        branches,
        transformers,
        interchanges,
        two_terminal_dc,
        vsc_dc,
        switched_shunts,
        impedance_corrections,
        multi_terminal_dc,
        multi_section_lines,
        zones,
        area_transfers,
        owners,
        facts,
    )
end

function parse_records!(rec::R, bytes, pos, len, options)::Tuple{R, Int} where {R <: Records}
    # Records terminated by specifying a bus number of zero or `Q`.
    while !(
        eof(bytes, pos, len) ||
        peekbyte(bytes, pos) == UInt8('0') ||
        peekbyte(bytes, pos) == UInt8(' ') && !eof(bytes, pos+1, len) && peekbyte(bytes, pos+1) == UInt8('0') ||
        peekbyte(bytes, pos) == UInt8('Q')
    )
        _, pos = parse_row!(rec, bytes, pos, len, options)
    end
    pos = next_line(bytes, pos, len)  # Move past a "0 bus" line.
    @debug 1 "Parsed $R: nrows = $(length(rec)), pos = $pos"
    return rec, pos
end

# Taken from `Parsers.checkcmtemptylines`
# TODO: move to Parsers.jl?
function next_line(bytes, pos, len)
    eof(bytes, pos, len) && return pos
    b = peekbyte(bytes, pos)
    while b !== UInt8('\n') && b !== UInt8('\r')
        pos += 1
        eof(bytes, pos, len) && break
        b = peekbyte(bytes, pos)
    end
    # Move forward to be past the `\r` or `\n` byte.
    pos += 1
    # if line ends `\r\n`, then we're at `\n`and need to move forward again.
    if b === UInt8('\r') && !eof(bytes, pos, len) && peekbyte(bytes, pos) === UInt8('\n')
        pos += 1
    end
    return pos
end

function parse_value(::Type{T}, bytes, pos, len, options) where {T}
    res = xparse(T, bytes, pos, len, options)
    code = res.code
    if invalid(code)
        if !(newline(code) && invaliddelimiter(code))  # not due to end-of-line comments
            @warn codes(res.code) pos
        end
    end
    pos += res.tlen
    return res.val, pos, res.code
end

function parse_value!(rec, col::Int, ::Type{T}, bytes, pos, len, options) where {T}
    val, pos, code = parse_value(nonmissingtype(T), bytes, pos, len, options)
    push!(getfield(rec, col)::Vector{T}, val)
    return rec, pos, code
end

@generated function parse_row!(rec::R, bytes, pos, len, options) where {R <: Records}
    block = Expr(:block)
    for col in 1:fieldcount(R)
        T = eltype(fieldtype(R, col))
        push!(block.args, quote
            rec, pos, code = parse_value!(rec, $col, $T, bytes, pos, len, options)
        end)
    end
    # @show block
    return block
end

###
### transformers
###

function _setmissing(a::Int, b::Int)
    exprs = Expr[]
    for col in a:b
        push!(exprs, :(push!(getfield(rec, $col), missing)))
    end
    return exprs
end

function _parse_values(::Type{R}, a::Int, b::Int) where {R <: Records}
    exprs = Expr[]
    for col in a:b
        T = eltype(fieldtype(R, col))
        push!(exprs, :((rec, pos, code) = parse_value!(rec, $col, $T, bytes, pos, len, options)))
    end
    return exprs
end

function _parse_maybemissing(R, col)
    T = eltype(fieldtype(R, col))
    return quote
        if newline(code)
            push!(getfield(rec, $col), missing)
        else
            (rec, pos, code) = parse_value!(rec, $col, $T, bytes, pos, len, options)
        end
    end
end

function _parse_maybemissing(R, col1, col2)
    T1 = eltype(fieldtype(R, col1))
    T2 = eltype(fieldtype(R, col2))
    return quote
        if newline(code)
            push!(getfield(rec, $col1), missing)
            push!(getfield(rec, $col2), missing)
        else
            (rec, pos, code) = parse_value!(rec, $col1, $T1, bytes, pos, len, options)
            (rec, pos, code) = parse_value!(rec, $col2, $T2, bytes, pos, len, options)
        end
    end
end

function _parse_maybezero(R, col1, col2)
    T1 = eltype(fieldtype(R, col1))
    T2 = eltype(fieldtype(R, col2))
    return quote
        if newline(code)
            push!(getfield(rec, $col1), zero($T1))
            push!(getfield(rec, $col2), zero($T2))
        else
            (rec, pos, code) = parse_value!(rec, $col1, $T1, bytes, pos, len, options)
            (rec, pos, code) = parse_value!(rec, $col2, $T2, bytes, pos, len, options)
        end
    end
end

function _parse_t2()
    block = Expr(:block)
    append!(block.args, _setmissing(EOL_COLS[1]+1+T2_COLS[2], EOL_COLS[2]))
    append!(block.args, _parse_values(Transformers, EOL_COLS[2]+1, EOL_COLS[3]-1))
    push!(block.args, _parse_maybemissing(Transformers, EOL_COLS[3]))
    append!(block.args, _parse_values(Transformers, EOL_COLS[3]+1, EOL_COLS[3]+T2_COLS[4]))
    append!(block.args, _setmissing(EOL_COLS[3]+1+T2_COLS[4], EOL_COLS[5]))
    return block
end

function _parse_t3()
    block = Expr(:block)
    append!(block.args, _parse_values(Transformers, EOL_COLS[1]+1+T2_COLS[2], EOL_COLS[3]-1))
    push!(block.args, _parse_maybemissing(Transformers, EOL_COLS[3]))
    append!(block.args, _parse_values(Transformers, EOL_COLS[3]+1, EOL_COLS[4]-1))
    push!(block.args, _parse_maybemissing(Transformers, EOL_COLS[4]))
    append!(block.args, _parse_values(Transformers, EOL_COLS[4]+1, EOL_COLS[5]-1))
    push!(block.args, _parse_maybemissing(Transformers, EOL_COLS[5]))
    return block
end

# 2-winding transformers (T2) have a subset of the data for 3-winding transformers (T3).
# T2 has data over 4 lines, T3 has data for 5 lines:
# - Line 1 is the same for both
# - Line 2 has only 3 entries for T2 data, and more for T3
# - Line 3 is the same for both
# - Line 4 has only 2 entries for T2 data, and more for T3 (similar to line 2)
# - Line 5 only exists for T3 data
# We determine data is T2 if there is a newline after 3 entries of line 2, else it's T3.
# This means T2 data with a comment after the last entry on line 2 will fool us.
@generated function parse_row!(rec::R, bytes, pos, len, options) where {R <: Transformers}
    block = Expr(:block)
    # parse row 1
    append!(block.args, _parse_values(R, 1, EOL_COLS[1]-7))
    # parse possibly missing o2,f2,o3,f3,o4,f4
    push!(block.args, _parse_maybemissing(Transformers, EOL_COLS[1]-6, EOL_COLS[1]-5))
    push!(block.args, _parse_maybemissing(Transformers, EOL_COLS[1]-4, EOL_COLS[1]-3))
    push!(block.args, _parse_maybemissing(Transformers, EOL_COLS[1]-2, EOL_COLS[1]-1))
    push!(block.args, _parse_maybemissing(R, EOL_COLS[1]))  # last col only in v33 data (not v30)
    # parse first part of row 2
    append!(block.args, _parse_values(R, EOL_COLS[1]+1, EOL_COLS[1]+T2_COLS[2]))
    # now we can detect if it is two-winding or three-winding data
    push!(block.args, :(newline(code) ? $(_parse_t2()) : $(_parse_t3())))
    push!(block.args, :(return rec, pos))
    # @show block
    return block
end

###
### SwitchedShunts, ImpedanceCorrections
###

const N_SPECIAL = IdDict(
    # SwitchedShunts can have anywhere between 1 - 8 `N` and `B` values in the data itself,
    # if n2, b2, ..., n8, b8 are not present, we set them to zero.
    # i.e. the last 14 = 7(n) + 7(b) columns reqire special handling.
    SwitchedShunts30 => 14,
    SwitchedShunts33 => 14,
    # ImpedanceCorrections can have anywhere between 2 - 11 `T` and `F` values in the data itself,
    # if t3, f3, ..., t11, f11 are not present, we set them to zero.
    # i.e. the last 18 = 9(t) + 9(f) columns reqire special handling.
    ImpedanceCorrections => 18,
    # MultiSectionLineGroups can have between 1 - 9 `DUM_i` columns
    MultiSectionLineGroups30 => 8,
    MultiSectionLineGroups33 => 8,
    # Loads have 2 extra columns in v33 compared to v30
    Loads => 2,
    # Generators have 1 - 4 `Oi`, `Fi` values, plus 2 extra columns in v33 compared to v30
    Generators => 8, # 3*2 + 2
    # Branches have 1 - 4 `Oi`, `Fi` values
    Branches30 => 6, # 3*2
    Branches33 => 6, # 3*2
)

@generated function parse_row!(rec::R, bytes, pos, len, options) where {R <: Union{SwitchedShunts, ImpedanceCorrections, Branches}}
    block = Expr(:block)
    N = fieldcount(R) - N_SPECIAL[R]
    append!(block.args, _parse_values(R, 1, N))
    coln = N + 1
    colb = N + 2
    for _ in 1:(N_SPECIAL[R] ÷ 2)
        push!(block.args, _parse_maybezero(R, coln, colb))
        coln += 2
        colb += 2
    end
    push!(block.args, :(return rec, pos))
    # @show block
    return block
end

###
### Loads, Generators, MultiSectionLineGroups
###

@generated function parse_row!(
    rec::R, bytes, pos, len, options
) where {R <: Union{Loads,Generators,MultiSectionLineGroups}}
    block = Expr(:block)
    N = fieldcount(R) - N_SPECIAL[R]
    append!(block.args, _parse_values(R, 1, N))
    for col in (N + 1):fieldcount(R)
        push!(block.args, _parse_maybemissing(R, col))
    end
    push!(block.args, :(return rec, pos))
    # @show block
    return block
end

###
### MultiTerminalDCLines
###

function parse_row!(rec::R, bytes, pos, len, options) where {I, R <: MultiTerminalDCLines{I}}
    line_id, pos = parse_idrow(I, bytes, pos, len, options)

    nconv = line_id.nconv
    converters = ACConverters(nconv)
    for _ in 1:nconv
        converters, pos = parse_row!(converters, bytes, pos, len, options)
    end

    ndcbs = line_id.ndcbs
    dc_buses = DCBuses(ndcbs)
    for _ in 1:ndcbs
        dc_buses, pos = parse_row!(dc_buses, bytes, pos, len, options)
    end

    ndcln = line_id.ndcln
    dc_links = DCLinks(ndcln)
    for _ in 1:ndcln
        dc_links, pos = parse_row!(dc_links, bytes, pos, len, options)
    end
    line = MultiTerminalDCLine(line_id, converters, dc_buses, dc_links)
    push!(rec.lines, line)
    return rec, pos
end

###
### IDRow
###

function parse_value!(args, ::Type{T}, bytes, pos, len, options) where {T}
    res = xparse(T, bytes, pos, len, options)
    code = res.code
    val = ifelse(valueok(code), res.val, missing)
    push!(args, val)
    pos += res.tlen
    return pos, code
end

@generated function parse_idrow(::Type{R}, bytes, pos, len, options) where {R <: IDRow}
    block = Expr(:block)
    nfields = fieldcount(R)
    T = nonmissingtype(fieldtype(R, 1))
    push!(block.args, quote
        args = Any[]
        (pos, code) = parse_value!(args, $T, bytes, pos, len, options)
    end)
    for i in 2:nfields
        T = nonmissingtype(fieldtype(R, i))
        push!(block.args, quote
            if invalid(code) || newline(code)
                push!(args, missing)
            else
                (pos, code) = parse_value!(args, $T, bytes, pos, len, options)
            end
        end)
    end
    push!(block.args, quote
        if !newline(code)
            pos = next_line(bytes, pos, len)
        end
        return R(args...), pos
    end)
    # @show block
    return block
end


text:
module TestSamples

using FGenerators
using FLoops
using SplittablesBase
using SplittablesTesting
using Test
using Transducers: Map

@fgenerator noone() = nothing
@fgenerator oneone() = @yield 1
@fgenerator function onetwothree()
    @yield 1
    @yield 2
    @yield 3
end

@fgenerator function organpipe(n::Integer)
    i = 0
    while i != n
        i += 1
        @yield i
    end
    while true
        i -= 1
        i == 0 && return
        @yield i
    end
end

struct OrganPipe <: Foldable
    n::Int
end

@fgenerator(foldable::OrganPipe) do
    n = foldable.n
    @yieldfrom 1:n
    @yieldfrom n-1:-1:1
end

function SplittablesBase.halve(foldable::OrganPipe)
    n = foldable.n
    return (1:n, n-1:-1:1)
end

Base.length(foldable::OrganPipe) = 2 * foldable.n - 1

@fgenerator function flatten2(a, b)
    @yieldfrom a
    @yieldfrom b
end

struct Count <: Foldable
    start::Int
    stop::Int
end

Base.length(itr::Count) = max(0, itr.stop - itr.start + 1)
Base.eltype(::Type{<:Count}) = Int

@fgenerator(itr::Count) do
    i = itr.start
    i > itr.stop && return
    while true
        @yield i
        i == itr.stop && break
        i += 1
    end
end

repeat3 = @fgenerator function (x)
    @yield x
    @yield x
    @yield x
end

repeat3_arrow = @fgenerator x -> begin
    @yield x
    @yield x
    @yield x
end

repeat3_arrow2 = @fgenerator (x,) -> begin
    @yield x
    @yield x
    @yield x
end

@fgenerator function ffilter(f, xs)
    @floop for x in xs
        if f(x)
            @yield x
        end
    end
end

asval(x::Val) = x
asval(x) = Val(x)
const FlagType = Union{Val{true},Val{false},Bool}

@fgenerator function linesin(str::AbstractString; keep::FlagType = Val(false))
    keep = asval(keep)
    start = firstindex(str)
    for (i, c) in pairs(str)
        if c == '\n'
            if keep === Val(true)
                @yield SubString(str, start, i)
            else
                @yield SubString(str, start, prevind(str, i))
            end
            start = nextind(str, i)
        end
    end
    if start <= ncodeunits(str)
        @yield SubString(str, start)
    end
end

raw_testdata = raw"""
noone() == []
oneone() == [1]
onetwothree() == [1, 2, 3]
organpipe(2) == [1, 2, 1]
organpipe(3) == [1, 2, 3, 2, 1]
OrganPipe(2) == [1, 2, 1]
OrganPipe(3) == [1, 2, 3, 2, 1]
flatten2((1, 2), (3,)) == [1, 2, 3]
Count(0, 2) == [0, 1, 2]
Count(0, -1) == Int[]
repeat3(:a) == [:a, :a, :a]
repeat3_arrow(:a) == [:a, :a, :a]
repeat3_arrow2(:a) == [:a, :a, :a]
ffilter(isodd, noone()) == []
ffilter(isodd, oneone()) == [1]
ffilter(isodd, onetwothree()) == [1, 3]
ffilter(isodd, 1:5) == [1, 3, 5]
ffilter(isodd, organpipe(3)) == [1, 3, 1]
linesin("a\nbb\nccc") == ["a", "bb", "ccc"]
linesin("a\nbb\nccc"; keep=true) == ["a\n", "bb\n", "ccc"]
linesin("a\nbb\nccc\n") == ["a", "bb", "ccc"]
linesin("a\nbb\nccc\n"; keep=true) == ["a\n", "bb\n", "ccc\n"]
"""

args_and_kwargs(args...; kwargs...) = args, (; kwargs...)

# An array of `(label, (f, args, kwargs, comparison, desired))`
testdata = map(split(raw_testdata, "\n", keepempty = false)) do x
    comp_ex = Meta.parse(x)
    @assert comp_ex.head == :call
    @assert length(comp_ex.args) == 3
    comparison, ex, desired = comp_ex.args
    f = ex.args[1]
    ex.args[1] = args_and_kwargs

    label = strip(x[1:prevind(x, first(findlast(String(comparison), x)))])
    Meta.parse(label)  # validation

    @eval ($label, ($(Symbol(f)), $ex..., $comparison, $desired))
end

@testset "$label" for (label, (f, args, kwargs, comparison, desired)) in testdata
    ==′ = comparison
    @test collect(f(args...; kwargs...)) ==′ desired
    @test collect(Map(identity), f(args...; kwargs...)) ==′ desired
end

@testset "inference" begin
    @test @inferred(sum(ffilter(isodd, onetwothree()))) == 4
    @test @inferred(sum(ffilter(isodd, organpipe(3)))) == 5
    @test @inferred(sum(ffilter(isodd, OrganPipe(3)))) == 5
end

SplittablesTesting.test_ordered(Any[OrganPipe(n) for n in 1:10])

end  # module


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _xyz(jd::T, equinox::T) where {T<:AbstractFloat}
    t = (jd - 15020) / JULIANCENTURY # Reduced Julian century from 1900

    # NOTE: longitude arguments below are given in *equinox* of date.  Precess
    # these to equinox 1950 to give everything an even footing.  Compute
    # argument of precession from equinox of date back to 1950
    pp = (1.396041 + 0.000308*(t + 0.5))*(t-0.499998)
    # Compute mean solar longitude, precessed back to 1950
    el = @evalpoly t (279.696678 - pp) 36000.76892 0.000303

    # Compute Mean longitude of the Moon
    c = @evalpoly t (270.434164 - pp) (480960.0 + 307.883142) -0.001133

    # Compute longitude of Moon's ascending node
    n = @evalpoly t (259.183275 - pp) -(1800.0 + 134.142008) 0.002078

    # Compute mean solar anomaly
    g = @evalpoly t 358.475833 35999.04975 -0.00015

    # Compute the mean jupiter anomaly
    j = @evalpoly t 225.444651 2880.0 154.906654

    # Compute mean anomaly of Venus
    v = @evalpoly t 212.603219 (58320.0 + 197.803875) 0.001286

    # Compute mean anomaly of Mars
    m = @evalpoly t 319.529425 (19080.0 + 59.8585) 0.000181

    # Convert degrees to radians for trig functions
    el = deg2rad(el)
    g  = deg2rad(g)
    j  = deg2rad(j)
    c  = deg2rad(c)
    v  = deg2rad(v)
    n  = deg2rad(n)
    m  = deg2rad(m)

    # Calculate x, y, z using trigonometric series
    sin1,  cos1  = sincos(el)
    sin2,  cos2  = sincos(g - el)
    sin3,  cos3  = sincos(g + el)
    sin4,  cos4  = sincos(g + g + el)
    sin5,  cos5  = sincos(g + g - el)
    sin6,  cos6  = sincos(g - el - j)
    sin7,  cos7  = sincos(2 * g + el - 2 * v)
    sin8,  cos8  = sincos(c)
    sin9,  cos9  = sincos(c - 2 * el)
    sin10, cos10 = sincos(4 * g + el - 8 * m + 3 * j)
    sin11, cos11 = sincos(4 * g - el - 8 * m + 3 * j)
    sin12, cos12 = sincos(g + el - v)
    sin13, cos13 = sincos(2 * g - el - 2 * j)

    x = 0.999860 * cos1     -
        0.025127 * cos2     +
        0.008374 * cos3     +
        0.000105 * cos4     +
        0.000063 * t * cos2 +
        0.000035 * cos5     -
        0.000026 * sin6     -
        0.000021 * t * cos3 +
        0.000018 * sin7     +
        0.000017 * cos8     -
        0.000014 * cos9     +
        0.000012 * cos10    -
        0.000012 * cos11    -
        0.000012 * cos12    +
        0.000011 * cos7     +
        0.000011 * cos13

    y = 0.917308 * sin1     +
        0.023053 * sin2     +
        0.007683 * sin3     +
        0.000097 * sin4     -
        0.000057 * t * sin2 -
        0.000032 * sin5     -
        0.000024 * cos6     -
        0.000019 * t * sin3 -
        0.000017 * cos7     +
        0.000016 * sin8     +
        0.000013 * sin9     +
        0.000011 * sin10    +
        0.000011 * sin11    -
        0.000011 * sin12    +
        0.000010 * sin7     -
        0.000010 * sin13

    z = 0.397825 * sin1     +
        0.009998 * sin2     +
        0.003332 * sin3     +
        0.000042 * sin4     -
        0.000025 * t * sin2 -
        0.000014 * sin5     -
        0.000010 * cos6

    # Precess to new equator?  Avoid useless calculations.
    if isfinite(equinox) && equinox != 1950
        x, y, z = precess_xyz(x, y, z, 1950.0, equinox)
    end

    # Velocities
    xvel = -0.017200 * sin(el)   -
        0.000288*sin(g + el)     -
        0.000005*sin(2.0*g + el) -
        0.000004*sin(c)          +
        0.000003*sin(c - 2.0*el) +
        0.000001*t*sin(g + el)   -
        0.000001*sin(2.0*g - el)

    yvel = 0.015780*cos(el)       +
        0.000264*cos(g + el)      +
        0.000005*cos(2.0*g + el)  +
        0.000004*cos(c)           +
        0.000003*cos(c - 2.0*el)  -
        0.000001*t*cos(g + el)

    zvel = 0.006843*cos(el)      +
        0.000115*cos(g  + el)    +
        0.000002*cos(2.0*g + el) +
        0.000002*cos(c)          +
        0.000001*cos(c - 2.0*el)

    # Precess to new equator?  Avoid useless calculations.
    if isfinite(equinox) && equinox != 1950
        xvel, yvel, zvel = precess_xyz(xvel, yvel, zvel, 1950, equinox)
    end

    return x, y, z, xvel, yvel, zvel
end

"""
    xyz(jd[, equinox]) -> x, y, z, v_x, v_y, v_z

### Purpose ###

Calculate geocentric \$x\$, \$y\$, and \$z\$ and velocity coordinates of the
Sun.

### Explanation ###

Calculates geocentric \$x\$, \$y\$, and \$z\$ vectors and velocity coordinates
(\$dx\$, \$dy\$ and \$dz\$) of the Sun.  (The positive \$x\$ axis is directed
towards the equinox, the \$y\$-axis, towards the point on the equator at right
ascension 6h, and the \$z\$ axis toward the north pole of the equator).  Typical
position accuracy is \$<10^{-4}\$ AU (15000 km).

### Arguments ###

* `jd`: number of Reduced Julian Days for the wanted date.  It can be either a
  scalar or a vector.
* `equinox` (optional numeric argument): equinox of output. Default is 1950.

You can use `juldate` to get the number of Reduced Julian Days for the selected
dates.

### Output ###

The 6-tuple ``(x, y, z, v_x, v_y, v_z)``, where

* ``x, y, z``: scalars or vectors giving heliocentric rectangular
  coordinates (in AU) for each date supplied.  Note that ``\\sqrt{x^2 + y^2 + z^2}``
  gives the Earth-Sun distance for the given date.
* ``v_x, v_y, v_z``: velocity vectors corresponding to ``x, y``, and
  ``z``.

### Example ###

What were the rectangular coordinates and velocities of the Sun on
1999-01-22T00:00:00 (= JD 2451200.5) in J2000 coords?  Note: Astronomical
Almanac (AA) is in TDT, so add 64 seconds to UT to convert.

```jldoctest
julia> using AstroLib, Dates

julia> jd = juldate(DateTime(1999, 1, 22))
51200.5

julia> xyz(jd + 64/86400, 2000)
(0.514568709240398, -0.7696326261820209, -0.33376880143023935, 0.014947267514079971, 0.008314838205477328, 0.003606857607575486)
```

Compare to Astronomical Almanac (1999 page C20)

                x  (AU)        y  (AU)     z (AU)
    xyz:      0.51456871   -0.76963263  -0.33376880
    AA:       0.51453130   -0.7697110   -0.3337152
    abs(err): 0.00003739    0.00007839   0.00005360
    abs(err)
        (km):   5609          11759         8040

NOTE: Velocities in AA are for Earth/Moon barycenter
      (a very minor offset) see AA 1999 page E3

               x vel (AU/day) y vel (AU/day)   z vel (AU/day)
    xyz:      -0.014947268   -0.0083148382    -0.0036068576
    AA:       -0.01494574    -0.00831185      -0.00360365
    abs(err):  0.000001583    0.0000029886     0.0000032076
    abs(err)
     (km/sec): 0.00265        0.00519          0.00557

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
xyz(jd::Real, equinox::Real=NaN) = _xyz(promote(float(jd), float(equinox))...)

# Can't use @vectorize_1arg because of the optional keyword.
function xyz(jd::AbstractArray{J}, equinox::Real=NaN) where {J<:Real}
    typej = float(J)
    x     = similar(jd, typej)
    y     = similar(jd, typej)
    z     = similar(jd, typej)
    xvel  = similar(jd, typej)
    yvel  = similar(jd, typej)
    zvel  = similar(jd, typej)
    for i in eachindex(jd)
        x[i], y[i], z[i], xvel[i], yvel[i], zvel[i] = xyz(jd[i], equinox)
    end
    return x, y, z, xvel, yvel, zvel
end


text:
mutable struct MyMutableStruct end
struct MyStruct
    a::Int
    b::Ptr{Cvoid}
    c::NTuple{16,Cfloat}
end

function f(x::Int, y, z; p=5)
    println(50)
    println("finished")
end

g(x) = 50 * x

@enum my_enum a b c d e f

@enum my_other_enum begin
    a = 3
    b = 2
    c = 200
end

const myconst = 50