text:
# octSock5
# Copyright (c) 2017, Markus Nentwig
# All rights reserved.
module octsock5

export octsock5_cl, octsock5_new, octsock5_accept, octsock5_write, octsock5_read, octsock5_delete;

__precompile__()

# profiling note: global "const" refers to the type, not the value.
# A macro provides the compiler with a literal instead of a variable.
macro H0_FLOAT() 	return Int64(0x00000001); end
macro H0_INTEGER() 	return Int64(0x00000002); end
macro H0_SIGNED() 	return Int64(0x00000004); end
macro H0_COMPLEX() 	return Int64(0x00000008); end
macro H0_ARRAY() 	return Int64(0x00000010); end
macro H0_TUPLE() 	return Int64(0x00000020); end
macro H0_STRING() 	return Int64(0x00000040); end
macro H0_DICT() 	return Int64(0x00000080); end
macro H0_TERM() 	return Int64(0x00000100); end
macro HEADERTYPE() 	return Int64; end
macro HEADERLEN()	return 6; end
macro HEADERSIZE()	return 48; end

# unique number that identifies each scalar / array / complex / number type combination
function TYPEID(H0::Int64, nBytes::Int64)
    return (H0 + (nBytes << 16))::Int64;
end

type octsock5_cl
    header::Array{@HEADERTYPE};
    stringMem::Array{UInt8};
    stringMemSize::UInt64;
    stringMemPtr::Ptr{UInt8};
    server::Any;
    io::IO;
    octsock5_cl() = (x = new());
end

function octsock5_new(;isServer::Bool=false, portNum::Int64=-1)
    self::octsock5_cl = octsock5_cl();
    self.header = Array{@HEADERTYPE}(@HEADERLEN);
    # pre-allocate initial memory for inbound strings
    self.stringMemSize = 1000;
    self.stringMem = Array{UInt8}(self.stringMemSize);
    self.stringMemPtr = pointer(self.stringMem);

    if (portNum < 0)
        # === Windows named pipes ===
        # default implementation. Use this.
        if (false == is_windows())
            error("named pipes are windows-OS specific, others must use TCP/IP ports");
        end
        conn::Any = string("\\\\.\\pipe\\octsock5_", @sprintf("%i", portNum));
    else
        # === TCP/IP ===
        # TBD support connection to another host
        # on localhost, expect 2x latency, 1/3 throughput relative to windows named pipes.
        conn = portNum;
    end

    if isServer
        self.server = listen(conn);
    else
        self.server = Void;
        self.io = connect(conn); # TBD: Could connect to another host here.
    end
    # auto-close if object is abandoned or on shutdown
    finalizer(self, octsock5_delete);

    return self;
end

# octsock5_accept is split from octsock5_new, because the client can only start when the server is up
# but accept is a blocking function
# server: octsock5_new, acknowledge
# client: octsock5_new
# server: octsock5_accept
function octsock5_accept(self::octsock5_cl)
    if (self.server == Void)
        error("not in server mode");
    end
    self.io = accept(self.server);
end

function octsock5_delete(self::octsock5_cl)
    try
        close(self.io);
        if self.server!= Void
            close(self.server);
        end
    catch end
    nothing;
end

function writeHeader(self::octsock5_cl)
    unsafe_write(self.io, Ref(self.header), @HEADERSIZE);
    nothing;
end

function writePointer(self::octsock5_cl, ptr::Ptr{Void}, len::UInt64)
    unsafe_write(self.io, ptr, len);
    nothing;
end

function writeScalar{T}(self::octsock5_cl, arg::T)
    unsafe_write(self.io, Ref{T}(arg), sizeof(T));
    nothing;
end

function octsock5_write(self::octsock5_cl, arg::Dict)
    @inbounds self.header[1] = @H0_DICT;
    writeHeader(self);

    for item::Pair in arg
        # === write key ===
        octsock5_write(self, item[1]);
        # === write value ===
        octsock5_write(self, item[2]);
    end

    @inbounds self.header[1] = @H0_TERM;
    writeHeader(self);
    nothing;
end

function octsock5_write(self::octsock5_cl, arg::Tuple)
    tupLen::UInt64 = length(arg);

    @inbounds self.header[1] = @H0_TUPLE;
    @inbounds self.header[2] = tupLen;
    writeHeader(self);

    for ix::UInt64 = 1:tupLen
        octsock5_write(self, arg[ix]);
    end
    nothing;
end

function octsock5_write(self::octsock5_cl, arg::String)
    len::UInt64 = length(arg);
    @inbounds self.header[1] = @H0_STRING;
    @inbounds self.header[2] = len;
    writeHeader(self);
    writePointer(self, Ptr{Void}(pointer(arg)), len);
    nothing;
end

function octsock5_write(self::octsock5_cl, arg::T) where T <: Union{RowVector, UnitRange, LinSpace, StepRangeLen}
    octsock5_write(self, Array(arg));
    nothing;
end

function octsock5_write{T}(self::octsock5_cl, arg::Array{T})
    H0::Int64 = (T <: Complex) 		? (@H0_ARRAY) | (@H0_COMPLEX) 	: (@H0_ARRAY);
    H0 |= (real(T) <: AbstractFloat) 	? (@H0_FLOAT) 			: 0;
    H0 |= (real(T) <: Integer) 		? (@H0_INTEGER) 		: 0;
    H0 |= (real(T) <: Signed) 		? (@H0_SIGNED) 			: 0;

    nd::Int64 = ndims(arg);
    if (nd > 3) error("unsupported number of matrix dimensions"); end
    nBytes::UInt64 = sizeof(T);
    @inbounds self.header[1] = H0;
    @inbounds self.header[2] = nBytes;
    @inbounds self.header[3] = nd;

    for dim = 1 : nd
        @inbounds self.header[3+dim] = size(arg, dim);
    end

    writeHeader(self);
    writePointer(self, Ptr{Void}(pointer(arg)), UInt64(sizeof(arg)));
    nothing;
end

function octsock5_write{T}(self::octsock5_cl, arg::T)
    H0::Int64 = (T <: Complex) 		? (@H0_COMPLEX) : 0;
    H0 |= (real(T) <: AbstractFloat)	? (@H0_FLOAT) 	: 0;
    H0 |= (real(T) <: Integer) 		? (@H0_INTEGER) : 0;
    H0 |= (real(T) <: Signed) 		? (@H0_SIGNED) 	: 0;
    nBytes::UInt64 = sizeof(T);
    @inbounds self.header[1] = H0;
    @inbounds self.header[2] = nBytes;
    @inbounds self.header[3] = 1;

    writeHeader(self);
    writeScalar(self, arg);
    nothing;
end

function readScalar{T}(self::octsock5_cl, dummy::T)
    r::Ref{T} = Ref{T}(dummy);
    p::Ptr{T} = Base.unsafe_convert(Ptr{T}, r);

    unsafe_read(self.io, p, sizeof(dummy));
    return r.x::T;
end

function readArray{T}(self::octsock5_cl, dummy::T)
    @inbounds H3::Int64 = self.header[3];
    @inbounds H4::Int64 = self.header[4];
    @inbounds H5::Int64 = self.header[5];
    @inbounds H6::Int64 = self.header[6];
    if (H3 == 1)
        obj::Array{T} = Array{T, 1}(H4);
        nElem::UInt64 = H4;
        p1::Ptr{Void} = pointer(obj);
    elseif (H3 == 2)
        obj = Array{T, 2}(H4, H5);
        nElem = H4*H5;
        p1 = pointer(obj);
    elseif (H3 == 3)
        obj = Array{T, 3}(H4, H5, H6);
        nElem = H4*H5*H6;
        p1 = pointer(obj);
    else
        error("invalid number of dimensions");
    end

    nBytes::UInt64 = sizeof(T) * nElem;

    unsafe_read(self.io, p1, nBytes);
    return obj;
end

function octsock5_read(self::octsock5_cl)

    # === read header ===
    unsafe_read(self.io, Ref(self.header), @HEADERSIZE);
    #print("header 1: 0x", hex(self.header[1], 8), "\n");
    #print("header 2: 0x", hex(self.header[2], 8), "\n");
    #print("header 3: 0x", hex(self.header[3], 8), "\n");
    #print("header 4: 0x", hex(self.header[4], 8), "\n");
    #print("header 5: 0x", hex(self.header[5], 8), "\n");
    @inbounds H0::Int64 = self.header[1];

    # === handle terminator (dict) ===
    if ((H0 & @H0_TERM) != 0)
        return Void;
    end

    # === handle tuples ===
    if ((H0 & @H0_TUPLE) != 0)
        @inbounds local nElemTuple::UInt64 = self.header[2];
        local objTuple::Array{Any} = Array{Any}(nElemTuple);
        for ix::UInt64 = 1:nElemTuple
            @inbounds objTuple[ix] = octsock5_read(self);
        end
        return Tuple(objTuple);
    end

    # === handle string ===
    if ((H0 & @H0_STRING) != 0)
        @inbounds nBytesStr::Int64 = self.header[2];
        # allocate memory
        if self.stringMemSize < nBytesStr
            self.stringMemSize = nBytesStr;
            self.stringMem = Array{UInt8}(self.stringMemSize);
            self.stringMemPtr = pointer(self.stringMem);
        end
        unsafe_read(self.io, self.stringMemPtr, nBytesStr);
        return unsafe_string(self.stringMemPtr, nBytesStr);
    end

    # === handle Dict ===
    if ((H0 & @H0_DICT) != 0)
        local objDict::Dict = Dict();
        while true
            key = octsock5_read(self);
            if (key == Void)
                return objDict;
            end
            objDict[key] = octsock5_read(self);
        end
    end

    @inbounds nElemBytes::Int64 = self.header[2];
    isComplex::Bool = ((H0 & @H0_COMPLEX) != 0);
    isSigned::Bool = ((H0 & @H0_SIGNED) != 0);
    @inbounds tId::UInt64 = TYPEID(H0, self.header[2]);
    if ((H0 & @H0_ARRAY) == 0)
        # === handle array ===
        if ((H0 & @H0_COMPLEX) != 0)
            if (nElemBytes == 2)
                if 	(tId == TYPEID((@H0_INTEGER) | (@H0_COMPLEX), 2)) return readScalar(self, Complex{UInt8}(0));
                elseif	(tId == TYPEID((@H0_INTEGER) | (@H0_COMPLEX) | (@H0_SIGNED), 2)) return readScalar(self, Complex{Int8}(0));
                else error(); end
            elseif (nElemBytes == 4)
                if 	(tId == TYPEID((@H0_INTEGER) | (@H0_COMPLEX), 4)) return readScalar(self, Complex{UInt16}(0));
                elseif 	(tId == TYPEID((@H0_INTEGER) | (@H0_COMPLEX) | (@H0_SIGNED), 4)) return readScalar(self, Complex{Int16}(0));
                else error(); end
            elseif (nElemBytes == 8)
                if 	(tId == TYPEID((@H0_FLOAT)   | (@H0_COMPLEX), 8)) return readScalar(self, Complex{Float32}(0));
                elseif 	(tId == TYPEID((@H0_INTEGER) | (@H0_COMPLEX), 8)) return readScalar(self, Complex{UInt32}(0));
                elseif 	(tId == TYPEID((@H0_INTEGER) | (@H0_COMPLEX) | (@H0_SIGNED), 8)) return readScalar(self, Complex{Int32}(0));
                else error(); end
            elseif (nElemBytes == 16)
                if 	(tId == TYPEID((@H0_COMPLEX) | (@H0_FLOAT), 16)) return readScalar(self, Complex{Float64}(0));
                elseif 	(tId == TYPEID((@H0_COMPLEX) | (@H0_INTEGER), 16)) return readScalar(self, Complex{UInt64}(0));
                elseif 	(tId == TYPEID((@H0_COMPLEX) | (@H0_INTEGER) | (@H0_SIGNED), 16)) return readScalar(self, Complex{Int64}(0));
                else error() end
            else error(); end
        else
            if (nElemBytes == 1)
                if 	(tId == TYPEID((@H0_INTEGER), 1)) return readScalar(self, UInt8(0));
                elseif	(tId == TYPEID((@H0_INTEGER) | (@H0_SIGNED), 1)) return readScalar(self, Int8(0));
                else error(); end
            elseif (nElemBytes == 2)
                if 	(tId == TYPEID((@H0_INTEGER), 2)) return readScalar(self, UInt16(0));
                elseif 	(tId == TYPEID((@H0_INTEGER) | (@H0_SIGNED), 2)) return readScalar(self, Int16(0));
                else error(); end
            elseif (nElemBytes == 4)
                if 	(tId == TYPEID((@H0_FLOAT), 4)) return readScalar(self, Float32(0));
                elseif 	(tId == TYPEID((@H0_INTEGER), 4)) return readScalar(self, UInt32(0));
                elseif 	(tId == TYPEID((@H0_INTEGER) | (@H0_SIGNED), 4)) return readScalar(self, Int32(0));
                else error(); end
            elseif (nElemBytes == 8)
                if 	(tId == TYPEID((@H0_FLOAT), 8)) return readScalar(self, Float64(0));
                elseif 	(tId == TYPEID((@H0_INTEGER), 8)) return readScalar(self, UInt64(0));
                elseif 	(tId == TYPEID((@H0_INTEGER) | (@H0_SIGNED), 8)) return readScalar(self, Int64(0));
                else error(); end
            else error(); end
        end
    else
        # === handle scalar ===
        if ((H0 & @H0_COMPLEX) != 0)
            if (nElemBytes == 2)
                if 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_COMPLEX), 2)) return readArray(self, Complex{UInt8}(0));
                elseif	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_COMPLEX) | (@H0_SIGNED), 2)) return readArray(self, Complex{Int8}(0));
                else error(H0);
                end
            elseif (nElemBytes == 4)
                if 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_COMPLEX), 4)) return readArray(self, Complex{UInt16}(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_COMPLEX) | (@H0_SIGNED), 4)) return readArray(self, Complex{Int16}(0));
                else error() end
            elseif (nElemBytes == 8)
                if 	(tId == TYPEID((@H0_ARRAY) | (@H0_FLOAT)   | (@H0_COMPLEX), 8)) return readArray(self, Complex{Float32}(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_COMPLEX), 8)) return readArray(self, Complex{UInt32}(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_COMPLEX) | (@H0_SIGNED), 8)) return readArray(self, Complex{Int32}(0));
                else error() end
            elseif (nElemBytes == 16)
                if 	(tId == TYPEID((@H0_ARRAY) | (@H0_COMPLEX) | (@H0_FLOAT), 16)) return readArray(self, Complex{Float64}(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_COMPLEX) | (@H0_INTEGER), 16)) return readArray(self, Complex{UInt64}(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_COMPLEX) | (@H0_INTEGER) | (@H0_SIGNED), 16)) return readArray(self, Complex{Int64}(0));
                else error() end
            else error(); end
        else
            if (nElemBytes == 1)
                if 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER), 1)) return readArray(self, UInt8(0));
                elseif	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_SIGNED), 1)) return readArray(self, Int8(0));
                else error();
                end
            elseif (nElemBytes == 2)
                if 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER), 2)) return readArray(self, UInt16(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_SIGNED), 2)) return readArray(self, Int16(0));
                else error() end
            elseif (nElemBytes == 4)
                if 	(tId == TYPEID((@H0_ARRAY) | (@H0_FLOAT), 4)) return readArray(self, Float32(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER), 4)) return readArray(self, UInt32(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_SIGNED), 4)) return readArray(self, Int32(0));
                else error() end
            elseif (nElemBytes == 8)
                if 	(tId == TYPEID((@H0_ARRAY) | (@H0_FLOAT), 8)) return readArray(self, Float64(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER), 8)) return readArray(self, UInt64(0));
                elseif 	(tId == TYPEID((@H0_ARRAY) | (@H0_INTEGER) | (@H0_SIGNED), 8)) return readArray(self, Int64(0));
                else error() end
            else error(); end
        end
end
end
end # module


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Result (ModiaMath/Result/_module.jl)
#

#---------------------------- Acessing the time series results (various options)
appendUnit(name, unit) = unit == "" ? string(name) : string(name, " [", unit, "]")

function getStringDictResult(model, res::RawResult)
    nt     = res.nt
    data   = res.data
    result = Dict{AbstractString,Any}()
    for (i, name) in enumerate(res.names)
        result[name] = @view data[1:nt,i]
    end

    return result, nt
end

const StringDictAnyResult = Dict{AbstractString,Any}
resultHeading(result::StringDictAnyResult) = ""

function getSignal(seriesDict::StringDictAnyResult, name::AbstractString)
    if haskey(seriesDict, name)
        sig     = seriesDict[name]
        keyName = name
    else
        if name[end] == ']'
            @static if VERSION >= v"0.7.0-DEV.2005"
                indexRange = something(findlast("[", name), 0:-1)
            else
                indexRange = rsearch(name, "[")
            end
            i = indexRange[1]
            @assert(i >= 2)
            keyName = name[1:i-1]

            if haskey(seriesDict, keyName)
                sig2 = seriesDict[keyName]

                # Build AST and evaluate: sig = sig2[:,...]
                ex1 = Meta.parse(name[i:end])
                ex2 = :( $sig2[:] )
                append!(ex2.args, ex1.args)
                sig = @eval $ex2
                if ndims(sig) != 1
                    @static if VERSION >= v"0.7.0-DEV.2005"
                        @warn "ModiaMath.plot: argument name (= $name) does not characterize one array element\nIndex ranges are not yet supported."
                    else
                        warn("ModiaMath.plot: argument name (= $name) does not characterize one array element\nIndex ranges are not yet supported.")
                    end
                    sig = nothing
                end

            else
                @static if VERSION >= v"0.7.0-DEV.2005"
                    @warn "ModiaMath.plot: argument name (= $name) is not correct."
                else
                    warn("\nModiaMath.plot: argument name (= ", name, ") is not correct.")
                end
                sig = nothing
            end

        else
            @static if VERSION >= v"0.7.0-DEV.2005"
                @warn "ModiaMath.plot: argument name (= $name) is not correct or does not identify a signal in the result."
            else
                warn("\nModiaMath.plot: argument name (= ", name, ") is not correct or does not identify a signal in the result.")
            end
            sig = nothing
            keyName = name
        end
    end
    return (sig, keyName, name)
end

"""
    (xsig, xsigLegend, ysig, ysigLegend) =
           ModiaMath.resultTimeSeries(result, name, xLabel::Bool, xAxis)

For a desired `result` data structure, this function has to be provided
to return the x-vector (`xsig`), the y-vector (`ysig`) and the legend
of the x-vector (`xsigLegend`), and of the y-vector (`ysigLegend`) as Strings, given
the key of the y-vector (`name`) and the key of the x-vector (`xAxis`).
If `xLabel=false` the legend of the x-vector should be an empty string (`""`).
"""
function resultTimeSeries(result::StringDictAnyResult, name, xLabel::Bool, xAxis)
    # Get x-axis signal
    (xsig, xkeyName, xName) = getSignal(result, string(xAxis))
    if xsig == nothing
        return (nothing, nothing, nothing, nothing)
    end
    if ndims(xsig) != 1
        @static if VERSION >= v"0.7.0-DEV.2005"
            @warn "ModiaMath.plot: argument xAxis (= $xAxis) does not characterize a signal vector."
        else
            warn("\nModiaMath.plot: argument xAxis (= $xAxis) does not characterize a signal vector.")
        end
        return (nothing, nothing, nothing, nothing)
    end
    xsigLegend = xLabel ? appendUnit(string(xAxis), string(unit(xsig[1]))) : ""
    xsig       = ustrip.(xsig)

    # Get y-axis signals
    (ysig, ykeyName, yName) = getSignal(result, string(name))

    if ysig == nothing
        return (nothing, nothing, nothing, nothing)
    end

    if ndims(ysig) == 1
        # ysig is a vector
        ysigLegend = [appendUnit(yName, string(unit(ysig[1])))]
    elseif ndims(ysig) == 2
        # ysig is a matrix
        nLegend = size(ysig,2)
        @static if VERSION >= v"0.7.0-DEV.2005"
            ysigLegend = Array{String}(undef, nLegend)
        else
            ysigLegend = Array{String}(nLegend)
        end

        for i = 1:nLegend
            ysigLegend[i] = appendUnit(ykeyName*"["*string(i)*"]", string(unit(ysig[1,i])))
        end
    else
        @static if VERSION >= v"0.7.0-DEV.2005"
            @warn "ModiaMath.plot: Variable $ykeyName is not plotted because variables with more as one dimension not yet supported."
        else
            warn("\nModiaMath.plot: Variable $ykeyName is not plotted because variables with more as one dimension not yet supported.")
        end
        return (nothing, nothing, nothing, nothing)
    end

    ysig = ustrip.(ysig)

    return (xsig, xsigLegend, ysig, ysigLegend)
end

"""
    mutable struct ResultWithVariables

Struct that is generated as a result of ModiaMath.simulate!(...),
if the model is a Modia3D.AbstractComponentWithVariables struct.
"""
mutable struct ResultWithVariables
    "Dictionary of time series (elements can be Float64/Int64/Bool Vectors or Matrices)"
    series::StringDictAnyResult

    "Dictionary of variables (elements are <: ModiaMath.AbstractVariable with required fields value and unit)"
    var::Dict{Symbol,Any}

    "String used as optional heading of a plot window"
    resultHeading::String

    ResultWithVariables(series, var, resultHeading="") = new(series, var, resultHeading)
end

resultHeading(result::ResultWithVariables) = result.resultHeading

function resultTimeSeries(result::ResultWithVariables, name, xLabel::Bool, xAxis)
    seriesDict = result.series
    (ysig, ykeyName, yNameAsString) = getSignal(seriesDict, string(name))
    symbol_ykeyName = Symbol(ykeyName)

    if ysig == nothing
        return (nothing, nothing, nothing, nothing)
    end

    yvar = result.var[symbol_ykeyName]
    if ndims(ysig) == 1
        # ysig is a vector
        ysigLegend = [appendUnit(yNameAsString, yvar.unit)]
    else
        # sig has more as one dimension
        @static if VERSION >= v"0.7.0-DEV.2005"
            ysigLegend = Array{String}(undef, length(yvar.value))
        else
            ysigLegend = Array{String}(length(yvar.value))
        end

        for i in eachindex(yvar.value)
            ysigLegend[i] = appendUnit(ModiaMath.indexToString(symbol_ykeyName, yvar.value, i), yvar.unit)
        end
    end

    # Get x-axis signal
    (xsig, xkeyName, xName) = getSignal(seriesDict, string(xAxis))
    if xsig == nothing
        return (nothing, nothing, nothing, nothing)
    end
    if ndims(xsig) != 1
        @static if VERSION >= v"0.7.0-DEV.2005"
            @warn "ModiaMath.plot: argument xAxis (= $xAxis) does not characterize a signal vector."
        else
            warn("\nModiaMath.plot: argument xAxis (= $xAxis) does not characterize a signal vector.")
        end
        return (nothing, nothing, nothing, nothing)
    end
    xsigLegend = xLabel ? appendUnit(xName, result.var[Symbol(xkeyName)].unit) : ""
    return (xsig, xsigLegend, ysig, ysigLegend)
end

getHeading(result, heading) = heading != "" ? heading : resultHeading(result)

ext: jl


text:
# dummy type wrapping a Float64 used in tests
struct F64 <: AbstractFloat
    x::Float64
end
F64(x::F64) = x

# operations
for op in (:+, :-)
    @eval Base.$op(a::F64) = F64($op(a.x))
end
for op in (:+, :-, :*, :/)
    @eval Base.$op(a::F64, b::F64) = F64($op(a.x, b.x))
end
for op in(:zero, :one)
    @eval Base.$op(::Type{F64}) = F64($op(Float64))
end

using Random
Random.rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{F64}}) = F64(rand(rng, Float64))
Random.randn(rng::Random.AbstractRNG, ::Type{F64}) = F64(randn(rng, Float64))

Base.sqrt(a::F64) = F64(sqrt(a.x))

# comparison
Base.isapprox(a::F64, b::F64) = isapprox(a.x, b.x)
Base.:<(a::F64, b::F64) = a.x < b.x
Base.:<=(a::F64, b::F64) = a.x <= b.x
Base.eps(::Type{F64}) = eps(Float64)

# promotion
Base.promote_type(::Type{Float32}, ::Type{F64}) = Float64 # for eig
Base.promote_type(::Type{Float64}, ::Type{F64}) = Float64 # for vecnorm
Base.promote(a::F64, b::T) where {T <: Number} = a, F64(b)
Base.promote(a::T, b::F64) where {T <: Number} = F64(a), b
Base.convert(::Type{F64}, a::F64) = a
Base.convert(::Type{Float64}, a::F64) = a.x
Base.convert(::Type{F64}, a::T) where {T <: Number} = F64(a)

# for testing of eigen
Base.acos(a::F64) = F64(acos(a.x))
Base.cos(a::F64) = F64(cos(a.x))
Base.sin(a::F64) = F64(sin(a.x))


text:
@recipe function f(d::UvSpectrumDensity, interval::Tuple, npoints::Int = 100)
    xs = range(interval[1], stop = interval[2], length = npoints)
    ys = map(x -> evaluate(d, x), xs)
    xs, ys
end
@recipe function f(d::UvSpectrumDensity, interval::Interval, npoints::Int = 100)
    xs = range(interval.left, stop = interval.right, length = npoints)
    ys = map(x -> evaluate(d, x), xs)
    xs, ys
end
@recipe function f(d::UvSpectrumDensity, h::Histogram{<:Any,1}, npoints::Int = 5*length(h.edges[1]))
    @assert isa(h.edges[1], AbstractRange) "Plotting for uneven spaced histgrams is not yet implemented. "
    xs = range(h.edges[1][1], stop = h.edges[1][end], length = npoints)
    ys = map(x -> evaluate(d, x), xs)
    xs, ys .* step(h.edges[1])
end


text:
function solve(d)
    pos = [indexin(["|"], d[1])[1], 1]
    output = ""
    dir = [0, 1]
    steps = 0
    while true
        while d[pos[2]][pos[1]] != "+"
            val = d[pos[2]][pos[1]]
            if ismatch(r"[A-Z]", val)
                output *= val
            end
            pos += dir
            steps += 1
            if d[pos[2]][pos[1]] == " "
                return output, steps
            end
        end
        if abs.(dir) == [0, 1]
            if d[pos[2]][pos[1]-1] == "-"
                dir = [-1, 0]
            elseif d[pos[2]][pos[1]+1] == "-"
                dir = [1, 0]
            end
        else
            if d[pos[2]-1][pos[1]] == "|"
                dir = [0, -1]
            elseif d[pos[2]+1][pos[1]] == "|"
                dir = [0, 1]
            end
        end
        pos += dir
        steps += 1
    end
end

d = split.(readlines("input.txt"), "")
result = solve(d)
println(result[1])
println(result[2])


text:
include("common.jl")

const OUTPUT_DIR = joinpath(OUTPUT_ROOT_DIR, string(Dates.now()))
mkdir(OUTPUT_DIR)

const VERSION_PATH = joinpath(OUTPUT_DIR, "version.txt")
const RESULT_PATH = joinpath(OUTPUT_DIR, "result.csv")
const STDOUT_PATH = joinpath(OUTPUT_DIR, "stdout.txt")
const STDERR_PATH = joinpath(OUTPUT_DIR, "stderr.txt")

results_table = Justitia.CSVRecord(RESULT_PATH, NearestGoodBoxQPResult)

open(STDOUT_PATH, "w") do out
    open(STDERR_PATH, "w") do err
        redirect_stdout(out) do
            redirect_stderr(err) do
                Justitia.run_experiments!(
                    results_table,
                    INSTANCES,
                    APPROACHES,
                    NearestGoodBoxQPResult,
                    config = Dict("factory_config" => CONFIG),
                )
            end
        end
    end
end


text:
type ConstantLagDDEProblem{uType,tType,lType,isinplace,F,H,C,MM} <: AbstractConstantLagDDEProblem{uType,tType,lType,isinplace}
  f::F
  h::H
  u0::uType
  lags::lType
  tspan::Tuple{tType,tType}
  callback::C
  mass_matrix::MM
end

function ConstantLagDDEProblem(f,h,u0,lags,tspan; iip = isinplace(f,4), callback = nothing,mass_matrix=I)
  ConstantLagDDEProblem{typeof(u0),eltype(tspan),typeof(lags),iip,typeof(f),typeof(h),typeof(callback),typeof(mass_matrix)}(f,h,u0,lags,tspan,callback,mass_matrix)
end

type DDEProblem{uType,tType,lType,isinplace,F,H,C,MM} <: AbstractDDEProblem{uType,tType,lType,isinplace}
  f::F
  h::H
  u0::uType
  lags::lType
  tspan::Tuple{tType,tType}
  callback::C
  mass_matrix::MM
end

function DDEProblem(f,h,u0,lags,tspan; iip = isinplace(f,4), callback = nothing,mass_matrix=I)
  DDEProblem{typeof(u0),eltype(tspan),typeof(lags),iip,typeof(f),typeof(h),typeof(callback),typeof(mass_matrix)}(f,h,u0,lags,tspan,callback,mass_matrix)
end


text: println("go")

max_stars_repo_name: timm/misc

max_issues_repo_name: timm/misc

max_forks_repo_name: timm/misc


text:
using Documenter
using Krippendorff

push!(LOAD_PATH,"../src/")
makedocs(
    sitename = "Krippendorff.jl",
    format = Documenter.HTML(),
    pages = [
            "Krippendorffs Alpha" => "index.md",
            "Important Functions" => "importantFunctions.md",
            "About" => "about.md"
         ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/FPGro/Krippendorff.jl.git",
    devbranch = "main"
)

max_stars_repo_name: FPGro/Krippendorff

max_issues_repo_name: FPGro/Krippendorff

max_forks_repo_name: FPGro/Krippendorff


text:
"""
Julia program to find the N'th Lucas Number.
Lucas Numbers is a sequence similar to Fibonacci numbers and is defined
as the sum of its two immediately previous terms.
The first and second numbers in Lucas sequence are 2 and 1 respectively.
"""

function lucas_num(n)
    l1 = 2
    l2 = 1
    if(n == 0)
        return l1
    end

    for i in (2:n)
        next = l1 + l2
        l1 = l2
        l2 = next
    end
    return l2
end

print("Enter the value of n(where you need the nth lucas number): ")
n = readline()
n = parse(Int, n)
if(n < 0)
    print("Invalid Value of n !!!")
    exit()
end
res = lucas_num(n)
print("The $n'th Lucas Number is $res.")

"""
Time Complexity - O(n), where `n` is the given number.
Space Complexity - O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE I

Enter the value of n(where you need the nth lucas number): 25
The 25'th Lucas Number is 167761.
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
# using packages that provide functionality for testing
using Test
using Logging
using TestSetExtensions
# import this package and its dependencies
#import TimeSeries
#import DataFrames
import InfrastructureSystems
const IS = InfrastructureSystems
import WaterSystems
const WSY = WaterSystems

const LOG_FILE = "water-systems-test.log"

LOG_LEVELS = Dict(
    "Debug" => Logging.Debug,
    "Info" => Logging.Info,
    "Warn" => Logging.Warn,
    "Error" => Logging.Error,
)

function get_logging_level(env_name::String, default)
    level = get(ENV, env_name, default)
    log_level = get(LOG_LEVELS, level, nothing)
    if isnothing(log_level)
        error("Invalid log level $level: Supported levels: $(values(LOG_LEVELS))")
    end

    return log_level
end

function run_tests()
    console_level = get_logging_level("SYS_CONSOLE_LOG_LEVEL", "Error")
    console_logger = ConsoleLogger(stderr, console_level)
    file_level = get_logging_level("SYS_LOG_LEVEL", "Info")

    IS.open_file_logger(LOG_FILE, file_level) do file_logger
        multi_logger = IS.MultiLogger(
            [console_logger, file_logger],
            IS.LogEventTracker((Logging.Info, Logging.Warn, Logging.Error)),
        )
        global_logger(multi_logger)

        include("test_utils/get_test_data.jl")

        @time @testset "Begin WaterSystems tests" begin
            @includetests ARGS
        end

        # TODO: Enable this once all expected errors are not logged.
        #@test length(IS.get_log_events(multi_logger.tracker, Logging.Error)) == 0

        @info IS.report_log_summary(multi_logger)
    end
end

logger = global_logger()

try
    run_tests()
finally
    # Guarantee that the global logger is reset.
    global_logger(logger)
    nothing
end

## old tests -- probably do not work, JJS 12/19/19
# println("Read data from WNTR networks 1, 2, and 3.")
# #@time @test include(joinpath(TEST_DIR, "readnetworkdata.jl"))
# println("Test all constructors.")
# @time @test include(joinpath(TEST_DIR, "constructors.jl"))
# println("Test Network Struct.")
# @time @test include(joinpath(TEST_DIR, "network.jl"))
# println("Test WaterSystem struct.")
# @time @test include(joinpath(TEST_DIR, "watersystem.jl"))
# println("Test wntr parser.")
# @time @test include(joinpath(TEST_DIR, "test_dict_parser.jl"))
# @test true


text:
using WaterLily
using GLMakie
using FileIO
using MeshIO

path = abspath(joinpath(@__DIR__, "assets", "sleigh.obj"))
sleigh = load(path)

mesh(sleigh)


text:
#!/usr/bin/env julia

using SphericalFourierBesselDecompositions
SFB = SphericalFourierBesselDecompositions

using Test
using DelimitedFiles
using LinearAlgebra

using Profile

win_features = [(), (:separable,)]

@testset "Covariance Matrix f=$features" for features in win_features
    @testset "Full window" begin
        rmin = 500.0
        rmax = 1000.0
        nbar = 3e-4
        amodes = SFB.AnlmModes(2, 2, rmin, rmax)
        cmodes = SFB.ClnnModes(amodes, Δnmax=0)
        wmodes = SFB.ConfigurationSpaceModes(rmin, rmax, 1000, amodes.nside)
        win = SFB.make_window(wmodes, :radial, :ang_quarter, :rotate, features...)

        # calc shot noise
        w̃mat, vmat = SFB.bandpower_binning_weights(cmodes; Δℓ=1, Δn=1)
        bcmodes = SFB.ClnnBinnedModes(w̃mat, vmat, cmodes)
        bcmix = SFB.power_win_mix(win, w̃mat, vmat, wmodes, bcmodes)
        NW_th = SFB.win_lnn(win, wmodes, cmodes) ./ nbar
        N_th = bcmix \ NW_th

        # calc theoretical power
        pkdata, header = readdlm((@__DIR__)*"/data/pk_m.dat", header=true)
        pk = SFB.Spline1D(Float64.(pkdata[:,1]), Float64.(pkdata[:,2]),
                          extrapolation=SFB.Splines.powerlaw)
        C_th = SFB.gen_Clnn_theory(pk, cmodes)
        CN_th = C_th .+ N_th

        # calc covariance, chain
        @time VW_A1 = SFB.Covariance.calc_covariance_exact_A1(CN_th, win, wmodes, cmodes)
        @time VW_chain = SFB.calc_covariance_exact_chain(C_th, nbar, win, wmodes, cmodes)

        #Profile.clear()
        #@time @profile SFB.calc_covariance_exact_chain(CN_th, win, wmodes, cmodes)
        #Profile.print()

        @show size(VW_A1) size(VW_chain)
        @test size(VW_A1) == size(VW_chain)

        @test issymmetric(VW_A1)
        @test issymmetric(VW_chain)
        @test all(isfinite.(VW_A1))
        @test all(isfinite.(VW_chain))

        @show VW_chain
        @show VW_A1
        @show VW_chain ./ VW_A1
        @test_broken VW_chain ≈ VW_A1  # We don't really expect these to agree.
    end
end

# vim: set sw=4 et sts=4 :


text:
using IRTools
using SymbolicUtils: term
using MacroTools: isexpr
using SymbolicUtils
using SymbolicCodegen

struct SymCall{D}
    dict::D
end

function _symcall(s::SymCall, f, args...)
    dict = s.dict
    applicable(f, args...) || return symterm(dict, f, args...)
    return s(f, args...)
end

function symterm(dict, f, args...)
    constargs = (SymbolicCodegen.toconst(arg, dict) for arg in args)
    T = typeof(f(constargs...))
    term(f, args...; type=T)
end

IRTools.@dynamo function (s::SymCall)(a...)
    ir = IRTools.IR(a...)
    ir == nothing && return
    for (x, st) in ir
        isexpr(st.expr, :call) || continue
        ir[x] = IRTools.xcall(_symcall, IRTools.self, st.expr.args...)
    end
    return ir
end

# EXAMPLE

# f(x::Real) = x+1
# h(x) = x^2 + 2
# g(x) = h(x) + f(x)

# @syms x::Int

# dict = Dict(:x => 3)
# using SymbolicUtils: symtype
# s = SymCall(dict)(g,x)
# symtype(s)


text:
#=
Struct defs for using pycachesim
=#
function make_power_of_two(n)
    if n < 2
        return Culong(2);
    end
    if (n & (n - 1)) == 0
        return Culong(n);
    else
        # not a power of two, reduce to next lower power of two
        # n is a Cuint which is 32bit unsigned
        m = 2^31;
        while m>n
            m = m/2;
        end
        return Culong(m);
    end
end

function rep_policy_id(replacement_policy)
    if replacement_policy == "FIFO"
        return Cint(0);
    elseif replacement_policy == "LRU"
        return Cint(1);
    elseif replacement_policy == "MRU"
        return Cint(2);
    elseif replacement_policy == "RR"
        return Cint(3);
    else
        return Cint(1); # default to LRU
    end
end

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

struct cache_entry
    cl_id::Culong;

    dirty::Culong;  # if 0, content is in sync with main memory. if 1, it is not.
                        #used for write-back
    invalid::Culong;# denotes an entry which does not contain a valid cacheline.
                        # it is empty.
    cache_entry() = new(0,1,1);
end

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

struct addr_range
    # Address range used to communicate consecutive accesses
    # last addr of range is addr+length-1
    addr::Culong;
    length::Culong;
end

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

struct stats
    count::Culonglong;
    byte::Culonglong;
    #cl::Culong; # might be used later
end

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

struct Cache
    name::Cstring;
    sets::Culong;
    ways::Culong;
    cl_size::Culong;
    cl_bits::Culong;
    subblock_size::Culong;
    subblock_bits::Culong;
    replacement_policy_id::Cint;    #// 0 = FIFO, 1 = LRU, 2 = MRU, 3 = RR
                                    #// (state is kept in the ordering)
                                    #// for LFU an additional field would be required to capture state
    write_back::Cint;       # 1 = write-back
                            # 0 = write-through
    write_allocate::Cint;   # 1 = write-allocate,
                            # 0 = non-write-allocate
    write_combining::Cint;  # 1 = this is a write-combining cache
                            # 0 = regular cache_entry

    load_from::Ptr{Cache};
    store_to::Ptr{Cache};
    victims_to::Ptr{Cache};

    swap_on_load::Cint;

    placement::Ptr{cache_entry};
    subblock_bitfield::Ptr{UInt8};

    LOAD::stats;
    STORE::stats;
    HIT::stats;
    MISS::stats;
    EVICT::stats;

    verbosity::Cint;

    Cache(name, sets, ways, cl_size, replacement_policy, load_from = 0, store_to = 0) = (
        cname = Base.unsafe_convert(Cstring, name);
        replacement_policy_id = rep_policy_id(replacement_policy);
        cl_size2 = make_power_of_two(cl_size);
        load_from = load_from==0 ? Ptr{Cache}() : load_from;
        store_to = store_to==0 ? Ptr{Cache}() : store_to;

        new(cname, sets, ways, cl_size2, 0, cl_size2, 0, replacement_policy_id, 1, 1, 0,
            load_from, store_to, Ptr{Cache}(), 0, Ptr{cache_entry}(), Ptr{UInt8}(),
            stats(0,0),stats(0,0),stats(0,0),stats(0,0),stats(0,0), 0)
    );
end


text:
"""
We define the `StochasticMatrix` type and associated functionality herein. This script was
written as an exercise in interfacing with and extending the Base library. If one intends to
work with stochastic matrices in practice, the more idiomatic approach is to work with
matrices, but perform checks with the `isstochastic` function. Why does this make sense?
Core operations such as matrix multiplication, getindex, setindex!, show, etc. work the same
on matrices as they do on stochastic matrices. One would need to (trivially) extend all such
functions to operate on stochastic matrices. This is wasteful and more trouble than it is
worth. Note that it DOES make sense to implement e.g. `Diagonal` matrices seperately from
general matrices, for operations such as matrix multiplication can be much more efficient on
the former than on the latter.
"""

"""
    StochasticMatrix{T, S<:AbstractMatrix{T} <: AbstractMatrix{T}

A square `Matrix` whose elements are all between `0` and `1` and whose rows all sum to `1`.

Our inner constructor calls `isstochastic` for safety.

##### Fields
- `matrix::S` : the stochastic matrix.

##### Complexity (for default constructor)
- Time complexity: `O(n^2)`, where `n` denotes the number of rows (or columns) in a square
matrix.
"""
struct StochasticMatrix{T, S<:AbstractMatrix{T}} <: AbstractMatrix{T}
    matrix::S
    function StochasticMatrix{T, S}(matrix) where {T, S<:AbstractMatrix{T}}
        if !isstochastic(matrix)
            error("Oh no! The matrix that you passed is not stochastic.")
        end
        return new{T, S}(matrix)
    end
end

"""
    StochasticMatrix(matrix::AbstractMatrix{T}) where {T} -> StochasticMatrix

Construct a `StochasticMatrix`.
"""
function StochasticMatrix(matrix::AbstractMatrix{T}) where {T}
    return StochasticMatrix{T, typeof(matrix)}(matrix)
end

"""
    Base.show(io::IO, stochasticmatrix::StochasticMatrix)

Show a `StochasticMatrix` in a human-friendly manner.
"""
function Base.show(io::IO, stochasticmatrix::StochasticMatrix)
    println(typeof(stochasticmatrix), ":")
    println(stochasticmatrix.matrix)
end

"""
    Base.size(stochasticmatrix::StochasticMatrix) -> Tuple{Int, Int}

Return a 2-tuple `(nrows, ncolumns)`.
"""
function Base.size(stochasticmatrix::StochasticMatrix)
    return size(stochasticmatrix.matrix)
end

"""
    Base.getindex(stochasticmatrix::StochasticMatrix)
"""
function Base.getindex(stochasticmatrix::StochasticMatrix, i::Integer, j::Integer)
    return Base.getindex(stochasticmatrix.matrix, i, j)
end

"""
    Base.setindex!(stochasticmatrix::StochasticMatrix, i::Integer, j::Integer) -> nothing
"""
function Base.setindex!(stochasticmatrix::StochasticMatrix, i::Integer, j::Integer)
    Base.setindex!(stochasicmatrix.matrix, i, j)
    return nothing
end

"""
    isstochastic(matrix::AbstractMatrix, tol::Float64 = 1e-5) -> Bool

Decide if `matrix` is stochastic.

`tol` is the tolerance associated with the row sums. We consider the row sums to be valid if
and only if they are in the interval `[1-tol, 1+tol]`. This is important when the element
type is `Float64`, for example.

##### Complexity
- Time complexity: `O(n^2)`, where `n` denotes the number of rows (or columns) in a square
matrix.
"""
function isstochastic(matrix::AbstractMatrix, tol::Float64 = 1e-5)
    # if non-square
    if size(matrix)[1] != size(matrix)[2]
        return false
    end
    # check if elements are in [0,1]
    for element in matrix
        if element < 0 || element > 1
            return false
        end
    end
    # check row sums
    for row in eachrow(matrix)
        rowsum = sum(row)
        if rowsum < 1 - tol || rowsum > 1 + tol
            return false
        end
    end
    return true  # this line is reached only if all above checks pass
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
d1 = [1 => 4.2, 2 => 5.3]
d1 = {1 => 4.2, 2 => 5.3}
d2 = {"a"=>1, (2,3)=>true}

d3 = [:A => 100, :B => 200]
d3[:B] #> 200
# d3[:Z] #> ERROR: key not found: :Z
get(d3, :Z, 999) #> 999
in( (:Z, 999), d3) #> false
in( (:A, 100), d3) #> true
d3[:A] = 150 #> d3 is now [:A => 150, :B => 200]
d3[:C] = 300 #> d3 is now [:A => 150, :B => 200, :C => 300]
length(d3) #> 3
# d3["CVO"] = 500 #> ERROR: CVO is not a valid key for type Symbol#
# d3[:CVO] = "Julia" #> ERROR: `convert` has no method matching     # convert(::Type{Int64}, ::ASCIIString)

d4 = Dict()
d4["lang"] = "Julia"

d5 = Dict{Float64, Int64}()
# d5["c"] = 6 #> ERROR: `convert` has no method matching convert(::Type{Float64}, ::ASCIIString)

d3 = [:A => 100, :B => 200]
delete!(d3, :B) #> Dict{Symbol,Int64} with 1 entry:  :A => 100

d3 = [:A => 100, :B => 200]
ki = keys(d3) #> KeyIterator
for k in ki
   println(k)
end #> A   B
:A in ki #> true
:Z in ki #> false
haskey(d3, :A) #> true
collect(keys(d3)) #>
vi = values(d3) #> ValueIterator
for v in values(d3)
    println(v)
end #> 100  200

keys1 = ["J.S. Bach", "Woody Allan", "Barak Obama"]
values1 = [ 1685, 1935, 1961]
d5 = Dict(keys1, values1)
#>
d5 = [ "J.S. Bach" => 1685,
       "Woody Allan" => 1935,
       "Barak Obama" => 1961 ]
#> Dict{ASCIIString,Int64} with 3 entries:
for (k, v) in d5
  println("$k was born in $v")
end
for p in d5
  println("$(p[1]) was born in $(p[2])")
end
#J.S. Bach was born in 1685
#Barak Obama was born in 1961
#Woody Allan was born in 1935

dpairs = ["a", 1, "b", 2, "c", 3]
d6 = [dpairs[i]=>dpairs[i+1] for i in 1:2:length(dpairs)]
d6 = (String => Int64)[dpairs[i]=>dpairs[i+1] for i in 1:2:length(dpairs)]

# use of a dictionary:
function showfactor(n)
 d = factor(n)
 println("factors for $n")
 for (k, v) in d
     print("$k^$v\t")
 end
end

# neat tricks:
dict = d6
arrkey = [key for (key, value) in dict] #> 3-element Array{Any,1}: "c" "b" "a"
arrval = [value for (key, value) in dict] #> 3-element Array{Any,1}: 3 2 1

dict = d1
n = 1
arrn = [dict[i] for i = 1:n] #> 1-element Array{Any,1}: 4.2
arrn = map((i) -> dict[i], [1:n]) #> 1-element Array{Any,1}: 4.2

max_stars_repo_path: Module 1/Chapter05/Old/Old/dicts.jl

max_issues_repo_path: Module 1/Chapter05/Old/Old/dicts.jl

max_forks_repo_path: Module 1/Chapter05/Old/Old/dicts.jl


text:
module VSCodeLiveUnitTesting

include("../../OrderedCollections/src/OrderedCollections.jl")
include("../../CodeTracking/src/CodeTracking.jl")

module JuliaInterpreter
    using ..CodeTracking

    include("../../JuliaInterpreter/src/packagedef.jl")
end

module LoweredCodeUtils
    using ..JuliaInterpreter
    using ..JuliaInterpreter: SSAValue, SlotNumber, Frame
    using ..JuliaInterpreter: @lookup, moduleof, pc_expr, step_expr!, is_global_ref, whichtt,
                        next_until!, finish_and_return!, nstatements, codelocation,
                        is_return, lookup_return, is_GotoIfNot, is_ReturnNode

    include("../../LoweredCodeUtils/src/packagedef.jl")
end

module Revise
    using ..OrderedCollections
    using ..CodeTracking
    using ..JuliaInterpreter
    using ..LoweredCodeUtils

    using ..CodeTracking: PkgFiles, basedir, srcfiles, line_is_decl, basepath
    using ..JuliaInterpreter: whichtt, is_doc_expr, step_expr!, finish_and_return!, get_return,
                        @lookup, moduleof, scopeof, pc_expr, is_quotenode_egal,
                        linetable, codelocs, LineTypes, is_GotoIfNot, isassign, isidentical
    using ..LoweredCodeUtils: next_or_nothing!, trackedheads, structheads, callee_matches

    include("../../Revise/src/packagedef.jl")
end

include("copied_from_test_reports.jl")

function run_test_loop(test_file::AbstractString)
    Pkg.status()

    println()

    test_folder = dirname(test_file)

    try
        cd(test_folder)
        VSCodeLiveUnitTesting.Revise.track(test_file; mode=:eval, skip_include=false)
    catch err
        Base.display_error(err, catch_backtrace())
    end

    VSCodeLiveUnitTesting.Revise.entr([test_file]; all=true, postpone=true) do
        println()
        println("Rerunning tests...")
        println()

        try
            cd(test_folder)
            VSCodeLiveUnitTesting.Revise.include(test_file)
        catch err
            Base.display_error(err, catch_backtrace())
        end
    end
end

function live_unit_test(pkg_name::AbstractString, test_file::AbstractString)
    absolute_path_of_test_file = joinpath(pwd(), test_file)

    pkgspec = deepcopy(Pkg.PackageSpec(pkg_name))

    ctx = Pkg.API.Context()

    if !isinstalled!(ctx, pkgspec)
        error("Package not in here")
    end

    Pkg.instantiate(ctx)

    testfilepath = gettestfilepath(ctx, pkgspec)

    if !isfile(testfilepath)
        error("")
    end

    test_folder_has_project_file = isfile(joinpath(dirname(testfilepath), "Project.toml"))

    if VERSION >= v"1.4.0" || (VERSION >= v"1.2.0" && test_folder_has_project_file)
        # Operations.sandbox() has different arguments between versions
        sandbox_args = (ctx,
                            pkgspec,
                            pkgspec.path,
                            joinpath(pkgspec.path, "test"))

        if VERSION >= v"1.4.0"
            test_project_override = test_folder_has_project_file ?
                    nothing :
                    gen_target_project(ctx, pkgspec, pkgspec.path, "test")

            sandbox_args = (sandbox_args..., test_project_override)
        end

        sandbox(sandbox_args...) do
            run_test_loop(absolute_path_of_test_file)
        end
    else
        with_dependencies_loadable_at_toplevel(ctx, pkgspec; might_need_to_resolve=true) do localctx
            Pkg.activate(localctx.env.project_file)
            try
                run_test_loop(absolute_path_of_test_file)
            finally
                Pkg.activate(ctx.env.project_file)
            end
        end
    end
end

end

max_stars_repo_name: novitk/julia-vscode

max_issues_repo_name: novitk/julia-vscode

max_forks_repo_name: novitk/julia-vscode


text:
module PolicyViz2

using GridInterpolations, Interact, PGFPlots, Colors, ColorBrewer, HDF5

include("./viz_policy_constants.jl")
include("./viz_policy.jl")
include("./nnet_calculations.jl")
include("./policy_functions.jl")

end #module PolicyViz

max_stars_repo_name: alisamiii/VerticalCAS

max_issues_repo_name: alisamiii/VerticalCAS

max_forks_repo_name: bethelhall/CoCoCAS


text:
using Statistics
"""
    simu_qp(h2 = 0.8)
---
# Description
Simulation of QTL and phenotypes with the toy dataset comes with this package.
The data were randomly sampled on ID and SNP from a project I worked on.

# Results
- `toy.raw`: from `plink`.
- `toy.frq`: from `raw2gt`.
- `toy.gt`:  from `raw2gt`. genotype of `012` only, without spaces.
- `toy.ph`:  simulated phenotypes corresponding to ID in `toy.fam`.
- `toy.bv`: the true breeding values.
- `toy.qtl`: QTL location index and value.
- `toy.G`: `G` matrix with vanRaden method I, with diagonals added 1e-4.
- `toy.giv`: Above inversed.
"""
function simu_qp(h2 = 0.8)
    subtitle("Simulation QTL and phenotypes with toy data")
    item("Prepare genotypes")
    isdir(abg_rst) || mkdir(abg_rst)
    nlc, nid = countlines.(joinpath.(abg_dat, ["toy.bim", "toy.fam"]))

    run(`$plink --bfile $abg_dat/toy --recode A --out $abg_rst/toy`)
    run(pipeline("$abg_rst/toy.raw", `$abg_bin/raw2gt $abg_rst/toy.frq`, "$abg_rst/toy.gt"))
    _, twop, Z = read_gt_n_frq("$abg_rst/toy.gt", "$abg_rst/toy.frq")
    done()

    item("QTL, BV, and phenotype simulaiton")
    nq = nlc ÷ 20
    lq = sort(rand(1:nlc, nq))
    eq = randn(nq) .* sqrt(1. / nq)
    bv = Z[lq, :]'eq
    vg = var(bv)
    ve = vg / h2 * (1 - h2)
    er = randn(nid) .* sqrt(ve)
    ph = bv + er
    write(joinpath(abg_rst, "toy.ph"), join(ph, '\n'), '\n')
    write(joinpath(abg_rst, "toy.bv"), join(bv, '\n'), '\n')
    open(joinpath(abg_rst, "toy.qtl"), "w") do io
        for i in 1:nq
            println(io, lq[i], ' ', eq[i])
        end
    end
    done()

    item("G and its inverse")
    G = vanRaden(Z, twop, 1e-4)
    giv = inv(G)
    write(joinpath(abg_rst, "toy.G"), G)
    write(joinpath(abg_rst, "toy.giv"), giv)
    done()
end


text:
using Unitful
using DataStructures

t = range(0.0, stop=10.0, length=100)

result = OrderedDict{String,Any}()

result["time"] = t*u"s"
result["phi"]  = sin.(t)u"rad"
result["w"]    = cos.(t)u"rad/s"

plot(result, ("phi", "w"), prefix="Sim 1:", heading="Test reuse")

result["phi"]  = 1.2*sin.(t)u"rad"
result["w"]    = 0.8*cos.(t)u"rad/s"

plot(result, ("phi", "w"), prefix="Sim 2:", reuse=true)

println("\n... result info:")
printResultInfo(result)


text:
"""
Package: ForecastPlots

    dplot(x; dlabels = ["Data","Trend","Seasonal","Remainder"])

Plot a seasonal, trend and remainder decomposition of a time series.

# Arguments
- `x`: Matrix with three columns containing sesonality, trend and remainder decomposition
- `dlabels`: String array of length four containing labels for data, seasonality, trend and remainder.

# Returns

Sesonal decomposition plot with scale bars

# Examples
```julia-repl
dplot(randn(100,3))
dplot(randn(100,3); labels = ["數據","趨勢","季節性","剩餘"])
dplot(rand(100,3),now()-Day(99):Day(1):now())
```
"""
function dplot(x::Matrix{<:T}, t=1:size(x,1);
               dlabels = ["Data","Trend","Seasonality","Remainder"], kw...) where T<: Real

    @assert length(dlabels) == 4 "Four labels needed for 'Data','Trend','Seasonality' and 'Remainder'"
    @assert x isa Matrix{<:Real} "Numerical Matrix with columns expected"
    @assert size(x,2) == 3 "Numerical Matrix with columns expected"

    seasonality = x[:,1]
    trend = x[:,2]
    remainder = x[:,3]

    xtickfont = font(3, "Courier")
    ytickfont = font(3, "Courier")
    yguidefont = font(5, "Courier")

    # reference bar
    a,b = extrema(skipmissing(values(seasonality.+
                                     trend.+
                                     remainder))); hd = b-a
    a,b = extrema(skipmissing(values(seasonality))); hs = b-a
    a,b = extrema(skipmissing(values(trend))); ht = b-a
    a,b = extrema(skipmissing(values(remainder))); hr = b-a
    mh = min(hd,hs,ht,hr)

    inset_subplots  = [(1, bbox(-0.012, 0, 0.01, 1.0, :left)),
                       (2, bbox(-0.012, 0, 0.01, 1.0, :left)),
                       (3, bbox(-0.012, 0, 0.01, 1.0, :left)),
                       (4, bbox(-0.012, 0, 0.01, 1.0, :left))]

    layout = @layout [Data
                      Seasonal
                      Trend
                      Remainder]

    rect(w, h, x, y) = Shape(x .+ [0, w, w, 0, 0], y .+ [0, 0, h, h, 0])

    ttype = eltype(t) <: TimeType
    sbd = 10
    sbx = ttype ? Dates.value(t[1]) - (sbd+1)*Dates.value(t[2]-t[1])  : -sbd
    sbw = ttype ? sbd*Dates.value(t[2]-t[1])÷4 : 10÷4

    sbc = :gray

    p1 = plot(t,seasonality .+ trend .+ remainder;
              yguide = dlabels[1], xguide = "",
              ytickfont, yguidefont,
              bottom_margin = -5Plots.mm,
              legend = false, grid = false,
              kw...,
              xticks = false)

    p1 = plot!(rect(sbw, hr,sbx, minimum(seasonality .+ trend .+ remainder));
               seriescolor = :grey,
               linewidth = 0,
               label = nothing)

    p2 = plot(t,trend;
              yguide = dlabels[2], xguide = "", title = "", xticks = false,
              ytickfont, yguidefont,
              bottom_margin = -5Plots.mm,
              legend = false, grid = false)

    p2 = plot!(rect(sbw, hr, sbx, minimum(trend));
               seriescolor = :grey,
               linewidth = 0,
               label = nothing)

    p3 = plot(t,seasonality;
              yguide = dlabels[3], xguide = "", title = "", xticks = false,
              ytickfont, yguidefont, seriestype = :sticks,
              bottom_margin = -5Plots.mm,
              legend = false, grid = false)

    p3 = plot!(rect(sbw, hr, sbx, minimum(seasonality));
               seriescolor = :grey,
               linewidth = 0,
               label = nothing)

    p4 = plot(t,remainder;
              yguide = dlabels[4], xguide = "",
              xtickfont, ytickfont, yguidefont, seriestype = :sticks,
              legend = false, grid = false,
              kw...,
              title = "")

    p4 = plot!(rect(sbw , hr, sbx, minimum(remainder));
               seriescolor = :grey,
               linewidth = 0,
               label = nothing)

    plot(p1,p2,p3,p4; layout)

end


text:
#-----------------------------------------------------------------------
module squirrel     # BEGIN MODULE SQUIRREL
#-----------------------------------------------------------------------

using LinearAlgebra, Combinatorics, Statistics , DoubleFloats
using ForwardDiff, DiffResults

include("type.jl")
include("seval.jl")
include("geosol.jl")
include("broyden.jl")
include("outlier.jl")

include("metric.jl")
include("metrics/Minkowski.jl")
include("srl/FHC22.jl")
include("srl/RTC21.jl")
include("srl/mloc.jl")

#-----------------------------------------------------------------------
#
#   GEODESIC ENDPOINT AND JACOBIAN
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    gsolve( Xi::RealVec , Vi::RealVec , g::Function , tol::Real , integrator=AutoVern7(Rodas5()) )

The `gsolve` function takes an initial point `Xi` and four velocity `Vi`
and computes the endpoint of a future pointing null geodesic in the
metric func `g`. The variable `integrator` specifies the integration
scheme, and `tol` is the tolerance parameter.

"""
function gsolve( Xi::RealVec , Vi::RealVec , g::Function , tol::Real
                 , integrator=AutoVern7(Rodas5()) )
    tpfl=typeof(Xi[1])
    Z0 = zeros(tpfl,8)
    V0 = nullenforcerf( Vi , Xi , g )
    Z0 = vcat( Xi , V0 )
    return solveZ( Z0 , g , tol , tol , integrator , tol )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    V34( V3::RealVec )

The function `V34` takes a three-vector and constructs from it a
four-vector with vanishing time component.

"""
function V34( V3::RealVec )
    tpfl = typeof(V3[1])
    return [ tpfl(0) ; V3[1] ; V3[2] ; V3[3] ]
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    VidF( Zi::RealMtx )

The function `VidF` takes a matrix of initial data points and constructs
from it a single vector of initial 3-velocities for null vectors.

"""
function VidF( Zi::RealMtx )
    return vcat( Zi[6:8,1] , Zi[6:8,2] , Zi[6:8,3] , Zi[6:8,4] )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    zF( Vid::RealVec , Zi::RealMtx , g::Function , tol::Real )

The function `zF` returns differences between the endpoints of four
geodesics for the initial data encoded in `Vid` and `Zi`, and the metric
function `g`. The variable `tol` is the tolerance parameter for the
integration. This function vanishes when the four geodesics intersect.

"""
function zF( Vid::RealVec , Zi::RealMtx , g::Function , tol::Real )
    # Want to find roots of this wrt initial velocity
    tpfl=typeof(Vid[1])

    Xf = [zeros(tpfl,4) for _ in 1:4]

    Threads.@threads for i=1:4
        ( k1 , k2 ) = ( 1 + 3*(i-1) , 3 + 3*(i-1) )
        Xf[i] = gsolve( Zi[1:4,i] , V34( Vid[k1:k2] ) , g , tol )[1:4]
    end

    return vcat( Xf[1] - Xf[2] , Xf[1] - Xf[3] , Xf[1] - Xf[4] )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    gejac( Xi::RealVec , Vi::RealVec , g::Function , δ::Real )

The function `gejac` computes the endpoint of a geodesic and its
Jacobian. The variables have the same meaning as those in `gsolve`.

"""
function gejac( Xi::RealVec , Vi::RealVec , g::Function , tol::Real )
    v = Vi[2:4]
    result = DiffResults.JacobianResult(vcat(Xi,Vi),v)
    result = ForwardDiff.jacobian!(result,var->gsolve(Xi,V34(var),g,tol)
                                    , v )
    return ( DiffResults.value(result)
             , DiffResults.jacobian(result)[1:4,:] )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    geocJ( Zi::RealMtx , g::Function , δ::Real )

The function `geocJ` computes the Jacobian of the function `zF` from the
endpoint Jacobians computed by calling `gejac` function.

"""
function geocJ( Zi::RealMtx , g::Function , δ::Real )
    tpfl=typeof(Zi[1,1])
    v   = zeros(tpfl,3)
    J   = zeros(tpfl,12,12)
    dXV = [zeros(tpfl,4,3) for _ in 1:4]
    Zf  = copy(Zi)

    Threads.@threads for i=1:4
        res = gejac( Zi[1:4,i] , Zi[5:8,i] , g , δ )
        ( Zf[:,i] , dXV[i] ) = ( res[1] , res[2] )
    end

    J[1:4,1:3]  = dXV[1]
    J[5:8,1:3]  = dXV[1]
    J[9:12,1:3] = dXV[1]

    J[1:4,4:6]    = - dXV[2]
    J[5:8,7:9]    = - dXV[3]
    J[9:12,10:12] = - dXV[4]

    return (Zf,J)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#
#   SOLVING FOR INITIAL DATA
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    zFc( Zf::RealMtx )

The function `zFc` rearranges the endpoints `Zf` to match the output of
the function `zF`.

"""
function zFc( Zf::RealMtx )
    return vcat( Zf[1:4,1] - Zf[1:4,2] , Zf[1:4,1] - Zf[1:4,3] ,
                 Zf[1:4,1] - Zf[1:4,4] )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    idf( Zi::RealMtx , gfunc::Function , tol::Real , nb::Int )

The function `idf` takes as input `Zi`, a matrix formed from the
emission points and guesses for the initial four-velocities, and the
metric function `gfunc`. It outputs the corrected initial data for
intersecting null geodesics. The argument `tol` is the tolerance
parameter and `nb` is the Broyden termination limit. This function calls
`geocJ` to compute the initial Jacobian, then calls and applies the
Broyden root finding function `bsolve` to the function `zF` (also
called).

"""
function idf( Zi::RealMtx , gfunc::Function , tol::Real , nb::Int )
    tpfl=typeof(Zi[1,1])

    Zf = copy(Zi)
    V0 = VidF(Zi)
    res = geocJ( Zi , gfunc , tol )

    ( F0 , J ) = ( zFc(res[1]) , res[2] )

    V = bsolve( v->zF(v,Zi,gfunc,tol) , J , F0 , V0 , nb )

    Zf[5:8,1] = V34(  V[1:3]  )
    Zf[5:8,2] = V34(  V[4:6]  )
    Zf[5:8,3] = V34(  V[7:9]  )
    Zf[5:8,4] = V34( V[10:12] )

    return Zf
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#
#   RELATIVISTIC LOCATOR
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    locator4( X::RealMtx , Xc::RealVec , gfunc::Function , tol::Real , nb::Int=24 , idv::Bool=false , V::RealMtx=zeros(Float64,4,4) )

The function `locator4` computes the intersection point from a set of
four emission points `X` using the guess `Xc`. The intersection point is
computed by first finding the initial data using the function `idf`,
then integrating geodesics to obtain the result. If an improved guess
for the four-velocities is available, one can set `ìdv=true` and specify
the four-velocities as column vectors in the matrix `V`.

"""
function locator4( X::RealMtx , Xc::RealVec , gfunc::Function ,
                   tol::Real , nb::Int=24 , idv::Bool=false ,
                   V::RealMtx=zeros(Float64,4,4) )
    tpfl=typeof(X[1,1])

    Zi  = zeros(tpfl,8,4)
    Zf = copy(Zi)

    if idv
        for i=1:4
            Zi[:,i] = vcat( X[:,i] , V[:,i] )
        end
    else
        for i=1:4
            Zi[:,i] = vcat( X[:,i] , Xc - X[:,i] )
        end
    end

    Zi = idf( Zi , gfunc , tol , nb )
    Threads.@threads for i=1:4
        Zf[:,i] = gsolve( Zi[1:4,i] , Zi[5:8,i] , gfunc , tol )
    end
    return (Zf[1:4,1] + Zf[1:4,2] + Zf[1:4,3] + Zf[1:4,4])/4
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#
#   RELATIVISTIC LOCATION WITH >4 EMISSION POINTS
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    locator( X::RealMtx , gfunc::Function , δ::Real , ne::Int=5 , nb::Int=24 , outthresh::Real=2e1 , tpflc::DataType=Double64 )

The function `locator` computes the intersection point from a set of
`ne>4` emission points `X` by applying `locator4` to all combinations of
4 points out of `ne` in `X`. A basic outlier detection algorithm
(implemented in the function `odetc`) is applied to reduce errors.

"""
function locator(  X::RealMtx , gfunc::Function , δ::Real , ne::Int=5 ,
                   nb::Int=24 , outthresh::Real=2e1 ,
                   tpflc::DataType=Double64 )
    tpfl  = typeof(X[1,1])

    l = size(X)

    if gfunc == x->ημν(x)
        return mlocator( tpflc.(X) )
    elseif  l[2] < 4 || ne < 4
        print("Need more than four emission points.")
        return zeros(tpfl,4)
    elseif l[2] == 4 || ne == 4
        Xdual   = locator4FHC22( tpflc.(X) )
        X1 = locator4( X , Xdual[1] , gfunc , δ , nb , false )
        X2 = locator4( X , Xdual[2] , gfunc , δ , nb , false )
        return (X1,X2)
    elseif l[2] >= 5
        if  ne >= 5 && ne < size(X)[2]
            Y = X[:,1:ne]
        elseif size(X)[2] == 5 || ne >= size(X)[2]
            Y = X
        end

        Xc      = mlocator( tpflc.(Y) )
        W       = combX(Y)
        nr      = length(W)
        Xs      = [zeros(tpfl,4) for _ in 1:nr]

        for i=1:nr
            Xs[i]  = locator4( W[i] , tpfl.(Xc) , gfunc , δ , nb ,
                               false )
        end

        return mean(odetc( Xs , outthresh )[1])
    end
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
end     # END MODULE SQUIRREL
#-----------------------------------------------------------------------


text:
register(DataDep(
    "pudl-raw-eia860m",
    """
    Dataset: PUDL Raw EIA Form 860M
    Website: https://zenodo.org/record/4281337
    Author: Gosnell, Christina, Selvans, Zane
    Date of Publication: November 19, 2020
    License:

    Raw US Energy Information Administration (EIA) Form 860M data, archived from https://www.eia.gov/electricity/data/eia860m/

    EIA Form 860 M data for electric power plants with 1 megawatt or greater combined nameplate capacity.
    Preliminary Monthly Electric Generator Inventory (based on Form EIA-860M as a supplement to Form EIA-860).

 $PUDL_DESCRIPTION
    """,
    [
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/datapackage.json",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2015-07.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2015-08.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2015-09.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2015-10.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2015-11.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2015-12.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-01.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-02.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-03.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-04.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-05.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-06.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-07.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-08.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-09.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-10.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-11.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2016-12.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-01.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-02.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-03.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-04.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-05.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-06.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-07.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-08.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-09.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-10.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-11.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2017-12.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-01.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-02.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-03.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-04.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-05.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-06.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-07.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-08.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-09.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-10.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-11.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2018-12.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-01.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-02.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-03.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-04.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-05.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-06.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-07.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-08.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-09.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-10.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-11.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2019-12.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2020-01.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2020-02.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2020-03.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2020-04.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2020-05.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2020-06.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2020-07.xlsx",
        "https://zenodo.org/api/files/a23e66f5-0b22-473b-a11c-628c3dc9c0f9/eia860m-2020-08.xlsx",
    ],
))


text:
# specific to Julia, that is, related to the output file format (language).

#=
F: read the file containing constants
I:
O:
BA
=#
function jl_include_constants_from_literature(src_file_name, pad_string)
  # create src_buffer -
  src_buffer = ""
  # path to distrubtion -
  path_to_src_file = src_file_name
  open(path_to_src_file,"r") do src_file
    for line in eachline(src_file)
      src_buffer *= pad_string*"$line"
    end
  end
  return src_buffer
end

#=
F: generate kinetic buffer
I:
O: return monod affinity constant symbol, mRNA species array,
  protein species array, W_string_array, disassociation_const_string_array.
BA
=#
function jl_build_kinetics_buffer(all_species_dict::Dict{String, Int},
  all_rnx_list::Array, all_txtl_dict::Dict, sys2user::Dict)
  kinetics = "function calculate_kinetics(X, data_dictionary)" *
             "\n\n\t# load all species dictionary" *
             "\n\tall_species_dict = data_dictionary[\"all_species_dict\"]  # String-->Int"
  # for signaling
  kinetics *= "\n\n\t# generate kinetics for signaling, assume saturation kinetics"
  totalRnxNo = length(all_rnx_list)
  kinetics *= "\n\t# load Monod affinity constants" *
              "\n\tMonodK = data_dictionary[\"MonodAffinityConstantDict\"]  # String-->Float"
  MonodAffinityConstant_String_Array = Array{String,1}()  # collect monod affinity constant symbol
  kinetics *= "\n\t# load reaction kinetic constants" *
              "\n\tkcat = data_dictionary[\"kcat_signaling\"]  # in order of rnx" *
              "\n\n\t# write reaction rate equations" *
              "\n\trnx_rate_vector = zeros($totalRnxNo)"
  for (index, rnx) in enumerate(all_rnx_list)  # go thru every rnx
    tmp_rnx_rate = "kcat[$index]"
    if isdefined(rnx, :catalysts)  # enzyme
      for token in rnx.catalysts
        tmp_rnx_rate *= "*X[$(all_species_dict[token.oriBioName])]"
      end
    end
    if isdefined(rnx,:reactants)  # substrate
      for token in rnx.reactants
        tmp_rnx_rate *= "*X[$(all_species_dict[token.oriBioName])]/(MonodK[\"MonodK~rnx$(index)~$(token.oriBioName)\"]+X[$(all_species_dict[token.oriBioName])])"
        push!(MonodAffinityConstant_String_Array, "MonodK~rnx$(index)~$(token.oriBioName)")
      end
    end
    kinetics *= "\n\trnx_rate_vector[$index] = "*tmp_rnx_rate
  end

  # for TXTL
  W_string_array = Array{String,1}()  # W: weight of protein action on transcription
  disassociation_const_string_array = Array{String,1}()  # disassociation constant in tranfer function
  kinetics *= "\n\n\t# generate kinetics & control terms for TXTL" *
              "\n\t# load data from data dictionary" *
              "\n\tW_value_dict = data_dictionary[\"W_value_dict\"]  #String-->Float" *
              "\n\tcoop = data_dictionary[\"cooperativity\"]" *
              "\n\tdisassociation_constant_dict = data_dictionary[\"transferFunctionDisassociationConstantDict\"]" *
              "\n\tbackground_control_term = data_dictionary[\"backgroundGeneExpressionControlTermDict\"]" *
              "\n\tkcatTX = data_dictionary[\"kcatTranscription\"]" *
              "\n\tRNAPconc = data_dictionary[\"RNAPConcentration\"]" *
              "\n\tgeneConc = data_dictionary[\"avgGeneConcentration\"]" *
              "\n\tsaturationTX = data_dictionary[\"transcriptionSaturationConstant\"]" *
              "\n\tkcatTL = data_dictionary[\"kcatTranslation\"]" *
              "\n\tRIBOconc = data_dictionary[\"RIBOConcentration\"]" *
              "\n\tsaturationTL = data_dictionary[\"translationSaturationConstant\"]"
  totalTXTLNo = length(all_txtl_dict)
  kinetics *= "\n\n\t# write control terms, transcription rate equations, and translation rate equations" *
             "\n\tTX_control_term = Dict{String, Float64}()  # control term" *
             "\n\tTX_rate_vector = Dict{String, Float64}()  # transcription rate" *
             "\n\tTL_rate_vector = Dict{String, Float64}()  # translation rate"
  for (key, txtl) in all_txtl_dict # go thru every txtl
    tmp_protein_string = replace(key, sys2user["MRNA"] => sys2user["PROTEIN"], count=1)
    kinetics *= "\n\t# $key and $tmp_protein_string"
    up_factors_array = Array{String,1}()  # collection of upregulation factors name: targetedmRNA_factor(s)
    if !isempty(txtl.activationProtein)  # upregulation
      for token_array in txtl.activationProtein
        tmp_up = ""
        tmp_up_name = key*"_"
        for token in token_array
          tmp_W = "W~$(token.oriBioName)~$(key)"
          tmp_up *= "W_value_dict[\"$(tmp_W)\"] * X[$(all_species_dict[token.oriBioName])]^coop/((disassociation_constant_dict[\"KD~$(token.oriBioName)\"])^coop+X[$(all_species_dict[token.oriBioName])]^coop)*"
          push!(W_string_array, tmp_W)
          push!(disassociation_const_string_array, "KD~$(token.oriBioName)")
          tmp_up_name *= token.oriBioName*"_"
        end
        tmp_up = chop(tmp_up)
        tmp_up_name = chop(tmp_up_name)
        push!(up_factors_array, tmp_up_name)
        kinetics *= "\n\t$(tmp_up_name) = $(tmp_up)"
      end
    end
    down_factors_array = Array{String,1}()  # collection of downregulation factors name: targetedmRNA_factor(s)
    if !isempty(txtl.inhibitionProtein)  # downregulation
      for token_array in txtl.inhibitionProtein
        tmp_down = ""
        tmp_down_name = key*"_"
        for token in token_array
          tmp_W = "W~$(token.oriBioName)~$(key)"
          tmp_down *= "W_value_dict[\"$(tmp_W)\"] * X[$(all_species_dict[token.oriBioName])]^coop/((disassociation_constant_dict[\"KD~$(token.oriBioName)\"])^coop+X[$(all_species_dict[token.oriBioName])]^coop)*"
          push!(W_string_array, tmp_W)
          push!(disassociation_const_string_array, "KD~$(token.oriBioName)")
          tmp_down_name *= token.oriBioName*"_"
        end
        tmp_down = chop(tmp_down)
        tmp_down_name = chop(tmp_down_name)
        push!(down_factors_array, tmp_down_name)
        kinetics *= "\n\t$(tmp_down_name) = $(tmp_down)"
      end
    end
    # combine to form up and down, then merge up and down to form control term
    if !isempty(up_factors_array) && !isempty(down_factors_array)
      up_action = ""
      for up in up_factors_array
        up_action *= up*" +"
      end
      kinetics *= "\n\tup_action_on_$key = $(chop(up_action))"
      down_action = ""
      for down in down_factors_array
        down_action *= down*" +"
      end
      kinetics *= "\n\tdown_action_on_$key = $(chop(down_action))" *
                  "\n\tTX_control_term[\"$key\"] = (background_control_term[\"$key\"] + up_action_on_$key)/(1 + background_control_term[\"$key\"] + up_action_on_$key + down_action_on_$key)"
    elseif !isempty(up_factors_array)
      up_action = ""
      for up in up_factors_array
        up_action *= up*" +"
      end
      kinetics *= "\n\tup_action_on_$key = $(chop(up_action))" *
                  "\n\tTX_control_term[\"$key\"] = (background_control_term[\"$key\"] + up_action_on_$key)/(1 + background_control_term[\"$key\"] + up_action_on_$key)"
    elseif !isempty(down_factors_array)
      down_action = ""
      for down in down_factors_array
        down_action *= down*" +"
      end
      kinetics *= "\n\tdown_action_on_$key = $(chop(down_action))" *
                  "\n\tTX_control_term[\"$key\"] = background_control_term[\"$key\"]/(1 + background_control_term[\"$key\"] + down_action_on_$key)"
    else
      # seems impossible, a TXTL dict is constructed b/c something acts on the key
    end
    kinetics *= "\n\tTX_rate_vector[\"$key\"] = TX_control_term[\"$key\"]*kcatTX*RNAPconc*geneConc/(saturationTX+geneConc)" *
                "\n\tTL_rate_vector[\"$tmp_protein_string\"] = kcatTL*RIBOconc*X[$(all_species_dict[key])]/(saturationTL+X[$(all_species_dict[key])])"
  end

  kinetics *= "\n\n\treturn (rnx_rate_vector, TX_rate_vector, TL_rate_vector)" *
              "\nend"
  return (kinetics, MonodAffinityConstant_String_Array, W_string_array, disassociation_const_string_array)
end

function jl_build_data_dictionary_buffer(host_type::AbstractString, all_species_array::Array,
  all_species2index_dict::Dict,
  rnx_species_array::Array, all_rnx_list::Array, all_txtl_dict::Dict,
  Monod_affinity_constant_array::Array, W_string_array::Array,
  disassociation_const_array::Array, mRNA_species_array::Array, protein_species_array::Array)
  # data dictionary
  buffer = "function generate_model_parameters_dictionary()"
  # stoichiometry
  buffer *= "\n\t# load stoichiometry"
  buffer *= "\n\tstoichiometric_matrix = readdlm(\"./stoichiometry.dat\")"

  No_species = length(all_species_array)
  # initial concentration
  buffer *= "\n\n\t# initial species concentration" *
    "\n\tinitial_condition = [\n\t\t"
  buffer *= join(["0.0,  # $(i)  $(all_species_array[i])" for i = 1:(No_species - 1)], "\n\t\t")
  buffer *= "\n\t\t0.0   # $(No_species) $(all_species_array[No_species])"
  buffer *= "\n\t]"

  # all species dictionary
  buffer *= "\n\n\t# all species dictionary"
  buffer *= "\n\tall_species_dict = Dict{String, Int64}()  # String --> Int"
  #all_species_dict = Dict{String, Int64}()  # create all species dict: string --> int64
  for (st, index) in all_species2index_dict
    buffer *= "\n\tall_species_dict[\"$st\"] = $index"
    #all_species_dict[st] = index
  end
  # all species dictionary_reversed
  buffer *= "\n\n\t# all species dictionary_reversed"
  buffer *= "\n\tall_species_reversed_dict = Dict{Int64, String}()  # Int --> String"
  buffer *= "\n\tfor (key,val) in all_species_dict"
  buffer *= "\n\t\tall_species_reversed_dict[val] = key"
  buffer *= "\n\tend"
  # all rnx species array
  buffer *= "\n\n\t# all rnx species array"
  buffer *= "\n\trnx_species_array::Array{String,1} = ["
  for st in rnx_species_array
    buffer *= "\n\t\t\"$st\","
  end
  buffer = chop(buffer)*"\n\t]"
  # all mRNA species array
  buffer *= "\n\n\t# all mRNA species array"
  buffer *= "\n\tmRNA_species_array::Array{String,1} = ["
  for st in mRNA_species_array
    buffer *= "\n\t\t\"$st\","
  end
  buffer = chop(buffer)*"\n\t]"
  # all protein species array
  buffer *= "\n\n\t# all protein species array"
  buffer *= "\n\tprotein_species_array::Array{String,1} = ["
  for st in protein_species_array
    buffer *= "\n\t\t\"$st\","
  end
  buffer = chop(buffer)*"\n\t]"

  # signaling reaction kinetic constants, reaction name: reactant(s)<rnx type:catalyst(s)>product(s)
  buffer *= "\n\n\t# signaling reaction kinetic constants, reaction name: reactant(s)<rnx type:catalyst(s)>product(s)"
  buffer *= "\n\tkcat_signaling = ones($(length(all_rnx_list)))  # kcat[#reaction]: reaction name"
  for (id, rnx) in enumerate(all_rnx_list)
    buffer*= "\n\tkcat_signaling[$id] = 1e-4  # kcat: $(rnx.rnxName)"
  end
  # Monod affinity constant in signaling reaction, nomenclature: MonodK_#reaction_Substrate
  buffer *= "\n\n\t# Monod affinity constant in signaling reaction, nomenclature: MonodK_#reaction_Substrate"
  buffer *= "\n\tMonodAffinityConstantDict = Dict{String, Float64}()"
  for MK in Monod_affinity_constant_array
    # buffer *= "\n\t$MK = 1.0"
    # buffer *= "\n\tMonodAffinityConstant_dict[\"$MK\"] = $MK"
    buffer *= "\n\tMonodAffinityConstantDict[\"$MK\"] = 0.1"
  end

  # W_value in transcription control term, nomenclature: W_targetmRNA_actor
  buffer *= "\n\n\n\t# W_value in transcription control term, nomenclature: W_targetmRNA_actor"
  buffer *= "\n\tW_value_dict = Dict{String, Float64}()"
  for WS in W_string_array
    buffer *= "\n\tW_value_dict[\"$WS\"] = 1.0"
  end
  # disassociation constants in transfer function, nomenclature: KD_speciesName
  buffer *= "\n\n\t# disassociation constants in transfer function, nomenclature: KD_speciesName"
  buffer *= "\n\ttransferFunctionDisassociationConstantDict = Dict{String, Float64}()"
  for DC in disassociation_const_array
    buffer *= "\n\ttransferFunctionDisassociationConstantDict[\"$DC\"] = 1.0"
  end
  # transcription specific correction factor
  buffer *= "\n\n\t# transcription specific correction factor"
  buffer *= "\n\ttranscriptionSpecificCorrectionFactorDict = Dict{String, Float64}()"
  for mrna in mRNA_species_array
    buffer *= "\n\ttranscriptionSpecificCorrectionFactorDict[\"$mrna\"] = 1.0  # $mrna"
  end
  # background gene expression control term
  buffer *= "\n\n\t# background gene expression control term"
  buffer *= "\n\tbackgroundGeneExpressionControlTermDict = Dict{String, Float64}()"
  for mrna in mRNA_species_array
    buffer *= "\n\tbackgroundGeneExpressionControlTermDict[\"$mrna\"] = 0.001  # $mrna"
  end
  # translation specific correction factor
  buffer *= "\n\n\t# translation specific correction factor"
  buffer *= "\n\ttranslationSpecificCorrectionFactor = Dict{String, Float64}()"
  for pro in protein_species_array
    buffer *= "\n\ttranslationSpecificCorrectionFactor[\"$pro\"] = 1.0  # $pro"
  end
  # cooperativity number in transfer function
  buffer *= "\n\n\t# cooperativity number in transfer function"
  buffer *= "\n\tcooperativity = 1"
  #buffer *= "\n\n\tbackground_mRNA_synthesis_rate_vector = 0.01*ones($length_TXTL)"
  # load txtl constants buffer
  buffer *= "\n\n"
  if host_type == "bacteria"
    buffer *= jl_include_constants_from_literature(
              joinpath(Base.@__DIR__, "txtl_constants_ecoli.jl"),"\n\t")
  else
    buffer *= jl_include_constants_from_literature(
              joinpath(Base.@__DIR__, "txtl_constants_hl60.jl"), "\n\t")
  end

  #---------------------------------
  # put all stuff in a Dictionary
  buffer *= "\n\n\t####################################"
  buffer *= "\n\t# put all stuff in a Dictionary"
  buffer *= "\n\tdataDictionary = Dict{String, Any}()"
  buffer *= "\n\tdataDictionary[\"stoichiometric_matrix\"] = stoichiometric_matrix"
  buffer *= "\n\tdataDictionary[\"all_species_dict\"] = all_species_dict"
  buffer *= "\n\tdataDictionary[\"all_species_reversed_dict\"] = all_species_reversed_dict"
  buffer *= "\n\tdataDictionary[\"rnx_species_array\"] = rnx_species_array"
  buffer *= "\n\tdataDictionary[\"mRNA_species_array\"] = mRNA_species_array"
  buffer *= "\n\tdataDictionary[\"protein_species_array\"] = protein_species_array"
  buffer *= "\n\tdataDictionary[\"initial_condition\"] = initial_condition"

  buffer *= "\n\tdataDictionary[\"kcat_signaling\"] = kcat_signaling"
  buffer *= "\n\tdataDictionary[\"MonodAffinityConstantDict\"] = MonodAffinityConstantDict"
  buffer *= "\n\tdataDictionary[\"W_value_dict\"] = W_value_dict"
  buffer *= "\n\tdataDictionary[\"transferFunctionDisassociationConstantDict\"] = transferFunctionDisassociationConstantDict"
  buffer *= "\n\tdataDictionary[\"transcriptionSpecificCorrectionFactorDict\"] = transcriptionSpecificCorrectionFactorDict"
  buffer *= "\n\tdataDictionary[\"backgroundGeneExpressionControlTermDict\"] = backgroundGeneExpressionControlTermDict"
  buffer *= "\n\tdataDictionary[\"translationSpecificCorrectionFactor\"] = translationSpecificCorrectionFactor"
  buffer *= "\n\tdataDictionary[\"cooperativity\"] = cooperativity"
  buffer *= "\n\tdataDictionary[\"RNAPConcentration\"] = rnapII_concentration"
  buffer *= "\n\tdataDictionary[\"RIBOConcentration\"] = ribosome_concentration"
  buffer *= "\n\tdataDictionary[\"degradationConstantmRNA\"] = degradation_constant_mRNA"
  buffer *= "\n\tdataDictionary[\"degradationConstantProtein\"] = degradation_constant_protein"
  buffer *= "\n\tdataDictionary[\"kcatTranscription\"] = kcat_transcription"
  buffer *= "\n\tdataDictionary[\"kcatTranslation\"] = kcat_translation"
  buffer *= "\n\tdataDictionary[\"avgGeneConcentration\"] = avg_gene_concentration"
  buffer *= "\n\tdataDictionary[\"transcriptionSaturationConstant\"] = saturation_transcription"
  buffer *= "\n\tdataDictionary[\"translationSaturationConstant\"] = saturation_translation"
  buffer *= "\n\tdataDictionary[\"specificGrowthRate\"] = maximum_specific_growth_rate - death_rate_constant"

  buffer *= "\n\n\treturn dataDictionary"
  buffer *= "\nend"

  return buffer
end

function jl_build_simulation_buffer(NoExtracellularSpecies::Int64)
  # buffer = build_copyright_header_buffer()
  buffer = "\n# set up ODE, get the derivatives
function Balances(t,y,dataDictionary)
  # correct for negatives
  idx_small = findall(y.<0)
  y[idx_small] .= 0.0

  # load data from data dictionary
  rnx_species = dataDictionary[\"rnx_species_array\"]
  mRNA_species = dataDictionary[\"mRNA_species_array\"]
  protein_species = dataDictionary[\"protein_species_array\"]
  stoichiometry = dataDictionary[\"stoichiometric_matrix\"]
  all_species_dict = dataDictionary[\"all_species_dict\"]  # String --> number
  specific_growth_rate = dataDictionary[\"specificGrowthRate\"]
  transcription_correction = dataDictionary[\"transcriptionSpecificCorrectionFactorDict\"]
  translation_correction = dataDictionary[\"translationSpecificCorrectionFactor\"]
  degradation_constant_mRNA = dataDictionary[\"degradationConstantmRNA\"]
  degradation_constant_protein = dataDictionary[\"degradationConstantProtein\"]

  # kinetics
  (rnx_rate, transcription_rate, translation_rate) = calculate_kinetics(y, dataDictionary)

  # calculate the derivatives
  dydt = zeros(length(all_species_dict))

  # signaling reaction network
  X2 = zeros(length(rnx_species)-$(NoExtracellularSpecies))  # species inside the cell
  for (id, st) in enumerate(rnx_species[$(NoExtracellularSpecies)+1:end])  # get X2 value from y
    X2[id] = y[all_species_dict[st]]
  end
  dX1 = stoichiometry[1:$(NoExtracellularSpecies), :]*rnx_rate  # extracellular metabolites
  dX2 = stoichiometry[$(NoExtracellularSpecies)+1:end, :]*rnx_rate - specific_growth_rate*X2  # INTRA
  dX = [dX1; dX2]
  for (id, st) in enumerate(rnx_species)  # return derivatives back to dydt
    dydt[all_species_dict[st]] = dX[id]
  end
  dydt[1] = specific_growth_rate*y[1]  # for Biomass

  # TXTL network
  for i = 1:length(mRNA_species)
  # seems easy to add additional terms to account for the signaling effect
    mRNA_id = all_species_dict[mRNA_species[i]]
    p_id = all_species_dict[protein_species[i]]
    dydt[mRNA_id] = transcription_rate[mRNA_species[i]] - (specific_growth_rate + transcription_correction[mRNA_species[i]]*degradation_constant_mRNA) * y[mRNA_id]
    dydt[p_id] = translation_rate[protein_species[i]] - (specific_growth_rate + translation_correction[protein_species[i]]*degradation_constant_protein) * y[p_id]
  end

  return dydt
end"
  return buffer
end

function jl_build_solveODEBalances_buffer(all_species_array::Array,
  all_species_dict::Dict{String, Int}, mRNA_species::Array, protein_species::Array)

  buffer = "include(\"./include.jl\")" *
    "\n\nfunction SolveBalances(TStart, TStop, TStep, dataDictionary)" *
    "\n\tt = collect(TStart:TStep:TStop)" *
    "\n\tall_species_reversed_dict = dataDictionary[\"all_species_reversed_dict\"]"
  # initial concentration
  buffer *= "\n\t# initial species concentration"
  buffer *= "\n\ty0 = dataDictionary[\"initial_condition\"]"
  # run simulation
  buffer *= "\n\t# call ODE solver"
  buffer *= "\n\tf(t, y) = Balances(t, y, dataDictionary)" *
    "\n\tt, y = ode23s(f, y0, t; points=:specified)"
  # data transfer
  buffer *= "\n\n\t# data transfer" *
    "\n\trow = length(t)" *
    "\n\tcol = length(y0)" *
    "\n\tY = zeros(row, col)" *
    "\n\tforeach(x->(Y[x, :] = y[x]), collect(1:row))"
  # plotting
  buffer *= "\n\t# plot results" *
    "\n\tplt.figure(\"simulation \$TStart to \$TStop\")" *
    "\n\tsubplt_col = ceil(Int, col/4)" *
    "\n\tfor i = 1:col" *
    "\n\t\tplt.subplot(4, subplt_col, i)" *
    "\n\t\tplt.plot(t, Y[:,i])" *
    "\n\t\tplt.title(all_species_reversed_dict[i])" *
    "\n\tend" *
    "\n\tplt.show()"
  buffer *= "\n\n\treturn (t, Y)" *
    "\nend"

  # add a test
  buffer *= "\n\n\n"
  buffer *= "# # uncomment to test"
  buffer *= "\n# dataDictionary = generate_model_parameters_dictionary()"
  buffer *= "\n# SolveBalances(0, 10, 0.1, dataDictionary)"

  return buffer
end

# function build_copyright_header_buffer()
#   current_year = string(Dates.year(now()))
#   # Get comment data from
#   buffer = ""
#   buffer *=
# "# ----------------------------------------------------------------------------------- #
# # Copyright (c) $(current_year) Varnerlab
# # Robert Frederick Smith School of Chemical and Biomolecular Engineering
# # Cornell University, Ithaca NY 14850
# #
# # Permission is hereby granted, free of charge, to any person obtaining a copy
# # of this software and associated documentation files (the \"Software\"), to deal
# # in the Software without restriction, including without limitation the rights
# # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# # copies of the Software, and to permit persons to whom the Software is
# # furnished to do so, subject to the following conditions:
# #
# # The above copyright notice and this permission notice shall be included in
# # all copies or substantial portions of the Software.
# #
# # THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# # THE SOFTWARE.
# # ----------------------------------------------------------------------------------- #"
#
#   return buffer
# end

# FBA data dictionary generation
function jl_generate_FBA_data_dictionary(all_rnx_list::Array,
  rnx_species_array::Array, extra_species_num::Int)
  secrete_id = Set{Int64}()  # For initialize the coefficient array
  # data dictionary
  buffer = "function DataDictionary(time_start,time_stop,time_step)"
  # stoichiometry
  buffer *= "\n\n\t# load stoichiometry"
  buffer *= "\n\tstoichiometric_matrix = readdlm(\"./stoichiometry.dat\")"
  # Setup default flux bounds array
  buffer *= "\n\t# Setup default flux bounds array"
  buffer *= "\n\tdefault_bounds_array = ["
  for (id, rnx) in enumerate(all_rnx_list)
    buffer *= "\n\t\t0 100.0; # $(id) $(rnx.rnxName)"
  end
  buffer *= "\n\t]"
  # Setup default species bounds array
  buffer *= "\n\n\t# Setup default species bounds array"
  buffer *= "\n\tspecies_bounds_array = ["
  for (id, sp) in enumerate(rnx_species_array)
    if id > extra_species_num
      buffer *= "\n\t\t 0.0 0.0; # $(id) $(sp)"
    else
      buffer *= "\n\t\t-1.0 1.0; # $(id) $(sp)"
    end
  end
  buffer *= "\n\t]"
  # Setup the objective coefficient array
  buffer *= "\n\n\t# Setup the objective coefficient array"
  buffer *= "\n\tobjective_coefficient_array = ["
  for (id, rnx) in enumerate(all_rnx_list)
    buffer *= "\n\t\t0.0; # $(id) $(rnx.rnxName)"
    if rnx.rnxType == "secrete"
        push!(secrete_id, id)
    end
  end
  buffer *= "\n\t]"
  # Min/Max flag - default is minimum -
  buffer *= "\n\n\t# Min/Max flag - default is minimum -"
  buffer *= "\n\tis_minimum_flag = true"
  # List of reation strings - used to write flux report
  buffer *= "\n\n\t# List of reation strings - used to write flux report"
  buffer *= "\n\tlist_of_reaction_strings = ["
  for rnx in all_rnx_list
    buffer *= "\n\t\t\"$(rnx.rnxName)\""
  end
  buffer *= "\n\t]"
  # List of metabolite strings - used to write flux report
  buffer *= "\n\n\t# List of metabolite strings - used to write flux report"
  buffer *= "\n\tlist_of_metabolite_symbols = ["
  for sp in rnx_species_array
    buffer *= "\n\t\t\"$(sp)\""
  end
  buffer *= "\n\t]"
  # return block -
  buffer *= "\n\n"
  buffer *= "\t# =============================== DO NOT EDIT BELOW THIS LINE ============================== #\n"
  buffer *= "\tdata_dictionary = Dict{AbstractString,Any}()\n"
  buffer *= "\tdata_dictionary[\"stoichiometric_matrix\"] = stoichiometric_matrix\n"
  buffer *= "\tdata_dictionary[\"objective_coefficient_array\"] = objective_coefficient_array\n"
  buffer *= "\tdata_dictionary[\"default_flux_bounds_array\"] = default_bounds_array;\n"
  buffer *= "\tdata_dictionary[\"species_bounds_array\"] = species_bounds_array\n"
  buffer *= "\tdata_dictionary[\"list_of_reaction_strings\"] = list_of_reaction_strings\n"
  buffer *= "\tdata_dictionary[\"list_of_metabolite_symbols\"] = list_of_metabolite_symbols\n"
  buffer *= "\tdata_dictionary[\"is_minimum_flag\"] = is_minimum_flag\n"
  buffer *= "\t# =============================== DO NOT EDIT ABOVE THIS LINE ============================== #\n"
  buffer *= "\treturn data_dictionary\n"
  buffer *= "end\n"

  # maximize_product_dictionary
  buffer *= "\n\n\n# maximize_product_dictionary"
  buffer *= "\nfunction maximize_product_dictionary(time_start,time_stop,time_step)"
  # load the data dictionary -
  buffer *= "\n\t# load the data dictionary -"
  buffer *= "\n\tdata_dictionary = DataDictionary(time_start,time_stop,time_step)"
  # Modify the data dictionary -
  buffer *= "\n\t# Modify the data dictionary -"
  buffer *= "\n\tobjective_coefficient_array = data_dictionary[\"objective_coefficient_array\"]"
  secrete_id = sort(collect(secrete_id))
  for id in secrete_id
      buffer *= "\n\tobjective_coefficient_array[$(id)] = -1"
  end
  buffer *= "\n\t# return -"
  buffer *= "\n\treturn data_dictionary"
  buffer *= "\nend"

  return buffer
end


text:
#########################################################################
#
#   Derivation rule and type equivalence declaration functions / macros
#
#########################################################################

drules = Dict() # holds the derivation expressions

#### function derivation rules declaration functions/macros

# macro version
macro deriv_rule(func::Expr, dv::Symbol, diff)
    # func = :(  colon(x,y) ) ; dv = :x ; diff = 0.
    emod = current_module()

    ss = Symbol[]
    ts = Type[]
    for e in func.args[2:end]
        if isa(e, Symbol)
            push!(ss, e)
            push!(ts, Any)

        elseif isa(e, Expr) && e.head== :(::)
            push!(ss, e.args[1])
            push!(ts, emod.eval(e.args[2]))

        else
            error("[deriv_rule] cannot parse $e")
        end
    end

    deriv_rule(emod.eval(func.args[1]), collect(zip(ss, ts)), dv, diff)
end

function deriv_rule(func::Union{Function, Type},
    args::Vector, dv::Symbol, diff::Union{Expr, Symbol, Real})
    emod = current_module()

    sig = Tuple{ Type[ e[2] for e in args ]... }

    ss  = Symbol[ e[1] for e in args ]

    index = findfirst(dv .== ss)
    (index == 0) && error("[deriv_rule] cannot find $dv in function arguments")

    # non generic functions are matched by their name
    fidx = isa(func, Function) && isbuiltin(func) ?
                  builtin_name(func) : func

    haskey(drules, (fidx, index)) || (drules[(fidx, index)] = Dict())

    g = tograph(diff, emod)  # make the graph
    push!(ss, :ds)

    drules[(fidx, index)][sig] = (g, ss)
    nothing
end

isbuiltin(f) = VERSION >= v"0.5.0-" ? isa(f, Core.Builtin) : isa(f.env, Symbol)
builtin_name(f) = VERSION >= v"0.5.0-" ? Base.function_name(f) : f.env

function hasrule(f, pos)
    haskey(drules, (f, pos)) && return true
    # non generic functions are matched by their name
    if isa(f, Function) && isbuiltin(f)
        fname = builtin_name(f)
        return haskey(drules, (fname, pos))
    end
    false
end

function getrule(f, pos)
    if isa(f, Function) && isbuiltin(f) # non generic functions are matched by their name
        fname = builtin_name(f)
        haskey(drules, (fname, pos)) && return drules[(fname, pos)]
    else
        haskey(drules, (f, pos)) && return drules[(f, pos)]
    end
    error("no derivation rule for $(f) at arg #$(pos)")
end

#### Type tuple matching  (naive multiple dispatch)
function tmatch(sig, keys)
    keys2 = filter(k -> length(k.parameters) == length(sig.parameters), keys)
    sort!(keys2, lt=issubtype)
    idx = findfirst(s -> sig <: s, keys2)
    return idx==0 ? nothing : keys2[idx]
end