text:
###########################################################################
#
# This file contains following functions
#   glmxls - GLM regression models
#   bivariatexls - for two-way tables for discrete or continuous variables
#   univariatexls - for continuous variables
#   dfxls - for exporting a dataframe
#
###########################################################################

"""
"""
function hltest(glmout,q = 10)
    d = DataFrame(yhat = predict(glmout),y = glmout.model.rr.y)
    d[:group] = xtile(d[:yhat],nq = q)
    df3 = by(d,:group) do subdf
        n = size(subdf,1)
        o1 = sum(subdf[:y])
        e1 = sum(subdf[:yhat])
        DataFrame(
            observed1 = o1,
            observed0 = n - o1,
            expected1 = e1,
            expected0 = n - e1
        )
    end

    hlstat = sum((df3[:observed1] .- df3[:expected1]).^2 ./ df3[:expected1]
        .+ (df3[:observed0] .- df3[:expected0]).^2 ./ df3[:expected0])

    dof = q - 2

    pval = Distributions.ccdf(Distributions.Chisq(dof),hlstat)

    return (hlstat, dof, pval)
end

function nagelkerke(glmout)
    return StatsBase.r2(glmout,:Nagelkerke)
    # L = loglikelihood(glmout)
    # L0 = loglikelihood(glm(@eval(@formula($(glmout.mf.terms.eterms[1]) ~ 1)),glmout.mf.df,Bernoulli(),LogitLink()))
    # pow = 2/nobs(glmout)
    # return (1 - exp(-pow*(L - L0))) / (1 - exp(pow*L0))
end

function macfadden(glmout)
    return StatsBase.r2(glmout,:McFadden)
    # L = loglikelihood(glmout)
    # L0 = loglikelihood(glm(@eval(@formula($(glmout.mf.terms.eterms[1]) ~ 1)),glmout.mf.df,Bernoulli(),LogitLink()))
    # return 1 - L/L0
end

"""
    glmxls(glmout::DataFrames.DataFrameRegressionModel, workbook::PyObject, worksheet::AbstractString; labels::Union{Nothing,Label}=nothing,eform=false,ci=true, row = 0, col =0)

Outputs a GLM regression table to an excel spreadsheet.
To use this function, `PyCall` is required with a working version python and
a python package called `xlsxwriter` installed. If a label is found for a variable
or a value of a variable in a `Label`, the label will be output. Options are:

- `glmout`: returned value from a GLM regression model
- `workbook`: a returned value from xlsxwriter.Workbook() function (see an example below)
- `worksheet`: a string for the worksheet name
- `labels`: an option to specify a `label` dictionary (see an example below)
- `eform`: use `eform = true` to get exponentiated estimates, standard errors, or 95% confidence intervals
- `ci`: use `ci = true` (default) to get 95% confidence intervals. `ci = false` will produce standard errors and Z values instead.
- `row`: specify the row of the workbook to start the output table (default = 0 (for row 1))
- `col`: specify the column of the workbook to start the output table (default = 0 (for column A))

# Example 1
This example is useful when one wants to append a worksheet to an existing workbook.
It is responsibility of the user to open a workbook before the function call and close it
to actually create the physical file by close the workbook.

```
julia> using PyCall

julia> @pyimport xlsxwriter

julia> wb = xlsxwriter.Workbook("test_workbook.xlsx")
PyObject <xlsxwriter.workbook.Workbook object at 0x000000002A628E80>

julia> glmxls(ols1,wb,"OLS1",labels = label)

julia> bivairatexls(df,:incomecat,[:age,:race,:male,:bmicat],wb,"Bivariate",labels = label)

Julia> wb.close()
```

# Example 2
Alternatively, one can create a spreadsheet file directly. `PyCall` or `@pyimport`
does not need to be called before the function.

```
julia> glmxls(ols1,"test_workbook.xlsx","OLS1",labels = label)
```

"""
function glmxls(glmout,wbook::PyObject,wsheet::AbstractString;
    labels::Union{Nothing,Label} = nothing,
    eform::Bool = false, ci = true, row = 0, col = 0)

    if (typeof(glmout) <: StatsModels.RegressionModel) == false
        error("This is not a regression model output.")
    end

    if isa(glmout.model,GeneralizedLinearModel)
        distrib = glmout.model.rr.d
        linkfun = Link(glmout.model.rr)
    else
        distrib = nothing
        linkfun = nothing
    end

    # create a worksheet
    t = wbook.add_worksheet(wsheet)

    # attach formats to the workbook
    formats = ExcelTables.attach_formats(wbook)

    # starting location in the worksheet
    r = row
    c = col

    # set column widths
    t.set_column(c,c,40)
    t.set_column(c+1,c+4,9)

    # headings ---------------------------------------------------------
    # if ci == true, Estimate (95% CI) P-Value
    # if eform == true, Estimate is OR for logit, IRR for poisson
    otype = "Estimate"
    if eform == true
        otype = Stella.coeflab(distrib,linkfun)
    end

    t.write_string(r,c,"Variable",formats[:heading])
    if ci == true
        t.merge_range(r,c+1,r,c+3,string(otype," (95% CI)"),formats[:heading])
        t.write_string(r,c+4,"P-Value",formats[:heading])
    else
        t.write_string(r,c+1,otype,formats[:heading])
        t.write_string(r,c+2,"SE",formats[:heading])
        t.write_string(r,c+3,"Z Value",formats[:heading])
        t.write_string(r,c+4,"P-Value",formats[:heading])
    end

    #------------------------------------------------------------------
    r += 1
    c = col

    # go through each variable and construct variable name and value label arrays
    tdata = coeftable(glmout)
    nrows = length(tdata.rownms)
    varname = Vector{String}(undef,nrows)
    vals = Vector{String}(undef,nrows)
    nlev = zeros(Int,nrows)
    #nord = zeros(Int,nrow)
    for i = 1:nrows
    	# variable name
        # parse varname to separate variable name from value
        if occursin(":",tdata.rownms[i])
            (varname[i],vals[i]) = split(tdata.rownms[i],": ")
        elseif occursin(" - ",tdata.rownms[i])
            (varname[i],vals[i]) = split(tdata.rownms[i]," - ")
        else
            varname[i] = tdata.rownms[i]
            vals[i] = ""
        end

        # use labels if exist
        if labels != nothing
            if vals[i] != ""
                valn = vals[i] == "true" ? 1 : parse(Int,vals[i])
                vals[i] = vallab(labels,Symbol(varname[i]),valn)
            end
            varname[i] = varlab(labels,Symbol(varname[i]))
        end
    end
    for i = 1:nrows
        nlev[i] = countlev(varname[i],varname)
    end

    # write table
    tconfint = confint(glmout)
    lastvarname = ""

    for i = 1:nrows
        if otype == "OR" && varname[i] == "(Intercept)"
            continue
        end
        if varname[i] != lastvarname
            # output cell boundaries only and go to the next line
            if nlev[i] > 1
                t.write_string(r,c,varname[i],formats[:heading_left])

                if ci == true
                    t.write(r,c+1,"",formats[:empty_right])
                    t.write(r,c+2,"",formats[:empty_both])
                    t.write(r,c+3,"",formats[:empty_left])
                    t.write(r,c+4,"",formats[:p_fmt])
                else
                    t.write(r,c+1,"",formats[:empty_border])
                    t.write(r,c+2,"",formats[:empty_border])
                    t.write(r,c+3,"",formats[:empty_border])
                    t.write(r,c+4,"",formats[:p_fmt])
                end
                r += 1
                t.write_string(r,c,vals[i],formats[:varname_1indent])

            else
                if vals[i] != "" && vals[i] != "Yes"
                    t.write_string(r,c,string(varname[i]," - ",vals[i]),formats[:heading_left])
                else
                    t.write_string(r,c,varname[i],formats[:heading_left])
                end
            end
        else
            t.write_string(r,c,vals[i],formats[:varname_1indent])
        end

    	# estimates
        if eform == true
    	    t.write(r,c+1,exp(tdata.cols[1][i]),formats[:or_fmt])
        else
            t.write(r,c+1,tdata.cols[1][i],formats[:or_fmt])
        end

        if ci == true

            if eform == true
            	# 95% CI Lower
            	t.write(r,c+2,exp(tconfint[i,1]),formats[:cilb_fmt])

            	# 95% CI Upper
            	t.write(r,c+3,exp(tconfint[i,2]),formats[:ciub_fmt])
            else
                # 95% CI Lower
            	t.write(r,c+2,tconfint[i,1],formats[:cilb_fmt])

            	# 95% CI Upper
            	t.write(r,c+3,tconfint[i,2],formats[:ciub_fmt])
            end
        else
            # SE
            if eform == true
        	    t.write(r,c+2,exp(tdata.cols[1][i])*tdata.cols[2][i],formats[:or_fmt])
            else
                t.write(r,c+2,tdata.cols[1][i],formats[:or_fmt])
            end

            # Z value
            t.write(r,c+3,tdata.cols[3][i],formats[:or_fmt])

        end

        # P-Value
        t.write(r,c+4,tdata.cols[4][i] < 0.001 ? "< 0.001" : tdata.cols[4][i] ,formats[:p_fmt])

        lastvarname = varname[i]

        # update row
        r += 1
    end

    # Write model characteristics and goodness of fit statistics
    c = col

    # N
    t.write(r,c,"N",formats[:model_name])
    t.merge_range(r,c+1,r,c+4,nobs(glmout),formats[:n_fmt_center])

    # degress of freedom
    # r += 1
    # t.write(r,c,"DF",formats[:model_name])
    # t.merge_range(r,c+1,r,c+4,dof(glmout),formats[:n_fmt_center])

    # R² or pseudo R²
    r += 1
    if isa(linkfun,LogitLink)
        t.write(r,c,"Pseudo R² (MacFadden)",formats[:model_name])
        t.merge_range(r,c+1,r,c+4,macfadden(glmout),formats[:p_fmt_center])
        t.write(r+1,c,"Pseudo R² (Nagelkerke)",formats[:model_name])
        t.merge_range(r+1,c+1,r+1,c+4,nagelkerke(glmout),formats[:p_fmt_center])

        # -2 log-likelihood
        t.write(r+2,c,"-2 Log-Likelihood",formats[:model_name])
        t.merge_range(r+2,c+1,r+2,c+4,deviance(glmout),formats[:p_fmt_center])

        # Hosmer-Lemeshow GOF test
        t.write(r+3,c,"Hosmer-Lemeshow Chisq Test (df), p-value",formats[:model_name])
        hl = hltest(glmout)
        t.merge_range(r+3,c+1,r+3,c+4,string(round(hl[1],digits=4)," (",hl[2],"); p = ",round(hl[3],digits=4)),formats[:p_fmt_center])

        # ROC (c-statistic)
        t.write(r+4,c,"Area under the ROC Curve",formats[:model_name])
        roc = auc(roc(glmout.model.rr.y,predict(glmout)))
        n1 = sum(glmout.model.rr.y) # number of positive responses
        n2 = nobs(glmout) - n1
        q1 = roc / (2 - roc)
        q2 = (2*roc^2) / (1 + roc)
        rocse = sqrt((roc*(1-roc) + (n1-1)*(q1 - roc^2) + (n2 - 1)*(q2 - roc^2)) / (n1*n2))
        t.merge_range(r+4,c+1,r+4,c+4,string(round(roc,digits=4)," (95% CI, ", round(roc - 1.96*rocse,digits=4), " - ", round(roc + 1.96*rocse,digits=4),")"),formats[:p_fmt_center])

        r += 5
    else
        t.write(r,c,"R²",formats[:model_name])
        t.merge_range(r,c+1,r,c+4,r2(glmout),formats[:p_fmt_center])
        t.write(r+1,c,"Adjusted R²",formats[:model_name])
        t.merge_range(r+1,c+1,r+1,c+4,adjr2(glmout),formats[:p_fmt_center])

        r += 2
    end

    # AIC & BIC
    t.write(r,c,"AIC",formats[:model_name])
    t.merge_range(r,c+1,r,c+4,aic(glmout),formats[:p_fmt_center])

    r += 1
    t.write(r,c,"BIC",formats[:model_name])
    t.merge_range(r,c+1,r,c+4,bic(glmout),formats[:p_fmt_center])
end
function glmxls(glmout,
    wbook::AbstractString,
    wsheet::AbstractString;
    labels::Union{Nothing,Label} = nothing,
    eform::Bool = false,
    ci = true,
    row = 0,
    col = 0)

    xlsxwriter = pyimport("xlsxwriter")

    glmxls(glmout,xlsxwriter.Workbook(wbook),wsheet,labels=labels,eform=eform,ci=ci,row=row,col=col)
end

"""
    bivariatexls(df::DataFrame,colvar::Symbol,rowvars::Vector{Symbol},workbook::PyObject,worksheet::AbstractString; labels::Union{Nothing,Label}=nothing,row=0,col=0)

 Creates bivariate statistics and appends it in a nice tabular format to an existing workbook.
 To use this function, `PyCall` is required with a working version python and
 a python package called `xlsxwriter` installed.  If a label is found for a variable
 or a value of a variable in a `Label`, the label will be output. Options are:

- `df`: a DataFrame
- `colvar`: a categorical variable whose values will be displayed on the columns
- `rowvars`: a Vector of Symbols for variables to be displayed on the rows. Both continuous and categorical variables are allowed.
    For continuous variables, mean and standard deviations will be output and a p-value will be based on an ANOVA test. For categorical variables,
    a r x c table with cell counts and row percentages will be output with a p-value based on a chi-square test.
- `workbook`: a returned value from xlsxwriter.Workbook() function (see an example below)
- `worksheet`: a string for the worksheet name
- `labels`: an option to specify a `Label` object (see an example below)
- `row`: specify the row of the workbook to start the output table (default = 0 (for row 1))
- `col`: specify the column of the workbook to start the output table (default = 0 (for column A))
- `column_percent`: set this to `false` if you want row percentages in the output table (default = true)

### Example 1
This example is useful when one wants to append a worksheet to an existing workbook.
It is responsibility of the user to open a workbook before the function call and close it
to actually create the physical file by close the workbook.

```
julia> using PyCall

julia> xlsxwriter = pyimport("xlsxwriter")

julia> wb = xlsxwriter.Workbook("test_workbook.xlsx")
PyObject <xlsxwriter.workbook.Workbook object at 0x000000002A628E80>

julia> glmxls(ols1,wb,"OLS1",labels = label)

julia> bivairatexls(df,:incomecat,[:age,:race,:male,:bmicat],wb,"Bivariate",labels = label)

Julia> wb.close()
```

### Example 2
Alternatively, one can create a spreadsheet file directly. `PyCall` or `@pyimport`
does not need to be called before the function.

```
julia> bivariatexls(df,:incomecat,[:age,:race,:male,:bmicat],"test_workbook.xlsx","Bivariate",labels = label)
```
"""
function bivariatexls(df::DataFrame,
    colvar::Symbol,
    rowvars::Vector{Symbol},
    wbook::PyObject,
    wsheet::AbstractString;
    wt::Union{Nothing,Symbol} = nothing,
    labels::Union{Nothing,Label} = nothing,
    row::Int = 0,
    col::Int = 0,
    column_percent::Bool = true,
    verbose::Bool = false)

    # check data

    # colvar has to be a CategoricalArray and must have 2 or more categories
    if isa(df[!,colvar], CategoricalArray) == false || length(levels(df[!,colvar])) < 2
        error("`",colvar,"` is not a CategoricalArray or does not have two or more levels")
    end

    # create a worksheet
    t = wbook.add_worksheet(wsheet)

    # attach formats to the workbook
    formats = ExcelTables.attach_formats(wbook)

    # starting row and column
    r = row
    c = col

    # drop NAs in colvar
    df2 = df[completecases(df[!,[colvar]]),:]

    # number of columns
    # column values
    if wt == nothing
        collev = freqtable(df2,colvar,skipmissing=true)
    else
        collev = freqtable(df2,colvar,skipmissing=true,weights=df2[wt])
    end
    nlev = length(collev.array)
    colnms = names(collev,1)
    coltot = sum(collev.array,dims=1)

    # set column widths
    t.set_column(c,c,40)
    t.set_column(c+1,c+(nlev+1)*2+1,9)

    # create heading
    # column variable name
    # It uses three rows
    t.merge_range(r,c,r+2,c,"Variable",formats[:heading])

    # 1st row = variable name
    colvname = String(colvar)
    if labels != nothing
        colvname = varlab(labels,colvar)
    end
    t.merge_range(r,c+1,r,c+(nlev+1)*2+1,colvname,formats[:heading])

    # 2nd and 3rd rows
    r += 1

    t.merge_range(r,1,r,2,"All",formats[:heading])
    t.write_string(r+1,1,"N",formats[:n_fmt_right])
    t.write_string(r+1,2,"(%)",formats[:pct_fmt_parens])

    c = 3
    for i = 1:nlev

        # value label
        vals = string(colnms[i])
        if labels != nothing
            vals = vallab(labels,colvar,colnms[i])
        end

        t.merge_range(r,c+(i-1)*2,r,c+(i-1)*2+1,vals,formats[:heading])
        t.write_string(r+1,c+(i-1)*2,"N",formats[:n_fmt_right])
        t.write_string(r+1,c+(i-1)*2+1,"(%)",formats[:pct_fmt_parens])
    end

    # P-value
    t.merge_range(r,c+nlev*2,r+1,c+nlev*2,"P-Value",formats[:heading])

    # total
    c = 0
    r += 2
    t.write_string(r,c,"All, n (Row %)",formats[:model_name])
    if wt == nothing
        x = freqtable(df2,colvar,skipmissing=true)
    else
        x = freqtable(df2,colvar,skipmissing=true,weights=df2[wt])
    end
    tot = sum(x)
    t.write(r,c+1,tot,formats[:n_fmt_right])
    t.write(r,c+2,1.0,formats[:pct_fmt_parens])
    for i = 1:nlev
        t.write(r,c+i*2+1,x.array[i],formats[:n_fmt_right])
        t.write(r,c+i*2+2,x.array[i]/tot,formats[:pct_fmt_parens])
    end
    t.write(r,c+(nlev+1)*2+1,"",formats[:empty_border])

    # covariates
    c = 0
    r += 1
    for varname in rowvars

        if verbose == true
            println("Processing ",varname)
        end

        # sum of rowvars must be non-zero
        # if sum(length(skipmissing(df[varname]))) == 0
        #     warn("`",varname,"` in `rowvars` is empty.")
        #     continue
        # end

        # print the variable name
        vars = string(varname)
        if labels != nothing
            vars = varlab(labels,varname)
        end

        # determine if varname is categorical or continuous
        if isa(df2[!,varname], CategoricalArray) || eltype(df2[!,varname]) == String

            # categorial
            df3=df2[completecases(df2[:,[varname]]),[varname,colvar]]
            if wt == nothing
                x = freqtable(df3,varname,colvar,skipmissing=true)
            else
                x = freqtable(df3,varname,colvar,skipmissing=true,weights=df3[wt])
            end
            rowval = names(x,1)
            rowtot = sum(x.array,dims=2)
            coltot = sum(x.array,dims=1)

            # variable name
            # if there only two levels and one of the values is 1 or true
            # and the other values is 0 or false,
            # just output the frequency and percentage of the 1/true row

            # variable name
            t.write_string(r,c,vars,formats[:model_name])

            # two levels with [0,1] or [false,true]
            if length(rowval) <= 2 && rowval in ([1],[true],["Yes"],[0,1],[false,true],["No","Yes"])

                nrow = length(rowval)
                # row total
                t.write(r,c+1,rowtot[nrow],formats[:n_fmt_right])
                t.write(r,c+2,rowtot[nrow]/tot,formats[:pct_fmt_parens])

                for j = 1:nlev
                    t.write(r,c+j*2+1,x.array[nrow,j],formats[:n_fmt_right])
                    if column_percent
                        t.write(r,c+j*2+2, coltot[j] > 0 ? x.array[nrow,j]/coltot[j] : "",formats[:pct_fmt_parens])
                    elseif rowtot[2] > 0
                        t.write(r,c+j*2+2, rowtot[nrow] > 0 ? x.array[nrow,j]/rowtot[nrow] : "",formats[:pct_fmt_parens])
                    end
                end
                pval = pvalue(ChisqTest(x.array))
                if isnan(pval) || isinf(pval)
                    pval = ""
                elseif pval < 0.001
                    pval = "< 0.001"
                end
                t.write(r,c+(nlev+1)*2+1, pval,formats[:p_fmt])
                r += 1
            else
                for i = 1:nlev+1
                    t.write_string(r,c+(i-1)*2+1,"",formats[:empty_right])
                    t.write_string(r,c+(i-1)*2+2,"",formats[:empty_left])
                end
                t.write_string(r,c+(nlev+1)*2+1,"",formats[:empty_border])

                r += 1
                for i = 1:length(rowval)
                    # row value
                    vals = string(rowval[i])

                    if labels != nothing
                        vals = vallab(labels,varname,rowval[i])
                    end
                    t.write_string(r,c,vals,formats[:varname_1indent])

                    # row total
                    t.write(r,c+1,rowtot[i],formats[:n_fmt_right])
                    t.write(r,c+2,rowtot[i]/tot,formats[:pct_fmt_parens])

                    for j = 1:nlev
                        t.write(r,c+j*2+1,x.array[i,j],formats[:n_fmt_right])
                        if column_percent
                            t.write(r,c+j*2+2,coltot[j] > 0 ? x.array[i,j]/coltot[j] : "",formats[:pct_fmt_parens])
                        else
                            t.write(r,c+j*2+2,rowtot[i] > 0 ? x.array[i,j]/rowtot[i] : "",formats[:pct_fmt_parens])
                        end
                    end
                    # p-value - output only once
                    pval = pvalue(ChisqTest(x.array))

                    if isnan(pval) || isinf(pval)
                        pval = ""
                    elseif pval < 0.001
                        pval = "< 0.001"
                    end
                    if length(rowval) == 1
                        t.write(r,c+(nlev+1)*2,pval,formats[:p_fmt])
                    elseif i == 1
                        t.merge_range(r,c+(nlev+1)*2+1,r+length(rowval)-1,c+(nlev+1)*2+1,pval,formats[:p_fmt])
                    end
                    r += 1
                end
            end
        else
            # continuous variable
            df3=df2[completecases(df2[!,[varname]]),[varname,colvar]]
            y = tabstat(df3,varname,colvar) #,wt=df3[wt])

            # variable name
            t.write_string(r,c,string(vars,", mean (SD)"),formats[:model_name])

            # All
            t.write(r,c+1,mean(df3[!,varname]),formats[:f_fmt_right])
            t.write(r,c+2,std(df3[!,varname]),formats[:f_fmt_left_parens])

            # colvar levels
            for i = 1:nlev
                if i <= size(y,1) && y[i,:N] > 1
                    t.write(r,c+i*2+1,y[i,:mean],formats[:f_fmt_right])
                    t.write(r,c+i*2+2,y[i,:sd],formats[:f_fmt_left_parens])
                else
                    t.write(r,c+i*2+1,"",formats[:f_fmt_right])
                    t.write(r,c+i*2+2,"",formats[:f_fmt_left_parens])
                end
            end
            if size(y,1) > 1
                pval = anova(df3,varname,colvar).p[1]
                if isnan(pval) || isinf(pval)
                    pval = ""
                elseif pval < 0.001
                    pval = "< 0.001"
                end
                t.write(r,c+(nlev+1)*2+1,pval,formats[:p_fmt])
            else
                t.write(r,c+(nlev+1)*2+1,"",formats[:p_fmt])
            end

            r += 1
        end
    end
end
function bivariatexls(df::DataFrame,
    colvar::Symbol,
    rowvars::Vector{Symbol},
    wbook::AbstractString,
    wsheet::AbstractString;
    labels::Union{Nothing,Label} = nothing,
    row::Int = 0,
    col::Int = 0)

    xlsxwriter = pyimport("xlsxwriter")

    wb = xlsxwriter.Workbook(wbook)

    bivariatexls(df,colvar,rowvars,wb,wsheet,labels=labels,row=row,col=col)

    wb.close()
end

"""
    univariatexls(df::DataFrame,contvars::Vector{Symbol},workbook::PyObject,worksheet::AbstractString; labels::Union{Nothing,Label}=nothing,row=0,col=0)

Creates univariate statistics for a vector of continuous variable and
appends it to an existing workbook.
To use this function, `PyCall` is required with a working version python and
a python package called `xlsxwriter` installed.  If a label is found for a variable
in a `Label` object, the label will be output. Options are:

- `df`: a DataFrame
- `contvars`: a vector of continuous variables
- `workbook`: a returned value from xlsxwriter.Workbook() function (see an example below)
- `worksheet`: a string for the worksheet name
- `labels`: an option to specify a `label` dictionary (see an example below)
- `row`: specify the row of the workbook to start the output table (default = 0 (for row 1))
- `col`: specify the column of the workbook to start the output table (default = 0 (for column A))

# Example 1
This example is useful when one wants to append a worksheet to an existing workbook.
It is responsibility of the user to open a workbook before the function call and close it
to actually create the physical file by close the workbook.

```
julia> using PyCall

julia> @pyimport xlsxwriter

julia> wb = xlsxwriter.Workbook("test_workbook.xlsx")
PyObject <xlsxwriter.workbook.Workbook object at 0x000000002A628E80>

julia> glmxls(ols1,wb,"OLS1",labels = label)

julia> bivairatexls(df,:incomecat,[:age,:race,:male,:bmicat],wb,"Bivariate",labels = label)

julia> univariatexls(df,[:age,:income_amt,:bmi],wb,"Univariate",labels = label)

Julia> wb.close()
```

# Example 2
Alternatively, one can create a spreadsheet file directly. `PyCall` or `@pyimport`
does not need to be called before the function.

```jldoctest
julia> univariatexls(df,[:age,:income_amt,:bmi],"test_workbook.xlsx","Bivariate",labels = label)
```

"""
function univariatexls(df::DataFrame,
    contvars::Vector{Symbol},
    wbook::PyObject,
    wsheet::AbstractString;
    wt::Union{Nothing,Symbol} = nothing,
    labels::Union{Nothing,Label}=nothing,
    row = 0,
    col = 0)

    # create a worksheet
    t = wbook.add_worksheet(wsheet)

    # attach formats to the workbook
    formats = ExcelTables.attach_formats(wbook)

    # starting row and column
    r = row
    c = col

    # column width
    t.set_column(0,0,20)
    t.set_column(1,length(contvars),12)

    # output the row names
    t.write_string(r,c,"Statistic",formats[:heading])
    t.write_string(r+1,c,"N Total",formats[:heading_left])
    t.write_string(r+2,c,"N Miss",formats[:heading_left])
    t.write_string(r+3,c,"N Used",formats[:heading_left])
    t.write_string(r+4,c,"Sum",formats[:heading_left])
    t.write_string(r+5,c,"Mean",formats[:heading_left])
    t.write_string(r+6,c,"SD",formats[:heading_left])
    t.write_string(r+7,c,"Variance",formats[:heading_left])
    t.write_string(r+8,c,"Minimum",formats[:heading_left])
    t.write_string(r+9,c,"P25",formats[:heading_left])
    t.write_string(r+10,c,"Median",formats[:heading_left])
    t.write_string(r+11,c,"P75",formats[:heading_left])
    t.write_string(r+12,c,"Maximum",formats[:heading_left])
    t.write_string(r+13,c,"Skewness",formats[:heading_left])
    t.write_string(r+14,c,"Kurtosis",formats[:heading_left])
    t.write_string(r+15,c,"Smallest",formats[:heading_left])
    t.write_string(r+16,c,"",formats[:heading_left])
    t.write_string(r+17,c,"",formats[:heading_left])
    t.write_string(r+18,c,"",formats[:heading_left])
    t.write_string(r+19,c,"",formats[:heading_left])
    t.write_string(r+20,c,"Largest",formats[:heading_left])
    t.write_string(r+21,c,"",formats[:heading_left])
    t.write_string(r+22,c,"",formats[:heading_left])
    t.write_string(r+23,c,"",formats[:heading_left])
    t.write_string(r+24,c,"",formats[:heading_left])

    col = 1
    for vsym in contvars

        # if there is a label dictionary, pick up the variable label
        varstr = string(vsym)
        if labels != nothing
            varstr = varlab(labels,vsym)
        end

        t.write_string(0,col,varstr,formats[:heading])
        u = univariate(df,vsym) #,wt=df[wt])
        for j = 1:14
            if j<4
                fmttype = :n_fmt
            else
                fmttype = :p_fmt
            end
            if isnan(u[j,:Value]) || isinf(u[j,:Value])
                t.write(j,col,"",formats[fmttype])
            else
                t.write(j,col,u[j,:Value],formats[fmttype])
            end
        end
        smallest=Stella.smallest(df[vsym])
        if eltype(df) <: Integer
            fmttype = :n_fmt
        else
            fmttype = :p_fmt
        end
        for j = 1:5
            t.write(j+14,col,smallest[j],formats[fmttype])
        end
        largest=Stella.largest(df[vsym])
        for j = 1:5
            t.write(j+19,col,largest[j],formats[fmttype])
        end
        col += 1
    end
end
function univariatexls(df::DataFrame,contvars::Vector{Symbol},wbook::AbstractString,wsheet::AbstractString;
    wt::Union{Nothing,Symbol} = nothing,labels::Union{Nothing,Label}=nothing, row = 0, col = 0)

    xlsxwriter=pyimport("xlsxwriter")

    wb = wlsxwriter.Workbook(wbook)

    univariatexls(df,contvars,wb,wsheet,wt=wt,labels=labels,row=row,col=col)

    wb.close()
end

"""
    dfxls(df::DataFrame, workbook::PyObject, worksheet::AbstractString; nrows = 500, start = 1, row=0, col=0)

 To use this function, `PyCall` is required with a working version python and
 a python package called `xlsxwriter` installed. Options are:

- `df`: a DataFrame
- `workbook`: a returned value from xlsxwriter.Workbook() function (see an example below)
- `worksheet`: a string for the worksheet name (default: "Data1")
- `start`: specify the row number from which data will be output (default: 1)
- `nrows`: specify the number of rows to output (default: 500). If nows = 0, the entire dataframe will be output.
- `row`: specify the row of the workbook to start the output table (default = 0 (for row 1))
- `col`: specify the column of the workbook to start the output table (default = 0 (for column A))

# Example 1
This example is useful when one wants to append a worksheet to an existing workbook.
It is responsibility of the user to open a workbook before the function call and close it
to actually create the physical file by close the workbook.

```
julia> using PyCall

julia> xlsxwriter = pyimport("xlsxwriter")

julia> wb = xlsxwriter.Workbook("test_workbook.xlsx")
PyObject <xlsxwriter.workbook.Workbook object at 0x000000002A628E80>

julia> glmxls(ols1,wb,"OLS1",labels = label)

julia> bivairatexls(df,:incomecat,[:age,:race,:male,:bmicat],wb,"Bivariate",labels = label)

julia> univariatexls(df,[:age,:income_amt,:bmi],wb,"Univariate",labels = label)

julia> dfxls(df,wb,"dataframe",nrows = 0)

Julia> wb.close()
```

# Example 2
Alternatively, one can create a spreadsheet file directly. `PyCall` or `@pyimport`
does not need to be called before the function.

```
julia> dfxls(df,"test_workbook.xlsx","df",nrows = 0)
```

"""
function dfxls(df::DataFrame,
    wbook::PyObject,
    worksheet::AbstractString;
    nrows::Int64 = 0, start::Int64 = 1, col::Int64 = 0, row::Int64 = 0)

    # create a worksheet
    t = wbook.add_worksheet(worksheet)

    # attach formats to the workbook
    formats = ExcelTables.attach_formats(wbook)

    # starting row and column
    c = col

    # Eltype
    typ = Vector{DataType}(undef,size(df,2))
    for i=1:size(df,2)
        if typeof(df[i]) <: CategoricalArray
            typ[i] = eltype(df[i].pool.index)
        else
            typ[i] = nonmissingtype(eltype(df[i]))
        end
    end
    varnames = names(df)

    # if nrows = 0, output the full data set
    if nrows == 0
        start = 1
        nrows = size(df,1)
    else
        nrows = (start+nrows-1) < size(df,1) ? nrows : (size(df,1)-start+1)
    end

    for i = 1:size(df,2)

        r = row
        t.set_column(c,c,10)
        t.write_string(r,c,string(varnames[i]),formats[:heading])
        r += 1

        for j in start:(start+nrows-1)
            # println("i = ",i,"; j = ",j,". value = ",df[j,i])

            if ismissing(df[j,i])
                t.write_string(r,c," ",formats[:n_fmt])
            elseif typ[i] <: AbstractString
                if df[j,i] == ""
                    t.write_string(r,c," ",formats[:text])
                else
                    t.write(r,c,df[j,i],formats[:text])
                end
            elseif typ[i] <: Number
                if isnan(df[j,i]) || isinf(df[j,i])
                    t.write_string(r,c," ",formats[:text])
                elseif typ[i] <: Integer
                    t.write(r,c,df[j,i],formats[:n_fmt])
                elseif typ[i] <: AbstractFloat
                    t.write(r,c,df[j,i],formats[:f_fmt])
                end
            elseif typ[i] <: Date
                t.write(r,c,Dates.value(df[j,i] - Date(1899,12,30)),formats[:f_date])
            elseif typ[i] <: DateTime
                t.write(r,c,(Dates.value(df[j,i] - DateTime(1899,12,30,0,0,0)))/86400000,formats[:f_datetime])
            elseif typ[i] == Symbol || typ[i] == DataType
                t.write(r,c,string(df[j,i]),formats[:text])
            else
                # skip
            end

            r += 1
        end
        c += 1
    end

end
function dfxls(df::DataFrame,
    wbook::AbstractString,
    worksheet::AbstractString;
    nrows::Int64 = 0, start::Int64 = 1, col::Int64 = 0, row::Int64 = 0)

    # import xlsxwriter
    XlsxWriter = pyimport("xlsxwriter")

    # create a workbook
    wb = XlsxWriter.Workbook(wbook)

    dfxls(df,wb,worksheet, nrows = nrows, start = start, col = col, row = row)

    wb.close()
end

function newfilename(filen::AbstractString)
    while (isfile(filen))
        # separate the name into three parts, basename (number).extension
        (basename,ext) = splitext(filen)
        m = match(r"(.*)(\(.*\))$",basename)
        if m == nothing
            filen = string(filen," (1)")
        else
            m2 = parse(Int,replace(m[2],r"\((.*)\)",s"\1"))
            filen = string(m[1]," (",m2+1,")")
        end
    end
    return filen
end

function countlev(str::AbstractString,sarray::Vector{String})
    k = 0
    for i = 1:length(sarray)
        if str == sarray[i]
            k += 1
        end
    end
    return k
end


text:
using StaticArrays
using MeshCore: nshapes, VecAttrib
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using MeshCore: attribute
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: IncRel

"""
    orient_surface_mesh(t2v)

Orient surface mesh.

Return
- `orientedt2v`: oriented incidence relation
- `orientable`: Can the surface mesh be oriented? Bool flag.
    This would be false for a Mobius strip, for instance.

This is a purely topological operation. Creases in the surface are not
recognized as edges. The operation proceeds by flooding the surface across
manifold edges, stopping at sheet or non-manifold edges.

The incidence relation `orientedt2v` includes an attribute for each triangle
that provides the surface id to which the triangle belongs
(`orientedt2v.left.attributes["surfid"]`).
"""
function orient_surface_mesh(t2v)
    orientable = true
    all_oriented = false
    surfid = fill(0, nshapes(t2v.left))
    e2v = ir_skeleton(t2v)
    t2e = ir_bbyfacets(t2v, e2v)
    ot2ev = [SVector(t2e[i]) for i in 1:nshapes(t2e.left)]
    e2t = ir_transpose(t2e)
    while !all_oriented
        all_oriented = true
        for i in 1:nshapes(t2v.left)
            if surfid[i] == 0
                orientable = orientable && _orient_surface!(ot2ev, surfid, t2v, i, t2e, e2t);
                all_oriented = false;
                break;
            end
        end
    end
    elements = ShapeColl(T3, nshapes(t2v.left))
    vrts = ShapeColl(P1, nshapes(t2v.right))
    oc = [_vertfromedges(ot2ev[i], e2v) for i in 1:length(ot2ev)]
    orientedt2v = IncRel(elements, vrts, oc)
    orientedt2v.left.attributes["surfid"] = VecAttrib(surfid)
    orientedt2v.right.attributes["geom"] = attribute(t2v.right, "geom")
    return orientedt2v, orientable
end

function _vertfromedges(e, e2v)
    e1id = abs(e[1])
    e2id = abs(e[2])
    e3id = abs(e[3])
    _commonv(ea, eb) = (ea[1] == eb[1] ? ea[1] : (ea[2] == eb[1] ? ea[2] : (ea[2] == eb[2] ? ea[2] : ea[1])))
    if e[1] > 0
        return SVector(e2v[e1id][1], e2v[e1id][2], _commonv(e2v[e2id], e2v[e3id]))
    else
        return SVector(e2v[e1id][2], e2v[e1id][1], _commonv(e2v[e2id], e2v[e3id]))
    end
end

function _iscompatible(je, kes)
    for e in 1:3
        if je == -kes[e]
            return true
        end
    end
    return false
end

function _deduce_edge_order(je, kes)
    if _iscompatible(je, kes)
        return kes
    else
        return SVector(-kes[3], -kes[2], -kes[1])
    end
end

function _orient_surface!(ot2ev, surfid, t2v, i, t2e, e2t)
    orientable = true
    currsurfid = maximum(surfid) + 1
    ot2ev[i] = t2e[i] # First Triangle: Accept its default orientation
    surfid[i] = currsurfid
    stack = Int[]
    push!(stack, i)
    while (!isempty(stack))
        j = pop!(stack)
        for e in 1:3
            eid = abs(ot2ev[j][e])
            if length(e2t[eid]) == 2 # is this a manifold edge?
                k = e2t[eid][1] == j ? e2t[eid][2] : e2t[eid][1]
                if surfid[k] == 0
                    ot2ev[k] = _deduce_edge_order(ot2ev[j][e], ot2ev[k])
                    surfid[k] = currsurfid
                    push!(stack, k)
                else
                    orientable = orientable && _iscompatible(ot2ev[j][e], ot2ev[k])
                end
            end
        end
    end
    return orientable
end

function _normal(c, geom)
    n = StaticArrays.cross(geom[c[2]] - geom[c[1]], geom[c[3]] - geom[c[1]])
    n /= StaticArrays.norm(n)
    return n
end

"""
    _classify_mesh_edges(t2v, e2t, crease_ang)

Classify mesh edges as on/not on topological edge.

- `e2t`: the attribute `isonte`
"""
function _classify_mesh_edges(t2v, e2t, crease_ang)
    cada = cos(crease_ang);
    geom = attribute(t2v.right, "geom")
    surfid = attribute(t2v.left, "surfid")
    # e2v = ir_skeleton(t2v)
    # t2e = ir_bbyfacets(t2v, e2v)
    # ot2ev = [SVector(t2e[i]) for i in 1:nshapes(t2e.left)]
    # e2t = ir_transpose(t2e)
    isonte = fill(false, nshapes(e2t.left))
    for e in 1:nshapes(e2t.left)
        numt_at_e = length(e2t[e])
        @assert numt_at_e >= 1
        if numt_at_e == 2 # manifold edge
            # if the triangles belong to different surfaces, this is a TE
            if surfid[e2t[e][1]] != surfid[e2t[e][2]]
                println("different surfaces, $(surfid[e2t[e][1]]) vs $(surfid[e2t[e][2]])")
                isonte[e] = true
            else
                # check if it is a crease
                n1 = _normal(t2v[e2t[e][1]], geom)
                n2 = _normal(t2v[e2t[e][2]], geom)
                if StaticArrays.dot(n1, n2) < cada
                    isonte[e] = true
                end
            end
        elseif numt_at_e == 1 # sheet edge
            isonte[e] = true
        else # non manifold edge
            isonte[e] = true
        end
    end
    e2t.left.attributes["isonte"] = VecAttrib(isonte)
    return e2t
end

"""
    make_topo_faces(t2v, e2t)

Make topological faces.

Returns
- `t2v`: the attribute `tfid` lists the numbers of the topological faces.
"""
function make_topo_faces(t2v, crease_ang = 30/180*pi)
    e2v = ir_skeleton(t2v)
    e2v.right.attributes["geom"] = t2v.right.attributes["geom"]
    t2e = ir_bbyfacets(t2v, e2v)
    e2t = ir_transpose(t2e)
    e2t = _classify_mesh_edges(t2v, e2t, crease_ang)
    all_done = false
    tfid = fill(0, nshapes(t2v.left))
    while !all_done
        all_done = true
        for i in 1:nshapes(t2v.left)
            if tfid[i] == 0
                _collect_tf!(tfid, t2v, i, t2e, e2t)
                all_done = false;
                break;
            end
        end
    end
    t2v.left.attributes["tfid"] = VecAttrib(tfid)
    return t2v
end

function _collect_tf!(tfid, t2v, i, t2e, e2t)
    isonte = attribute(e2t.left, "isonte")
    currtfid = maximum(tfid) + 1
    tfid[i] = currtfid
    stack = Int[]
    push!(stack, i)
    while (!isempty(stack))
        j = pop!(stack)
        for e in 1:3
            eid = abs(t2e[j][e])
            if !isonte[eid]
                k = e2t[eid][1] == j ? e2t[eid][2] : e2t[eid][1]
                if tfid[k] == 0
                    tfid[k] = currtfid
                    push!(stack, k)
                end
            end
        end
    end
    return tfid
end

ext: jl


text:
module PDDL

using Base: @kwdef
using Julog
using AutoHashEquals

# Logical formulae and terms
export Term, Compound, Var, Const, Clause
# Abstract data types
export Domain, Problem, State, Action
# Generic concrete data types
export GenericDomain, GenericProblem, GenericState, GenericAction
# Interface methods
export satisfy, satisfiers, evaluate
export initstate, goalstate, transition, transition!
export available, execute, execute!, relevant, regress, regress!
# Parsing and writing
export parse_domain, parse_problem, parse_pddl, @pddl, @pddl_str
export write_domain, write_problem, write_pddl
export load_domain, load_problem
export save_domain, save_problem
# Simulators
export Simulator, EndStateSimulator, StateRecorder
# Abstract interpretation
export abstracted, abstractstate, lub, glb, widen, widen!
export AbstractedDomain
# Domain compilation
export compiled, compilestate
export CompiledDomain, CompiledAction, CompiledState
# Domain and action grounding
export ground, groundargs, groundactions
export GroundDomain, GroundAction, GroundActionGroup
# Extension interfaces
export attach!, register!
# Analysis tools
export infer_affected_fluents, infer_static_fluents
# Utilities
export find_matches, effect_diff, precond_diff

# PDDL requirement definitions and dependencies
include("requirements.jl")
# Abstract interface for PDDL-based planners and applications
include("interface/interface.jl")
# Generic concrete representations of PDDL types
include("generic/generic.jl")
# Parser for PDDL files
include("parser/parser.jl")
# Writer for PDDL files
include("writer/writer.jl")
# Simulators for executing action sequences and recording outputs
include("simulators/simulators.jl")
# Built-in functions and operators recognized by interpreters / compilers
include("builtins.jl")
# Methods for extending built-in functions and operators
include("extensions.jl")
# Theories for fluents with non-standard data types (e.g. sets, arrays)
include("theories/theories.jl")
# Abstractions for fluents of various types
include("abstractions/abstractions.jl")
# Interpreter-based interface implementation
include("interpreter/interpreter.jl")
# Compiler-based interface implementations
include("compiler/compiler.jl")
# Domain and action grounding
include("grounding/grounding.jl")
# Tools for analyzing domains
include("analysis/analysis.jl")

using .Parser, .Writer

end # module


text:
@symbol_func function v_te(cur_reactor::AbstractReactor, cur_rho::AbstractSymbol)
  cur_T = T_rho(cur_reactor, cur_rho)

  cur_v = 1.88e7

  cur_v *= sqrt(cur_T)

  cur_v
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAppsV1DeploymentStrategy <: SwaggerModel
    rollingUpdate::Any # spec type: Union{ Nothing, IoK8sApiAppsV1RollingUpdateDeployment } # spec name: rollingUpdate
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiAppsV1DeploymentStrategy(;rollingUpdate=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1DeploymentStrategy, Symbol("rollingUpdate"), rollingUpdate)
        setfield!(o, Symbol("rollingUpdate"), rollingUpdate)
        validate_property(IoK8sApiAppsV1DeploymentStrategy, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiAppsV1DeploymentStrategy

const _property_map_IoK8sApiAppsV1DeploymentStrategy = Dict{Symbol,Symbol}(Symbol("rollingUpdate")=>Symbol("rollingUpdate"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiAppsV1DeploymentStrategy = Dict{Symbol,String}(Symbol("rollingUpdate")=>"IoK8sApiAppsV1RollingUpdateDeployment", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiAppsV1DeploymentStrategy }) = collect(keys(_property_map_IoK8sApiAppsV1DeploymentStrategy))
Swagger.property_type(::Type{ IoK8sApiAppsV1DeploymentStrategy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1DeploymentStrategy[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1DeploymentStrategy }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1DeploymentStrategy[property_name]

function check_required(o::IoK8sApiAppsV1DeploymentStrategy)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1DeploymentStrategy }, name::Symbol, val)
end


text:
##############
# CONVERSION #
##############

function convert(::Type{RevisedAffineForm{N}}, x::T) where {N,T<:Real}
    TS = float(T)
    return RevisedAffineForm(x, SVector{N,TS}(zeros(TS, N)), zero(TS))
end

function convert(::Type{RevisedAffineForm{N}}, x::Interval{T}) where {N,T<:AbstractFloat}
    c, r = midpoint_radius(x)
    return RevisedAffineForm(c, SVector{N,T}(zeros(T, N)), r)
end

function convert(::Type{RevisedAffineForm{N,T}}, x::Real) where {N,T<:AbstractFloat}
    return RevisedAffineForm(convert(T, x), SVector{N,T}(zeros(T, N)), zero(T))
end

function convert(::Type{RevisedAffineForm{N,T}}, x::Interval) where {N,T<:AbstractFloat}
    c, r = midpoint_radius(x)
    return RevisedAffineForm(convert(T, c), SVector{N,T}(zeros(T, N)), convert(T, r))
end

#############
# PROMOTION #
#############

function promote_rule(
    ::Type{RevisedAffineForm{N,T}}, ::Type{Interval{S}}
) where {N,T<:AbstractFloat,S<:Real}
    return RevisedAffineForm{N,promote_type(T, S)}
end

function promote_rule(
    ::Type{RevisedAffineForm{N,T}}, ::Type{S}
) where {N,T<:AbstractFloat,S<:Real}
    return RevisedAffineForm{N,promote_type(T, S)}
end

function promote_rule(
    ::Type{RevisedAffineForm{N,T}}, ::Type{RevisedAffineForm{N,S}}
) where {N,T<:AbstractFloat,S<:AbstractFloat}
    return RevisedAffineForm{N,promote_type(T, S)}
end


text:
##
using DelimitedFiles
input_path = joinpath(@__DIR__(), "input.txt")
flip(x, y) = (y, x)
flip((x,y)) = (y, x)
flip(c::CartesianIndex) = CartesianIndex(flip(c.I))

##
input_data = readdlm(input_path,',',String)
encode_move(::Val{'U'}) = 0x00
encode_move(::Val{'D'}) = 0x11
encode_move(::Val{'R'}) = 0x01
encode_move(::Val{'L'}) = 0x10
data = map(input_data) do v
    k = encode_move(Val(v[1]))
    l = parse(Int, v[2:end])
    k, l
end |> permutedims
paths = [eachcol(data)...]

##
function decode_move(k, l)
    k == 0x00 && return CartesianIndex(0,l)
    k == 0x11 && return CartesianIndex(0,-l)
    k == 0x01 && return CartesianIndex(l,0)
    k == 0x10 && return CartesianIndex(-l,0)
    CartesianIndex(0,0)
end
decode_move(t::Tuple) = decode_move(t...)
function follow_path(paths::AbstractVector{<:AbstractVector{Tuple{UInt8,Int}}})
    map(paths) do mvs
        accumulate(mvs; init=CartesianIndex(0,0)) do cur, (k, l)
            cur + decode_move(k,l)
        end
    end
end
function intersections(paths::AbstractVector{<:AbstractVector{CartesianIndex{2}}})
    intersect(paths...)
end
dist(xs) = sum(abs,xs)
dist(c::CartesianIndex) = dist(c.I)

## TESTS
test_cases = (
    "R8,U5,L5,D3",
    "U7,R6,D4,L4",
    "R75,D30,R83,U83,L12,D49,R71,U7,L72,U62,R66,U55,R34,D71,R55,D58,R83",
    "R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51,U98,R91,D20,R16,D67,R40,U7,R15,U6,R7",
)
res = map(test_cases) do x
    a = split(x,",")
    path = map(a) do v
        k = encode_move(Val(v[1]))
        l = parse(Int, v[2:end])
        k, l
    end
    follow_path([path])
end

##
test_dist = (
    "R8,U5,L5,D3",
    "U7,R6,D4,L4",
)
paths = map(test_dist) do x
    a = split(x,",")
    path = map(a) do v
        k = encode_move(Val(v[1]))
        l = parse(Int, v[2:end])
        k, l
    end
end
res = follow_path(paths|>collect)
display(res)

##
points = map(res) do ps
    prec = CartesianIndex(0,0)
    map(ps) do p
        x, y = p.I
        px, py = prec.I
        dx, dy = sign(x-px), sign(y-py)
        prec = p
        dx == 0 && return [CartesianIndex(px, j) for j in py+dy:dy:y]
        [CartesianIndex(j, py) for j in px+dx:dx:x]
    end
end
points = map(x -> reduce(vcat, x), points)
map(length, points)

##
res = follow_path(paths)
nothing

##
minimum(dist.(intersections(points)))

##
inters = intersections(points)

##
# steps2inters = map(points) do ps
#     findall(ps) do p
#         any(==(p), inters)
#     end
# end
steps2inters = map(inters) do i
    reduce(vcat, map(points) do ps
        findfirst(==(i), ps)
    end)
end

##
minimum(sum, steps2inters)


text:
module Models

using Parameters
using ACEtb.SlaterKoster: TwoCentreModel, SKH
using JuLIP.Potentials: fcut

import ACEtb.SlaterKoster: eval_bonds!, eval_diag!, getskh
import JuLIP: cutoff

export KwonTB

"""
`KwonTB <: SKModel`

Hamiltonian for an orthogonal sp TB model of Si developed by Kwon et al [1].

This implementation deviates  from [1] in how the cut-off is applied:
instead of "patching" a cubic spline between r1 and rcut, we simply multiply
with a quintic spline on the interval [0.5 (rcut + r0), rcut].

[1] I. Kwon, R. Biswas, C. Z. Wang, K. M. Ho and C. M. Soukoulis.
Transferable tight-binding models for silicon.
Phys Rev B 49 (11), 1994.
"""
@with_kw struct KwonTB <: TwoCentreModel
   skh::SKH = SKH("sp")
   # -------------------------------
   r0::Float64 = 2.360352   # Å
   Es::Float64 = -5.25      # eV
   Ep::Float64 = 1.2        # eV
   E0::Float64 = 8.7393204  # eV   8.7393204  is the original value; but it is just a constant shift
   # ------------------------------- Electronic Parameters
   # α   1     2    3     4
   #    ssσ   spσ  ppσ   ppπ
   hr0::NTuple{4, Float64} = (-2.038, 1.745, 2.75, -1.075)
   nc::NTuple{4, Float64} = (9.5, 8.5, 7.5, 7.5)
   rc::NTuple{4, Float64} = (3.4, 3.55, 3.7, 3.7)
   # ------------------------------- 2-body parameters
   m::Float64 = 6.8755
   mc::Float64 = 13.017
   dc::Float64 = 3.66995     # Å
   C::NTuple{4, Float64} = (2.1604385, -0.1384393, 5.8398423e-3, -8.0263577e-5)
   # --------------------------------
   r1::Float64 = 3.260176     # Å  (start of cut-off) >>> different r1 from Kwon paper
   rcut::Float64 = 4.16       # Å  (end of cut-off)
end

cutoff(model::KwonTB) = model.rcut

getskh(model::KwonTB) = model.skh

# isorthogonal(H::KwonTB) = true

kwon_bond(m::KwonTB, r::Real, α) = (
          m.hr0[α] * (m.r0 / r)^2 *
               exp( - 2 * (r/m.rc[α])^m.nc[α] + 2 * (m.r0/m.rc[α])^m.nc[α] ) )

function eval_bonds!(V, model::KwonTB, r)
   fc = fcut(r, model.r1, model.rcut)
   for α = 1:4
      V[α] = kwon_bond(model, r, α) * fc
   end
   return V
end

function eval_diag!(D, model::KwonTB, r, R)
   D[1] = model.Es
   D[2:4] .= model.Ep
   return D
end

end


text:
# # Fractional Zener
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/model_zener.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/model_zener.ipynb)

using RHEOS
## include a helper function for plotting
include("assets/plothelper.jl");
#-

Fract_Zener

# #### Constitutive Equation
# ```math
# \sigma(t) + \frac{c_\alpha}{c_\beta} \frac{d^{\alpha-\beta} \sigma(t)}{dt^{\alpha-\beta}}= c_{\alpha} \frac{d^\alpha \epsilon(t)}{dt^\alpha} + c_\gamma \frac{d^\gamma \epsilon(t)}{dt^\gamma}+ \frac{c_\alpha c_\gamma}{c_\beta}\frac{d^{\alpha+\gamma-\beta} \epsilon(t)}{dt^{\alpha+\gamma-\beta}}
# ```

# ```math
# \text{for}\; \ 0 \leq \beta \leq \alpha \leq 1
# ```

# #### Relaxation Modulus
# ```math
# G(t) = c_\beta t^{-\beta} E_{\alpha-\beta,1-\beta}\left(-\frac{c_\beta}{c_\alpha} t^{\alpha-\beta}\right) + \frac{c_{\gamma} }{\Gamma(1-\gamma)} t^{-\gamma}
# ```

# #### Creep Modulus
# ```math
# \tilde{J}(s)= \frac{1}{s}\frac{c_\alpha s^{\alpha}+c_\beta s^{\beta}}{c_\alpha s^\alpha c_\beta s^{\beta} + c_\gamma s^\gamma (c_\alpha s^\alpha+c_\beta s^\beta)}
# ```

# #### Storage Modulus
# ```math
# G^{\prime}(\omega) = \frac{\left(c_\beta \omega^\beta\right)^2 \cdot c_\alpha \omega^\alpha \cos(\alpha \frac{\pi}{2}) + \left(c_\alpha \omega^\alpha\right)^2 \cdot c_\beta \omega^\beta \cos(\beta \frac{\pi}{2})}{\left(c_\alpha \omega^\alpha\right)^2+\left(c_\beta \omega^\beta\right)^2+2c_\alpha \omega^\alpha \cdot c_\beta \omega^\beta \cos((\alpha-\beta)\frac{\pi}{2})} + c_\gamma \omega^\gamma \cos\left(\gamma \frac{\pi}{2}\right)
# ```

# #### Loss Modulus
# ```math
# G^{\prime\prime}(\omega) = \frac{\left(c_\beta \omega^\beta\right)^2 \cdot c_\alpha \omega^\alpha \sin(\alpha \frac{\pi}{2}) + \left(c_\alpha \omega^\alpha\right)^2 \cdot c_\beta \omega^\beta \sin(\beta \frac{\pi}{2})}{\left(c_\alpha \omega^\alpha\right)^2+\left(c_\beta \omega^\beta\right)^2+2c_\alpha \omega^\alpha \cdot c_\beta \omega^\beta \cos((\alpha-\beta)\frac{\pi}{2})}+ c_\gamma \omega^\gamma \sin\left(\gamma \frac{\pi}{2}\right)
# ```

# # Fractional SLS (Zener)

FractSLS_Zener
#-

models = Vector{RheoModel}()

## plot moduli for varying α
for alpha in [0.1, 0.25, 0.5, 0.74, 0.9]

    push!(models, RheoModel(FractSLS_Zener, (cₐ = 1, a = alpha, kᵦ = 1, kᵧ = 1)))

end

plotmodel(models)

# # Fractional Jeffreys (Zener)

FractJeffreys_Zener
#-

models = Vector{RheoModel}()

## plot moduli for varying β
for beta in [0.1, 0.25, 0.5, 0.74, 0.9]

    push!(models, RheoModel(FractJeffreys_Zener, (ηₐ = 1, cᵦ = 1, β = beta, ηᵧ = 1)))

end

plotmodel(models, ymaxG = 2.0)

# # Standard Linear Solid (Zener)

SLS_Zener
#-

models = Vector{RheoModel}()

## plot moduli for varying kᵦ
for k in [1.0, 3.0, 5.0]

    push!(models, RheoModel(SLS_Zener, (η = 1, kᵦ = k, kᵧ = 1)))

end

plotmodel(models)

# # Jeffreys (Zener)

Jeffreys_Zener
#-

models = Vector{RheoModel}()

## plot moduli for varying ηₐ
for eta in [1.0, 5.0, 8.0]

    push!(models, RheoModel(Jeffreys_Zener, (ηₐ = eta, k = 3, ηᵧ = 1)))

end

plotmodel(models)

# # Fractional Solid

FractSolid
#-

models = Vector{RheoModel}()

## plot moduli for varying β
for beta in [0.1, 0.25, 0.5, 0.74, 0.9]

    push!(models, RheoModel(FractSolid, (η = 1, cᵦ = 1, β = beta, k = 0.5)))

end

plotmodel(models, ymaxG = 2)


text:
using Test, GenericLinearAlgebra, LinearAlgebra

@testset "The General eigenvalue problem" begin

cplxord = t -> (real(t), imag(t))

@testset "General eigen problem with n=$n and element type=$T" for
    n in (10, 23, 100),
    T in (Float64, Complex{Float64})

    A = randn(T, n, n)
    vGLA    = GenericLinearAlgebra._eigvals!(copy(A))
    vLAPACK = eigvals(A)
    vBig    = eigvals(big.(A)) # not defined in LinearAlgebra so will dispatch to the version in GenericLinearAlgebra
    @test sort(vGLA, by = cplxord) ≈ sort(vLAPACK, by = cplxord)
    @test sort(vGLA, by = cplxord) ≈ sort(complex(eltype(A)).(vBig), by = cplxord)
    if VERSION > v"1.2.0-DEV.0"
        @test issorted(vBig, by = cplxord)
    end
end

@testset "make sure that solver doesn't hang" begin
    for i in 1:1000
        A = randn(8, 8)
        sort(abs.(GenericLinearAlgebra._eigvals!(copy(A)))) ≈ sort(abs.(eigvals(A)))
    end
end

@testset "Convergence in corner cases. Issue 29." begin
    function H(n::Int)
        H = zeros(2n, 2n)
        for i = 1 : 2 : 2n
            H[i, i+1] = 1
            H[i+1, i] = 1
        end
        H
    end

    function E(n::Int)
        E = zeros(2n, 2n)
        for i = 1 : (n - 1)
            E[2i + 1, 2i] = 1
        end
        E[1, 2n] = 1
        E
    end

    my_matrix(n::Int, η::Float64 = 1e-9) = H(n) .+ η .* E(n)

    A = my_matrix(4, 1e-3);
    @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(A), by = t -> (real(t), imag(t)))
end

@testset "Convergence in with 0s Issue 58." begin
    A = [0.0 1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 0.0]
    @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(A))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(A), by = t -> (real(t), imag(t)))
    B = [0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 -1.0 0.0]
    @test sort(GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(B))), by = t -> (real(t), imag(t))) ≈ sort(eigvals(B), by = t -> (real(t), imag(t)))
end

@testset "Extract Schur factor" begin
    A = randn(5, 5)
    @test sum(eigvals(schur(A).T))     ≈ sum(eigvals(Float64.(schur(big.(A)).T)))
    @test sum(eigvals(schur(A).Schur)) ≈ sum(eigvals(Float64.(schur(big.(A)).Schur)))
end

@testset "Issue 63" begin
    A = [1 2 1 1 1 1
         0 1 0 1 0 1
         1 1 2 0 1 0
         0 1 0 2 1 1
         1 1 1 0 2 0
         0 1 1 1 0 2]

    z1 = [complex(1, sqrt(big(3))), complex(1, -sqrt(big(3)))]
    z2 = [
        4*2^(big(2)/3)/complex(5, 3*sqrt(big(111)))^(big(1)/3),
        complex(5, 3*sqrt(big(111)))^(big(1)/3)/2^(big(2)/3)
        ]

    truevals = real.([
        (7 - transpose(z1)*z2),
        3,
        3,
        3,
        (7 - adjoint(z1)*z2),
        (7 + [2, 2]'*z2)])/3

    if VERSION > v"1.2.0-DEV.0"
        @test eigvals(big.(A)) ≈ truevals
    else
        @test sort(eigvals(big.(A)), by=real) ≈ truevals
    end
end

Demmel(η) = [0  1 0 0
             1  0 η 0
             0 -η 0 1
             0  0 1 0]

@testset "Demmel matrix" for t in (1e-10, 1e-9, 1e-8)
    # See "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it"
    A = Demmel(t)
    vals = GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(A, maxiter=35))
    @test abs.(vals) ≈ ones(4)
end

function Hevil2(θ, κ, α, γ)
    # Eq (13) and (14)
    β = ω = 0.0
    ν = cos(θ)*cos(2γ) + cos(α + β + ω)*sin(2γ)*κ/2
    σ = 1 + κ*sin(2γ)*cos(α + β + ω - θ) + κ^2*sin(γ)^2
    μ = -sin(θ)*cos(2γ) - sin(α + β + ω)*sin(2γ)*κ/2
    ρ = sqrt(σ - ν^2)

    return [ν    (cos(2θ) - ν^2)/ρ              μ/ρ*(cos(2θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2)      (-2*μ*ν - sin(2*θ))/sqrt(ρ^2 - μ^2)
            ρ   -ν - sin(2*θ)*μ/ρ^2            -μ/ρ^2*(μ*sin(2*θ) + 2*ν*ρ^2)/sqrt(ρ^2 - μ^2)    -μ/ρ*(cos(2*θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2)
            0    sin(2*θ)*sqrt(ρ^2 - μ^2)/ρ^2   ν + sin(2*θ)*μ/ρ^2                              (cos(2*θ) - ν^2)/ρ
            0    0                              ρ                                               -ν]
end

@testset "Complicate matrix from Sandia technical report" begin
    H = Hevil2(0.111866322512629152, 1.08867072154101741, 0.338146383137297168, -0.313987810419091240)

    @test Float64.(abs.(eigvals(big.(H)))) ≈ ones(4)
end

@testset "Issue 67" for (A, λs) in (
    ([1 -2 1 -1 -1 0; 0 1 0 1 0 1; 1 -1 2 0 -1 0; 0 1 0 2 1 1; 1 0 1 0 0 0; 0 -1 1 -1 -2 0]    ,
        [0.0, 0.0, (3 - √big(3)*im)/2, (3 + √big(3)*im)/2, (3 - √big(3)im)/2, (3 + √big(3)im)/2]),
    ([1 0 -1 0 0 0; 0 1 1 1 -1 0; 1 0 -1 0 0 0; 1 -1 0 -1 1 1; 0 0 1 0 0 0; -1 0 1 0 0 0]      ,
        zeros(6)),
    ([1 -2 -1 1 -1 0; 0 1 0 -1 0 1; -1 1 2 0 1 0; 0 -1 -2 2 -1 -1; 1 0 -1 0 0 0; 0 -1 -1 1 0 0],
        [(1 - √big(3)*im)/2, (1 + √big(3)*im)/2, (1 - √big(3)im)/2, (1 + √big(3)im)/2, 2, 2]),
    ([2 0 -1 0 0 1; 0 2 0 -1 -1 0; -1 0 1 0 0 -1; 0 -1 0 1 1 0; 0 1 0 -1 0 0; -1 0 1 0 0 0]    ,
        ones(6)),
    ([1 0 1 0 -1 0; -2 1 2 -1 1 1; -1 0 -1 0 1 0; 2 1 2 -1 -2 1; 1 0 1 0 -1 0; 1 -1 -2 1 0 -1] ,
        [-1, -1, 0, 0, 0, 0]))

    @testset "shouldn't need excessive iterations (30*n) in the Float64 case" begin
        GenericLinearAlgebra._schur!(float(A))
    end

    # For BigFloats, many iterations are required for convergence
    # Improving this is probably a hard problem
    vals = eigvals(big.(A), maxiter=1500)

    # It's hard to compare complex conjugate pairs so we compare the real and imaginary parts separately
    @test sort(real(vals)) ≈ sort(real(λs)) atol=1e-25
    @test sort(imag(vals)) ≈ sort(imag(λs)) atol=1e-25
end

end


text:
struct Trajectory{T}
    q::Vector{Vector{T}} # generalized coordinates
    v::Vector{Vector{T}} # generalized velocities
    u::Vector{Vector{T}} # control inputs
    γ::Vector{Vector{T}} # impact impulses
    b::Vector{Vector{T}} # friction impulses
    w::Vector{Vector{T}} # disturbances
end

function Trajectory(model, T; nv=model.nq, nc=model.nc, nb=model.nc)
    nq = model.nq
    nu = model.nu
    nw = model.nw

    q = [zeros(nq) for t = 1:T+2]
    v = [zeros(nv) for t = 1:T+1] # midpoint velocities; TODO: knotpoint velocities
    u = [zeros(nu) for t = 1:T]
    γ = [zeros(nc) for t = 1:T]
    b = [zeros(nb) for t = 1:T]
    w = [zeros(nw) for t = 1:T]
    Trajectory(q, v, u, γ, b, w)
end

function reset!(traj::Trajectory)
    T = length(traj.u)
    for t = 1:T
        fill!(traj.q[t], 0.0)
        fill!(traj.v[t], 0.0)
        fill!(traj.u[t], 0.0)
        fill!(traj.γ[t], 0.0)
        fill!(traj.b[t], 0.0)
        fill!(traj.w[t], 0.0)
    end
    fill!(traj.q[T+1], 0.0)
    fill!(traj.q[T+2], 0.0)
    fill!(traj.v[T+1], 0.0)
    return nothing
end

struct GradientTrajectory{T}
    ∂q3∂q1::Vector{Matrix{T}}
    ∂q3∂q2::Vector{Matrix{T}}
    ∂q3∂u1::Vector{Matrix{T}}
    ∂γ1∂q1::Vector{Matrix{T}}
    ∂γ1∂q2::Vector{Matrix{T}}
    ∂γ1∂u1::Vector{Matrix{T}}
    ∂b1∂q1::Vector{Matrix{T}}
    ∂b1∂q2::Vector{Matrix{T}}
    ∂b1∂u1::Vector{Matrix{T}}
end

function GradientTrajectory(model, T; nc=model.nc, nb=model.nc)
    nq = model.nq
    nu = model.nu

    ∂q3∂q1 = [zeros(nq, nq) for t = 1:T]
    ∂q3∂q2 = [zeros(nq, nq) for t = 1:T]
    ∂q3∂u1 = [zeros(nq, nu) for t = 1:T]
    ∂γ1∂q1 = [zeros(nc, nq) for t = 1:T]
    ∂γ1∂q2 = [zeros(nc, nq) for t = 1:T]
    ∂γ1∂u1 = [zeros(nc, nu) for t = 1:T]
    ∂b1∂q1 = [zeros(nb, nq) for t = 1:T]
    ∂b1∂q2 = [zeros(nb, nq) for t = 1:T]
    ∂b1∂u1 = [zeros(nb, nu) for t = 1:T]

    GradientTrajectory(∂q3∂q1, ∂q3∂q2, ∂q3∂u1,
                       ∂γ1∂q1, ∂γ1∂q2, ∂γ1∂u1,
                       ∂b1∂q1, ∂b1∂q2, ∂b1∂u1)
end

function reset!(traj::GradientTrajectory)
    T = length(traj.∂q3∂q1)
    for t = 1:T
        fill!(traj.∂q3∂q1[t], 0.0)
        fill!(traj.∂q3∂q2[t], 0.0)
        fill!(traj.∂q3∂u1[t], 0.0)
        fill!(traj.∂γ1∂q1[t], 0.0)
        fill!(traj.∂γ1∂q2[t], 0.0)
        fill!(traj.∂γ1∂u1[t], 0.0)
        fill!(traj.∂b1∂q1[t], 0.0)
        fill!(traj.∂b1∂q2[t], 0.0)
        fill!(traj.∂b1∂u1[t], 0.0)
    end
    return nothing
end


text:
module Grid

function grid_xp(n::Int64; pad=false::Bool)
        if pad
                x = zeros(n + 1)
        else
                x = zeros(n)
        end
        h = 1.0 / (n - 1)
        for i in 1:n
                x[i] = (i - 1) * h
        end

        return x, h
end

function grid_xm(n::Int64)
        x = zeros(n)
        h = 1.0 / (n - 2)
        for i in 1:n
                x[i] = (i - 1 - 0.5) * h
        end
        x[1] = 0.0
        x[end] = 1.0

        return x, h
end

"""
Construct a 2D grid vector from a 1D grid vector `x1`. The resulting vector
repeats all x-values: x1[1], x1[1], ... , x1[1], x1[2], x1[2] ... x1[2] `ny`
times.
"""
function grid_2d_x(x1::AbstractArray, ny::Int64)
        nx = size(x1, 1)
        x = zeros(nx * ny)
        for i=1:nx
                for j=1:ny
                        @inbounds x[j + (i - 1) * ny] = x1[i]
                end
        end
        return x
end

"""
Construct a 2D grid vector from a 1D grid vector `y1`. The resulting vector
repeats all y-values: y1[1], y1[1], ... , y1[ny], y1[1], y1[2] ... y1[ny] `nx`
times.
"""
function grid_2d_y(y1::AbstractArray, nx::Int64)
        ny = size(y1, 1)
        y = zeros(nx * ny)
        for i=1:nx
                for j=1:ny
                        @inbounds y[j + (i - 1) * ny] = y1[j]
                end
        end
        return y
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

module FMICore

# check float size (32 or 64 bits)
juliaArch = Sys.WORD_SIZE
@assert (juliaArch == 64 || juliaArch == 32) "FMICore: Unknown Julia Architecture with $(juliaArch)-bit, must be 64- or 32-bit."
Creal = Cdouble
if juliaArch == 32
    Creal = Cfloat
end

"""
The mutable struct representing an abstract (version unknown) FMU.
"""
abstract type FMU end
abstract type FMUSolution end

include("FMI2_c.jl")
include("FMI2.jl")

include("FMI3_c.jl")
include("FMI3.jl")

### EXPORTING LISTS START ###

export FMU
export FMU2ExecutionConfiguration, FMU_EXECUTION_CONFIGURATION_RESET, FMU_EXECUTION_CONFIGURATION_NO_RESET, FMU_EXECUTION_CONFIGURATION_NO_FREEING

# FMI2.jl
export FMU2, FMU2Component, FMU2ComponentEnvironment, FMU2Solution, FMU2Event
export fmi2StatusToString
export fmi2CausalityToString, fmi2StringToCausality
export fmi2VariabilityToString, fmi2StringToVariability
export fmi2InitialToString, fmi2StringToInitial
export fmi2DependencyKindToString, fmi2StringToDependencyKind

# FMI2_c.jl

# data types
export fmi2ValueReference, fmi2Component
export fmi2Char, fmi2String, fmi2Boolean, fmi2Real, fmi2Integer, fmi2Byte, fmi2FMUstate, fmi2ComponentEnvironment, fmi2Enum
export fmi2Status, fmi2StatusOK, fmi2StatusWarning, fmi2StatusDiscard, fmi2StatusError, fmi2StatusFatal, fmi2StatusPending
export fmi2CallbackFunctions
export fmi2Annotation
export fmi2Causality, fmi2CausalityParameter, fmi2CausalityCalculatedParameter, fmi2CausalityInput, fmi2CausalityOutput, fmi2CausalityLocal, fmi2CausalityIndependent
export fmi2Variability, fmi2VariabilityConstant, fmi2VariabilityFixed, fmi2VariabilityTunable, fmi2VariabilityDiscrete, fmi2VariabilityContinuous
export fmi2Initial, fmi2InitialExact, fmi2InitialApprox, fmi2InitialCalculated
export fmi2True, fmi2False
export fmi2Type, fmi2TypeModelExchange, fmi2TypeCoSimulation
export fmi2DatatypeVariable
export fmi2ScalarVariable
export fmi2StatusKind, fmi2StatusKindDoStepStatus, fmi2StatusKindPendingStatus, fmi2StatusKindLastSuccessfulTime, fmi2StatusKindTerminated
export fmi2EventInfo
export fmi2Unit, fmi2SimpleType
export fmi2DependencyKind, fmi2DependencyKindDependent, fmi2DependencyKindConstant, fmi2DependencyKindFixed, fmi2DependencyKindTunable, fmi2DependencyKindUnknown
export fmi2VariableNamingConvention, fmi2VariableNamingConventionFlat, fmi2VariableNamingConventionStructured
export fmi2VariableDependency
export fmi2DefaultExperiment, fmi2Unknown, fmi2ModelStructure
export fmi2ModelDescription
export fmi2ComponentState, fmi2ComponentStateInstantiated, fmi2ComponentStateInitializationMode, fmi2ComponentStateEventMode, fmi2ComponentStateContinuousTimeMode, fmi2ComponentStateTerminated, fmi2ComponentStateError, fmi2ComponentStateFatal

# functions
export fmi2CallbackLogger, fmi2CallbackAllocateMemory, fmi2CallbackFreeMemory, fmi2CallbackStepFinished
export fmi2Instantiate, fmi2FreeInstance!, fmi2GetTypesPlatform, fmi2GetVersion
export fmi2SetDebugLogging, fmi2SetupExperiment, fmi2EnterInitializationMode, fmi2ExitInitializationMode, fmi2Terminate, fmi2Reset
export fmi2GetReal!, fmi2SetReal, fmi2GetInteger!, fmi2SetInteger, fmi2GetBoolean!, fmi2SetBoolean, fmi2GetString!, fmi2SetString
export fmi2GetFMUstate!, fmi2SetFMUstate, fmi2FreeFMUstate!, fmi2SerializedFMUstateSize!, fmi2SerializeFMUstate!, fmi2DeSerializeFMUstate!
export fmi2GetDirectionalDerivative!, fmi2SetRealInputDerivatives, fmi2GetRealOutputDerivatives!
export fmi2DoStep, fmi2CancelStep, fmi2GetStatus!, fmi2GetRealStatus!, fmi2GetIntegerStatus!, fmi2GetBooleanStatus!, fmi2GetStringStatus!
export fmi2SetTime, fmi2SetContinuousStates, fmi2EnterEventMode, fmi2NewDiscreteStates!, fmi2EnterContinuousTimeMode, fmi2CompletedIntegratorStep!
export fmi2GetDerivatives!, fmi2GetEventIndicators!, fmi2GetContinuousStates!, fmi2GetNominalsOfContinuousStates!

# FMI3.jl
export FMU3, FMU3Component
export fmi3StatusToString

# FMI3_c.jl

# data types
export fmi3ValueReference, fmi3ModelDescription, fmi3Instance, fmi3ModelVariable, fmi3DatatypeVariable
export fmi3False, fmi3True
export fmi3Status, fmi3StatusOK, fmi3StatusWarning, fmi3StatusDiscard, fmi3StatusError, fmi3StatusFatal
export fmi3Causality, fmi3CausalityParameter, fmi3CausalityCalculatedParameter, fmi3CausalityInput, fmi3CausalityOutput, fmi3CausalityLocal, fmi3CausalityIndependent, fmi3CausalityStructuralParameter
export fmi3Initial, fmi3InitialExact, fmi3InitialApprox, fmi3InitialCalculated
export fmi3Variability, fmi3VariabilityConstant, fmi3VariabilityFixed, fmi3VariabilityTunable, fmi3VariabilityDiscrete, fmi3VariabilityContinuous
export fmi3Type, fmi3TypeModelExchange, fmi3TypeCoSimulation, fmi3TypeScheduledExecution
export fmi3Clock
export fmi3IntervalQualifier, fmi3IntervalQualifierIntervalChanged, fmi3IntervalQualifierIntervalNotYetKnown, fmi3IntervalQualifierIntervalUnchanged
export fmi3Char, fmi3String, fmi3Boolean, fmi3Binary, fmi3Float32, fmi3Float64, fmi3Int8, fmi3UInt8, fmi3Int16, fmi3UInt16, fmi3Int32, fmi3UInt32, fmi3Int64, fmi3UInt64, fmi3Byte, fmi3FMUState, fmi3InstanceEnvironment, fmi3Enum
export fmi3DependencyKind, fmi3DependencyKindDependent, fmi3DependencyKindIndependent, fmi3DependencyKindConstant, fmi3DependencyKindTunable, fmi3DependencyKindDiscrete, fmi3DependencyKindFixed, fmi3DependencyKindGetVariableDependencies

# functions
# ToDo

### EXPORTING LISTS END ###

end # module


text:
"""
Gives a reward of +1 if sp is the target lane
"""
type TargetLaneReward <: AbstractMLRewardModel
    target_lane::Int
end
function reward{D<:AbstractMLDynamicsModel}(mdp::MLMDP{MLState, MLAction, D, TargetLaneReward},
          s::MLState,
          ::MLAction,
          sp::MLState)
    return isnull(s.terminal) && sp.cars[1].y == mdp.rmodel.target_lane
end

function reward{D<:AbstractMLDynamicsModel}(mdp::MLPOMDP{MLState, MLAction, MLPhysicalState, D, TargetLaneReward},
          s::MLState,
          ::MLAction,
          sp::MLState)
    return isnull(s.terminal) && sp.cars[1].y == mdp.rmodel.target_lane
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
using WaterLily
using BenchmarkTools
using LinearAlgebra: norm2
using Profile

radius = 8; Re = 250
body = AutoBody((x,t)-> norm2(x .- 2radius) - radius)
body_sim = ()->Simulation((6radius+2,4radius+2),[1.,0.],radius; body, ν=radius/Re)

L = 2^5
function uλ(i,vx)
    x,y,z = @. (vx-1.5)*π/L              # scaled coordinates
    i==1 && return -sin(x)*cos(y)*cos(z) # u_x
    i==2 && return  cos(x)*sin(y)*cos(z) # u_y
    return 0.                            # u_z
end
vort_sim = ()->Simulation((L+2,L+2,L+2),zeros(3),L;uλ,ν=L/Re,U=1.)

function benchmark_step()
    @benchmark mom_step!($sim.flow,$sim.pois)
end

benchmark_vort() = @benchmark sim_step!(sim,10) setup=(sim=vort_sim()) seconds=30
benchmark_body() = @benchmark sim_step!(sim,10) setup=(sim=body_sim())


text:
# ----- Number -----
typealias AbstractUniform Union{DiscreteUniform, Uniform}

immutable NumberStrategy{D <: AbstractUniform, T} <: AbstractStrategy
    dist::D
    specialvalues::Vector{T}
    p::Bernoulli
    shrinker::NoShrinker
end
Base.minimum(strategy::NumberStrategy) = dist.params(strategy.dist)[1]
shrinker(strategy::NumberStrategy) = strategy.shrinker
returntype{D, T}(::Type{NumberStrategy{D, T}}) = T

# ----- Float -----
function _specialvalues{T <: AbstractFloat}(::Type{T}, min::T, max::T)
    specialvalues = [convert(T, NaN)]
    append!(specialvalues, [min, max])
    if (min < zero(T)) && (max > zero(T))
        push!(specialvalues, zero(T))
    end
    if (min < one(T)) && (max > one(T))
        push!(specialvalues, one(T))
    end
    if (min < -one(T)) && (max > -one(T))
        push!(specialvalues, -one(T))
    end
    return specialvalues
end

"""
Generate arbitrary float values between a desired range
"""
function floats{T <: AbstractFloat}(
    min::T, max::T;
    specialvalues::Vector{T}=_specialvalues(T, min, max),
    p::Bernoulli=Bernoulli(0.1),
    shrinker=NoShrinker())

    NumberStrategy(Uniform(min, max), specialvalues, p, shrinker)
end

# ----- Int -----
function _specialvalues{T <: Signed}(::Type{T}, min::T, max::T)
    specialvalues = [min, max]
    if (min < zero(T)) && (max > zero(T))
        push!(specialvalues, zero(T))
    end
    if (min < one(T)) && (max > one(T))
        push!(specialvalues, one(T))
    end
    if (min < -one(T)) && (max > -one(T))
        push!(specialvalues, -one(T))
    end
    return specialvalues
end
"""
Generate arbitrary signed values between a desired range
"""
function signeds{T <: Signed}(
    min::T, max::T;
    specialvalues::Vector{T}=_specialvalues(T, min, max),
    p::Bernoulli=Bernoulli(0.1),
    shrinker=NoShrinker())

    NumberStrategy(DiscreteUniform(min, max), specialvalues, p, shrinker)
end

function _specialvalues{T <: Unsigned}(::Type{T}, min::T, max::T)
    specialvalues = [min, max]
    if (min < one(T)) && (max > one(T))
        push!(specialvalues, one(T))
    end
    return specialvalues
end
"""
Generate arbitrary unsigned values between a desired range
"""
function unsigneds{T <: Unsigned}(
    min::T, max::T;
    specialvalues::Vector{T}=_specialvalues(T, min, max),
    p::Bernoulli=Bernoulli(0.1),
    shrinker=NoShrinker())

    NumberStrategy(DiscreteUniform(min, max), specialvalues, p, shrinker)
end

function arbitrary{D, T}(strategy::NumberStrategy{D, T})
    if Bool(rand(strategy.p))
        if length(strategy.specialvalues) > 0
            return convert(T, rand(strategy.specialvalues))
        end
    end
    return convert(T, rand(strategy.dist))
end
shrink{T <: Number}(x::T, strategy::NumberStrategy) = Task(() -> nothing)

# ----- Booleans -----
immutable BoolStrategy <: AbstractStrategy
    bernoulli::Bernoulli
end
BoolStrategy(p::Float64) = BoolStrategy(Bernoulli(p))

bools(; p=0.5) = BoolStrategy(p)

arbitrary(strategy::BoolStrategy) = rand(strategy.bernoulli)
shrink(x::Bool, strategy::BoolStrategy) = Bool[]

max_stars_repo_name: cpritcha/PropertyCheck

max_issues_repo_name: cpritcha/PropertyCheck

max_forks_repo_name: cpritcha/PropertyCheck


text:
using Distributed
if nprocs() == 1
    addprocs(3)
end

mat = Dict()
println("----------------alskdjflaskfjasklf laksjflaskj-=-------")
mat["5x5x5"] = [1,2,3]
mat["15x15x15"] = [5,6]

matice = @spawn mat;
# mat[1] = 7
# matice = @spawn mat;

@everywhere function used_in_smallfunc(par1, par2)
  println(par1, par2, " ")
  return par1 + par2
end

function bigfunc(iterator)
  localvar = 2.5
  pars = (1, 15)

  smallfunc = function (elm)
    # hypothetical (expensive) operation with elm
    # (in my code this is solving a linear system)
    vysledek = localvar * elm
    global matice

    # println("$localvar * $elm = $vysledek ")
    # matice
    # print("pred fetch $matice")
    # fetch(matice);
    # print("po fetch $matice")
    local_matrix_storage = fetch(matice)
    print(local_matrix_storage)
    if local_matrix_storage["5x5x5"][1] == 1
      local_matrix_storage["5x5x5"][1] = 2
      matice = @spawn local_matrix_storage
    else
      print(" else...")
    end
    # used_in_smallfunc(pars...)  # LoadError
    used_in_smallfunc(pars[1], pars[2])
  end

  pmap(smallfunc, iterator)
end

bigfunc(1:10)


text:
include("../src/DataFile.jl")
include("../src/SolveBalances.jl")
using Debug

# =============================================================================== #
# Solve.jl
# Function library to solve the EMT model for different experimental/test cases
# Author: J.Varner
#
# Version: 1.0
# Date: 12/17/15
# =============================================================================== #

# addition indices
addition_start,addition_end = 85,100

# mRNA / transcriptional indices
mdeg_start,mdeg_end = 101,123
script_start,script_end = 124,146

# protein / translational indices
pdeg_start,pdeg_end = 39,78
lation_start,lation_end = 147,169

# Transcriptional mult
tscript_mult = 0.2 # 0.033

#@debug
# Run the model to steady state -
 function EstimateSteadyState(current_data_dictionary)
# WARNING UNLESS DEEPCOPY IS USED, CURRENT_DATA_DICTIONARY IS THE SAME OBJECT AS DATA_DICIONARY!
  # Do we have a data_dictionary (if not load a new one)
  if (isempty(current_data_dictionary))
    # DB 02/27/2016
    TSTART = 0.0
    Ts=1.0
    TSTOP=100
    # DB
    current_data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  #test = data_dictionary["INITIAL_CONDITION_ARRAY"]
  #println(data_dictionary)
  include("../src/edit_file_global.jl")
  UpdateArray(current_data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(current_data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)

  rate_constant_array = copy(current_data_dictionary["RATE_CONSTANT_ARRAY"]);
  rate_constant_array[script_start:script_end]=tscript_mult*rate_constant_array[script_start:script_end]*const_transcription
  rate_constant_array[lation_start:lation_end]=rate_constant_array[lation_start:lation_end]*const_translation
  #@show rate_constant_array[140:162]
 # Just try to adjust dosage level by similiar amount as protein (?)
  rate_constant_array[addition_start:addition_end]=10*rate_constant_array[addition_start:addition_end]#*const_translation

# factor of 10 just corrects for a 0.1 "gamma/delta" value in the data file
  rate_constant_array[mdeg_start:mdeg_end]=10*rate_constant_array[mdeg_start:mdeg_end]*degradation_constant_mRNA
  rate_constant_array[pdeg_start:pdeg_end]=10*rate_constant_array[pdeg_start:pdeg_end]*degradation_constant_protein
  current_data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array

# 137 - 139
# rate_constant_array[137:139]=30*rate_constant_array[137:139]*degradation_constant_mRNA

  #UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)
  # Setup loop -
  #TSIM_all = []
  #X_all    = Array(Float64,length(data_dictionary["INITIAL_CONDITION_ARRAY"]))
  # tmpI = 0 # save
  EPSILON = 1e-3;
  TSTART = 0.0;
  Ts = 1.0;
  TSTOP = 100;
  did_reach_steady_state = false
  #@show current_data_dictionary["RATE_CONSTANT_ARRAY"][141] #here

  while (!did_reach_steady_state)
    # println("ITERATION",)
    # solve the balances -
  #  @show current_data_dictionary["RATE_CONSTANT_ARRAY"][141]

    (TSIM,X1) = SolveBalances(TSTART,TSTOP,Ts,current_data_dictionary);
  #  @show current_data_dictionary["RATE_CONSTANT_ARRAY"][141]
    ###save
    # if tmpI == 0
    #   X_all = X1
    # else
    # X_all=vcat(X_all,X1) # Save
    # end
    # tmpI = tmpI+1
    # Take a few additional steps -
    TNEXT_START = TSTOP+Ts;
    TNEXT_STOP = TNEXT_START+1.0;
    Ts = 0.1;

    # solve the balances again 0

    initial_condition_array = vec(X1[end,:])
    current_data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

    (TSIM,X2) = SolveBalances(TNEXT_START,TNEXT_STOP,Ts,current_data_dictionary);
    # X_all=vcat(X_all,X2) # Save
    # view simulation trajectory

    # Find the difference -
    DIFF = norm((X2[end,:] - X1[end,:]));
    # Take the end, the last step is the one we compare
    # Should we stop -or- go around again?

# println(size(X1))
# println(size(X2))
#@bp

    if (DIFF<EPSILON)
      did_reach_steady_state = true;
      # savepath = "../output/all_XSS.dat"
      # println("Save",size(X_all))
      # writedlm(savepath,X_all)# save
      return (X2[end,:]);

    else

      # No, we did *not* reach steady state ....
      TSTART = TSTOP+Ts
      TSTOP = 1000 + TSTART;
      Ts = 100.0;

      initial_condition_array = vec(X2[end,:])
      current_data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
      # println("else")
    end
  end
  # return
  # println("final return")
  return XSS;

end

@debug function SolveModelGeneric_ram(TSTART,TSTOP,Ts,data_dictionary_active,sim_RCA)

  include("../src/edit_file_global.jl")
  UpdateArray(data_dictionary_active,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)

  ### THIS WAS MISSING BECAUSE IT WAS REPLACED BY SIM_RCA... BUT NEEDS TO COME BEFORE IT
  UpdateArray(data_dictionary_active,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  ### IT HAD BEEN COMPESATED FOR BY THE MEMORY LEAK FROM THE STEADYSTATE FUNCTION ..

  UpdateArray(data_dictionary_active,"RATE_CONSTANT_ARRAY",sim_RCA)
  # @show sim_RCA
  # @show data_dictionary_active["RATE_CONSTANT_ARRAY"]
  # @bp
  # Run the simulation -
  rate_constant_array = data_dictionary_active["RATE_CONSTANT_ARRAY"];
  rate_constant_array[script_start:script_end]=tscript_mult*rate_constant_array[script_start:script_end]*const_transcription
  rate_constant_array[lation_start:lation_end]=rate_constant_array[lation_start:lation_end]*const_translation
  #@show rate_constant_array[140:162]
 # Just try to adjust dosage level by similiar amount as protein (?)
  rate_constant_array[addition_start:addition_end]=10*rate_constant_array[addition_start:addition_end]#*const_translation

# factor of 10 just corrects for a 0.1 "gamma/delta" value in the data file
  rate_constant_array[mdeg_start:mdeg_end]=10*rate_constant_array[mdeg_start:mdeg_end]*degradation_constant_mRNA
  rate_constant_array[pdeg_start:pdeg_end]=10*rate_constant_array[pdeg_start:pdeg_end]*degradation_constant_protein
  #@show rate_constant_array # inspect scaled array
  data_dictionary_active["RATE_CONSTANT_ARRAY"] = rate_constant_array

# 137 - 139
# rate_constant_array[137:139]=30*rate_constant_array[137:139]*degradation_constant_mRNA

#

  #@bp
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary_active);
  #@show data_dictionary_active["RATE_CONSTANT_ARRAY"][141]
  #  @show data_dictionary_active["RATE_CONSTANT_ARRAY"][182]

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

function SolveModelGenericPOETS(TSTART,TSTOP,Ts,local_data_dictionary,simulation)

  include(simulation)
  # !!!! I think th initial condition array update is not required
  #UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(local_data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)

  # Run the simulation -
  # println(TSTOP)
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,local_data_dictionary);

  # return tuple -
  return (TSIM,X)

end

#  #
#  #
#  #
#  #

function SolveModelGenericPOETS_ram(TSTART,TSTOP,Ts,local_data_dictionary,sim_RCA)
  # Set simulation rate constants
  UpdateArray(local_data_dictionary,"RATE_CONSTANT_ARRAY",sim_RCA)
  # Scale parameters
  rate_constant_array = local_data_dictionary["RATE_CONSTANT_ARRAY"];
  rate_constant_array[script_start:script_end]=tscript_mult*rate_constant_array[script_start:script_end]*const_transcription
  rate_constant_array[lation_start:lation_end]=rate_constant_array[lation_start:lation_end]*const_translation
  #@show rate_constant_array[140:162]
 # Just try to adjust dosage level by similiar amount as protein (?)
  rate_constant_array[addition_start:addition_end]=10*rate_constant_array[addition_start:addition_end]#*const_translation

# factor of 10 just corrects for a 0.1 "gamma/delta" value in the data file
  rate_constant_array[mdeg_start:mdeg_end]=10*rate_constant_array[mdeg_start:mdeg_end]*degradation_constant_mRNA
  rate_constant_array[pdeg_start:pdeg_end]=10*rate_constant_array[pdeg_start:pdeg_end]*degradation_constant_protein
  local_data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_array
  # Run the simulation -
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,local_data_dictionary);

  # return tuple -
  return (TSIM,X)

end

# Simulate the addition of TGFB -
function SolveModelAddTGFB(TSTART,TSTOP,Ts,data_dictionary)

  # Do we have a data_dictionary (if not load a new one)
  if (isempty(data_dictionary))
    data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  # Setup simulation time -
  # TSIM = collect(TSTART:Ts:TSTOP); # redundant
  include("../src/edit_file_global.jl")

  UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  #UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)

  #println(data_dictionary)
  # Run the model to steady-state *before* we add the TGFB -
  XSS = EstimateSteadyState(data_dictionary);

  # ok, we established a steady state -
  # reset the initial conditions -
  initial_condition_array = vec(XSS);
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  # Add TGFB1/2 -
  # RESUME TGFB signaling for simulation
  include("/Users/dbassen/work/KwateeServer-v1.0/src/edit_file_ss.jl")
  UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  # Run the simulation -
  #println(TSTOP)
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary);

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

# Simulate the addition of VEGF -
function SolveModelAddVEGF(TSTART,TSTOP,Ts,data_dictionary)

  # Do we have a data_dictionary (if not load a new one)
  if (isempty(data_dictionary))
    data_dictionary = DataFile(TSTART,TSTOP,Ts)
  end

  # Setup simulation time -
  # TSIM = collect(TSTART:Ts:TSTOP); # redundant

  # Modify rate constant array

  include("/Users/dbassen/work/KwateeServer-v1.0/src/edit_file_global.jl")
  UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  #UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)

  # Run the model to steady-state *before* we add the TGFB -
  XSS = EstimateSteadyState(data_dictionary);

  # ok, we established a steady state -
  # reset the initial conditions -
  initial_condition_array = vec(XSS);
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  # Add VEGF -
  include("/Users/dbassen/work/KwateeServer-v1.0/src/edit_file_vegf.jl")
  UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)

    # Run the simulation -
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary);

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

function SolveModelGeneric(TSTART,TSTOP,Ts,data_dictionary_active,simulation)

  # Do we have a data_dictionary (if not load a new one)
  # if (isempty(data_dictionary_active))
  #   data_dictionary_active = DataFile(TSTART,TSTOP,Ts)
  # end

  # Setup simulation time -
  # TSIM = collect(TSTART:Ts:TSTOP); # redundant

  include("../src/edit_file_global.jl")
  UpdateArray(data_dictionary_active,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  #UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  #UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)

  #println(data_dictionary)
  # Run the model to steady-state *before* we add the TGFB -
  #XSS = EstimateSteadyState(data_dictionary);
#  println(data_dictionary["INITIAL_CONDITION_ARRAY"])
  # ok, we established a steady state -
  # reset the initial conditions -
  #initial_condition_array = vec(XSS);
  #data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;

  # Add TGFB1/2 -
  # RESUME TGFB signaling for simulation
  include(simulation)
  #UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
  UpdateArray(data_dictionary_active,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
  #println(data_dictionary["INITIAL_CONDITION_ARRAY"]) ###
  # Run the simulation -
  # println(TSTOP)
  (TSIM,X) = SolveBalances(TSTART,TSTOP,Ts,data_dictionary_active);

  # return tuple -
  return (TSIM,X,initial_condition_array)

end

#