text:
# -- PRIVATE METHODS BELOW ------------------------------------------------------------------------------ #
function build_function_header_buffer(comment_dictionary)

    # initialize -
    buffer = ""

    # get some data from the comment_dictionary -
    function_name = comment_dictionary["function_name"]
    function_description = comment_dictionary["function_description"]
    input_arg_array = comment_dictionary["input_args"]
    output_arg_array = comment_dictionary["output_args"]

    buffer *= "# ----------------------------------------------------------------------------------- #\n"
    buffer *= "# Function: $(function_name)\n"
    buffer *= "# Description: $(function_description)\n"
    buffer *= "# Generated on: $(now())\n"
    buffer *= "#\n"
    buffer *= "# Input arguments:\n"

    for argument_dictionary in input_arg_array

      arg_symbol = argument_dictionary["symbol"]
      arg_description = argument_dictionary["description"]

      # write the buffer -
      buffer *= "# $(arg_symbol) => $(arg_description) \n"
    end

    buffer *= "#\n"
    buffer *= "# Output arguments:\n"
    for argument_dictionary in output_arg_array

      arg_symbol = argument_dictionary["symbol"]
      arg_description = argument_dictionary["description"]

      # write the buffer -
      buffer *= "# $(arg_symbol) => $(arg_description) \n"
    end
    buffer *= "# ----------------------------------------------------------------------------------- #\n"

    # return the buffer -
    return buffer
end

function build_copyright_header_buffer(problem_dictionary::Dict{String,Any})

    # What is the current year?
    current_year = string(Dates.year(now()))

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

    # return -
    return buffer
end

function build_transcription_rate_buffer(problem_dictionary::Dict{String,Any})::Array{String,1}

  # initialize -
  buffer = Array{String,1}()

  # get the comment buffer -
  comment_header_dictionary = problem_dictionary["configuration_dictionary"]["function_comment_dictionary"]["calculate_transcription_rates"]
  function_comment_buffer = build_function_header_buffer(comment_header_dictionary)
  +(buffer, function_comment_buffer)

  # function -
  path_to_impl = "$(path_to_package)/distribution/julia/include/compute_transcription_rates.jl"
  buffer = include_function(path_to_impl)

  # return -
  return buffer
end

function build_translation_rate_buffer(problem_dictionary::Dict{String,Any})::Array{String,1}

  # initialize -
  buffer = Array{String,1}()

  # get the comment buffer -
  comment_header_dictionary = problem_dictionary["configuration_dictionary"]["function_comment_dictionary"]["calculate_translation_rates"]
  function_comment_buffer = build_function_header_buffer(comment_header_dictionary)
  +(buffer, function_comment_buffer)

  # function -
  path_to_impl = "$(path_to_package)/distribution/julia/include/compute_translation_rates.jl"
  buffer = include_function(path_to_impl)

  # return -
  return buffer
end

function build_translation_control_buffer(problem_dictionary::Dict{String,Any})::Array{String,1}

  # initialize -
  buffer = Array{String,1}()

  # get the comment buffer -
  comment_header_dictionary = problem_dictionary["configuration_dictionary"]["function_comment_dictionary"]["translation_control_function"]
  function_comment_buffer = build_function_header_buffer(comment_header_dictionary)
  +(buffer, function_comment_buffer)

  # function -
  +(buffer,"function calculate_translation_control(state::Array{Float64,1}, data_dictionary::Dict{String,Any},vargs...)\n")
  +(buffer,"\n")

  # how many proteins do we have?
  list_of_seq_objects::Array{VLSequenceRecord} = problem_dictionary["sequence_object_array"]
  idx_proteins = findall(x->(x.operationType == :L), list_of_seq_objects)

  # initialize -
  +(buffer,"\t# initialize -\n")
  +(buffer, "\tw = ones($(length(idx_proteins)))\n")
  +(buffer, "\n")

  # comment -
  +(buffer,"\t# User: override the default translation control impl here - \n")
  +(buffer,"\t# TODO: override ... \n")

  # # list of proteins -
  # list_of_protein_objects = list_of_species_objects[idx_proteins]
  # for (index,species_object) in enumerate(list_of_protein_objects)

  #   +(buffer,"\tw[$(index)] = 1.0")
  #   +(buffer,"\t# $(index) TL control variable: $(species_object.species_symbol)\n")

  # end

  +(buffer,"\n")
  +(buffer,"\t# return - \n")
  +(buffer,"\treturn w\n")
  +(buffer,"end\n")

  # return -
  return buffer
end

function does_set_contain_species(species_set::Array{VLSpeciesSymbol},species_object::VLSpeciesSymbol)

  # initilize -
  does_set_contains_species_flag = false
  for local_species_object in species_set

    # get type and name -
    local_species_type = local_species_object.species_type
    local_species_symbol = "$(string(local_species_type))_"*local_species_object.species_symbol

    # get the test symbol and type -
    test_species_type = species_object.species_type
    test_species_symbol = species_object.species_symbol

    if (local_species_type == test_species_type && test_species_symbol == local_species_symbol)
      does_set_contains_species_flag = true
      return does_set_contains_species_flag
    end
  end

  return does_set_contains_species_flag
end

function is_species_a_target_in_connection_list(list_of_connections::Array{VLGRNConnectionObject}, target_species::VLSpeciesSymbol, connection_type::Symbol)

  target_connection_list = VLGRNConnectionObject[]
  for connection_object in list_of_connections

    # get targets -
    local_target_set = connection_object.connection_target_set
    local_connection_type = connection_object.connection_type
    if (local_connection_type == connection_type && does_set_contain_species(local_target_set,target_species) == true)
      push!(target_connection_list,connection_object)
    end
  end

  return target_connection_list
end

function iterate_grn_control_control_connection(gene_object::VLSpeciesSymbol, list_of_connections::Array{VLGRNConnectionObject})

  # What is my gene_symbol -
  gene_symbol = gene_object.species_symbol

  tmp_buffer = ""
  for connection_object in list_of_connections

    # actor -
    connection_symbol = connection_object.connection_symbol
    tmp_buffer *= "\tW_$(gene_symbol)_$(connection_symbol) = parameter_value_dictionary[\"W_$(gene_symbol)_$(connection_symbol)\"]\n"
  end

  return tmp_buffer
end

function iterate_grn_binding_control_connection(gene_object::VLSpeciesSymbol,list_of_connections::Array{VLGRNConnectionObject})

  # What is my gene_symbol -
  gene_symbol = gene_object.species_symbol

  tmp_buffer = ""
  for connection_object in list_of_connections

    # actor -
    connection_symbol = connection_object.connection_symbol

    tmp_buffer *= "\tn_$(gene_symbol)_$(connection_symbol) = parameter_value_dictionary[\"n_$(gene_symbol)_$(connection_symbol)\"]\n"
    tmp_buffer *= "\tK_$(gene_symbol)_$(connection_symbol) = parameter_value_dictionary[\"K_$(gene_symbol)_$(connection_symbol)\"]\n"

  end

  return tmp_buffer
end

function build_transcription_control_buffer(problem_dictionary::Dict{String,Any})::Array{String,1}

  # initialize -
  buffer = Array{String,1}()

  # get the comment buffer -
  comment_header_dictionary = problem_dictionary["configuration_dictionary"]["function_comment_dictionary"]["transcription_control_function"]
  function_comment_buffer = build_function_header_buffer(comment_header_dictionary)
  +(buffer, function_comment_buffer)

  # get list of connections and species -
  list_of_grn_connections::Array{VLGRNConnectionObject} = problem_dictionary["list_of_grn_connections"]
  list_of_grn_species::Array{VLSpeciesSymbol} = problem_dictionary["list_of_species"]
  list_of_seq_objects::Array{VLSequenceRecord} = problem_dictionary["sequence_object_array"]

  # get the list of genes -
  idx_genes = findall(x->(x.species_type == :GENE), list_of_grn_species)
  list_of_genes = list_of_grn_species[idx_genes]

  # function -
  +(buffer,"function calculate_transcription_control(state::Array{Float64,1}, data_dictionary::Dict{String,Any}, vargs...)\n")
  +(buffer, "\n")
  +(buffer, "\t# initialize the control_array - \n")
  +(buffer, "\tcontrol_array = zeros($(length(list_of_genes)))\n")
  +(buffer, "\n")

  # alias the species list -
  +(buffer, "\t# Alias the species - \n")
  for (index,species_object) in enumerate(list_of_grn_species)

    # Grab the symbol -
    species_symbol = species_object.species_symbol

    # write the record -
    +(buffer, "\t$(species_symbol) = state[$(index)]\n")
  end
  +(buffer, "\n")

  # alias the binding parameters -
  +(buffer, "\t# Alias the model parameters - \n")
  +(buffer, "\tparameter_value_dictionary = data_dictionary[\"parameter_dictionary\"][\"parameter_value_dictionary\"]\n")

  # alias the parameter values -
  runtime_model_parameter_name_array = problem_dictionary["runtime_model_parameter_name_array"]
  for parameter_name in runtime_model_parameter_name_array
    +(buffer,"\t$(parameter_name)=parameter_value_dictionary[\"$(parameter_name)\"]\n")
  end
  +(buffer,"\n")

  # get list of genes -
  for (gene_index,gene_object) in enumerate(list_of_genes)

    # get gene symbol -
    gene_symbol = gene_object.species_symbol

    # for this gene, what is my enzyme symbol?
    enzyme_symbol = get_enzyme_symbol(gene_symbol, list_of_seq_objects)

    # connections -
    activating_connections = is_species_a_target_in_connection_list(list_of_grn_connections,gene_object,:activate)
    inhibiting_connections = is_species_a_target_in_connection_list(list_of_grn_connections,gene_object,:inhibit)

    # generate the binding functions -
    list_of_all_connections = VLGRNConnectionObject[]
    append!(list_of_all_connections, activating_connections)
    append!(list_of_all_connections, inhibiting_connections)
    for (index,connection_object) in enumerate(list_of_all_connections)

      # actor -
      actor_list = connection_object.connection_actor_set
      connection_symbol = connection_object.connection_symbol

      +(buffer, "\t# Transfer function target:$(gene_symbol) actor:$(connection_symbol)\n")
      +(buffer, "\tactor_set_$(gene_symbol)_$(connection_symbol) = [\n")
      for actor_object in actor_list

        actor_symbol = actor_object.species_symbol
        actor_type = actor_object.species_type
        +(buffer, "\t\tPROTEIN_$(actor_symbol)\n")
      end

      +(buffer, "\t]\n")
      +(buffer, "\tactor = prod(actor_set_$(gene_symbol)_$(connection_symbol))\n")
      +(buffer, "\tb_$(gene_symbol)_$(connection_symbol) = (actor^(n_$(gene_symbol)_$(connection_symbol)))/")
      +(buffer, "(K_$(gene_symbol)_$(connection_symbol)^(n_$(gene_symbol)_$(connection_symbol))+actor^(n_$(gene_symbol)_$(connection_symbol)))\n")
      +(buffer, "\n")
    end

    +(buffer, "\t# Control function for $(gene_symbol) - \n")
    +(buffer,"\tcontrol_array[$(gene_index)] = (")
    numerator = ""
    if (isempty(activating_connections) == true)
      +(buffer, "W_$(gene_symbol)_$(enzyme_symbol)")
    else

      numerator *= "W_$(gene_symbol)_$(enzyme_symbol)+"
      for connection_object in activating_connections
        # actor -
        connection_symbol = connection_object.connection_symbol
        numerator *= "W_$(gene_symbol)_$(connection_symbol)*b_$(gene_symbol)_$(connection_symbol)+"
      end
      +(buffer, numerator[1:end-1])
    end

    +(buffer, ")/(1+W_$(gene_symbol)_$(enzyme_symbol)")

    if (isempty(activating_connections) == false)
      demoninator = ""
      for connection_object in activating_connections
        # actor -
        connection_symbol = connection_object.connection_symbol
        demoninator *= "+W_$(gene_symbol)_$(connection_symbol)*b_$(gene_symbol)_$(connection_symbol)"
      end

      +(buffer, demoninator[1:end])
    end

    # ok - do we have inhibitory statements?
    if (isempty(inhibiting_connections) == true)
      +(buffer, ")\n")
      +(buffer, "\n")
    else

      demoninator = ""
      for connection_object in inhibiting_connections
        # actor -
        connection_symbol = connection_object.connection_symbol
        demoninator *= "+W_$(gene_symbol)_$(connection_symbol)*b_$(gene_symbol)_$(connection_symbol)"
      end

      +(buffer, demoninator[1:end])
      +(buffer, ")\n")
      +(buffer, "\n")
    end
  end

  +(buffer, "\t# return - \n")
  +(buffer, "\treturn control_array\n")
  +(buffer,"end\n")

  # return -
  return buffer
end

function build_allosteric_control_buffer(problem_dictionary::Dict{String,Any})::Array{String,1}

  # initialize -
  buffer = Array{String,1}()

  # get the comment buffer -
  comment_header_dictionary = problem_dictionary["configuration_dictionary"]["function_comment_dictionary"]["transcription_control_function"]
  function_comment_buffer = build_function_header_buffer(comment_header_dictionary)
  +(buffer, function_comment_buffer)

  # get a list of reaction names -
  list_of_reactions = problem_dictionary["metabolic_reaction_order_array"]
  number_of_reactions = length(list_of_reactions)

  # function -
  +(buffer,"function calculate_allosteric_control(state::Array{Float64,1}, data_dictionary::Dict{String,Any}, vargs...)\n")
  +(buffer,"\n")

  +(buffer,"\t# initialize the v-vector - \n")
  +(buffer,"\tv = ones($(number_of_reactions))\n");

  +(buffer, "\n")
  +(buffer,"\t# User: override the default allosteric control impl here - \n")
  +(buffer,"\t# TODO: override ... \n")

  +(buffer,"\n")
  +(buffer,"\t# return - \n")
  +(buffer,"\treturn v\n")
  +(buffer,"end\n")

  # return -
  return buffer
end

# -- PRIVATE METHODS ABOVE ------------------------------------------------------------------------------ #

# -- PUBLIC METHODS BELOW ------------------------------------------------------------------------------- #

# function to build Control.jl
function build_control_program_component(problem_dictionary::Dict{String,Any})::VLProgramComponent

  # initialize -
  filename = "Control.jl"
  program_component = VLProgramComponent()
  buffer = Array{String,1}()

  # build the header -
  header_buffer = build_copyright_header_buffer(problem_dictionary)
  +(buffer, header_buffer)
  +(buffer, "\n")

  # build the TX buffer -
  tx_control_buffer = build_transcription_control_buffer(problem_dictionary)
  +(buffer, tx_control_buffer)
  +(buffer, "\n")

  # build the TL buffer -
  tl_control_buffer = build_translation_control_buffer(problem_dictionary)
  +(buffer, tl_control_buffer)

  # build the TL buffer -
  allosteric_control_buffer = build_allosteric_control_buffer(problem_dictionary)
  +(buffer, allosteric_control_buffer)

  # collapse -
  flat_buffer = ""
  [flat_buffer *= line for line in buffer]

  # add data to program_component -
  program_component.filename = filename;
  program_component.buffer = flat_buffer
  program_component.type = :buffer

  # return -
  return program_component
end

# function to build Kinetics.jl
function build_kinetics_program_component(problem_dictionary::Dict{String,Any})::VLProgramComponent

  # initialize -
  filename = "Kinetics.jl"
  program_component = VLProgramComponent()
  buffer = Array{String,1}()

  # build the header -
  header_buffer = build_copyright_header_buffer(problem_dictionary)
  +(buffer, header_buffer)
  +(buffer, "\n")

  # build the TX buffer -
  tx_buffer = build_transcription_rate_buffer(problem_dictionary)
  +(buffer, tx_buffer)
  +(buffer, "\n")

  # build the TL buffer -
  tl_buffer = build_translation_rate_buffer(problem_dictionary)
  +(buffer, tl_buffer)

  # collapse -
  flat_buffer = ""
  [flat_buffer *= line for line in buffer]

  # add data to program_component -
  program_component.filename = filename;
  program_component.buffer = flat_buffer
  program_component.type = :buffer

  # return -
  return program_component
end

# function to build the Data.jl -
function build_data_dictionary_program_component(problem_dictionary::Dict{String,Any})::VLProgramComponent

  # initialize -
  filename = "Data.jl"
  program_component = VLProgramComponent()
  buffer = Array{String,1}()

  # grab the default value dictionary, the get default values that we put in the code -
  default_parameter_dictionary = problem_dictionary["default_parameter_dictionary"]
  lower_flux_bound_value = default_parameter_dictionary["default_flux_lower_bound"]
  upper_flux_bound_value = default_parameter_dictionary["default_flux_upper_bound"]
  lower_species_bound_value = default_parameter_dictionary["default_species_lower_bound"]
  upper_species_bound_value = default_parameter_dictionary["default_species_upper_bound"]
  default_transcription_time_constant_modifier = default_parameter_dictionary["default_transcription_time_constant_modifier"]
  default_translation_time_constant_modifier = default_parameter_dictionary["default_translation_time_constant_modifier"]

  # build the header -
  header_buffer = build_copyright_header_buffer(problem_dictionary)
  +(buffer, header_buffer)

  # get the comment buffer -
  comment_header_dictionary = problem_dictionary["configuration_dictionary"]["function_comment_dictionary"]["data_dictionary_function"]
  function_comment_buffer = build_function_header_buffer(comment_header_dictionary)
  +(buffer, function_comment_buffer)

  # before we get too far, get some stuff from the problem_dictionary -
  list_of_species = problem_dictionary["list_of_species"]
  list_of_sequence_objects = problem_dictionary["sequence_object_array"]

  # get the order of reactions -
  metabolic_reaction_name_array = problem_dictionary["metabolic_reaction_order_array"]
  metabolic_reaction_object_dictionary = problem_dictionary["metabolic_reaction_dictionary"]

  # function -
  +(buffer,"function generate_default_data_dictionary(")
  +(buffer,"path_to_biophysical_constants_file::String = \"./src/config/Biophysics.toml\",\n");
  +(buffer,"\tpath_to_ec_file::String = \"./src/config/EC.toml\", path_to_control_constants_file::String=\"./src/config/Control.toml\",\n");
  +(buffer,"\tbiophysical_constants_file_parser::Function=build_biophysical_dictionary, ec_constants_file_parser::Function=build_ec_data_dictionary,\n\tcontrol_constants_parser::Function=build_control_constants_dictionary, vargs...)\n")
  +(buffer,"\n")
  +(buffer,"\t# initialize - \n")
  +(buffer,"\tdata_dictionary = Dict{String,Any}()\n")
  +(buffer,"\n")
  +(buffer,"\t# load the stoichiometric array - \n")
  +(buffer,"\tstoichiometric_matrix = readdlm(\"./src/Network.dat\")")
  +(buffer,"\t# WARNING: loads default stm from the expected location. If you've changed the location (or stm) => change the path.\n")
  +(buffer,"\t(number_of_species, number_of_reactions) = size(stoichiometric_matrix)\n")
  +(buffer, "\n")
  +(buffer, "\t# Load the biophysical dictionary - \n")
  +(buffer, "\tbiophysical_dictionary = biophysical_constants_file_parser(path_to_biophysical_constants_file)\n")
  +(buffer, "\n")
  +(buffer, "\t# Load the ec data dictionary - \n")
  +(buffer, "\tec_data_dictionary = ec_constants_file_parser(path_to_ec_file)\n")
  +(buffer, "\n")
  +(buffer, "\t# Load the parameter dictionary - \n")
  +(buffer, "\tparameter_dictionary = control_constants_parser(path_to_control_constants_file)\n")
  +(buffer, "\n")

  # generate default value for is_minimum_flag -
  +(buffer, "\n")
  +(buffer, "\t# The default behavior for the solver is to minimize - \n")
  +(buffer, "\tis_minimum_flag = true\n")
  +(buffer, "\n")

  # generate default object coeff array -
  +(buffer, "\n")
  +(buffer, "\t# Initialize an empty objective coefficient array - \n")
  +(buffer, "\tobjective_coefficient_array = zeros(number_of_reactions)\n")
  +(buffer, "\n")

  # generate species_symbol_type_array -
  +(buffer, "\t# List of the species types - \n")
  +(buffer, "\tspecies_symbol_type_array = [\n")
  for (index,species_object) in enumerate(list_of_species)
    +(buffer, "\t\t:$(species_object.species_type)\t;\t")
    +(buffer, "# $(index) $(species_object.species_symbol)\n")
  end
  +(buffer,"\t]\n")

  # generate gene_coding_length_array -
  +(buffer, "\n")
  +(buffer, "\t# Gene coding length array holds the length (nt) for model genes - \n")
  +(buffer, "\tgene_coding_length_array = [\n")
  idx_gene = findall(x->(x.operationType == :X), list_of_sequence_objects)
  list_of_genes = list_of_sequence_objects[idx_gene]
  for (index, seq_object) in enumerate(list_of_genes)

    +(buffer, "\t\t$(seq_object.length)\t;\t")
    +(buffer, "# $(index) GENE $(seq_object.modelSpecies)\n")
  end
  +(buffer,"\t]\n")

  # generate protein_coding_length_array
  +(buffer, "\n")
  +(buffer, "\t# Protein coding length array holds the length (aa) for model proteins - \n")
  +(buffer, "\tprotein_coding_length_array = [\n")
  idx_proteins = findall(x->(x.operationType == :L), list_of_sequence_objects)
  list_of_proteins = list_of_sequence_objects[idx_proteins]
  for (index, seq_object) in enumerate(list_of_proteins)

    +(buffer, "\t\t$(seq_object.length)\t;\t")
    +(buffer, "# $(index) PROTEIN $(seq_object.modelSpecies)\n")
  end
  +(buffer,"\t]\n")

  # generate time_constant_modifier_array -
  +(buffer, "\n")
  +(buffer, "\t# Time constant modifier array holds a modifier for the TX and TL time constants - \n")
  +(buffer, "\ttime_constant_modifier_array = [\n")
  ordered_seq_array = [list_of_genes list_of_proteins]
  for (index, seq_object) in enumerate(ordered_seq_array)

    # what type of operation -
    if (seq_object.operationType == :X)
      +(buffer, "\t\t$(default_transcription_time_constant_modifier)\t;\t")
    else
      +(buffer, "\t\t$(default_translation_time_constant_modifier)\t;\t")
    end

    +(buffer, "# $(index) $(seq_object.operationType) $(seq_object.modelSpecies)\n")
  end
  +(buffer,"\t]\n")

  # generate a default initial state array -
  +(buffer, "\n")
  +(buffer, "\t# Default: Initial concentrations for all model species are assumed to be 0 \n")
  +(buffer, "\t# Default: You can override this setting in the update_species_concentration_array function \n")
  +(buffer, "\tspecies_concentration_array = zeros(number_of_species)\n")

  # generate dummy impl for species bounds array -
  +(buffer, "\n")
  +(buffer, "\t# Default: all species are assumed to be bounded by [0,0] \n")
  +(buffer, "\t# Default: You can override this setting in the update_species_bounds_array function \n")
  +(buffer, "\tLSB = $(lower_species_bound_value)\n")
  +(buffer, "\tUSB = $(lower_species_bound_value)\n")
  +(buffer, "\tspecies_bounds_array = [LSB*ones(number_of_species,1) USB*ones(number_of_species,1)]\n")

  # generate dummy impl for flux bounds array -
  +(buffer, "\n")
  +(buffer, "\t# Default: all fluxes are assumed to assumed to be bounded by [0,U] where U is some uper bound \n")
  +(buffer, "\t# Default: You can override this choice in the update_flux_bounds_array function \n")
  +(buffer, "\tLFB = $(lower_flux_bound_value)\n")
  +(buffer, "\tUFB = $(upper_flux_bound_value)\n")
  +(buffer, "\tflux_bounds_array = [LFB*ones(number_of_reactions,1) UFB*ones(number_of_reactions,1)]\n")

  # generate list of species symbols -
  +(buffer,"\n")
  +(buffer, "\t# List of the species symbol names - \n")
  +(buffer, "\tspecies_symbol_array = [\n")
  for (index,species_object) in enumerate(list_of_species)
    +(buffer, "\t\t\"$(species_object.species_symbol)\"\t;\t")
    +(buffer, "# $(index) $(species_object.species_type)\n")
  end
  +(buffer,"\t]\n")

  # generate species-order dictionary -
  +(buffer,"\n")
  +(buffer,"\t# Forward species order dictionary -\n")
  +(buffer,"\tspecies_order_dictionary = Dict{String,Int}(species_symbol_array[i]=>i for i=1:number_of_species)\n")

  # generate inverse species-order dictionary -
  +(buffer,"\n")
  +(buffer,"\t# Inverse species order dictionary -\n")
  +(buffer,"\tinverse_species_order_dictionary = Dict{Int,String}(i=>species_symbol_array[i] for i=1:number_of_species)\n")

  # generate list of reaction names -
  +(buffer,"\n")
  +(buffer, "\t# List of the reaction names - \n")
  +(buffer, "\treaction_name_array = [\n")
  for (index,reaction_tag) in enumerate(metabolic_reaction_name_array)

    # lookup the reaction object -
    reaction_object = metabolic_reaction_object_dictionary[reaction_tag]

    # get stuff -
    +(buffer, "\t\t\"$(reaction_tag)\"\t;\t")
    +(buffer, "# $(index) $(reaction_object.left_phrase)=$(reaction_object.right_phrase) $(reaction_object.reversible)\n")
  end
  +(buffer,"\t]\n")

  # generate species-order dictionary -
  +(buffer,"\n")
  +(buffer,"\t# Forward reaction order dictionary - \n")
  +(buffer,"\treaction_order_dictionary = Dict{String,Int}(reaction_name_array[i]=>i for i=1:number_of_reactions)\n")

  # generate inverse species-order dictionary -
  +(buffer,"\n")
  +(buffer,"\t# Inverse reaction order dictionary - \n")
  +(buffer,"\tinverse_reaction_order_dictionary = Dict{Int,String}(i=>reaction_name_array[i] for i=1:number_of_reactions)\n")

  +(buffer, "\n")
  +(buffer,"\t# -- DO NOT EDIT BELOW THIS LINE ------------------------------------------------------------------------------------ # \n")
  +(buffer,"\tdata_dictionary[\"stoichiometric_matrix\"] = stoichiometric_matrix\n")
  +(buffer,"\tdata_dictionary[\"number_of_species\"] = number_of_species\n")
  +(buffer,"\tdata_dictionary[\"number_of_reactions\"] = number_of_reactions\n")
  +(buffer,"\tdata_dictionary[\"is_minimum_flag\"] = is_minimum_flag\n")
  +(buffer,"\tdata_dictionary[\"objective_coefficient_array\"] = objective_coefficient_array\n")
  +(buffer,"\tdata_dictionary[\"species_concentration_array\"] = species_concentration_array\n")
  +(buffer,"\tdata_dictionary[\"species_bounds_array\"] = species_bounds_array\n")
  +(buffer,"\tdata_dictionary[\"flux_bounds_array\"] = flux_bounds_array\n")
  +(buffer,"\tdata_dictionary[\"species_symbol_type_array\"] = species_symbol_type_array\n")
  +(buffer,"\tdata_dictionary[\"gene_coding_length_array\"] = gene_coding_length_array\n")
  +(buffer,"\tdata_dictionary[\"protein_coding_length_array\"] = protein_coding_length_array\n")
  +(buffer,"\tdata_dictionary[\"time_constant_modifier_array\"] = time_constant_modifier_array\n")
  +(buffer,"\tdata_dictionary[\"biophysical_dictionary\"] = biophysical_dictionary\n")
  +(buffer,"\tdata_dictionary[\"ec_data_dictionary\"] = ec_data_dictionary\n")
  +(buffer,"\tdata_dictionary[\"parameter_dictionary\"] = parameter_dictionary\n")

  +(buffer, "\n")
  +(buffer, "\t# Order information - \n")
  +(buffer,"\tdata_dictionary[\"reaction_name_array\"] = reaction_name_array\n")
  +(buffer,"\tdata_dictionary[\"species_symbol_array\"] = species_symbol_array\n")
  +(buffer,"\tdata_dictionary[\"species_order_dictionary\"] = species_order_dictionary\n")
  +(buffer,"\tdata_dictionary[\"inverse_species_order_dictionary\"] = inverse_species_order_dictionary\n")
  +(buffer,"\tdata_dictionary[\"reaction_order_dictionary\"] = reaction_order_dictionary\n")
  +(buffer,"\tdata_dictionary[\"inverse_reaction_order_dictionary\"] = inverse_reaction_order_dictionary\n")
  +(buffer, "\n")
  +(buffer,"\treturn data_dictionary\n")
  +(buffer,"\t# -- DO NOT EDIT ABOVE THIS LINE ------------------------------------------------------------------------------------ # \n")
  +(buffer,"end\n")

  # collapse -
  flat_buffer = ""
  [flat_buffer *= line for line in buffer]

  # add data to program_component -
  program_component.filename = filename;
  program_component.buffer = flat_buffer
  program_component.type = :buffer

  # return -
  return program_component
end
# -- PUBLIC METHODS ABOVE ------------------------------------------------------------------------------- #"

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
# Residual Sum of Squares
function compute_objval_rss(sim_data::Vector{Float64}, exp_data::Vector{Float64})::Float64
    error::Float64 = 0.0

    @simd for i in eachindex(exp_data)
        @inbounds error += (sim_data[i] - exp_data[i])^2
    end

    return error
end

# Cosine similarity
function compute_objval_cos(sim_data::Vector{Float64}, exp_data::Vector{Float64})::Float64

    error::Float64 = 1.0 - dot(sim_data,exp_data)/(norm(sim_data)*norm(exp_data))

    return error
end

function conditions_index(condition_name::String)::Int

    return findfirst(isequal(condition_name),Sim.conditions)
end

function diff_sim_and_exp(
        sim_matrix::Matrix{Float64},
        exp_dict::Dict{String,Array{Float64,1}},
        exp_timepoint::Vector{Float64},
        conditions::Vector{String};
        sim_norm_max::Float64)
    sim_result::Vector{Float64} = []
    exp_result::Vector{Float64} = []

    for (idx,condition) in enumerate(conditions)
        if condition in keys(exp_dict)
            append!(sim_result,sim_matrix[Int.(exp_timepoint.+1),idx])
            append!(exp_result,exp_dict[condition])
        end
    end

    return (sim_result./sim_norm_max, exp_result)
end

# Define an objective function to be minimized.
function objective(indiv_gene::Vector{Float64})::Float64
    indiv::Vector{Float64} = decode_gene2val(indiv_gene)

    (p,u0) = update_param(indiv)

    if Sim.simulate!(p,u0) isa Nothing
        error::Vector{Float64} = zeros(length(observables))
        for (i,obs_name) in enumerate(observables)
            if isassigned(Exp.experiments,i)
                error[i] = compute_objval_rss(
                    diff_sim_and_exp(
                        Sim.simulations[i,:,:],
                        Exp.experiments[i],
                        Exp.get_timepoint(obs_name),
                        Sim.conditions,
                        sim_norm_max=ifelse(
                            Sim.normalization,
                            maximum(Sim.simulations[i,:,:]),
                            1.0
                        ),
                    )...
                )
            end
        end
        #=
        error::Vector{Float64} = zeros(16)

        sim_norm_max = maximum(Sim.simulations[observables_index("Phosphorylated_MEKc"),:,:])
        exp_t = Exp.get_timepoint(observables_index("Phosphorylated_MEKc"))
        error[1] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_MEKc"), Int.(exp_t.+1), conditions_index("EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_MEKc")]["EGF"]
        )
        error[2] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_MEKc"), Int.(exp_t.+1), conditions_index("HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_MEKc")]["HRG"]
        )

        sim_norm_max = maximum(Sim.simulations[observables_index("Phosphorylated_ERKc"),:,:])
        exp_t = Exp.get_timepoint(observables_index("Phosphorylated_ERKc"))
        error[3] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_ERKc"), Int.(exp_t.+1), conditions_index("EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_ERKc")]["EGF"]
        )
        error[4] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_ERKc"), Int.(exp_t.+1), conditions_index("HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_ERKc")]["HRG"]
        )

        sim_norm_max = maximum(Sim.simulations[observables_index("Phosphorylated_RSKw"),:,:])
        exp_t = Exp.get_timepoint(observables_index("Phosphorylated_RSKw"))
        error[5] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_RSKw"), Int.(exp_t.+1), conditions_index("EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_RSKw")]["EGF"]
        )
        error[6] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_RSKw"), Int.(exp_t.+1), conditions_index("HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_RSKw")]["HRG"]
        )

        sim_norm_max = maximum(Sim.simulations[observables_index("Phosphorylated_CREBw"),:,:])
        exp_t = Exp.get_timepoint(observables_index("Phosphorylated_CREBw"))
        error[7] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_CREBw"), Int.(exp_t.+1), conditions_index("EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_CREBw")]["EGF"]
        )
        error[8] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_CREBw"), Int.(exp_t.+1), conditions_index("HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_CREBw")]["HRG"]
        )

        sim_norm_max = maximum(Sim.simulations[observables_index("dusp_mRNA"),:,:])
        exp_t = Exp.get_timepoint(observables_index("dusp_mRNA"))
        error[9] = compute_objval_rss(
            Sim.simulations[observables_index("dusp_mRNA"), Int.(exp_t.+1), conditions_index("EGF")]./sim_norm_max,
            Exp.experiments[observables_index("dusp_mRNA")]["EGF"]
        )
        error[10] = compute_objval_rss(
            Sim.simulations[observables_index("dusp_mRNA"), Int.(exp_t.+1), conditions_index("HRG")]./sim_norm_max,
            Exp.experiments[observables_index("dusp_mRNA")]["HRG"]
        )

        sim_norm_max = maximum(Sim.simulations[observables_index("cfos_mRNA"),:,:])
        exp_t = Exp.get_timepoint(observables_index("cfos_mRNA"))
        error[11] = compute_objval_rss(
            Sim.simulations[observables_index("cfos_mRNA"), Int.(exp_t.+1), conditions_index("EGF")]./sim_norm_max,
            Exp.experiments[observables_index("cfos_mRNA")]["EGF"]
        )
        error[12] = compute_objval_rss(
            Sim.simulations[observables_index("cfos_mRNA"), Int.(exp_t.+1), conditions_index("HRG")]./sim_norm_max,
            Exp.experiments[observables_index("cfos_mRNA")]["HRG"]
        )

        sim_norm_max = maximum(Sim.simulations[observables_index("cFos_Protein"),:,:])
        exp_t = Exp.get_timepoint(observables_index("cFos_Protein"))
        error[13] = compute_objval_rss(
            Sim.simulations[observables_index("cFos_Protein"), Int.(exp_t.+1), conditions_index("EGF")]./sim_norm_max,
            Exp.experiments[observables_index("cFos_Protein")]["EGF"]
        )
        error[14] = compute_objval_rss(
            Sim.simulations[observables_index("cFos_Protein"), Int.(exp_t.+1), conditions_index("HRG")]./sim_norm_max,
            Exp.experiments[observables_index("cFos_Protein")]["HRG"]
        )

        sim_norm_max = maximum(Sim.simulations[observables_index("Phosphorylated_cFos"),:,:])
        exp_t = Exp.get_timepoint(observables_index("Phosphorylated_cFos"))
        error[15] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_cFos"), Int.(exp_t.+1), conditions_index("EGF")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_cFos")]["EGF"]
        )
        error[16] = compute_objval_rss(
            Sim.simulations[observables_index("Phosphorylated_cFos"), Int.(exp_t.+1), conditions_index("HRG")]./sim_norm_max,
            Exp.experiments[observables_index("Phosphorylated_cFos")]["HRG"]
        )
        =#
        return sum(error)
    else
        return Inf
    end
end

ext: jl

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim

max_forks_repo_name: himoto/paramEstim


text:
function appendsymbols(symbols::Vector{Symbol}, N::Int)
    nSymbols = length(symbols)
    @assert nSymbols <= N
    @assert nSymbols > 0
    return vcat(symbols, [ Symbol(:var, n) for n in 1:(N-nSymbols) ])
end


text:
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
using Bootstrap #External packages, must be added
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Esben's path
#cd("$(homedir())/Documents/GitHub/Thesis/Data")
#path = "$(homedir())/Documents/GitHub/Thesis/Data"

#Skipper's path
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data"
#HPC path
#path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data"
mainData = loadIndexDataNoDur(path)
fileName = path*"/Results/IndexData/IndexData"

"""
LASSO WITH WEIGHTED ERRORS
"""
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
using Bootstrap #External packages, must be added
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Esben's path
cd("$(homedir())/Documents/GitHub/Thesis/Data")
path = "$(homedir())/Documents/GitHub/Thesis/Data"
mainData = loadIndexDataNoDur(path)
fileName = path*"/Results/IndexData/IndexData"
mainDataArr = Array(mainData)

colNames = names(mainData)

nRows = size(mainDataArr)[1]
nCols = size(mainDataArr)[2]

mainXarr = mainDataArr[:,1:nCols-1]

# Transform #
mainXarr = expandWithTransformations(mainXarr)

# Standardize
standX = zScoreByColumn(mainXarr)
standY = zScoreByColumn(mainDataArr[:, nCols:nCols])
allData = hcat(standX, standY)
bCols = size(standX)[2]
nRows = size(standX)[1]
println(" \n \n \n solving Convex \n \n \n")

function solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma, errorWeights)
	bCols = size(Xtrain)[2]
	M = JuMP.Model(solver = GurobiSolver())
	@variables M begin
	        b[1:bCols]
	        t[1:bCols]
	        w
	end
	@objective(M,Min,0.5*w+gamma*ones(bCols)'*t)
	@constraint(M, soc, norm( [1-w;2*(errorWeights.*(Xtrain*b-Ytrain))] ) <= 1+w)
	@constraint(M,  b .<= t)
	@constraint(M, -t .<= b)

	solve(M)
	bSolved = getvalue(b)

	#In-Sample R-squared value
	errors = (Ytrain-Xtrain*bSolved)
	errorTotal = sum(errors[i]^2 for i=1:length(errors))
	errorsMean = Ytrain-mean(Ytrain)
	errorMeanTotal = sum(errorsMean[i]^2 for i=1:length(errorsMean))
	ISRsquared = 1-(errorTotal/errorMeanTotal)

	#OOS R-squared value
	oosErrors = Ypred - Xpred*bSolved
	oosErrorTotal = sum(oosErrors[i]^2 for i=1:length(oosErrors))
	oosErrorsMean = Ypred - mean(Ypred)
	oosErrorsMeanTotal = sum(oosErrorsMean[i]^2 for i=1:length(oosErrorsMean))
	OOSRsquared = 1-(oosErrorTotal/oosErrorsMeanTotal)

	return ISRsquared, OOSRsquared
end
trainingSize = 10
errorWeights = linspace(1,0, trainingSize)
testRuns = Int64(floor(1000/trainingSize))
nGammas = 500
predictions = 5
ISR = zeros(nGammas, testRuns)
OOSR = zeros(nGammas, testRuns)
gammaArray = logspace(0, 7, nGammas)
for r = 1:testRuns
	Xtrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), 1:bCols]
	Ytrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), bCols+1]
	Xpred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), 1:bCols]
	Ypred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), bCols+1]
	for g = 1:nGammas
		gamma = gammaArray[g]
		ISR[g, r], OOSR[g, r] = solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma, errorWeights)
	end
end

combinedArray = hcat(round.(gammaArray,3), ISR)
runCounter = collect(0:testRuns)
ISRcomb = vcat(runCounter', combinedArray)
writedlm("ISRsquaredweighted105.CSV", ISRcomb,",")
combinedArray = hcat(round.(gammaArray,3), OOSR)
OOSRcomb = vcat(runCounter', combinedArray)
writedlm("OOSRsquaredweighted105.CSV", OOSRcomb,",")

"""
LASSO
"""
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
using Bootstrap #External packages, must be added
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Esben's path
#cd("$(homedir())/Documents/GitHub/Thesis/Data")
#path = "$(homedir())/Documents/GitHub/Thesis/Data"

#Skipper's path
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data"

mainData = loadIndexDataNoDur(path)
<<<<<<< HEAD
=======
cd("$(homedir())/Documents/GitHub/Thesis/Data")
path = "$(homedir())/Documents/GitHub/Thesis/Data"
mainData = loadIndexDataNoDur(path)
fileName = path*"/Results/IndexData/IndexData"
>>>>>>> 4fe6210679f2fdd18244c00e7397df164bf0740d
=======
>>>>>>> parent of be6f44e... Merge branch 'master' of https://github.com/skipperafrosenborg/Thesis
mainDataArr = Array(mainData)

colNames = names(mainData)

nRows = size(mainDataArr)[1]
nCols = size(mainDataArr)[2]

mainXarr = mainDataArr[:,1:nCols-1]
mainYarr = mainDataArr[:, nCols:nCols]

# Transform with time elements
mainXarr = expandWithTime3612(mainXarr)

# Transform #
mainXarr = expandWithTransformations(mainXarr)

mainXarr = expandWithMAandMomentum(mainXarr, mainYarr, (nCols-1))

# Standardize
standX = zScoreByColumn(mainXarr)
standY = mainDataArr[:, nCols:nCols]
#standY = zScoreByColumn(mainDataArr[:, nCols:nCols])
allData = hcat(standX, standY)
bCols = size(standX)[2]
nRows = size(standX)[1]
println(" \n \n \n solving Convex \n \n \n")

function solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma)
	bCols = size(Xtrain)[2]
	M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
	@variables M begin
	        b[1:bCols]
	        t[1:bCols]
	        w
	end
	@objective(M,Min,0.5*w+gamma*ones(bCols)'*t)
	@constraint(M, soc, norm( [1-w;2*(Xtrain*b-Ytrain)] ) <= 1+w)
	@constraint(M,  b .<= t)
	@constraint(M, -t .<= b)

	solve(M)
	bSolved = getvalue(b)

	#In-Sample R-squared value
	errors = (Ytrain-Xtrain*bSolved)
	errorTotal = sum(errors[i]^2 for i=1:length(errors))
	errorsMean = Ytrain-mean(Ytrain)
	errorMeanTotal = sum(errorsMean[i]^2 for i=1:length(errorsMean))
	ISRsquared = 1-(errorTotal/errorMeanTotal)
	#=
	#OOS R-squared value (for more predictions)
	oosErrors = Ypred - Xpred*bSolved
	oosErrorTotal = sum(oosErrors[i]^2 for i=1:length(oosErrors))
	oosErrorsMean = Ypred - mean(Ypred)
	oosErrorsMeanTotal = sum(oosErrorsMean[i]^2 for i=1:length(oosErrorsMean))
	OOSRsquared = 1-(oosErrorTotal/oosErrorsMeanTotal)
	=#

	#OOS R-squared value (for single prediction)
	oosErrors = Ypred - Xpred*bSolved
	oosErrorTotal = sum(oosErrors[i]^2 for i=1:length(oosErrors))
	oosSize = sum(Ypred[i]^2 for i=1:length(Ypred))
	OOSRsquared = sqrt(oosErrorTotal) #/ sqrt(oosSize)

	#Indicator Results
	YpredValue = Ypred[1]
	Yestimate = Xpred*bSolved
	YestimateValue = Yestimate[1]
	if YpredValue >= 0 && YestimateValue >= 0
		Indicator = 1
	elseif YpredValue < 0 && YestimateValue < 0
		Indicator = 1
	else
		Indicator = 0
	end

	if YestimateValue >= 0
		Indicator2 = 1
	elseif YestimateValue < 0
		Indicator2 = 0
	else
		Indicator2 = 0
	end

	return ISRsquared, OOSRsquared, Indicator, Indicator2
end

<<<<<<< HEAD
<<<<<<< HEAD
<<<<<<< HEAD
<<<<<<< HEAD
=======
>>>>>>> parent of be6f44e... Merge branch 'master' of https://github.com/skipperafrosenborg/Thesis
nGammas = 10
trainingSize = 12
predictions = 1
testRuns = nRows-trainingSize-predictions
ISR = zeros(testRuns, nGammas)
OOSR = zeros(testRuns, nGammas)
Indi = zeros(testRuns, nGammas)
gammaArray = logspace(0, 3, nGammas)
<<<<<<< HEAD
=======
nGammas = 5
trainingSize = 48
=======
=======
>>>>>>> bee04ad7b02cee3642894bc88ca8f415908cae9a
nGammas = 2
trainingSize = 12
>>>>>>> bee04ad7b02cee3642894bc88ca8f415908cae9a
predictions = 1
testRuns = nRows-trainingSize-predictions
ISR = zeros(nGammas, testRuns)
OOSR = zeros(nGammas, testRuns)
Indi = zeros(nGammas, testRuns)
Indi2 = zeros(nGammas, testRuns)
gammaArray = logspace(0, 1, nGammas)
<<<<<<< HEAD
>>>>>>> 4fe6210679f2fdd18244c00e7397df164bf0740d
=======
>>>>>>> parent of be6f44e... Merge branch 'master' of https://github.com/skipperafrosenborg/Thesis
=======
>>>>>>> bee04ad7b02cee3642894bc88ca8f415908cae9a

for r = 1:(nRows-trainingSize-predictions)
	Xtrain = allData[r:(trainingSize+(r-1)), 1:bCols]
	Ytrain = allData[r:(trainingSize+(r-1)), bCols+1]
	Xpred  = allData[(trainingSize+(r-1)+1):(trainingSize+(r-1)+predictions), 1:bCols]
	Ypred  = allData[(trainingSize+(r-1)+1):(trainingSize+(r-1)+predictions), bCols+1]
	for g = 1:nGammas
		gamma = gammaArray[g]
<<<<<<< HEAD
<<<<<<< HEAD
		ISR[r, g], OOSR[r, g], Indi[r, g] = solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma)
=======
		ISR[g, r], OOSR[g, r], Indi[g, r], Indi2[g, r] = solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma)
>>>>>>> bee04ad7b02cee3642894bc88ca8f415908cae9a
=======
		ISR[g, r], OOSR[g, r], Indi[g, r], Indi2[g, r] = solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma)
>>>>>>> bee04ad7b02cee3642894bc88ca8f415908cae9a
	end
	println("Row $r/",(nRows-trainingSize-predictions))
end

#=
for r = 1:testRuns
	Xtrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), 1:bCols]
	Ytrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), bCols+1]
	Xpred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), 1:bCols]
	Ypred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), bCols+1]
	for g = 1:nGammas
		gamma = gammaArray[g]
		ISR[g, r], OOSR[g, r] = solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma)
	end
end
=#
fileName = path*"/Results/IndexData/LassoTests/121TimeExpTA"
combinedArray = vcat((round.(gammaArray,3))', ISR)

runCounter = collect(0:testRuns)
ISRcomb = hcat(runCounter, combinedArray)
writedlm(fileName*"ISRsquared.CSV", ISRcomb,",")
combinedArray = vcat((round.(gammaArray,3))', OOSR)
OOSRcomb = hcat(runCounter, combinedArray)
writedlm(fileName*"OSSRsquared.CSV", OOSRcomb,",")

combinedArray = vcat((round.(gammaArray,3))', Indi)
Indicomb = hcat(runCounter, combinedArray)
writedlm(fileName*"Indi.CSV", Indicomb,",")
<<<<<<< HEAD
=======
combinedArray = hcat(round.(gammaArray,3), ISR)
runCounter = collect(0:testRuns)'
ISRcomb = vcat(runCounter, combinedArray)
writedlm("ISRTestNoDur121IndiPortfolio.CSV", ISRcomb,",")
combinedArray = hcat(round.(gammaArray,3), OOSR)
OOSRcomb = vcat(runCounter, combinedArray)
writedlm("OOSRTestNoDur121Portfolio.CSV", OOSRcomb,",")

combinedArray = hcat(round.(gammaArray,3), Indi)
Indicomb = vcat(runCounter, combinedArray)
<<<<<<< HEAD
<<<<<<< HEAD
writedlm("Indicator481TimeTAOtherNotStandY.CSV", Indicomb,",")
>>>>>>> 4fe6210679f2fdd18244c00e7397df164bf0740d
=======
writedlm("IndicatorTestNoDur121IndiPortfolio.CSV", Indicomb,",")

>>>>>>> bee04ad7b02cee3642894bc88ca8f415908cae9a
=======
>>>>>>> parent of be6f44e... Merge branch 'master' of https://github.com/skipperafrosenborg/Thesis

=======
writedlm("IndicatorTestNoDur121IndiPortfolio.CSV", Indicomb,",")

>>>>>>> bee04ad7b02cee3642894bc88ca8f415908cae9a
combinedArray = hcat(round.(gammaArray,3), Indi2)
Indicomb = vcat(runCounter, combinedArray)
writedlm("Indicator2TestNoDur121IndiPortfolio.CSV", Indicomb,",")

"""
LASSO with daily returns
"""
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
using Bootstrap #External packages, must be added
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Esben's path
cd("$(homedir())/Documents/GitHub/Thesis/Data")
path = "$(homedir())/Documents/GitHub/Thesis/Data/DailyReturns"

"""
REMEMBER: Select index to predict. Then change dataloading. Then decide
observations and predictions to include - change fileName accordingly.
"""

mainData = loadIndexDataDailyOther(path)
fileName = path*"/Results/IndexData/DailyReturnsXobservationsYpredictions"
mainDataArr = Array(mainData)

colNames = names(mainData)

nRows = size(mainDataArr)[1]
nCols = size(mainDataArr)[2]

mainXarr = mainDataArr[:,1:nCols-1]

# Transform #
mainXarr = expandWithTransformations(mainXarr)

# Standardize
standX = zScoreByColumn(mainXarr)
standY = zScoreByColumn(mainDataArr[:, nCols:nCols])
allData = hcat(standX, standY)
bCols = size(standX)[2]
nRows = size(standX)[1]
println(" \n \n \n solving Convex \n \n \n")

function solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma)
	bCols = size(Xtrain)[2]
	M = JuMP.Model(solver = GurobiSolver())
	@variables M begin
	        b[1:bCols]
	        t[1:bCols]
	        w
	end
	@objective(M,Min,0.5*w+gamma*ones(bCols)'*t)
	@constraint(M, soc, norm( [1-w;2*(Xtrain*b-Ytrain)] ) <= 1+w)
	@constraint(M,  b .<= t)
	@constraint(M, -t .<= b)

	solve(M)
	bSolved = getvalue(b)

	#In-Sample R-squared value
	errors = (Ytrain-Xtrain*bSolved)
	errorTotal = sum(errors[i]^2 for i=1:length(errors))
	isValue = sum(Ytrain[i]^2 for i=1:length(Ytrain))
	isErrorPercentage = (errorTotal / isValue) * 100
	errorsMean = Ytrain-mean(Ytrain)
	errorMeanTotal = sum(errorsMean[i]^2 for i=1:length(errorsMean))
	ISRsquared = 1-(errorTotal/errorMeanTotal)

	#OOS R-squared value
	oosErrors = Ypred - Xpred*bSolved
	oosErrorTotal = sum(oosErrors[i]^2 for i=1:length(oosErrors))
	oosValue = sum(Ypred[i]^2 for i=1:length(Ypred))
	oosErrorPercentage = (oosErrorTotal / oosValue)*100
	oosErrorsMean = Ypred - mean(Ypred)
	oosErrorsMeanTotal = sum(oosErrorsMean[i]^2 for i=1:length(oosErrorsMean))
	OOSRsquared = 1-(oosErrorTotal/oosErrorsMeanTotal)

	return isErrorPercentage, ISRsquared, oosErrorPercentage, OOSRsquared
end

testRuns = 300
nGammas = 50
trainingSize = 50
predictions = 5
ISR = zeros(testRuns, nGammas)
ISP = zeros(testRuns, nGammas)
OOSR = zeros(testRuns, nGammas)
OOSP = zeros(testRuns, nGammas)
varianceTrainArray = zeros(testRuns)
gammaArray = logspace(0, 2, nGammas)
for r = 1:testRuns
	Xtrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), 1:bCols]
	Ytrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), bCols+1]
	varianceTrainArray[r] = var(Ytrain)
	Xpred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), 1:bCols]
	Ypred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), bCols+1]
	for g = 1:nGammas
		gamma = gammaArray[g]
		ISP[r, g], ISR[r, g], OOSP[r, g], OOSR[r, g] = solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma)
	end
end

runCounter = collect(0:testRuns)
varArr = vcat("variance",varianceTrainArray)
ISR
gammaArray'
combinedArray = vcat(round.(gammaArray',3), ISR)
ISRcomb = hcat(runCounter, varArr, combinedArray)
writedlm("ISRsquaredDailyOther505.CSV", ISRcomb,",")

combinedArray = vcat(round.(gammaArray',3), ISP)
ISPcomb = hcat(runCounter, varArr, combinedArray)
writedlm("ISPDailyOther505.CSV", ISPcomb,",")

combinedArray = vcat(round.(gammaArray',3), OOSR)
OOSRcomb = hcat(runCounter, varArr, combinedArray)
writedlm("OOSRsquaredDailyOther505.CSV", OOSRcomb,",")

combinedArray = vcat(round.(gammaArray',3), OOSP)
OOSPcomb = hcat(runCounter, varArr, combinedArray)
writedlm("OOSPDailyOther505.CSV", OOSPcomb,",")

"""
LASSO with daily returns and weighted
"""
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
using Bootstrap #External packages, must be added
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Esben's path
cd("$(homedir())/Documents/GitHub/Thesis/Data")
path = "$(homedir())/Documents/GitHub/Thesis/Data/DailyReturns"

"""
REMEMBER: Select index to predict. Then change dataloading. Then decide
observations and predictions to include - change fileName accordingly.
"""

mainData = loadIndexDataDailyOther(path)
fileName = path*"/Results/IndexData/DailyReturnsXobservationsYpredictions"
mainDataArr = Array(mainData)

colNames = names(mainData)

nRows = size(mainDataArr)[1]
nCols = size(mainDataArr)[2]

mainXarr = mainDataArr[:,1:nCols-1]

# Transform #
mainXarr = expandWithTransformations(mainXarr)

# Standardize
standX = zScoreByColumn(mainXarr)
standY = zScoreByColumn(mainDataArr[:, nCols:nCols])
allData = hcat(standX, standY)
bCols = size(standX)[2]
nRows = size(standX)[1]
println(" \n \n \n solving Convex \n \n \n")

function solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma, errorWeights)
	bCols = size(Xtrain)[2]
	M = JuMP.Model(solver = GurobiSolver())
	@variables M begin
	        b[1:bCols]
	        t[1:bCols]
	        w
	end
	@objective(M,Min,0.5*w+gamma*ones(bCols)'*t)
	@constraint(M, soc, norm( [1-w;2*(errorWeights.*(Xtrain*b-Ytrain))] ) <= 1+w)
	@constraint(M,  b .<= t)
	@constraint(M, -t .<= b)

	solve(M)
	bSolved = getvalue(b)

	#In-Sample R-squared value
	errors = (Ytrain-Xtrain*bSolved)
	errorTotal = sum(errors[i]^2 for i=1:length(errors))
	isValue = sum(Ytrain[i]^2 for i=1:length(Ytrain))
	isErrorPercentage = (errorTotal / isValue) * 100
	errorsMean = Ytrain-mean(Ytrain)
	errorMeanTotal = sum(errorsMean[i]^2 for i=1:length(errorsMean))
	ISRsquared = 1-(errorTotal/errorMeanTotal)

	#OOS R-squared value
	oosErrors = Ypred - Xpred*bSolved
	oosErrorTotal = sum(oosErrors[i]^2 for i=1:length(oosErrors))
	oosValue = sum(Ypred[i]^2 for i=1:length(Ypred))
	oosErrorPercentage = (oosErrorTotal / oosValue)*100
	oosErrorsMean = Ypred - mean(Ypred)
	oosErrorsMeanTotal = sum(oosErrorsMean[i]^2 for i=1:length(oosErrorsMean))
	OOSRsquared = 1-(oosErrorTotal/oosErrorsMeanTotal)

	return isErrorPercentage, ISRsquared, oosErrorPercentage, OOSRsquared
end

testRuns = 300
nGammas = 50
trainingSize = 50
errorWeights = linspace(0, 1, trainingSize)
predictions = 5
ISR = zeros(testRuns, nGammas)
ISP = zeros(testRuns, nGammas)
OOSR = zeros(testRuns, nGammas)
OOSP = zeros(testRuns, nGammas)
varianceTrainArray = zeros(testRuns)
gammaArray = logspace(0, 2, nGammas)
for r = 1:testRuns
	Xtrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), 1:bCols]
	Ytrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), bCols+1]
	varianceTrainArray[r] = var(Ytrain)
	Xpred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), 1:bCols]
	Ypred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), bCols+1]
	for g = 1:nGammas
		gamma = gammaArray[g]
		ISP[r, g], ISR[r, g], OOSP[r, g], OOSR[r, g] = solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma, errorWeights)
	end
end

runCounter = collect(0:testRuns)
varArr = vcat("variance",varianceTrainArray)
ISR
gammaArray'
combinedArray = vcat(round.(gammaArray',3), ISR)
ISRcomb = hcat(runCounter, varArr, combinedArray)
writedlm("ISRsquaredDailyOtherWeighted2505.CSV", ISRcomb,",")

combinedArray = vcat(round.(gammaArray',3), ISP)
ISPcomb = hcat(runCounter, varArr, combinedArray)
writedlm("ISPDailyOtherWeighted2505.CSV", ISPcomb,",")

combinedArray = vcat(round.(gammaArray',3), OOSR)
OOSRcomb = hcat(runCounter, varArr, combinedArray)
writedlm("OOSRsquaredDailyOtherWeighted2505.CSV", OOSRcomb,",")

combinedArray = vcat(round.(gammaArray',3), OOSP)
OOSPcomb = hcat(runCounter, varArr, combinedArray)
writedlm("OOSPDailyOtherWeighted2505.CSV", OOSPcomb,",")

"""
LASSO with monthly returns
"""
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
using Bootstrap #External packages, must be added
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Esben's path
cd("$(homedir())/Documents/GitHub/Thesis/Data")
path = "$(homedir())/Documents/GitHub/Thesis/Data"

mainData =loadIndexDataNoDur(path)
mainDataArr = Array(mainData)

colNames = names(mainData)

nRows = size(mainDataArr)[1]
nCols = size(mainDataArr)[2]

mainXarr = mainDataArr[:,1:nCols-1]

# Transform #
mainXarr = expandWithTransformations(mainXarr)

# Standardize
standX = zScoreByColumn(mainXarr)
standY = zScoreByColumn(mainDataArr[:, nCols:nCols])
allData = hcat(standX, standY)
bCols = size(standX)[2]
nRows = size(standX)[1]
println(" \n \n \n solving Convex \n \n \n")

function solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma)
	bCols = size(Xtrain)[2]
	M = JuMP.Model(solver = GurobiSolver())
	@variables M begin
	        b[1:bCols]
	        t[1:bCols]
	        w
	end
	@objective(M,Min,0.5*w+gamma*ones(bCols)'*t)
	@constraint(M, soc, norm( [1-w;2*(Xtrain*b-Ytrain)] ) <= 1+w)
	@constraint(M,  b .<= t)
	@constraint(M, -t .<= b)

	solve(M)
	bSolved = getvalue(b)

	#In-Sample R-squared value
	errors = (Ytrain-Xtrain*bSolved)
	errorTotal = sum(errors[i]^2 for i=1:length(errors))
	isValue = sum(Ytrain[i]^2 for i=1:length(Ytrain))
	isErrorPercentage = (errorTotal / isValue) * 100
	errorsMean = Ytrain-mean(Ytrain)
	errorMeanTotal = sum(errorsMean[i]^2 for i=1:length(errorsMean))
	ISRsquared = 1-(errorTotal/errorMeanTotal)

	#OOS R-squared value
	oosErrors = Ypred - Xpred*bSolved
	oosErrorTotal = sum(oosErrors[i]^2 for i=1:length(oosErrors))
	oosValue = sum(Ypred[i]^2 for i=1:length(Ypred))
	oosErrorPercentage = (oosErrorTotal / oosValue)*100
	oosErrorsMean = Ypred - mean(Ypred)
	oosErrorsMeanTotal = sum(oosErrorsMean[i]^2 for i=1:length(oosErrorsMean))
	OOSRsquared = 1-(oosErrorTotal/oosErrorsMeanTotal)

	return isErrorPercentage, ISRsquared, oosErrorPercentage, OOSRsquared
end

testRuns = 100
nGammas = 50
trainingSize = 10

predictions = 5
ISR = zeros(testRuns, nGammas)
ISP = zeros(testRuns, nGammas)
OOSR = zeros(testRuns, nGammas)
OOSP = zeros(testRuns, nGammas)
varianceTrainArray = zeros(testRuns)
gammaArray = logspace(0, 2, nGammas)
for r = 1:testRuns
	Xtrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), 1:bCols]
	Ytrain = allData[(1+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize), bCols+1]
	varianceTrainArray[r] = var(Ytrain)
	Xpred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), 1:bCols]
	Ypred  = allData[(trainingSize+(r-1)*trainingSize):(trainingSize+(r-1)*trainingSize+(predictions-1)), bCols+1]
	for g = 1:nGammas
		gamma = gammaArray[g]
		ISP[r, g], ISR[r, g], OOSP[r, g], OOSR[r, g] = solveLasso(Xtrain, Ytrain, Xpred, Ypred, gamma)
	end
end

runCounter = collect(0:testRuns)
varArr = vcat("variance",varianceTrainArray)
ISR
gammaArray'
combinedArray = vcat(round.(gammaArray',3), ISR)
ISRcomb = hcat(runCounter, varArr, combinedArray)
writedlm("ISRsquaredMonthlyOther105.CSV", ISRcomb,",")

combinedArray = vcat(round.(gammaArray',3), ISP)
ISPcomb = hcat(runCounter, varArr, combinedArray)
writedlm("ISPMonthlyOther105.CSV", ISPcomb,",")

combinedArray = vcat(round.(gammaArray',3), OOSR)
OOSRcomb = hcat(runCounter, varArr, combinedArray)
writedlm("OOSRsquaredMonthlyOther105.CSV", OOSRcomb,",")

combinedArray = vcat(round.(gammaArray',3), OOSP)
OOSPcomb = hcat(runCounter, varArr, combinedArray)
writedlm("OOSPMonthlyOther105.CSV", OOSPcomb,",")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
# # Custom Title 3-1

# This is the content

# # [Custom Title 3-2](@id custom_id_3-2)

# This is the content


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Here every gf is a gfvec, i.e every green-function is a vector of green function.
To lighten the writing we supress the "vec". Maybe consider putting everything as in the README.txt.
"""

module Green

using Base.Test

function readgreen_c(fin_gf_to::String; zn_col::Int = 1)

    gf_to = readdlm(fin_gf_to)

    # Get the Matsubara frequencies on imaginary axis
    zn_vec = deepcopy(gf_to[:, zn_col])

    # construct complex gf_t (each gf is formed by the consecutif real part
    # (pair columns) and imaginary part (odd columns))
    gf_c = to_to_c(zn_vec, gf_to)
    return (zn_vec, gf_c)

end

function to_to_t(gf_to::Array{Float64, 2})
    assert(size(gf_to)[2] == 11)
    gf_to = deepcopy(gf_to)
    return (1.0*gf_to[:, 2:2:end] + 1.0im*gf_to[:, 3:2:end])
end

function t_to_to(zn_vec::Array{Float64, 1}, gf_t::Array{Complex{Float64}, 2})

    zn_vec = deepcopy(zn_vec) ; gf_t = deepcopy(gf_t)
    assert(size(zn_vec)[1] == size(gf_t)[1])
    assert(size(gf_t)[2] == 5)

    gf_to = zeros(Float64, (size(zn_vec)[1], 11))

    gf_to[:, 1] = zn_vec
    gf_to[:, 2:2:end] =  real(gf_t)
    gf_to[:, 3:2:end] =  imag(gf_t)

    return gf_to
end

function t_to_c(zn_vec::Array{Float64, 1}, gf_t::Array{Complex{Float64}, 2})
    """convert from tabular form to cluster form """
    zn_vec = deepcopy(zn_vec) ; gf_t = deepcopy(gf_t)
    #assert(size(zn_vec)[1] == size(gf_t)[1])
    #assert(size(gf_t) == 5)

    gf_c = zeros(Complex{Float64}, (size(zn_vec)[1], 4, 4))

    gf_c[:, 1, 1] = gf_t[:, 1] ; gf_c[:, 1, 2] = gf_t[:, 3] ; gf_c[:, 1, 3] = gf_t[:, 4] ; gf_c[:, 1, 4] = gf_t[:, 3]
    gf_c[:, 2, 1] = gf_t[:, 3] ; gf_c[:, 2, 2] = gf_t[:, 2] ; gf_c[:, 2, 3] = gf_t[:, 3] ; gf_c[:, 2, 4] = gf_t[:, 5]
    gf_c[:, 3, 1] = gf_t[:, 4] ; gf_c[:, 3, 2] = gf_t[:, 3] ; gf_c[:, 3, 3] = gf_t[:, 1] ; gf_c[:, 3, 4] = gf_t[:, 3]
    gf_c[:, 4, 1] = gf_t[:, 3] ; gf_c[:, 4, 2] = gf_t[:, 5] ; gf_c[:, 4, 3] = gf_t[:, 3] ; gf_c[:, 4, 4] = gf_t[:, 2]

    return gf_c
end

function c_to_t(zn_vec::Array{Float64, 1}, gf_c::Array{Complex{Float64}, 3})
    """ """
    zn_vec = deepcopy(zn_vec) ; gf_c = deepcopy(gf_c)
    gf_to = c_to_to(zn_vec, gf_c)
    gf_t = to_to_t(gf_to)
    return gf_t
end

function to_to_c(zn_vec::Array{Float64, 1}, gf_to::Array{Float64, 2})
    """convert from tabular-out form to cluster form"""
    zn_vec = deepcopy(zn_vec) ; gf_to = deepcopy(gf_to)

    # println("size(zn_vec) =  ", size(zn_vec))
    # println("size(gf_to) = ", size(gf_to))

    assert(size(zn_vec)[1] == size(gf_to)[1])
    assert(size(gf_to)[2] == 11)

    # the following three lines could be replaced by: gf_t = to_to_t(gf_to)
    # but I do this for tests purposes (np.testing)
    gf_t = zeros(Complex{Float64}, (size(gf_to)[1], 5))
    gf_t = 1.0*gf_to[:, 2:2:end] + 1.0im*gf_to[:, 3:2:end]
    @test isapprox(gf_t, to_to_t(gf_to))

    gf_c = zeros(Complex{Float64}, (size(zn_vec)[1], 4, 4))

    gf_c[:, 1, 1] = gf_t[:, 1] ; gf_c[:, 1, 2] = gf_t[:, 3] ; gf_c[:, 1, 3] = gf_t[:, 4] ; gf_c[:, 1, 4] = gf_t[:, 3]
    gf_c[:, 2, 1] = gf_t[:, 3] ; gf_c[:, 2, 2] = gf_t[:, 2] ; gf_c[:, 2, 3] = gf_t[:, 3] ; gf_c[:, 2, 4] = gf_t[:, 5]
    gf_c[:, 3, 1] = gf_t[:, 4] ; gf_c[:, 3, 2] = gf_t[:, 3] ; gf_c[:, 3, 3] = gf_t[:, 1] ; gf_c[:, 3, 4] = gf_t[:, 3]
    gf_c[:, 4, 1] = gf_t[:, 3] ; gf_c[:, 4, 2] = gf_t[:, 5] ; gf_c[:, 4, 3] = gf_t[:, 3] ; gf_c[:, 4, 4] = gf_t[:, 2]

    @test isapprox(gf_c, t_to_c(zn_vec, gf_t))
    #println("gf_c = ")
    #println(gf_c)
    #println("gf_c = ")
    #println(t_to_c(zn_vec, gf_t))
    return gf_c
end

function c_to_to(zn_vec::Array{Float64}, gf_c::Array{Complex{Float64}, 3})
    """convert from cluster form to tabular-out form """
    zn_vec = deepcopy(zn_vec) ; gf_c = deepcopy(gf_c)
    assert(size(zn_vec)[1] == size(gf_c)[1])
    assert(size(gf_c)[2] == size(gf_c)[2]) ; assert(size(gf_c)[2] ==  size(gf_c)[3] == 4)

    gf_to = zeros(Float64, (size(zn_vec)[1], 11))

    gf_to[:, 1] = zn_vec
    gf_to[:, 2] = real(gf_c[:, 1, 1]) ; gf_to[:, 3] = imag(gf_c[:, 1, 1])
    gf_to[:, 4] = real(gf_c[:, 2, 2]) ; gf_to[:, 5] = imag(gf_c[:, 2, 2])
    gf_to[:, 6] = real(gf_c[:, 1, 2]) ; gf_to[:, 7] = imag(gf_c[:, 1, 2])
    gf_to[:, 8] = real(gf_c[:, 1, 3]) ; gf_to[:, 9] = imag(gf_c[:, 1, 3])
    gf_to[:, 10] = real(gf_c[:, 2, 4]) ; gf_to[:, 11] = imag(gf_c[:, 2, 4])

    return gf_to
end

end


text:
##  Copyright 2021 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

module RankedChoices

include("structures.jl")
include("sampler.jl")
include("hamiltonian-sampler.jl")
include("clustering.jl")
include("simulator.jl")
include("analysis.jl")
include("other-methods.jl")
include("parser.jl")
include("scrambler.jl")
include("embedding.jl")

export VoterRealization, VoterCohort, VoterMixture, Membership, Utility, RankedChoice
export IssueVote, MultiIssueVote
export validate_choices, sample_mixture_shares, sample_mixture_posterior
export obeys_ranking, shoehorn_ranking
export rejection_sample_utilities!, hamiltonian_sample_utilities!
export Prior, make_impartial_prior, simulate, HamiltonianSim, RejectionSim
export normalize_utility, estimate_quantiles, mutual_information, contest_probability, bootstrap_statistic
export count_simple_plurality, find_instant_runoff_winner, find_condorcet_winner
export read_prm_file, read_csv, read_xlsx, parse_matrix, sanitize
export scramble_latent_utilities, sample_possible_rankings
export optimize_von_neumann_entropy, embed_diffusion

end # module


text:
"""
Gage the geometric center (mean value) of an input array of arrays.

Arguments
---------
:Type{GeometricCenter}
    a Julia Type object and must be ``GeometricCenter``

array::AbstractArray
    an array of floating point arrays
}
"""
function gage_array_of_arrays(::Type{GeometricCenter}, array::AbstractArray)
    if length(array) == 0
        return Array{AbstractFloat,1}[]
    else
        return sum(array)/length(array)
    end
end


text:
function is_valid_iban(iban::String)
    ## convert to common format
    iban = unify_iban(iban)
    m = match(iban_global_re, iban)
    if m == nothing
        return false
    end
    country, check_digits, bban = m.captures
    ## supported country
    if !haskey(iban_country_length, country)
        return false
    end
    ## check total length
    if length(iban) != iban_country_length[country]
        return false
    end
    ## reformat iban
    if !is_iban_base97(iban)
        return false
    end
    return true
end

function split_iban(iban::String)
    iban = unify_iban(iban)
    m = match(iban_global_re, iban)
    if m == nothing
        error("IBAN '$(iban)' not according to format")
    end
    return m.captures
end

### iban validation utils

function unify_iban(iban::String)
    iban = replace(iban, " ", "")
    iban = replace(iban, ",", "")
    iban = uppercase(iban)
    return iban
end

function is_iban_base97(iban::String)
    ## reformat
    rban = *(iban[5:end], iban[1:4]) ## preallocate and fill instead
    ## base97 check
    digit_ban = [string(pos2number(x)) for x in rban]
    digit_int = parseint(BigInt, join(digit_ban))
    d, rest = divrem(digit_int, 97)
    return rest == 1
end

function pos2number(char::Char)
    if isdigit(char)
        ind = convert(Int, char) - 48
    elseif isalpha(char)
        ind = convert(Int, char) - 55
    else
        error("Not a valid IBAN entry: '$(char)'")
    end
    return ind
end

function construct_bban_re(pattern)
    re = join([*(re_format[f], "{", string(n), "}") for (n, f) in pattern])
    re = Regex(re)
    return re
end

bban_length(pattern) = sum([n for (n, f) in pattern])
iban_length(pattern) = bban_length(pattern) + 4


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE
using Test, Tensors
dtime = 0.25
parameters = ChabocheParameterState(E = 200.0e3,
                                    nu = 0.3,
                                    R0 = 100.0,
                                    Kn = 100.0,
                                    nn = 10.0,
                                    C1 = 10000.0,
                                    D1 = 100.0,
                                    C2 = 50000.0,
                                    D2 = 1000.0,
                                    Q = 50.0,
                                    b = 0.1)
mat = Chaboche(parameters = parameters)
times = [mat.drivers.time]
stresses = [copy(tovoigt(mat.variables.stress))]
dstrain11 = 1e-3*dtime
dstrain12 = 1e-3*dtime

dtimes = [dtime, dtime, dtime, dtime, 1.0]
dstrains11 = [dstrain11, dstrain11, dstrain11, -dstrain11, -4*dstrain11]
dstrains12 = [dstrain12, dstrain12, dstrain12, -dstrain12, -4*dstrain12]
plasticity_test = zeros(length(dstrains11)-1)
for i in 1:length(dtimes)
    dstrain11 = dstrains11[i]
    dtime = dtimes[i]
    biaxial_increment!(mat, dstrain11, dstrain12, dtime)
    update_material!(mat)
    push!(stresses, copy(tovoigt(mat.variables.stress)))
    if i > 1
        plasticity_test[i-1] = mat.variables.cumeq > 0.0
    end
    @test !iszero(tovoigt(mat.variables.stress)[1]) && !iszero(tovoigt(mat.variables.stress)[end])
    @test isapprox(tovoigt(mat.variables.stress; offdiagscale=2.0)[2:5],zeros(4); atol=1e-5)
end
@test any(x->x==1,plasticity_test)


text:
#########################################################################################################
#=
Glyph type for Text rendering.
It doesn't offer any functionality, and is only used for multiple dispatch.
=#
immutable GLGlyph{T} <: AbstractFixedVector{T, 4}
  glyph::T
  line::T
  row::T
  style_group::T
end

function GLGlyph(glyph::Integer, line::Integer, row::Integer, style_group::Integer)
  if !isascii(char(glyph))
    glyph = char('1')
  end
  GLGlyph{Uint16}(uint16(glyph), uint16(line), uint16(row), uint16(style_group))
end
function GLGlyph(glyph::Char, line::Integer, row::Integer, style_group::Integer)
  if !isascii(glyph)
    glyph = char('1')
  end
  GLGlyph{Uint16}(uint16(glyph), uint16(line), uint16(row), uint16(style_group))
end

GLGlyph() = GLGlyph(' ', typemax(Uint16), typemax(Uint16), 0)

Base.length{T}(::GLGlyph{T})                   = 4
Base.length{T}(::Type{GLGlyph{T}})             = 4
Base.eltype{T}(::GLGlyph{T})                   = T
Base.eltype{T}(::Type{GLGlyph{T}})             = T
Base.size{T}(::GLGlyph{T})                     = (4,)

Base.start{T}(::GLGlyph{T})                    = 1
Base.next{T}(x::GLGlyph{T}, state::Integer)    = (getfield(x, state), state+1)
Base.done{T}(x::GLGlyph{T}, state::Integer)    = state > 4

import Base: (+)

function (+){T}(a::Array{GLGlyph{T}, 1}, b::GLGlyph{T})
  for i=1:length(a)
    a[i] = a[i] + b
  end
end
function (+){T}(a::GLGlyph{T}, b::GLGlyph{T})
  GLGlyph{T}(a.glyph + b.glyph, a.line + b.line, a.row + b.row, a.style_group + b.style_group)
end
#########################################################################################################


text:
# Tested in Julia 1.5.1
using DelimitedFiles, Combinatorics
data = readdlm("./data/adventInput09.txt", Int64)
function ValidOr(preamble, nextNum; d = 2)
    pairs = permutations(preamble, d)
    c = 0
    for p in pairs
        c += 1
        if p[1] != p[2] && sum(p) == nextNum
            return true
            break
        elseif c == length(pairs)
            return false
        end
    end
end
function findNotValid(data; idx = 1, lastStep = 25)
    for i in 1:length(data)
        preamble = data[idx:lastStep]
        nextNum = data[lastStep + 1]
        if ValidOr(preamble, nextNum; d = 2) == false
            return nextNum
            break
        end
        lastStep += 1
        idx += 1
    end
end
function encryptionWeakness(data)
    notvalid = findNotValid(data)
    for i in 1:length(data)
        sumInvalid = data[i]
        for j in i+1:length(data)
            sumInvalid += data[j]
            if sumInvalid == notvalid
                return minimum(data[i:j]) + maximum(data[i:j])
                break
            end
        end
    end
end
encryptionWeakness(data)


text:
# Common functions for dynamics and rendering of pendulums
function pendulum_dynamics(env, s, a, x = isnothing(env.px) ? 0 : rand(env.px); rng::AbstractRNG = Random.GLOBAL_RNG)
    # Deal with terminal states
    # println("failur thresh: ", env.failure_thresh, " val: ", abs(s[1]))
    if (isnothing(env.failure_thresh) ?  false : abs(s[1]) > env.failure_thresh)
        # println("here")
        return fill(-100, size(s)), 0
    elseif env.include_time_in_state && (s[3] > env.maxT || s[3] ≈ env.maxT)
        return fill(100, size(s)), 0
    end

    θ, ω = s[1], s[2]
    dt, g, m, l = env.dt, env.g, env.m, env.l

    a = a[1]
    a = clamp(a, -env.max_torque, env.max_torque)
    costs = angle_normalize(θ)^2 + 0.1f0 * ω^2 + 0.001f0 * a^2

    a = a + env.ashift
    a = a + x

    ω = ω + (-3. * g / (2 * l) * sin(θ + π) + 3. * a / (m * l^2)) * dt
    θ = angle_normalize(θ + ω * dt)
    ω = clamp(ω, -env.max_speed, env.max_speed)

    if env.include_time_in_state
        sp = Float32.([θ, ω, s[3] + dt])
    else
        sp = Float32.([θ, ω])
    end
    r = Float32(env.Rstep - env.λcost*costs)
    return sp, r
end

angle_normalize(x) = mod((x+π), (2*π)) - π

## POMDP implementation -- Allows for pixel observation
@with_kw struct PendulumPOMDP <: POMDP{Array{Float32}, Float64, Array{Float32}}
    failure_thresh::Union{Nothing, Float64} = nothing # if set, defines the operating range fo the pendulum. Episode terminates if abs(θ) is larger than this.
    θ0 = Distributions.Uniform(-π, π) # Distribution to sample initial angular position
    ω0 = Distributions.Uniform(-1., 1.) # Distribution to sample initial angular velocity
    Rstep = 0 # Reward earned on each step of the simulation
    λcost = 1 # Coefficient to the traditional OpenAIGym Reward
    max_speed::Float64 = 8.
    max_torque::Float64 = 2.
    dt::Float64 = .05
    g::Float64 = 10.
    m::Float64 = 1.
    l::Float64 = 1.
    γ::Float32 = 0.99
    ashift::Float64 = 0.0
    actions::Vector{Float64} = [-1., 1.]
    observation_fn = (s) -> s
    include_time_in_state = false
    maxT = 99*dt
    px = nothing # Distribution over disturbances
end

InvertedPendulumPOMDP(failure_thresh = deg2rad(20),
                    θ0 = Distributions.Uniform(-failure_thresh/2., failure_thresh/2.),
                    ω0 = Distributions.Uniform(-.1, .1),
                    Rstep = 1,
                    λcost = 1;
                    kwargs...) = PendulumPOMDP(failure_thresh = failure_thresh, θ0 = θ0, ω0 = ω0, Rstep = Rstep, λcost = λcost; kwargs...)

ImageInvertedPendulum(;dt = .05, kwargs...) = InvertedPendulumPOMDP(observation_fn=(s)->simple_render_pendulum(s, dt=dt); dt=dt, kwargs...)

function POMDPs.gen(mdp::PendulumPOMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, r = pendulum_dynamics(mdp, s, a, rng=rng)
    (sp = sp, o=rand(rng, observation(mdp, sp)), r = r)
end

function POMDPs.gen(mdp::PendulumPOMDP, s, a, x, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, r = pendulum_dynamics(mdp, s, a, x, rng=rng)
    (sp = sp, o=rand(rng, observation(mdp, sp)), r = r)
end

function POMDPs.observation(mdp::PendulumPOMDP, s)
    o = mdp.observation_fn(s)
    Deterministic(Float32.(o))
end

function POMDPs.initialstate(mdp::PendulumPOMDP)
    ImplicitDistribution((rng) -> Float32.([rand(rng, mdp.θ0), rand(rng, mdp.ω0)]))
end

POMDPs.initialobs(mdp::PendulumPOMDP, s) = observation(mdp, s)

POMDPs.actions(mdp::PendulumPOMDP) = mdp.actions
isfailure(mdp::PendulumPOMDP, s) = s[1] == -100
# Terminate if there is a failure or if we exceed the maximum time (in the case that we are tracking such things)
POMDPs.isterminal(mdp::PendulumPOMDP, s) = isfailure(mdp, s) || s[1] == 100
POMDPs.discount(mdp::PendulumPOMDP) = mdp.γ

render(mdp::PendulumPOMDP, s, a::AbstractArray) = render(mdp, s, a...)
render(mdp::PendulumPOMDP, s, a = 0) = render_pendulum(mdp, s, a)

## MDP formulation
@with_kw struct PendulumMDP <: MDP{Array{Float32}, Array{Float32}}
    failure_thresh::Union{Nothing, Float64} = nothing # if set, defines the operating range fo the pendulum. Episode terminates if abs(θ) is larger than this.
    θ0 = Distributions.Uniform(-π, π) # Distribution to sample initial angular position
    ω0 = Distributions.Uniform(-1., 1.) # Distribution to sample initial angular velocity
    Rstep = 0 # Reward earned on each step of the simulation
    λcost = 1 # Coefficient to the traditional OpenAIGym Reward
    max_speed::Float64 = 8.
    max_torque::Float64 = 2.
    ashift::Float64 = 0.0
    dt::Float64 = .05
    g::Float64 = 10.
    m::Float64 = 1.
    l::Float64 = 1.
    γ::Float32 = 0.99
    actions::Vector{Float64} = [-1., 1.]
    include_time_in_state = false
    maxT = 99*dt
    px = nothing # Distribution over disturbances
end

InvertedPendulumMDP(failure_thresh = deg2rad(20),
                    θ0 = Distributions.Uniform(-failure_thresh/2., failure_thresh/2.),
                    ω0 = Distributions.Uniform(-.1, .1),
                    Rstep = 1,
                    λcost = 1;
                    kwargs...) = PendulumMDP(failure_thresh = failure_thresh, θ0 = θ0, ω0 = ω0, Rstep = Rstep, λcost = λcost; kwargs...)

function POMDPs.gen(mdp::PendulumMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, r = pendulum_dynamics(mdp, s, a, rng=rng)
    (sp = sp, r = r)
end

# Adversarial version
function POMDPs.gen(mdp::PendulumMDP, s, a, x, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, r = pendulum_dynamics(mdp, s, a, x, rng=rng)
    (sp = sp, r = r)
end

function POMDPs.initialstate(mdp::PendulumMDP)
    if mdp.include_time_in_state
        ImplicitDistribution((rng) -> Float32.([rand(rng, mdp.θ0), rand(rng, mdp.ω0), 0f0]))
    else
        ImplicitDistribution((rng) -> Float32.([rand(rng, mdp.θ0), rand(rng, mdp.ω0)]))
    end
end

POMDPs.actions(mdp::PendulumMDP) = mdp.actions

isfailure(mdp::PendulumMDP, s) = s[1] == -100
# Terminate if there is a failure or if we exceed the maximum time (in the case that we are tracking such things)
POMDPs.isterminal(mdp::PendulumMDP, s) = isfailure(mdp, s) || s[1] == 100
POMDPs.discount(mdp::PendulumMDP) = mdp.γ

render(mdp::PendulumMDP, s, a::AbstractArray) = render(mdp, s, a...)
render(mdp::PendulumMDP, s, a = 0) = render_pendulum(mdp, s, a)

## Rendering
function render_pendulum(env, s, a)
    θ = s[1] + π/2.
    point_array = [(0.5,0.5), (0.5 + 0.3*cos(θ), 0.5 - 0.3*sin(θ))]

    a_rad = abs(a)/10.
    if a < 0
        θstart = -3π/4
        θend = -θstart
        θarr = θend
    else
        θstart = π/4
        θend = -θstart
        θarr = θstart
    end

    # Draw the arrow
    endpt = (0.5, 0.5) .+ a_rad.*(cos(θarr), sin(θarr))
    uparr = endpt .+ 0.1*a_rad.*(cos(θarr)-sign(a)*sin(θarr), sign(a)*cos(θarr)+sin(θarr))
    dwnarr = endpt .+ 0.1*a_rad.*(-cos(θarr)-sign(a)sin(θarr), sign(a)*cos(θarr)-sin(θarr))
    arrow_array = [[endpt, uparr], [endpt, dwnarr]]

    img = compose(context(),
        (context(), line(arrow_array), arc(0.5, 0.5, a_rad, θstart, θend),  linewidth(0.5mm), fillopacity(0.), stroke("red")),
        (context(), circle(0.5, 0.5, 0.01), fill("blue"), stroke("black")),
        (context(), line(point_array), stroke("black"), linewidth(1mm)),
        (context(), rectangle(), fill("white"))
    )
    tmpfilename = tempname()
    img |> PNG(tmpfilename, 10cm, 10cm)
    load(tmpfilename)
end

function simple_render_pendulum(state; show_prev=true, dt=1, stride=4, noise=nothing)
    θ_curr = state[1]

    curr_frame = downsample(full_resolution_pendulum(θ_curr); stride)
    if show_prev
        θ_prev = θ_curr - state[2]*dt
        prev_frame = downsample(full_resolution_pendulum(θ_prev); stride)
        curr_frame = cat(curr_frame, prev_frame, dims=3)
    end

    if !isnothing(noise)
        curr_frame .+= rand(noise, size(curr_frame))
        curr_frame = clamp.(curr_frame, 0, 1)
    end

    # return  Gray.(reverse(curr_frame, dims=2))'
    return Gray.(curr_frame)
end

function full_resolution_pendulum(θ; intensity = 1)
    rot_mat_cw = @SArray [cos(θ) -sin(θ); sin(θ) cos(θ)]
    height = 70
    width = 10
    # center = height + width
    # shift = @SArray [center, width]

    # Empty matrix to store full size image before cropping and downsampling
    pic = ones(4*width, height)
    shift = [Int(size(pic, 1) / 2), 0]

    # Loop through all pendulum indices
    for point in Iterators.product(-width/2:width/2, 0:height)
        # Rotate to correct angle
        rotated_point = round.(Int, rot_mat_cw * collect(point))
        # Shift to location in pic array (can't have negative indices)
        shifted_point = rotated_point .+ shift
        # Make that index black
        if shifted_point[1] > 1 && shifted_point[1] <= size(pic, 1) &&
            shifted_point[2] > 1 && shifted_point[2] <= size(pic, 2)
            pic[shifted_point[1], shifted_point[2]] = 1 - intensity
        end
    end

    # Crop to 400x600
    # cropped_pic = pic[center - 200:center + 199, height - 599:height]

    return reverse(pic, dims=2)'
end

function downsample(pic; stride = 50, avg = true)
    height = convert(Int64, ceil(size(pic, 1) / stride))
    width = convert(Int64, ceil(size(pic, 2) / stride))

    downsampled_pic = zeros(height, width)

    for i = 1:width
        colmin = stride * (i - 1) + 1
        colmax = min(stride * i, size(pic, 2))
        for j = 1:height
            rowmin = stride * (j - 1) + 1
            rowmax = min(stride * j, size(pic, 1))
            if avg
                downsampled_pic[j, i] = mean(pic[rowmin:rowmax, colmin:colmax])
            else
                downsampled_pic[j, i] = minimum(pic[rowmin:rowmax, colmin:colmax])
            end
        end
    end

    return downsampled_pic
end


text:
struct Trace
  ir::IR
  stack::Vector{Any}
  primitives
  nodes::IdDict{Union{Partial,Shape},Variable}
end

Trace(P) = Trace(IR(), [], P, IdDict())

function node!(tr::Trace, T::Union{Partial,Shape}, v)
  tr.nodes[T] = v
  return
end

struct TraceError
  error
  stack
end

function Base.showerror(io::IO, e::TraceError)
  print(io, "Tracing Error: ")
  showerror(io, e.error)
  for Ts in e.stack
    print(io, "\nin ", Tuple{widen.(Ts)...})
  end
end

rename(env, ex) = IRTools.prewalk(
  x -> x isa GlobalRef ? getfield(x.mod, x.name) :
       x isa Variable ? env[x] : x, ex)

returntype(ir) = exprtype(ir, returnvalue(IRTools.blocks(ir)[end]))

function unapply!(tr, Ts, args)
  if VERSION > v"1.4-" && Ts[1] isa AType{typeof(Core._apply_iterate)}
    Ts = [Const(Core._apply), Ts[3:end]...]
    args = [Core._apply, args[3:end]...]
  end
  Ts[1] isa AType{typeof(Core._apply)} || return Ts, args
  Ts′ = Any[Ts[2]]
  args′ = Any[args[2]]
  for (x, T) in zip(args[3:end], Ts[3:end])
    len = partial(tr.primitives, Const(length), T).value
    for i = 1:len
      t = partial(tr.primitives, Const(getindex), T, Const(i))
      push!(Ts′, t)
      ex = haskey(tr.nodes, t) ? tr.nodes[t] :
            push!(tr.ir, stmt(xcall(getindex, x, i), type = t))
      push!(args′, ex)
    end
  end
  return unapply!(tr, Ts′, args′)
end

function copyblock!(ir::IR, b)
  c = IRTools.block!(ir)
  env = Dict()
  for (arg, T) in zip(arguments(b), argtypes(b))
    env[arg] = argument!(c, nothing, T, insert = false)
  end
  for (v, st) in b
    env[v] = push!(c, rename(env, st))
  end
  for br in IRTools.branches(b)
    push!(branches(c), rename(env, br))
  end
  return c
end

function abstract_group(cfg, i)
  j = IRTools.Inner.idoms(cfg', entry = length(cfg))[i]
  group = Int[]
  group!(b) =
    (b == j || b in group) || (push!(group, b); foreach(group!, cfg[b]))
  foreach(group!, cfg[i])
  return sort!(group), j
end

function extract_group(bl)
  out = IR()
  args = filter(x -> x isa Variable, union(arguments.(branches(block(bl.ir, bl.id)))...))
  args = union(args, filter(x -> x != nothing, map(x -> x.condition, branches(block(bl.ir, bl.id)))))
  innerargs = Dict(a => argument!(out, Union{}) for a in args)
  for br in branches(block(bl.ir, bl.id))
    push!(branches(block(out, 1)), rename(innerargs, br))
  end
  bs, after = abstract_group(CFG(bl.ir), bl.id)
  bmap = Dict()
  for b in bs
    bmap[b] = copyblock!(out, block(bl.ir, b)).id
  end
  ret = IRTools.block!(out)
  bmap[after] = ret.id
  for (arg, T) in zip(arguments(block(bl.ir, after)), argtypes(block(bl.ir, after)))
    argument!(ret, nothing, T, insert = false)
  end
  for b in blocks(out), (i, br) in enumerate(branches(b))
    branches(b)[i] = IRTools.Branch(br, block = bmap[br.block])
  end
  IRTools.domorder!(out)
  return out, args, after
end

function inline!(out, ir, args)
  env = Dict(zip(arguments(ir), args))
  before = blocks(out)[end]
  for br in branches(block(ir, 1))
    push!(branches(before), rename(env, br))
  end
  bmap = Dict()
  for bl in blocks(ir)[2:end]
    bmap[bl.id] = copyblock!(out, bl).id
  end
  for b in before.id:blocks(out)[end].id, (i, br) in enumerate(branches(block(out, b)))
    branches(block(out, b))[i] = IRTools.Branch(br, block = bmap[br.block])
  end
end

function abstract!(tr, bl, env)
  ir, args, after = extract_group(bl)
  for i = 1:length(args)
    argtypes(ir)[i] = exprtype(tr.ir, rename(env, args[i]))
  end
  infer!(Inference(Frame(ir), tr.primitives))
  inline!(tr.ir, ir, rename.((env,), args))
  for (k, v) in zip(arguments(block(bl.ir, after)), arguments(blocks(tr.ir)[end]))
    env[k] = v
  end
  return after
end

function openbranches(out, env, bl)
  brs = []
  for br in branches(bl)
    br.condition == nothing && (push!(brs, br); break)
    cond = exprtype(out, env[br.condition])
    cond == Const(true) && continue
    cond == Const(false) && (push!(brs, br); break)
    push!(brs, br)
  end
  return brs
end

_nodetype(x, T::AType) = T
_nodetype(x, T::Type) = Node{T}(x)
nodetype(ir::IR, x) = _nodetype(x, exprtype(ir, x))

function replacement(P, args, Ts)
  r = instead(P, args, Ts...)
  r === nothing && return (args, Ts)
  return r
end

function traceblock!(tr::Trace, env, bl)
  for (k, v) in bl
    ex = v.expr
    if isexpr(ex, :call)
      Ts = map(v -> nodetype(tr.ir, rename(env, v)), ex.args)
      Ts, args = unapply!(tr, Ts, rename(env, ex).args)
      Ts[1] === Const(Base.not_int) && (Ts[1] = Const(!))
      args, Ts = replacement(tr.primitives, args, Ts)
      if (T = partial(tr.primitives, Ts...)) != nothing
        if T isa Node && !effectful(Ts...)
          env[k] = T.value
        elseif haskey(tr.nodes, T) && !effectful(Ts...)
          env[k] = tr.nodes[T]
        else
          env[k] = push!(tr.ir, stmt(Expr(:call, args...), type = T))
          T isa Union{Partial,Shape} && (tr.nodes[T] = env[k])
        end
      else
        env[k] = tracecall!(tr, args, Ts...)
      end
    elseif isexpr(ex, :meta)
    elseif isexpr(ex, :boundscheck)
      env[k] = true
    elseif isexpr(ex)
      error("Can't trace through $(ex.head) expression")
    else
      env[k] = rename(env, ex)
    end
  end
end

function trace!(tr::Trace, ir, args)
  if ir.meta.method.isva
    nargs = ir.meta.method.nargs
    splat = args[nargs:end]
    splat = push!(tr.ir, stmt(xcall(tuple, splat...),
                            type = ptuple(nodetype.((tr.ir,), splat)...)))
    args = [args[1:nargs-1]..., splat]
  end
  env = Dict{Any,Any}(zip(arguments(ir), args))
  bl = 1
  while true
    traceblock!(tr, env, block(ir, bl))
    brs = openbranches(tr.ir, env, block(ir, bl))
    if length(brs) == 1
      isreturn(brs[1]) && return rename(env, returnvalue(brs[1]))
      bl = brs[1].block
      foreach((a, b) -> env[a] = rename(env, b), arguments(block(ir, bl)), arguments(brs[1]))
    else
      bl = abstract!(tr, block(ir, bl), env)
    end
  end
end

function tracecall!(tr::Trace, args, Ts...)
  push!(tr.stack, Ts)
  ir = IR(widen.(Ts)...)
  ir == nothing && error("No IR for $(Tuple{widen.(Ts)...})")
  ir = ir |> merge_returns! |> prepare_ir!
  result = trace!(tr, ir, args)
  pop!(tr.stack)
  return result
end

"""
    trace(P, Ts...)

Trace the method signature `Ts` using the primtive set `P`. e.g.

    julia> trace(Defaults(), typeof(+), Int, Int)
    1: (%1 :: typeof(+), %2 :: Int64, %3 :: Int64)
      %4 = (%1)(%2, %3) :: Int64
      return %4

    julia> trace(Defaults(), typeof(+), Const(2), Const(2))
    1: (%1 :: typeof(+), %2 :: const(2), %3 :: const(2))
      return 4
"""
function trace(P, Ts...)
  tr = Trace(P)
  try
    argnames = [argument!(tr.ir, T) for T in Ts]
    for (T, x) in zip(Ts, argnames)
      T isa Union{Partial,Shape} && node!(tr, T, x)
    end
    args = [T isa Const ? T.value : arg for (T, arg) in zip(Ts, argnames)]
    args, Ts = replacement(P, args, Ts)
    if (T = partial(tr.primitives, Ts...)) != nothing
      return!(tr.ir, push!(tr.ir, stmt(Expr(:call, args...), type = T)))
    else
      return!(tr.ir, tracecall!(tr, args, Ts...))
    end
    return cleanup!(tr.ir)
  catch e
    throw(TraceError(e, tr.stack))
  end
end

atype(T::AType) = T
atype(x) = Const(x)

function tracem(P, ex)
  @capture(ex, f_(args__)) || error("@trace f(args...)")
  :(trace($(esc(P)), atype.(($(esc(f)), $(esc.(args)...)))...))
end

"""
    @trace f(args...)
    @trace P f(args...)

Get a typed trace for `f`, analagous to `@code_typed`. Note that unlike
`@code_typed`, you probably want to pass types rather than values, e.g.

    julia> @trace Int+Int
    1: (%1 :: const(+), %2 :: Int64, %3 :: Int64)
      %4 = (+)(%2, %3) :: Int64
      return %4

If you instead pass actual integers, Mjolnir will aggressively
constant-propagate them, resulting in a trivial trace.

    julia> @trace 2+2
    1: (%1 :: const(+), %2 :: const(2), %3 :: const(2))
      return 4

`P` is the primitive set used, which is `Mjolnir.Defaults()` by default.
"""
macro trace(P, ex)
  tracem(P, ex)
end

macro trace(ex)
  tracem(Defaults(), ex)
end


text:
@testset "14.longest-common-prefix.jl" begin
    @test longest_common_prefix(["flower", "flow", "flight"]) == "fl"
    @test longest_common_prefix(["dog", "racecar", "car"]) == ""
    @test longest_common_prefix(["reflower", "flow", "flight"]) == ""
end


text:
# Make sure the ncvreg R library is loaded
load_R_library("ncvreg")

const NcvregDefaultParams = Dict(
    :family => "gaussian", # "binomial", "posisson"
    :alpha => 1.0,
    :nlambda => 100,
    :eps => 0.01,
    :max_iter => 1000,
    :convex => true,
    :warn => true,
    :returnX => false,

    # The following params are set based on size of input array
    :lambda_min => 0.001,  # 0.001 if n > p, 0.05 otherwise

    # The following params are set based on which type of regression it is
    :penalty => "SCAD",    # "MCP", "lasso"
    :gamma => 3.7         # 3.7 when SCAD and 3.0 otherwise
)

# Both SCAD and MCP is implemented via ncvreg R library
abstract NcvRegression <: RegressionLearner

function learn(r::NcvRegression, dataset::DataSet)
    X = indepvars_as_row_major_matrix(dataset)
    y = depvar_as_vector(dataset)

    # Create params for this specific call based on default ones plus the ones for this
    # regression datum.
    default_params = copy(NcvregDefaultParams)
    set_parameters!(r, default_params, dataset)
    params = merge(default_params, r.params)

    # Now call ncvreg and copy back the result as a Dict
    retval = R"""
      p <- $params
      ncvreg($X, $y, family=p$family, penalty=p$penalty,
        gamma=p$gamma,
        alpha=p$alpha, lambda.min=p$lambda_min, nlambda=p$nlambda,
        lambda=p$lambda, eps=p$eps, max.iter=p$max_iter,
        convex=p$convex, warn=p$warn, returnX=p$returnX)
    """
    rcopy(retval)
end

function coefficients(r::NcvRegression, dataset::DataSet; includeIntercept = true)
    coeffs = learn(r, dataset)[:beta]
    includeIntercept ? coeffs : coeffs[2:end]
end

type ScadRegression <: NcvRegression
    params::Dict{Symbol, Any}
    ScadRegression(; params...) = new(Dict(params))
end

function set_parameters!(s::ScadRegression, params, ds::DataSet)
    params[:penalty] = "SCAD"
    params[:gamma] = 3.7
    set_lambda_min_param!(params, num_indepvars(ds), num_cases(ds))
    params
end

function set_lambda_min_param!(params, p::Int, n::Int)
    if n > p
        params[:lambda_min] = 0.001
    else
        params[:lambda_min] = 0.05
    end
end

type McpRegression <: NcvRegression
    params::Dict{Symbol, Any}
    McpRegression(; params...) = new(Dict(params))
end

function set_parameters!(s::McpRegression, params, ds::DataSet)
    params[:penalty] = "MCP"
    params[:gamma] = 3.0
    set_lambda_min_param!(params, num_indepvars(ds), num_cases(ds))
    params
end


text:
d1 = float([ 0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
d2 = sparse(float([ 0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
@test a_star(g3, 1, 4, d1) ==
    a_star(g4, 1, 4, d1) ==
    a_star(g3, 1, 4, d2)
@test a_star(g4, 4, 1) == nothing


text:
using GeodesicTracer
using ComputedGeodesicEquations
using GeodesicBase
using CarterBoyerLindquist

using StaticArrays
using Optim

# include("/Users/lucyackland-snow/Desktop/DevEnv1/dev/AccretionFormulae/src/johannsen.jl")
include("hunting-circular-orbits.jl")
include("alt-emissivity-deviations.jl")

using Plots

# import GeodesicBase, GeodesicTracer
# function GeodesicBase.metric(m::GeodesicTracer.AbstractAutoDiffStaticAxisSymmetricParams{T}, u) where {T}
#     rθ = @SVector [u[2], u[3]]
#     comps = GeodesicTracer.metric_components(m, rθ)
#     @SMatrix [
#     comps[1] 0 0 comps[5]
#     0 comps[2] 0 0
#     0 0 comps[3] 0
#     comps[5] 0 0 comps[4]
#     ]
# end

a=0.0
M=1.0
# R_isco = CarterBoyerLindquist.rms(M, a)
r_range = 4.0:0.1:8.0
# m = BoyerLindquist(M=M, a=a)

e_lz = find_orbit_range(;m=m, r_range = r_range, a = a, upper = 0.1)

E_vals = []

for i in 1:length(r_range)
        E = e_lz[i, 1]
        push!(E_vals, E)
end

plt = plot(r_range, E_vals, xlabel="r (M)", ylabel="E")
display(plt)

r_isco_num = r_range[findmin(E_vals)[2]]
r_isco_true = AccretionFormulae.r_isco(a, M)

@show(r_isco_num)
@show(r_isco_true)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

julia_version = v"1.6.0"

name = "OpenSpiel"
version = v"0.3.6"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/deepmind/open_spiel/archive/v0.3.1.tar.gz", "3d2a7d1c2fb29e3d0df6f70c8692cb6b922a840bd58a86a3e7ef25e509cdbccb"),
    ArchiveSource("https://github.com/findmyway/dds/archive/v0.1.1.tar.gz", "fd13ee77feb5b5c3dfcc3333a0523266beb2a3d27715703cf508313af25306e5"),
    ArchiveSource("https://github.com/abseil/abseil-cpp/archive/20200923.1.tar.gz", "808350c4d7238315717749bab0067a1acd208023d41eaf0c7360f29cc8bc8f21"),
    ArchiveSource("https://github.com/findmyway/hanabi-learning-environment/archive/v0.1.0.tar.gz", "6126936fd13a95f8cadeacaa69dfb38a960eaf3bd588aacc8893a6e07e4791a3"),
    ArchiveSource("https://github.com/findmyway/project_acpc_server/archive/v0.1.0.tar.gz", "e29f969dd62ba354b7019cae3f7f1dbfbd9a744687ea4a8f7494c2bb1ee87382"),
]

# Bash recipe for building across all platforms
script = raw"""

mv open_spiel-0.3.1 open_spiel
mv abseil-cpp-20200923.1 open_spiel/open_spiel/abseil-cpp
mv dds-0.1.1 open_spiel/open_spiel/games/bridge/double_dummy_solver
mv hanabi-learning-environment-0.1.0 open_spiel/open_spiel/games/hanabi/hanabi-learning-environment
mv project_acpc_server-0.1.0 open_spiel/open_spiel/games/universal_poker/acpc

mkdir open_spiel/build
cd open_spiel/build
export BUILD_WITH_PYTHON=OFF BUILD_WITH_JULIA=ON BUILD_WITH_HANABI=ON BUILD_WITH_ACPC=OFF
cmake \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DJulia_PREFIX=${prefix} \
    ../open_spiel/
make -j${nproc}
make install
install_license ${WORKSPACE}/srcdir/open_spiel/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
]
platforms = expand_cxxstring_abis(platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
    Dependency("libcxxwrap_julia_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"9",
    julia_compat="$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Test
using Diban
using DataFrames

tbl = DataFrame(
    int32 = Int32[-1, 0, 1],
    int64 = Int64[-10, 0, 10],
    float32 = Float32[-0.5, 0, 0.5],
    float64 = Float64[-0.5, 0, 0.5],
    bool = [true, false, true],
    string = ["abc", "def", "ghi"],
    int32m = Union{Missing, Int32}[-1, missing, 1],
    int64m = Union{Missing, Int64}[-10, missing, 10],
    float32m = Union{Missing, Float32}[-0.5, missing, 0.5],
    float64m = Union{Missing, Float64}[-0.5, missing, 0.5],
    boolm = Union{Missing, Bool}[true, missing, false],
    stringm = Union{Missing, String}["abc", missing, "ghi"],
)

path = "tmp.parquet"
#path = "c:/scratch\\tmp.parquet"
# write_parquet(path, tbl; compression_codec = "lz4")
# read_parquet(path)

write_parquet(path, tbl; compression_codec = "snappy")
read_parquet(path)

write_parquet(path, tbl; compression_codec = "uncompressed")
read_parquet(path)

write_parquet(path, tbl; compression_codec = "zstd")
read_parquet(path)

write_parquet(path, tbl; compression_codec = "gzip")
read_parquet(path)

try
    rm(path)
catch
    # do nothing
    sleep(8)
    rm(path)
end


text:
using MetabolicEP_ToyImplementation
using Test

@testset "MetabolicEP_ToyImplementation.jl" begin
    # Write your tests here.
end