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

#=
F: read the file containing constants
I:
O:
BA
=#
function py2_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 py2_build_kinetics_buffer(all_species_dict::Dict{String, Int},
  all_rnx_list::Array, all_txtl_dict::Dict, sys2user::Dict)
  kinetics = "from __future__ import division" *
             "\nimport numpy as np" *
             "\n\ndef 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 = np.zeros($totalRnxNo)"
  for (index, rnx) in enumerate(all_rnx_list)  # go thru every rnx
    tmp_rnx_rate = "kcat[$(index - 1)]"
    if isdefined(rnx, :catalysts)  # enzyme
      for token in rnx.catalysts
        tmp_rnx_rate *= "*X[$(all_species_dict[token.oriBioName] - 1)]"
      end
    end
    if isdefined(rnx,:reactants)  # substrate
      for token in rnx.reactants
        tmp_rnx_rate *= "*X[$(all_species_dict[token.oriBioName] - 1)]/(MonodK[\"MonodK~rnx$(index)~$(token.oriBioName)\"]+X[$(all_species_dict[token.oriBioName] - 1)])"
        push!(MonodAffinityConstant_String_Array, "MonodK~rnx$(index)~$(token.oriBioName)")
      end
    end
    kinetics *= "\n\trnx_rate_vector[$(index - 1)] = "*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 = {}  # control term: string -> float" *
             "\n\tTX_rate_vector = {}  # transcription rate: string -> float" *
             "\n\tTL_rate_vector = {}  # translation rate: string -> float"
  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] - 1)]**coop / ((disassociation_constant_dict[\"KD~$(token.oriBioName)\"])**coop+X[$(all_species_dict[token.oriBioName] - 1)]**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] - 1)]**coop/((disassociation_constant_dict[\"KD~$(token.oriBioName)\"])**coop+X[$(all_species_dict[token.oriBioName] - 1)]**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] - 1)]/(saturationTL+X[$(all_species_dict[key] - 1)])"
  end

  kinetics *= "\n\n\treturn rnx_rate_vector, TX_rate_vector, TL_rate_vector"

  return (kinetics, MonodAffinityConstant_String_Array, W_string_array, disassociation_const_string_array)
end

function py2_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 = "from __future__ import division" *
    "\nimport numpy as np" *
    "\nimport math"
  buffer *= "\n\ndef generate_model_parameters_dictionary():"
  # stoichiometry
  buffer *= "\n\t# load stoichiometry" *
    "\n\twith open(\"./stoichiometry.dat\", \"r\") as file:" *
    "\n\t\tmatrix = [[float(s) for s in line.strip().split()] for line in file]" *
    "\n\t\tprint(np.shape(matrix))" *
    "\n\tstoichiometric_matrix = np.matrix(matrix)" *
    "\n\tprint(np.shape(matrix))"

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

  # all species dictionary
  buffer *= "\n\n\t# all species dictionary"
  buffer *= "\n\tall_species_dict = {}  # 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 - 1)"
    #all_species_dict[st] = index
  end
  # all species dictionary_reversed
  buffer *= "\n\n\t# all species dictionary_reversed"
  buffer *= "\n\tall_species_reversed_dict = {}  # Int --> String"
  buffer *= "\n\tfor key, val in all_species_dict.items():"
  buffer *= "\n\t\tall_species_reversed_dict[val] = key"
  # all rnx species array
  buffer *= "\n\n\t# all rnx species array"
  buffer *= "\n\trnx_species_array = ["
  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 = ["
  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 = ["
  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 = 1e-3*np.ones($(length(all_rnx_list)))  # kcat[#reaction]: reaction name"
  for (id, rnx) in enumerate(all_rnx_list)
    buffer*= "\n\tkcat_signaling[$(id - 1)] = 1.0  # 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 = {}"
  for MK in Monod_affinity_constant_array
    buffer *= "\n\tMonodAffinityConstantDict[\"$MK\"] = 1.0"
  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 = {}"
  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 = {}"

  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 = {}"
  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 = {}"
  for mrna in mRNA_species_array
    buffer *= "\n\tbackgroundGeneExpressionControlTermDict[\"$mrna\"] = 0.1  # $mrna"
  end
  # translation specific correction factor
  buffer *= "\n\n\t# translation specific correction factor"
  buffer *= "\n\ttranslationSpecificCorrectionFactor = {}"
  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 \n"
  if host_type == "bacteria"
    buffer *= replace(replace(py2_include_constants_from_literature(
      joinpath(Base.@__DIR__, "txtl_constants_ecoli.jl"),"\n\t"),
      ")^3" => ")**3"), "log(" => "math.log(")
  else
    buffer *= replace(replace(py2_include_constants_from_literature(
      joinpath(Base.@__DIR__, "txtl_constants_hl60.jl"),"\n\t"),
      ")^3" => ")**3"), "log(" => "math.log(")
  end

  #---------------------------------
  # put all stuff in a Dictionary
  buffer *= "\n\n\t####################################"
  buffer *= "\n\t# put all stuff in a Dictionary"
  buffer *= "\n\tdataDictionary = {}  # 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"

  return buffer
end

function py2_build_simulation_buffer(NoExtracellularSpecies::Int64)
  # buffer = build_copyright_header_buffer()
  buffer = "from __future__ import division" *
    "\nfrom Kinetics import calculate_kinetics" *
    "\nimport numpy as np" *
  "\n# set up ODE, get the derivatives
def Balances(y, t, dataDictionary):
  # correct for negatives
  y = [x if x > 0 else 0 for x in y]

  # 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 = np.zeros(len(all_species_dict))

  # signaling reaction network
  X2 = np.zeros(len(rnx_species)-$(NoExtracellularSpecies))  # species inside the cell
  for id, st in enumerate(rnx_species[$(NoExtracellularSpecies):]):
    # get X2 value from y
    X2[id] = y[all_species_dict[st]]

  dX1 = np.dot(stoichiometry[0:$(NoExtracellularSpecies), :], rnx_rate)  # extracellular metabolites
  dX2 = np.dot(stoichiometry[$(NoExtracellularSpecies): , :], rnx_rate) - specific_growth_rate * X2  # INTRA
  # print(np.shape(dX1), np.shape(dX2))
  dX = np.concatenate((dX1, dX2), axis=1)
  # print(np.shape(dX))
  for id, st in enumerate(rnx_species):
    # return derivatives back to dydt
    dydt[all_species_dict[st]] = dX[0, id]

  dydt[0] = specific_growth_rate * y[0]  # for Biomass

  # TXTL network
  for i in range(len(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]

  return dydt"

  return buffer
end

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

  buffer = "from __future__ import division" *
    "\nimport scipy.integrate as integrate" *
    "\nimport matplotlib.pyplot as plt" *
    "\nimport numpy as np" *
    "\nimport math" *
    "\nfrom DataDictionary import generate_model_parameters_dictionary" *
    "\nfrom Kinetics import calculate_kinetics" *
    "\nfrom Balances import Balances" *
    "\n" *
    "\ndef SolveBalances(TStart, TStop, TStep, dataDictionary):" *
    "\n\ttimeNum = int((TStop - TStart)/TStep)" *
    "\n\tt = np.linspace(TStart, TStop, timeNum)"

  # data dataDictionary
  buffer *= "\n\tall_species_reversed_dict = dataDictionary[\"all_species_reversed_dict\"]"
  buffer *= "\n\tinitial_condition = dataDictionary[\"initial_condition\"]"
  buffer *= "\n\ty0 = dataDictionary[\"initial_condition\"]"
  # run simulation
  buffer *= "\n\t# Calling the solver"
  buffer *= "\n\ty = integrate.odeint(Balances, y0, t, args=(dataDictionary,))"
  # plotting
  buffer *= "\n\n\t# plotting results" *
    "\n\tcol = len(y0)" *
    "\n\tsubplt_col = math.ceil(col/4)" *
    "\n\tfor i in range(col):" *
    "\n\t\tplt.subplot(4, subplt_col, i+1)" *
    "\n\t\tplt.plot(t, y[:, i])" *
    "\n\t\tplt.title(all_species_reversed_dict[i])" *
    "\n\tplt.show()"
  buffer *= "\n\treturn t, y"

  # 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

# FBA data dictionary generation
function py2_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 = "from __future__ import division" *
    "\nimport numpy as np" *
    "\nimport math"
  buffer *= "\n\ndef DataDictionary(time_start,time_stop,time_step):"
  # stoichiometry
  buffer *= "\n\n\t# load stoichiometry" *
    "\n\twith open(\"./stoichiometry.dat\", \"r\") as file:" *
    "\n\t\tmatrix = [[float(s) for s in line.strip().split()] for line in file]" *
    # "\n\t\tprint(np.shape(matrix))" *
    "\n\tstoichiometric_matrix = np.matrix(matrix)"
    # "\n\tprint(np.shape(matrix))"

  # Setup default flux bounds array
  buffer *= "\n\t# Setup default flux bounds array"
  buffer *= "\n\tdefault_flux_bounds_array = ("
  for (id, rnx) in enumerate(all_rnx_list)
    buffer *= "\n\t\t(0, 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-1))
    end
  end
  buffer *= "\n\t]"
  # 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 = {}\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_flux_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[\"extra_species_num\"] = $(extra_species_num)\n"
  buffer *= "\t# =============================== DO NOT EDIT ABOVE THIS LINE ============================== #\n"
  buffer *= "\treturn data_dictionary\n"

  # maximize_product_dictionary
  buffer *= "\n\n\n# maximize_product_dictionary"
  buffer *= "\ndef 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"

  return buffer
end


text:
using Soss

simple_poisson = @model traps begin
    α ~ Normal(log(4), .1)
    β ~ Normal(-0.25, .1)

    complaints ~ For(eachindex(traps)) do i
        λ = exp(α + β * traps[i])
        Poisson(λ)
    end
end

alt_poisson = @model (N, mean_traps) begin
    traps ~ For(1:N) do i
        traps[i] ~ Poisson(mean_traps)
    end

    complaints ~ For(eachindex(traps)) do i
        λ = exp(α + β * traps[i])
        Poisson(λ)
    end
end

ext: jl


text:
function writeNodePopData( popReport::DataFrame, mpSim::MPsim,
    timeFactor::Float64, reportGenerationTime::Float64, filename::String,
    overwrite::Bool )

    tStart = now()
    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )

    XLSX.openxlsx( filename, mode = overwrite ? "w" : "rw" ) do xf
        dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
        ws = xf[1]
        nReports = 1

        # Prep the new sheet.
        if overwrite
            XLSX.rename!( ws, string( "Report 1  ", dateStr ) )
        else
            nReports = length( filter( name -> startswith( name, "Report" ),
                XLSX.sheetnames( xf ) ) ) + 1
            ws = XLSX.addsheet!( xf, string( "Report ", nReports, "  ",
                dateStr ) )
        end  # if overwrite

        # Sheet header
        ws["A1"] = "Simulation length"
        ws["B1"] = mpSim.simLength / timeFactor
        ws["C1"] = "years"
        ws["A2"] = "Data generation time"
        ws["B2"] = reportGenerationTime
        ws["C2"] = "seconds"
        ws["A3"] = "Excel generation time"
        ws["C3"] = "seconds"

        popReport[:, :timePoint] /= timeFactor
        nrows, ncols = size( popReport )
        ii = 0

        for name in names( popReport )
            ii += 1
            ws[5, ii] = string( name )

            for jj in 1:nrows
                ws[jj + 5, ii] = popReport[jj, ii]
            end  # for jj in 1:nrows
        end  # for ii in 1:ncols

        excelGenerationTime = ( now() - tStart ).value / 1000.0
        ws["B3"] = excelGenerationTime
    end  # XLSX.openxlsx( ... ) do xf

end  # writeNodePopData( popReport, mpSim, timeFactor, reportGenerationTime,
     #   filename, overwrite )

function writeNodeEvolutionData(
    popReport::Tuple{DataFrame,Dict{String,Tuple}}, mpSim::MPsim,
    timeFactor::Float64, reportGenerationTime::Float64, filename::String,
    overwrite::Bool )

    tStart = now()
    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )

    XLSX.openxlsx( filename, mode = overwrite ? "w" : "rw" ) do xf
        dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
        ws = xf[1]
        nReports = 1

        # Prep the new sheet.
        if overwrite
            XLSX.rename!( ws, string( "Report 1  ", dateStr ) )
        else
            nReports = length( filter( name -> startswith( name, "Report" ),
                XLSX.sheetnames( xf ) ) ) + 1
            ws = XLSX.addsheet!( xf, string( "Report ", nReports, "  ",
                dateStr ) )
        end  # if overwrite

        # Sheet header
        ws["A1"] = "Simulation length"
        ws["B1"] = mpSim.simLength / timeFactor
        ws["C1"] = "years"
        ws["A2"] = "Data generation time"
        ws["B2"] = reportGenerationTime
        ws["C2"] = "seconds"
        ws["A3"] = "Excel generation time"
        ws["C3"] = "seconds"

        report = popReport[1]
        report[:, :timePoint] /= timeFactor
        fluxReports = popReport[2]
        nrows, ncols = size( report )
        ii = 0

        for name in string.( names( report ) )
            ii += 1
            ws[5, ii] = ii == 1 ? "Time point" : name

            for jj in 1:nrows
                ws[jj + 5, ii] = report[jj, ii]
            end  # for jj in 1:nrows

            if ii > 1
                rws = XLSX.addsheet!( xf,
                    string( name, " flux  Report ", nReports ) )
                fluxReport = fluxReports[name]
                fluxReport[1][:, :timeStart] /= timeFactor
                fluxReport[1][:, :timePoint] /= timeFactor
                ncols = size( fluxReport[1], 2 )
                jj = 1

                # In flux breakdown.
                for name in names( fluxReport[1] )
                    rws[1, jj] = jj == 1 ? "Period start" :
                        ( jj == 2 ? "Period end" : string( name ) )

                    for kk in 1:nrows
                        rws[kk + 1, jj] = fluxReport[1][kk, name]
                    end  # for kk in 1:nrows

                    jj += 1
                end  # for jj in 1:ncols

                jj += 1

                # Out flux breakdown.
                for name in names( fluxReport[2] )[3:end]
                    rws[1, jj] = string( name )

                    for kk in 1:nrows
                        rws[kk + 1, jj] = fluxReport[2][kk, name]
                    end  # for kk in 1:nrows

                    jj += 1
                end  # for name in names( fluxReport[2] )[3:end]

                if !haskey( mpSim.compoundNodeList, name )
                    continue
                end  # if !haskey( mpSim.compoundNodeList, name )

                jj += 1

                # Within flux breakdown.
                for name in names( fluxReport[3] )[3:end]
                    jj += 1
                    rws[1, jj] = string( name )

                    for kk in 1:nrows
                        rws[kk + 1, jj] = fluxReport[3][kk, name]
                    end  # for kk in 1:nrows
                end  # for name in names( fluxReport[3] )[3:end]

                # Node composition report.
                rws = XLSX.addsheet!( xf,
                    string( name, " composition  Report ", nReports ) )
                compReport = fluxReport[4]
                compReport[:, :timePoint] /= timeFactor
                ncols = size( fluxReport[1], 2 )

                jj = 1

                for name in names( compReport )
                    rws[1, jj] = jj == 1 ? "Time point" : string( name )

                    for kk in 1:nrows
                        rws[kk + 1, jj] = compReport[kk, name]
                    end  # for kk in 1:nrows

                    jj += 1
                end  # for jj in 1:ncols
            end  # if ii > 1
        end  # for ii in 1:ncols

        excelGenerationTime = ( now() - tStart ).value / 1000.0
        ws["B3"] = excelGenerationTime
    end  # XLSX.openxlsx( ... ) do xf

end  # writeNodeEvolutionData( popReport, mpSim, timeFactor,
     #   reportGenerationTime, filename, overwrite )


text:
using Test
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil, @is_initialized, @get_package, @get_numbertype, @get_ndims, SUPPORTED_PACKAGES, PKG_CUDA, PKG_NONE, NUMBERTYPE_NONE, NDIMS_NONE
import ParallelStencil: @require, @symbols
import ParallelStencil: extract_posargs_init, extract_kwargs_init, check_already_initialized, set_initialized, is_initialized, check_initialized, set_package, set_numbertype, set_ndims
using ParallelStencil.Exceptions
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. initialization of ParallelStencil" begin
            @require !@is_initialized()
            @init_parallel_stencil($package, ComplexF32, 3)
            @testset "initialized" begin
                @test @is_initialized()
                @test @get_package() == $package
                @test @get_numbertype() == ComplexF32
                @test @get_ndims() == 3
            end;
            @testset "Data" begin
                @test @isdefined(Data)
                @test length(@symbols($(@__MODULE__), Data)) > 1
                @test Symbol("Number") in @symbols($(@__MODULE__), Data)
                @test Symbol("Array") in @symbols($(@__MODULE__), Data)
                @test Symbol("DeviceArray") in @symbols($(@__MODULE__), Data)
            end;
            @reset_parallel_stencil()
        end;
        @testset "2. initialization of ParallelStencil without numbertype" begin
            @require !@is_initialized()
            @init_parallel_stencil(package = $package, ndims = 3)
            @testset "initialized" begin
                @test @is_initialized()
                @test @get_package() == $package
                @test @get_numbertype() == NUMBERTYPE_NONE
                @test @get_ndims() == 3
            end;
            @testset "Data" begin
                @test @isdefined(Data)
                @test length(@symbols($(@__MODULE__), Data)) > 1
                @test !(Symbol("Number") in @symbols($(@__MODULE__), Data))
                @test Symbol("Array") in @symbols($(@__MODULE__), Data)
                @test Symbol("DeviceArray") in @symbols($(@__MODULE__), Data)
            end;
            @reset_parallel_stencil()
        end;
        @testset "3. Exceptions" begin
            @testset "already initialized" begin
                set_initialized(true)
                set_package(:CUDA)
                set_numbertype(Float64)
                set_ndims(3)
                @require is_initialized()
                @test_throws IncoherentCallError check_already_initialized(:Threads, Float64, 3)
                @test_throws IncoherentCallError check_already_initialized(:CUDA, Float32, 3)
                @test_throws IncoherentCallError check_already_initialized(:CUDA, Float64, 2)
                set_initialized(false)
                set_package(PKG_NONE)
                set_numbertype(NUMBERTYPE_NONE)
                set_ndims(NDIMS_NONE)
            end;
            @testset "arguments" begin
                @test_throws ArgumentError extract_posargs_init($(@__MODULE__), 99, :Float64, 3)
                @test_throws ArgumentError extract_posargs_init($(@__MODULE__), nameof($package), :Char, 3)
                @test_throws ArgumentError extract_posargs_init($(@__MODULE__), nameof($package), :Float64, 77)
                @test_throws ArgumentEvaluationError extract_posargs_init($(@__MODULE__), nameof($package), :MyType, 3)
                @test_throws ArgumentEvaluationError extract_posargs_init($(@__MODULE__), nameof($package), :Float64, :myndims)
                @test_throws ArgumentError extract_kwargs_init($(@__MODULE__), Dict(:package => 99, :numbertype => :Float64, :ndims => 3))
                @test_throws ArgumentError extract_kwargs_init($(@__MODULE__), Dict(:package => nameof($package), :numbertype => :Char, :ndims => 3))
                @test_throws ArgumentError extract_kwargs_init($(@__MODULE__), Dict(:package => nameof($package), :numbertype => :Float64, :ndims => 77))
                @test_throws ArgumentEvaluationError extract_kwargs_init($(@__MODULE__), Dict(:package => nameof($package), :numbertype => :MyType, :ndims => 3))
                @test_throws ArgumentEvaluationError extract_kwargs_init($(@__MODULE__), Dict(:package => nameof($package), :numbertype => :MyType, :ndims => :myndims))
                @test_throws ArgumentError extract_kwargs_init($(@__MODULE__), Dict(:ndims => 77))
                @test_throws ArgumentEvaluationError extract_kwargs_init($(@__MODULE__), Dict(:ndims => :myndims))
            end;
            @testset "check_initialized" begin
                @require !is_initialized()
                @test_throws NotInitializedError check_initialized()
            end;
            @reset_parallel_stencil()
         end;
    end;
)) end == nothing || true;


text:
module AVDevice
  include(joinpath(dirname(@__FILE__),"..","..","..","init.jl"))
  w(f) = joinpath(avdevice_dir, f)

  using AVUtil
  using AVFormat
  using AVCodecs
  include(w("LIBAVDEVICE.jl"))
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
using RayTraceEllipsoid, UnitfulAngles, LinearAlgebra, Test

r = Ray()
@test r.dir == Vec(1,0,0)

r = Ray(Vec(1, 2, 10), Vec(0, 0, -1))
s = Ellipsoid(Vec(1, 2, 2), Vec(5, 4, 2), Vec(0, 0, 1), 90u"°")
a = advance!(r, s)
@test norm(s.center(r.orig)) ≈ s.r[3]

r = Ray(Vec(1, 2, 10), Vec(0, 0, 1))
a = advance!(r, s)
@test a

r = Ray(Vec(1, 2, 10), Vec(0, 0, -1))
s = Ellipsoid(Vec(1, 2, 2), Vec(5, 4, 2), Vec(0, 0, -1), 90u"°")
a = advance!(r, s)
@test s.center(r.orig)[3] ≈ -s.r[3]

r = Ray(Vec(1, 2, 10), Vec(-.3, .4, -1))
dir_org = deepcopy(r.dir)
s = Ellipsoid(Vec(1, 2, 2), Vec(5, 4, 2), Vec(0, 0, -1), 90u"°")
a = advance!(r, s)
ou = OpticUnit(s, true, 1., true, "a")
bend!(r, ou.interface)
@test dir_org ≈ r.dir

r = Ray(Vec(1, 2, 10), Vec(-.3, .4, -1))
dir_org = deepcopy(r.dir)
s = Ellipsoid(Vec(1, 2, 2), Vec(5, 4, 2), Vec(0, 0, -1), 90u"°")
a = advance!(r, s)
ou = OpticUnit(s, true, 1.5, true, "a")
bend!(r, ou.interface)
@test dir_org ≠ r.dir

r = Ray(Vec(1, 2, 10), Vec(0, 0, -1))
dir_org = deepcopy(r.dir)
s = Ellipsoid(Vec(1, 2, 2), Vec(5, 4, 2), Vec(0, 0, 1), 90u"°")
ou = OpticUnit(s, false, 1., true, "a")
a = raytrace!(r, ou)
@test norm(s.center(r.orig)) ≈ s.r[3]
@test dir_org ≈ r.dir


text:
function validate_constraints(
        source::UVT1,
        target::UVT2,
        resampling::ConstrainedResampling{N}) where {
            UVT1<:Union{AbstractUncertainValueDataset, Vector{<:AbstractUncertainValue}},
            UVT2<:Union{AbstractUncertainValueDataset, Vector{<:AbstractUncertainValue}},
            N}
    if N != 2
        msg = "Both `source` and `target` are uncertain, so `resampling` must contain exactly two separate sets of sampling constraints"
        throw(ArgumentError(msg))
    end
end

function validate_constraints(
        source::Vector{<:Real},
        target::UVT,
        resampling::ConstrainedResampling{N}) where {
            UVT<:Union{AbstractUncertainValueDataset, Vector{<:AbstractUncertainValue}}, N}

    if N != 1
        msg = "Only `target` is uncertain, but $N sets of constraints were provided. Only one should be given."
        throw(ArgumentError(msg))
    end
end

function validate_constraints(
        source::UVT,
        target::Vector{<:Real},
        resampling::ConstrainedResampling{N}) where {
            UVT<:Union{AbstractUncertainValueDataset, Vector{<:AbstractUncertainValue}}, N}

    if N != 1
        msg = "Only `source` is uncertain, but $N sets of constraints were provided. Only one should be given."
        throw(ArgumentError(msg))
    end
end

function validate_constraints(
        source::Vector{<:Real},
        target::Vector{<:Real},
        resampling::ConstrainedResampling{N}) where {
            UVT<:AbstractUncertainValue, N}
    @warn "Neither `source` nor `target` are uncertain. Not constraining them."
end


text:
istainted(ctx::SparsityContext, x) = ismetatype(x, ctx, ProvinanceSet)

@inline Cassette.overdub(ctx::SparsityContext, f::typeof(istainted), x) = istainted(ctx, x)
@inline Cassette.overdub(ctx::SparsityContext, f::typeof(this_here_predicate!)) =
    this_here_predicate!(ctx)

# Must return 7 exprs
function rewrite_branch(ctx, stmt, extraslot, i)
    # turn
    #   gotoifnot %p #g
    # into
    #   %t = istainted(%p)
    #   gotoifnot %t #orig
    #   %rec = this_here_predicate!(path)
    #   gotoifnot %rec #orig+1 (the next statement after gotoifnot)

    exprs = Any[]
    cond = stmt.args[1]        # already an SSAValue

    # insert a check to see if SSAValue(i) isa Tainted
    istainted_ssa = Core.SSAValue(i)
    push!(exprs, :($istainted($cond)))

    # not tainted? jump to the penultimate statement
    push!(exprs, Expr(:gotoifnot, istainted_ssa, i+5))

    # tainted? then use this_here_predicate!(SSAValue(1))
    current_pred = i+2
    push!(exprs, :($this_here_predicate!()))

    # Store the interpreter-provided predicate in the slot
    push!(exprs, Expr(:(=), extraslot, SSAValue(i+2)))

    push!(exprs, Core.GotoNode(i+6))

    push!(exprs, Expr(:(=), extraslot, cond))

    # here we put in the original code
    stmt1 = copy(stmt)
    stmt.args[1] = extraslot
    push!(exprs, stmt)

    exprs
end

function rewrite_ir(ctx, ir)
    # turn
    #   <val> ? t : f
    # into
    #   istainted(<val>) ? this_here_predicate!(p) : <val> ? t : f

    ir = copy(ir)

    extraslot = gensym("tmp")
    push!(ir.slotnames, extraslot)
    push!(ir.slotflags, 0x00)
    extraslot = Core.SlotNumber(length(ir.slotnames))

    Cassette.insert_statements!(ir.code, ir.codelocs,
        (stmt, i) -> Base.Meta.isexpr(stmt, :gotoifnot) ? 7 : nothing,
        (stmt, i) -> rewrite_branch(ctx, stmt, extraslot, i))

    ir.ssavaluetypes = length(ir.code)
    # Core.Compiler.validate_code(ir)
    return ir
end

function rewrite_tainted_branches(ctx, ref)
    rewrite_ir(ctx, ref.code_info)
end

const BranchesPass = Cassette.@pass rewrite_tainted_branches


text:
@testset "action probabilities" begin
    let royalty_rates = [0.25,],
        geology_types = [3.5,],
        nroy = length(royalty_rates),
        ngeo = length(geology_types),
        roy = royalty_rates[1],
        geoid = geology_types[1],
        itype = (geoid, roy,),
        pids = [1,],
        dmax = ShaleDrillingModel._dmax(wp),
        rngs = (zspace...,                                        vspace, vspace, 0:dmax, 1:length(wp), geology_types, royalty_rates,),
        idxs = (2:2:length(zspace[1])-1, 2:3:length(zspace[2])-1, 1:5:nv, 1:5:nv, 0:dmax, 1:length(wp), 1:ngeo,        1:nroy,),
        idx_sz = length.((θfull, idxs...)),
        prim = dcdp_primitives(flowfuncname, β, wp, zspace, Πp, ψspace),
        tmpv = dcdp_tmpvars(prim),
        sev = SharedEV(pids, prim, rngs[end-1:end]...),
        isev = ItpSharedEV(sev, prim, σv),
        T = eltype(θfull),
        tmp = Vector{T}(undef, dmax+1),
        θ1 = similar(θfull),
        θ2 = similar(θfull),
        θttmp = Vector{Float64}(undef, prim.nθt),
        grad   = zeros(T, idx_sz),
        fdgrad = zeros(T, idx_sz)

        println("Testing logP")

        dograd = true

        zero!(sev)
        CR = CartesianIndices(length.(idxs))

        # test that things just solve...
        let uv = (1.,1.), z = getindex.(zspace, 7), dp1 = 1, s_idx = 1, itypidx = (1, 1,), tmpgrad = Vector{Float64}(undef, idx_sz[1])
            θtvw = _θt(θfull, prim.nθt)
            σ0 = _σv(θfull)
            @test idx_sz[1] == length(θtvw)+1
            @test length(tmpgrad) == idx_sz[1]
            @test length(tmpgrad) == length(θtvw)+1
            logP!(tmpgrad, tmp, θtvw, σ0, prim, isev, uv, z, dp1, s_idx, itypidx, true)
        end
        # -----------------------

        println("solving all")

        dograd = true

        serial_solve_vf_all!(sev, tmpv, prim, θfull, dograd; maxit0=40, maxit1=20, vftol=1e-10)
        println("solved round 1. doing logP")
        grad .= 0.0
        for CI in CR
            zpi, zvi, ui, vi, di, si, gi, ri = CI.I
            zp , zv , u , v , d , s , g , r  = getindex.(rngs, CI.I)

            θtvw = _θt(θfull, prim)
            σ0 = _σv(θfull)

            if d ∈ ShaleDrillingModel.actionspace(prim.wp, s) # && !(s ∈ ShaleDrillingModel.ind_lrn(prim.wp.endpts))
                @views lp = logP!(grad[:,CI], tmp, θtvw, σ0, prim, isev,  (u,v), (zp,zv,), d, s, (gi, ri,), dograd)
            end
        end

        dograd = false
        for k in 1:length(θfull)
            print("θfull[$k]...")
            θ1 .= θfull
            θ2 .= θfull
            h = peturb(θfull[k])
            θ1[k] -= h
            θ2[k] += h
            hh = θ2[k] - θ1[k]

            σ1 = _σv(θ1)
            σ2 = _σv(θ2)

            serial_solve_vf_all!(sev, tmpv, prim, θ1, dograd; maxit0=40, maxit1=20, vftol=1e-10)
            print("logp: θ[$k]-h\t")
            for CI in CR
                zpi, zvi, ui, vi, di, si, gi, ri = CI.I
                zp , zv , u , v , d , s , g , r  = getindex.(rngs, CI.I)
                θtvw1 = _θt(θ1, prim)
                if d ∈ ShaleDrillingModel.actionspace(prim.wp, s) # && !(s ∈ ShaleDrillingModel.ind_lrn(prim.wp.endpts))
                    fdgrad[k,CI] -= logP!(Vector{T}(), tmp, θtvw1, σ1, prim, isev, (u,v), (zp,zv,), d, s, (gi,ri), dograd)
                end
            end

            println("solving again for logp: θ[$k]+h....")
            serial_solve_vf_all!(sev, tmpv, prim, θ2, dograd; maxit0=40, maxit1=20, vftol=1e-10)
            println("updating fdgrad")
            for CI in CR
                zpi, zvi, ui, vi, di, si, gi, ri = CI.I
                zp , zv , u , v , d , s , g , r  = getindex.(rngs, CI.I)
                θtvw2 = _θt(θ2, prim)
                if d ∈ ShaleDrillingModel.actionspace(prim.wp, s) # && !(s ∈ ShaleDrillingModel.ind_lrn(prim.wp.endpts))
                    fdgrad[k,CI] += logP!(Vector{T}(), tmp, θtvw2, σ2, prim, isev, (u,v), (zp,zv,), d, s, (gi,ri,), dograd)
                    fdgrad[k,CI] /= hh
                end
            end
        end

        println("checking gradient")

        @test !all(fdgrad .== 0.0)
        @test !all(grad .== 0.0)
        @test all(isfinite.(fdgrad))
        @test all(isfinite.(grad))

        @test fdgrad ≈ grad
        maxv, idx =  findmax(abs.(fdgrad .- grad))
        println("worst value is $maxv at $(CartesianIndices(fdgrad)[idx]) for dlogP")
        @test maxv < 1.5e-7
    end
end

# println("With σ")
# maxv, idx =  findmax(abs.(fdgrad_vw1 .- grad_vw1))
# mae = mean(abs.(fdgrad_vw1 .- grad_vw1))
# mse = var(fdgrad_vw1 .- grad_vw1)
# med = median(abs.(fdgrad_vw1 .- grad_vw1))
# q90 = quantile(vec(abs.(fdgrad_vw1 .- grad_vw1)), 0.9)
# sub = ind2sub(grad_vw1, idx)
# vals = getindex.(rngs, sub)
# println("worst value is $maxv at $sub for dlogP. This has characteristics $vals")
# println("MAE = $mae. MSE = $mse. Median abs error = $med. 90pctile = $q90")
# @test 0.0 < maxv < 0.1

# end


text:
# Structs and functions related to the C file pointer type as used in eccodes

export GribFile, destroy

mutable struct File
end

"""
Represents a grib file. Implements many of the Julia I/O functions.
"""
mutable struct GribFile
    ptr::Ptr{File}
    filename::String
    mode::String
    pos::Int
end

function handle_fopen_errors(ptr, fname)
    if ptr == C_NULL
        # Automatically reads errno
        throw(SystemError("Failed to open $fname"))
    end
end

"""
    GribFile(filename::AbstractString, mode="r")

Open a grib file. `mode` is a mode as described by `Base.open`.
"""
function GribFile(filename::AbstractString; mode="r")
    f = ccall(:fopen, Ptr{File}, (Cstring, Cstring), filename, mode)
    handle_fopen_errors(f, filename)
    nref = Ref(Int32(0))

    GribFile(f, filename, mode, 0)
end

"""
    GribFile(f::Function, filename::AbstractString, mode="r")

Open a grib file and automatically close after exucuting `f`.
`mode` is a mode as described by `Base.open`.

# Example
```julia
GribFile(filename) do f
    # do things in read mode
end
```
"""
function GribFile(f::Function, filename::AbstractString; mode="r")
    file = GribFile(filename, mode=mode)
    try
        f(file)
    finally
        destroy(file)
    end
end

function Base.iterate(f::GribFile, state=())
    next = Message(f)
    if isnothing(next)
        nothing
    else
        (next, ())
    end
end

Base.eltype(::Type{GribFile}) = Message

Base.IteratorSize(::Type{GribFile}) = Base.SizeUnknown()

function Base.seekstart(f::GribFile)
    close(f.ptr)
    fptr = ccall(:fopen, Ptr{File}, (Cstring, Cstring), f.filename, f.mode)
    handle_fopen_errors(fptr, f.filename)
    f.ptr = fptr
    f.pos = 0
end

""" Read without allocating a return vector """
function readnoreturn(f::GribFile, nm::Integer)
    for i in 1:nm
        Message(f)
    end
end

"""
    read(f::GribFile[, nm::Integer])

Read `nm` messages from `f` and return as vector. Default is 1.
"""
function Base.read(f::GribFile, nm::Integer)
    vec = Vector{Message}(undef, nm)
    valid_count = 0
    for i in 1:nm
        msg = Message(f)
        if !isnothing(msg)
            vec[i] = msg
            valid_count += 1
        end
    end
    resize!(vec, valid_count)
    vec
end
# 0.4 TODO: This should return a Vector. Fixing it will be breaking.
Base.read(f::GribFile) = Message(f)

"""
    position(f::GribFile)

Get the current position of the file.
"""
Base.position(f::GribFile) = f.pos

"""
    seek(f::GribFile, n::Integer)

Seek the file to the given position `n`. Throws a `DomainError`
if `n` is negative.
"""
function Base.seek(f::GribFile, n::Integer)
    if n < 0
        throw(DomainError(n, "Cannot seek to before the start of the file"))
    elseif n < f.pos
        seekstart(f)
        readnoreturn(f, n)
    elseif n > f.pos
        readnoreturn(f, n-f.pos)
    end
end

"""
    skip(f::GribFile, offset::Integer)

Seek the file relative to the current position. Throws a `DomainError`
if `offset` brings the file position to before the start of the file.
"""
function Base.skip(f::GribFile, offset::Integer)
    if f.pos + offset < 0
        throw(DomainError(offset, "Cannot skip to before the start of the file"))
    elseif offset < 0
        oldpos = f.pos
        seekstart(f)
        readnoreturn(f, oldpos+offset)
    else
        readnoreturn(f, offset)
    end
end

"""
    destroy(f::GribFile)

Safely close the file.
"""
function destroy(f::GribFile)
    err = ccall(:fclose, Cint, (Ptr{File},), f.ptr)
    if err != 0
        # 0.4 TODO: Change this to SystemError
        throw(ErrorException("GribFile closed with errorcode $(LibC.errno())"))
    end
end

function close(f::Ptr{File})
    err = ccall(:fclose, Cint, (Ptr{File},), f)
    if err != 0
        # 0.4 TODO: Change this to SystemError
        throw(ErrorException("GribFile closed with errorcode $(LibC.errno())"))
    end
end

# Functions related to printing
function Base.show(io::IO, mime::MIME"text/plain", f::GribFile)
    str = "GribFile $(f.filename) at position $(f.pos) in mode $(f.mode)"
    write(io, str)
end


text:
n, m, l = parse.(Int, split(readline()))
A = zeros(Int, n, m)
B = zeros(Int, m, l)
for i in 1:n
    A[i,:] = parse.(Int, split(readline()))
end

for i in 1:m
    B[i,:] = parse.(Int, split(readline()))
end

C = A * B
for i in 1:n
    for j in 1:l
        if j == 1
            print(C[i,j])
        else
            print(" ", C[i,j])
        end
    end
    println()
end


text:
include("../src/delta.jl")
include("../src/vm_basis.jl")
include("../src/hecke.jl")
include("../src/prime_range.jl")

function test_hecke_on_qexp()
   print("hecke.qexp...")

   delta = delta_qexp(100)
   h1 = hecke_operator_on_qexp(delta, 1, 12, 30)
   @test h1 - delta_qexp(30) == 0
   h2 = truncate(hecke_operator_on_qexp(delta, 3, 12, 10), 7)
   S = parent(h2)
   q = gen(S)
   @test h2 - (252q - 6048q^2 + 63504q^3 - 370944q^4 +
               1217160q^5 - 1524096q^6 + O(q^7)) == 0

   println("PASS")
end

function test_hecke_on_basis()
   print("hecke.basis...")

   B = victor_miller_basis(32, 10, false)
   H_basis = hecke_operator_on_basis(B, 3, 32)
   b0, b1, b2 = B
   h0 = hecke_operator_on_qexp(b0, 3, 32)
   R = base_ring(h0)
   S = MatrixSpace(R, 3, 3)
   matrix = S([617673396283948 2611200 5615943999897600; 0 50220 965671206912; 0 432 17312940])

   println("PASS")
end

function test_hecke_vm()
   print("hecke.victormiller...")

   B = victor_miller_basis(24, 20, true)
   b1, b2 = B
   h1 = hecke_operator_on_qexp(b1, 3, 24, 7)
   S = parent(h1)
   q = gen(S)
   @test h1 == 195660*q - 982499328*q^2 + 85442803344*q^3 + 1302498570240*q^4 +
               23514204375720*q^5 - 333538871869440*q^6 + O(q^7)
   h2 = hecke_operator_on_qexp(b2, 3, 24, 7)
   S = parent(h2)
   q = gen(S)
   @test h2 == -48*q + 143820*q^2 - 16295040*q^3 - 424520544*q^4 -
               4306546080*q^5 + 67844160144*q^6 + O(q^7)

   #test hecke_operator_on_basis
   H_basis = hecke_operator_on_basis(B, 3, 24)
   R = base_ring(h1)
   S = MatrixSpace(R, 2, 2)
   matrix = S([195660 -982499328; -48 143820])

   println("PASS")
end

function test_hecke()
   test_hecke_on_qexp()
   test_hecke_on_basis()
   test_hecke_vm()

   println("")
end


text:
# -----------------------------------------------------------
#
# Methods for DG interpolation using readable
# dictionary-style structs
#
# -----------------------------------------------------------

using Cubature

const REL_TOL = 1.0e-8
const ABS_TOL = 1.0e-12
const MAX_EVALS = 1000

# -----------------------------------------------------
# Here we have all our methods for
# taking a function and returning an appropriate list
# of DG coefficients.
#
# -----------------------------------------------------

# Efficiency criticality: HIGH

# -----------------------------------------------------
# 1-D Basis
# -----------------------------------------------------

function v{T<:Real}(k::Int, level::Int, cell::Int, mode::Int, x::T)
	if level==0 # At the base level, we are not discontinuous, and we simply
				# use Legendre polynomials up to degree k-1 as a basis
		return LegendreP(mode-1,2*x-1)*sqrt(2.0)
	else
		return h(k, mode, (1<<level)*x - (2*cell-1)) * sqrt(one(T)*(1<<level))
		# Otherwise we use an appropriately shifted, scaled, and normalized
		# DG function
	end
end

function v(k::Int, level::Int, cell::Int, mode::Int)
	return (xs-> v(k,level,cell,mode,xs))
end

# -----------------------------------------------------
# Tensor Product Construction
# -----------------------------------------------------

# Returns the value of the function at x
function V{D,T<:Real}(k::Int, level::NTuple{D,Int},
	cell::CartesianIndex{D}, mode::CartesianIndex{D}, xs::Array{T,1})
	ans = one(T)
	for i = 1:D
		ans *= v(k, level[i], cell[i], mode[i], xs[i])
	end
	return ans
end
# Is there any fast way to precompute the ones I care about??

# Returns a function
function V{D}(k, level::NTuple{D,Int},
			cell::CartesianIndex{D}, mode::CartesianIndex{D})
	return (xs-> V(k, level, cell, mode, xs))
end

# -----------------------------------------------------
# Methods for obtaining the coefficients
# -----------------------------------------------------

# Given a 1-D position and level, this tells us which cell
# that position belongs to, at that level resolution
function cell_index(x::Real,l::Int)
	if l <= 1
		return 1
	end
	if x >= 1
		return 2^(l-1)
	else
		return 1 + floor(Int, 2^(l-1) * x)
	end
end

# This takes an inner product, but since for higher levels our inner product
# is only concerned with a specific region in the grid, we restrict to that
# appropriately, depending on the level
function inner_product{D}(f::Function, g::Function, lvl::NTuple{D,Int}, cell::CartesianIndex{D};
								rel_tol = REL_TOL, abs_tol = ABS_TOL, max_evals=MAX_EVALS)

	h = (x-> f(x)*g(x))
	xmin = ntuple(i-> (cell[i]-1)/(1<<max(0,lvl[i]-1)), D)
	xmax = ntuple(i-> (cell[i])/(1<<max(0,lvl[i]-1)), D)
	val = hcubature(h, xmin, xmax; reltol=rel_tol, abstol=abs_tol, maxevals=max_evals)[1]
	return val
end

# We obtain coefficients simply by doing inner products, it's easy :)
# Only hard part is inner product integrations can be slower than we want :(
function get_coefficient_DG{D}(k::Int, lvl::NTuple{D,Int},
									cell::CartesianIndex{D},
									mode::CartesianIndex{D},
									f::Function;
									rel_tol = REL_TOL,
									abs_tol = ABS_TOL,
									max_evals=MAX_EVALS)

	return inner_product(f, V(k,lvl,cell,mode),lvl,cell;
							rel_tol = rel_tol, abs_tol = abs_tol, max_evals=max_evals)
end

# -----------------------------------------------------
# Full or Sparse Galerkin Coefficients in n-D
# -----------------------------------------------------
function coeffs_DG(D::Int, k::Int, n::Int, f::Function;
										rel_tol = REL_TOL, abs_tol = ABS_TOL,
										max_evals=MAX_EVALS,
										scheme="sparse")
	cutoff	= get_cutoff(scheme, D, n)
	coeffs	= Dict{CartesianIndex{D}, Array{Array{Float64,D},D}}()
	modes	= ntuple(i-> k, D)
	ls		= ntuple(i->(n+1),D)

	for level in CartesianRange(ls) #This really goes from 0 to l_i for each i
		cutoff(level) && continue

		cells = ntuple(i -> 1<<max(0, level[i]-2), D)
		level_coeffs = Array{Array{Real,D}}(cells)
		lvl = ntuple(i -> level[i]-1,D)
		for cell in CartesianRange(cells)
			cell_coeffs = Array{Real}(modes)
			for mode in CartesianRange(modes)
				cell_coeffs[mode] = get_coefficient_DG(k, lvl, cell, mode, f)
			end
			level_coeffs[cell]=cell_coeffs
		end
		coeffs[level] = level_coeffs
	end
	return coeffs
end

# -----------------------------------------------------------
# Reconstruction (full and sparse) in n-D from a Dict of
# coefficients
# -----------------------------------------------------------

function reconstruct_DG{D,T1<:Real,T2<:Real}(coeffs::Dict{CartesianIndex{D}, Array{Array{T1,D},D}},
									xs::Array{T2,1})

	value	= zero(T2)
	k		= size(first(values(coeffs))[1])[1]
	modes	= ntuple(i-> k ,D)

	for key in keys(coeffs)
		level = ntuple(i->key[i]-1,D)
		cell = CartesianIndex{D}(ntuple(i->cell_index(xs[i],level[i]),D))
		coeff = coeffs[key][cell]::Array{T1,D}
		@inbounds for mode in CartesianRange(modes)
			value += coeff[mode]*V(k, level, cell, mode, xs)
		end
	end
	return value
end


text:
const default_warn_itensor_order = 14

"""
    get_warn_itensor_order()
Return the threshold for the order of an ITensor above which
ITensors will emit a warning.
You can set the threshold with the function `set_warn_itensor_order!(N::Int)`.
"""
get_warn_itensor_order() = ITensors.GLOBAL_PARAMS["WarnTensorOrder"]

"""
    set_warn_itensor_order!(N::Int)
After this is called, ITensor will warn about ITensor contractions
that result in ITensors above the order `N`.
This function returns the initial warning threshold (what it was
set to before this function was called).
You can get the current threshold with the function `get_warn_itensor_order(N::Int)`. You can reset to the default value with
`reset_warn_itensor_order!()`.
"""
function set_warn_itensor_order!(N::Union{Int, Nothing})
  N_init = get_warn_itensor_order()
  ITensors.GLOBAL_PARAMS["WarnTensorOrder"] = N
  return N_init
end

"""
    reset_warn_itensor_order!()
After this is called, ITensor will warn about ITensor contractions
that result in ITensors above the default order
$default_warn_itensor_order.
This function returns the initial warning threshold (what it was
set to before this function was called).
"""
reset_warn_itensor_order!() =
  set_warn_itensor_order!(default_warn_itensor_order)

"""
    disable_warn_itensor_order!()
After this is called, ITensor will not warn about ITensor
contractions that result in large ITensor orders.
This function returns the initial warning threshold (what it was
set to before this function was called).
"""
disable_warn_itensor_order!() =
  set_warn_itensor_order!(nothing)


text:
import KernelDensity.UnivariateKDE

######################################################################
# Resampling datasets consisting of uncertain values furnished by
# theoretical distributions
######################################################################
tol = 1e-7

uv = UncertainValue(Normal, 1, 0.5)

uvc_lq = constrain(uv, TruncateLowerQuantile(0.2))
uvc_uq = constrain(uv, TruncateUpperQuantile(0.8))
uvc_q = constrain(uv, TruncateQuantiles(0.2, 0.8))
uvc_min = constrain(uv, TruncateMinimum(0.5))
uvc_max = constrain(uv, TruncateMaximum(1.5))
uvc_range = constrain(uv, TruncateRange(0.5, 1.5))

@test minimum(resample(uvc_lq, 1000)) >= quantile(uv, 0.2) - tol
@test maximum(resample(uvc_uq, 1000)) <= quantile(uv, 0.8) + tol
@test all(quantile(uv, 0.2) - tol .<= resample(uvc_q, 1000) .<= quantile(uv, 0.8) + tol)
@test minimum(resample(uvc_min, 1000)) >= 0.5 - tol
@test maximum(resample(uvc_max, 1000)) <= 1.5 + tol
@test all(0.5 - tol .<= resample(uvc_range, 1000) .<= 1.5 + tol)

##########################################################################
# Resampling datasets consisting of uncertain values furnished by
# theoretical distributions with parameters estimated from empirical data
##########################################################################
tol = 1e-7

uv = UncertainValue(Normal, rand(Normal(-1, 0.2), 1000))

uvc_lq = constrain(uv, TruncateLowerQuantile(0.2))
uvc_uq = constrain(uv, TruncateUpperQuantile(0.8))
uvc_q = constrain(uv, TruncateQuantiles(0.2, 0.8))
uvc_min = constrain(uv, TruncateMinimum(0.5))
uvc_max = constrain(uv, TruncateMaximum(1.5))
uvc_range = constrain(uv, TruncateRange(0.5, 1.5))

@test minimum(resample(uvc_lq, 1000)) >= quantile(uv, 0.2) - tol
@test maximum(resample(uvc_uq, 1000)) <= quantile(uv, 0.8) + tol
@test all(quantile(uv, 0.2) - tol .<= resample(uvc_q, 1000) .<= quantile(uv, 0.8) + tol)
@test minimum(resample(uvc_min, 1000)) >= 0.5 - tol
@test maximum(resample(uvc_max, 1000)) <= 1.5 + tol
@test all(0.5 - tol .<= resample(uvc_range, 1000) .<= 1.5 + tol)

##########################################################################
# Resampling datasets consisting of uncertain values furnished by
# kernel density estimations to the distributions.
##########################################################################

# quantile estimates for `uv` is not precise for KDE estimates,
# so we need to lower the tolerance.
tol = 1e-2
uv = UncertainValue(UnivariateKDE, rand(Uniform(10, 15), 1000))
#uv = UncertainValue(rand(Normal(15, 2), 1000)) # does the same

# Verify that we get errors if trying to sample outside the support of the
# distributon furnishing the data point
@test_throws ArgumentError constrain(uv, TruncateMaximum(-100))
@test_throws ArgumentError constrain(uv, TruncateMinimum(100))
@test_throws DomainError constrain(uv, TruncateRange(100, -100))

uvc_lq = constrain(uv, TruncateLowerQuantile(0.2))
uvc_uq = constrain(uv, TruncateUpperQuantile(0.8))
uvc_q = constrain(uv, TruncateQuantiles(0.2, 0.8))
uvc_min = constrain(uv, TruncateMinimum(13))
uvc_max = constrain(uv, TruncateMaximum(13))
uvc_range = constrain(uv, TruncateRange(11, 12))

@test minimum(resample(uvc_lq, 1000)) >= quantile(uv, 0.2) - tol
@test maximum(resample(uvc_uq, 1000)) <= quantile(uv, 0.8) + tol
@test all(quantile(uv, 0.2) - tol .<= resample(uvc_q, 1000) .<= quantile(uv, 0.8) + tol)
@test minimum(resample(uvc_min, 1000)) >= 13 - tol
@test maximum(resample(uvc_max, 1000)) <= 13 + tol
@test all(11 - tol .<= resample(uvc_range, 1000) .<= 12 + tol)


text:
module PriceTools
using eve.CREST
using eve.SDE
using eve.MarketDatas

export ItemSummery, summerizeItem

immutable ItemSummery
  itemID :: Int64
  source :: DataSource
  timestamp :: DateTime
  minsell :: Float64
  maxbuy :: Float64
  margin :: Float64
  markup :: Float64
end

function summerizeItem(data :: MarketData)
  sells = sort(data.sells, cols = (:Price), rev = false)
  buys = sort(data.buys, cols = (:Price), rev = true)
  maxbuy = buys[:Price][1]
  minsell = sells[:Price][1]
  margin = (minsell - maxbuy) / minsell
  markup = (minsell - maxbuy) / maxbuy
  ItemSummery(data.itemID,
              data.source,
              data.timeFetched,
              minsell,
              maxbuy,
              margin,
              markup)
end

end


text:
function exists(dfg::LightDFG{P,V,F}, node::V) where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
    return haskey(dfg.g.variables, node.label)
end

function exists(dfg::LightDFG{P,V,F}, node::F) where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
    return haskey(dfg.g.factors, node.label)
end

exists(dfg::LightDFG, nId::Symbol) = haskey(dfg.g.labels, nId)

exists(dfg::LightDFG, node::DFGNode) = exists(dfg, node.label)

function isVariable(dfg::LightDFG{P,V,F}, sym::Symbol) where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
    return haskey(dfg.g.variables, sym)
end

function isFactor(dfg::LightDFG{P,V,F}, sym::Symbol) where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
    return haskey(dfg.g.factors, sym)
end

function addVariable!(dfg::LightDFG{<:AbstractParams, V, <:AbstractDFGFactor}, variable::V)::V where V <: AbstractDFGVariable
    #TODO should this be an error
    if haskey(dfg.g.variables, variable.label)
        error("Variable '$(variable.label)' already exists in the factor graph")
    end

    FactorGraphs.addVariable!(dfg.g, variable) || return false

    # Track insertion
    push!(dfg.addHistory, variable.label)

    return variable
end

function addVariable!(dfg::LightDFG{<:AbstractParams, VD, <:AbstractDFGFactor},
                      variable::AbstractDFGVariable)::VD where VD <: AbstractDFGVariable
    return addVariable!(dfg, VD(variable))
end

#moved to abstract
# function addFactor!(dfg::LightDFG{<:AbstractParams, V, F}, variables::Vector{<:V}, factor::F)::F where {V <: AbstractDFGVariable, F <: AbstractDFGFactor}
#
#     #TODO should this be an error
#     if haskey(dfg.g.factors, factor.label)
#         error("Factor '$(factor.label)' already exists in the factor graph")
#     end
#     # for v in variables
#     #     if !(v.label in keys(dfg.g.metaindex[:label]))
#     #         error("Variable '$(v.label)' not found in graph when creating Factor '$(factor.label)'")
#     #     end
#     # end
#
#     variableLabels = map(v->v.label, variables)
#
#     resize!(factor._variableOrderSymbols, length(variableLabels))
#     factor._variableOrderSymbols .= variableLabels
#     # factor._variableOrderSymbols = copy(variableLabels)
#
#     @assert FactorGraphs.addFactor!(dfg.g, variableLabels, factor)
#     return factor
# end
#
# function addFactor!(dfg::LightDFG{<:AbstractParams, <:AbstractDFGVariable, F}, variableLabels::Vector{Symbol}, factor::F)::F where F <: AbstractDFGFactor
#     #TODO should this be an error
#     if haskey(dfg.g.factors, factor.label)
#         error("Factor '$(factor.label)' already exists in the factor graph")
#     end
#
#     resize!(factor._variableOrderSymbols, length(variableLabels))
#     factor._variableOrderSymbols .= variableLabels
#
#     @assert FactorGraphs.addFactor!(dfg.g, variableLabels, factor)
#
#     return factor
# end

function addFactor!(dfg::LightDFG{<:AbstractParams, <:AbstractDFGVariable, F}, factor::F) where F <: AbstractDFGFactor
    if haskey(dfg.g.factors, factor.label)
        error("Factor '$(factor.label)' already exists in the factor graph")
    end
    # TODO
    # @assert FactorGraphs.addFactor!(dfg.g, getVariableOrder(factor), factor)
    @assert FactorGraphs.addFactor!(dfg.g, Symbol[factor._variableOrderSymbols...], factor)
    return factor
end

function addFactor!(dfg::LightDFG{<:AbstractParams, <:AbstractDFGVariable, F},
                      factor::AbstractDFGFactor) where F <: AbstractDFGFactor
    return addFactor!(dfg, F(factor))
end

function getVariable(dfg::LightDFG, label::Symbol)
    if !haskey(dfg.g.variables, label)
        error("Variable label '$(label)' does not exist in the factor graph")
    end

    return dfg.g.variables[label]
end

function getFactor(dfg::LightDFG, label::Symbol)
    if !haskey(dfg.g.factors, label)
        error("Factor label '$(label)' does not exist in the factor graph")
    end
    return dfg.g.factors[label]
end

function updateVariable!(dfg::LightDFG, variable::AbstractDFGVariable; warn_if_absent::Bool=true)
    if !haskey(dfg.g.variables, variable.label)
        warn_if_absent && @warn "Variable label '$(variable.label)' does not exist in the factor graph, adding"
        return addVariable!(dfg, variable)
    end
    dfg.g.variables[variable.label] = variable
    return variable
end

function updateFactor!(dfg::LightDFG, factor::AbstractDFGFactor; warn_if_absent::Bool=true)
    if !haskey(dfg.g.factors, factor.label)
        warn_if_absent && @warn "Factor label '$(factor.label)' does not exist in the factor graph, adding"
        return addFactor!(dfg, factor)
    end

    # Confirm that we're not updating the neighbors
    dfg.g.factors[factor.label]._variableOrderSymbols != factor._variableOrderSymbols && error("Cannot update the factor, the neighbors are not the same.")

    dfg.g.factors[factor.label] = factor
    return factor
end

function deleteVariable!(dfg::LightDFG, label::Symbol)#::Tuple{AbstractDFGVariable, Vector{<:AbstractDFGFactor}}
    if !haskey(dfg.g.variables, label)
        error("Variable label '$(label)' does not exist in the factor graph")
    end

    deleteNeighbors = true # reserved, orphaned factors are not supported at this time
    if deleteNeighbors
        neigfacs = map(l->deleteFactor!(dfg, l), getNeighbors(dfg, label))
    end
    variable = dfg.g.variables[label]
    rem_vertex!(dfg.g, dfg.g.labels[label])

    return variable, neigfacs
end

function deleteFactor!(dfg::LightDFG, label::Symbol; suppressGetFactor::Bool=false)::AbstractDFGFactor
    if !haskey(dfg.g.factors, label)
        error("Factor label '$(label)' does not exist in the factor graph")
    end
    factor = dfg.g.factors[label]
    variable = rem_vertex!(dfg.g,  dfg.g.labels[label])
    return factor
end

function getVariables(dfg::LightDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)

    # variables = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGVariable, vertices(dfg.g)))
    variables = collect(values(dfg.g.variables))
    if regexFilter !== nothing
        variables = filter(v -> occursin(regexFilter, String(v.label)), variables)
    end
    if solvable != 0
        variables = filter(v -> _isSolvable(dfg, v.label, solvable), variables)
    end
    if length(tags) > 0
        mask = map(v -> length(intersect(v.tags, tags)) > 0, variables )
        return variables[mask]
    end
    return variables
end

function listVariables(dfg::LightDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)

    # variables = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGVariable, vertices(dfg.g)))
    if length(tags) > 0
        return map(v -> v.label, getVariables(dfg, regexFilter, tags=tags, solvable=solvable))
    else
        variables = collect(keys(dfg.g.variables))
        regexFilter !== nothing && (variables = filter(v -> occursin(regexFilter, String(v)), variables))
        solvable != 0 && (variables = filter(vId -> _isSolvable(dfg, vId, solvable), variables))
        return variables::Vector{Symbol}
    end
end

function getFactors(dfg::LightDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)
    # factors = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGFactor, vertices(dfg.g)))
    factors = collect(values(dfg.g.factors))
    if regexFilter != nothing
        factors = filter(f -> occursin(regexFilter, String(f.label)), factors)
    end
    if solvable != 0
        factors = filter(f -> _isSolvable(dfg, f.label, solvable), factors)
    end
    if length(tags) > 0
        mask = map(v -> length(intersect(v.tags, tags)) > 0, factors )
        return factors[mask]
    end
    return factors
end

function listFactors(dfg::LightDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)
    # factors = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGFactor, vertices(dfg.g)))
    if length(tags) > 0
        return map(v -> v.label, getFactors(dfg, regexFilter, tags=tags, solvable=solvable))
    end
    factors = collect(keys(dfg.g.factors))
    if regexFilter != nothing
        factors = filter(f -> occursin(regexFilter, String(f)), factors)
    end
    if solvable != 0
        factors = filter(fId -> _isSolvable(dfg, fId, solvable), factors)
    end
    return factors::Vector{Symbol}
end

function isConnected(dfg::LightDFG)::Bool
    return LightGraphs.is_connected(dfg.g)
    # return length(LightGraphs.connected_components(dfg.g)) == 1
end

function _isSolvable(dfg::LightDFG, label::Symbol, ready::Int)::Bool

    haskey(dfg.g.variables, label) && (return dfg.g.variables[label].solvable >= ready)
    haskey(dfg.g.factors, label) && (return dfg.g.factors[label].solvable >= ready)

    #TODO should this be a breaking error?
    @error "Node not in factor or variable"
    return false
end

function getNeighbors(dfg::LightDFG, node::DFGNode; solvable::Int=0)
    label = node.label
    if !exists(dfg, label)
        error("Variable/factor with label '$(node.label)' does not exist in the factor graph")
    end

    neighbors_il =  FactorGraphs.outneighbors(dfg.g, dfg.g.labels[label])
    neighbors_ll = [dfg.g.labels[i] for i in neighbors_il]
    # Additional filtering
    solvable != 0 && filter!(lbl -> _isSolvable(dfg, lbl, solvable), neighbors_ll)

    # Variable sorting (order is important)
    if typeof(node) <: AbstractDFGFactor
        order = intersect(node._variableOrderSymbols, neighbors_ll)#map(v->v.dfgNode.label, neighbors))
        return order
    end

    return neighbors_ll::Vector{Symbol}
end

function getNeighbors(dfg::LightDFG, label::Symbol; solvable::Int=0)::Vector{Symbol}
    if !exists(dfg, label)
        error("Variable/factor with label '$(label)' does not exist in the factor graph")
    end

    neighbors_il =  FactorGraphs.outneighbors(dfg.g, dfg.g.labels[label])
    neighbors_ll = [dfg.g.labels[i] for i in neighbors_il]
    # Additional filtering
    solvable != 0 && filter!(lbl -> _isSolvable(dfg, lbl, solvable), neighbors_ll)

    # Variable sorting (order is important)
    if haskey(dfg.g.factors, label)
        order = intersect(dfg.g.factors[label]._variableOrderSymbols, neighbors_ll)#map(v->v.dfgNode.label, neighbors))
        return order
    end

    return neighbors_ll

end

function getNeighborhood(dfg::LightDFG, variableFactorLabels::Vector{Symbol}, distance::Int; solvable::Int=0)::Vector{Symbol}
    # find neighbors at distance to add
    nbhood = Int[]

    for l in variableFactorLabels
        union!(nbhood, neighborhood(dfg.g, dfg.g.labels[l], distance))
    end

    allvarfacs = [dfg.g.labels[id] for id in nbhood]

    solvable != 0 && filter!(nlbl -> (getSolvable(dfg, nlbl) >= solvable), allvarfacs)

    return allvarfacs

end

# TODO copy LightDFG to LightDFG overwrite
# function copyGraph!(destDFG::LightDFG,
#                     sourceDFG::LightDFG,
#                     variableFactorLabels::Vector{Symbol};
#                     copyGraphMetadata::Bool=false,
#                     overwriteDest::Bool=false,
#                     deepcopyNodes::Bool=false,
#                     verbose::Bool = true)

#  Biadjacency Matrix https://en.wikipedia.org/wiki/Adjacency_matrix#Of_a_bipartite_graph
function getBiadjacencyMatrix(dfg::LightDFG; solvable::Int=0)::NamedTuple{(:B, :varLabels, :facLabels),Tuple{LightGraphs.SparseMatrixCSC,Vector{Symbol}, Vector{Symbol}}}
    varLabels = listVariables(dfg, solvable=solvable)
    factLabels = listFactors(dfg, solvable=solvable)
    varIndex = [dfg.g.labels[s] for s in varLabels]
    factIndex = [dfg.g.labels[s] for s in factLabels]

    adj = adjacency_matrix(dfg.g)

    adjvf = adj[factIndex, varIndex]
    return (B=adjvf, varLabels=varLabels, facLabels=factLabels)
end

"""
    $(SIGNATURES)
Gets an empty and unique LightDFG derived from an existing DFG.
"""
function _getDuplicatedEmptyDFG(dfg::LightDFG{P,V,F})::LightDFG where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
    newDfg = LightDFG{P,V,F}(;
        userId=dfg.userId, robotId=dfg.robotId, sessionId=dfg.sessionId,
        solverParams=deepcopy(dfg.solverParams))
    newDfg.description ="(Copy of) $(dfg.description)"
    return newDfg
end

#TODO JT test.
"""
    $(SIGNATURES)
A replacement for to_dot that saves only hardcoded factor graph plotting attributes.
"""
function savedot_attributes(io::IO, dfg::LightDFG)
    write(io, "graph G {\n")

    for vl in listVariables(dfg)
        write(io, "$vl [color=red, shape=ellipse];\n")
    end
    for fl in listFactors(dfg)
        write(io, "$fl [color=blue, shape=box, fontsize=8, fixedsize=false, height=0.1, width=0.1];\n")
    end

    for e in edges(dfg.g)
        write(io, "$(dfg.g.labels[src(e)]) -- $(dfg.g.labels[dst(e)])\n")
    end
    write(io, "}\n")
end

function toDotFile(dfg::LightDFG, fileName::String="/tmp/dfg.dot")::Nothing
    open(fileName, "w") do fid
        savedot_attributes(fid, dfg)
    end
    return nothing
end

function toDot(dfg::LightDFG)::String
    m = PipeBuffer()
    savedot_attributes(m, dfg)
    data = take!(m)
    close(m)
    return String(data)
end

"""
    $SIGNATURES

Speciallized function available to only LightDFG at this time.

Notes
- Has option for various types of filters (increases memory usage)

Example
```julia
using IncrementalInference

# canonical example graph as example
fg = generateCanonicalFG_Kaess()

@show path = findShortestPathDijkstra(fg, :x1, :x3)
@show isVariable.(fg, path)
@show isFactor.(fg, path)
```

DevNotes
- TODO expand to other AbstractDFG entities.
- TODO use of filter resource consumption can be improved.

Related

[`findFactorsBetweenNaive`](@ref), `LightGraphs.dijkstra_shortest_paths`
"""
function findShortestPathDijkstra(  dfg::LightDFG,
                                    from::Symbol,
                                    to::Symbol;
                                    regexVariables::Union{Nothing, Regex}=nothing,
                                    regexFactors::Union{Nothing, Regex}=nothing,
                                    tagsVariables::Vector{Symbol}=Symbol[],
                                    tagsFactors::Vector{Symbol}=Symbol[],
                                    typeVariables::Union{Nothing, <:AbstractVector}=nothing,
                                    typeFactors::Union{Nothing, <:AbstractVector}=nothing,
                                    solvable::Int=0,
                                    initialized::Union{Nothing,Bool}=nothing  )
    #
    # helper function to filter on vector of types
    function _filterTypeList(thelist::Vector{Symbol}, typeList, listfnc=x->ls(dfg, x) )
        thelist_ = Symbol[]
        for type_ in typeList
            union!(thelist_, listfnc(type_))
        end
        intersect( thelist, thelist_ )
    end

    #
    duplicate = regexVariables !== nothing ||
                regexFactors !== nothing   ||
                0 < length(tagsVariables)  ||
                0 < length(tagsFactors)    ||
                typeVariables !== nothing  ||
                typeFactors !== nothing    ||
                initialized !== nothing    ||
                solvable != 0
    #
    dfg_ = if duplicate
        # use copy if filter is being applied
        varList = ls(dfg, regexVariables, tags=tagsVariables, solvable=solvable)
        fctList = lsf(dfg, regexFactors, tags=tagsFactors, solvable=solvable)
        varList = typeVariables !== nothing ? _filterTypeList(varList, typeVariables) : varList
        fctList = typeFactors !== nothing ? _filterTypeList(fctList, typeFactors, x->lsf(dfg, x)) : fctList
        varList = if initialized !== nothing
            initmask = isInitialized.(dfg, varList) .== initialized
            varList[initmask]
        else
            varList
        end
        deepcopyGraph(typeof(dfg), dfg, varList, fctList)
    else
        # no filter can be used directly
        dfg
    end

    if !exists(dfg_, from) || !exists(dfg_, to)
        # assume filters excluded either `to` or `from` and hence no shortest path
        return Symbol[]
    end
    # LightDFG internally uses Integers
    frI = dfg_.g.labels[from]
    toI = dfg_.g.labels[to]

    # get shortest path from graph provider
    path_state = LightGraphs.dijkstra_shortest_paths( dfg_.g.graph, [frI;] )
    path = LightGraphs.enumerate_paths( path_state, toI )
    dijkpath = map( x->dfg_.g.labels[x], path )

    # return the list of symbols
    return dijkpath
end


text:
#using Revise
using LinearAlgebra
using CartesianFDM

n = (17, 16)
top = (nonperiodic(), periodic())

ctx = cartesianfdmcontext(top, n)

# square domain (0, 1) x (0, 1)
h = spacing.(top, n)
x = coordinate.(top, n)

# cylinder centered at (0.5, 0.5) with radius 0.25
cylinder(x, y, z...) = 0.25 - (2x - 1) ^ 2 - (2y - 1) ^ 2

# capacities
using Vofinit

h₀ = Cdouble.((h..., sum(h) / 2))
xₑ = Cdouble[0, 0, 0, 0]

using Base.Iterators

# volume
V = map(zip(CartesianIndices(n), product(x...))) do (ind, el)
    x₀ = Cdouble.((el..., 0.))
    prod(h) * getcc(cylinder, x₀, h₀, xₑ)
end

V = map(mask(ctx, reshape(V, :))) do val
    iszero(val) ? eps(val) : val
end

A₁ = map(zip(CartesianIndices(n), product(x...))) do (ind, el)
    x₀ = Cdouble.((el..., 0.))
    h[2] * getcc(x₀, h₀, xₑ) do y...
        cylinder(x₀[1], y[2])
    end
end

A₂ = map(zip(CartesianIndices(n), product(x...))) do (ind, el)
    x₀ = Cdouble.((el..., 0.))
    h[1] * getcc(x₀, h₀, xₑ) do y...
        cylinder(y[1], x₀[2])
    end
end

A = mask(ctx, [reshape(A₁, :), reshape(A₂, :)])

P = scalar(:P, n)
U = vector(:U, n)
N = reshape(map(CartesianIndices(n)) do ind
                ind[1] ≤ 4 || ind[1] ≥ n[1] - 4 ? 1. : 0.
            end, :)

# divergence and gradient
cont = divergence(neu, ctx, A, U, N)
grad = -gradient(neu, ctx, A, 0, P, 0)

@assert begin
    mapreduce(*, eachindex(U)) do i
        D = Symbolics.jacobian(cont, U[i])
        G = Symbolics.jacobian(grad[i], P)
        all(iszero, D' - G)
    end
end

W = [reshape(map(CartesianIndices(n)) do ind
                ind[1] ≤ 4 || ind[1] ≥ n[1] - 4 ? 1. : 0.
            end, :), zeros(prod(n))]

#conv =
E = strainrate(dir, ctx, A, V, U, W)
visc = divergence(dir, ctx, A, E)

# diagonal
@assert begin
    mapreduce(*, eachindex(U)) do i
        D = Symbolics.jacobian(visc[i], U[i])
        all(≤(100eps(1.)), D-D')
    end
end

# off-diagonal
@assert begin
    B₁ = Symbolics.jacobian(visc[1], U[2])
    B₂ = Symbolics.jacobian(visc[2], U[1])
    all(≤(100eps(1.)), B₁-B₂')
end

Re = 10.

mom = grad ./ 2 .+ visc ./ 8Re
rhs = vcat(mom..., cont ./ 2)
var = vcat(U..., P)

filename = "stokes.jl"
(fun, fun!) = build_function(rhs, var)
open(filename, "w") do file
    write(file, string(fun))
    write(file, "\n")
    write(file, string(fun!))
end
#include("model.jl")

(σ⁻, _) = getproperty(ctx, :σ)

Mᵥ = map(σ⁻) do σ
    σ * V / 2
end

reshape(Mᵥ[2], n...)[n[1], :] .= prod(h)

Mₚ = map(cont) do el
    iszero(el) ? prod(h) : zero(first(h))
end

M = Diagonal(vcat(Mᵥ..., Mₚ))

diffvar = (!iszero).(vcat(Mᵥ..., Mₚ))

using DifferentialEquations

u₀ = zeros((length(n)+1)*prod(n))
du₀ = zero(u₀)
tspan = (0.0, 1.0)

### option 2
#=
function dae(out, du, u, M, t)
    stokes!(out, u, nothing, t)
    out .= out - M * du
    nothing
end

prob = DAEProblem(dae, du₀, u₀, tspan, M, differential_vars=diffvar)

using Sundials
sol = solve(prob,IDA())
=#

### option 1
#=
f = ODEFunction(stokes!, mass_matrix=M)
prob = ODEProblem(f, x₀, tspan)

sol = solve(prob, Rodas5(), reltol=1e-8, abstol=1e-8)
=#

# test divergence
#=
V̂ = map(mask(ctx, fill(prod(h), prod(n)))) do x
    iszero(x) ? eps(x) : x
end
=#

#=
val = map(eachindex(h)) do i
    prod(h[j] for j in eachindex(h) if j ≠ i)
end

Â = mask(ctx, [fill(val[i], prod(n)) for i in eachindex(n)])

Ŵ = mask(ctx, ones(prod(n)))
Û = [ones(prod(n)), zeros(prod(n))]

@assert all(iszero, divergence(neu, ctx, Â, Û, Ŵ))
foo = divergence(neu, ctx, mask(ctx, vector(:A, n)), Û, 0)
bar = divergence(neu, ctx, mask(ctx, vector(:A, n)), 0, Ŵ)
=#
#=
W =
=#