text:
# --- PRIVATE METHODS ----------------------------------------------------------- #
function parse_vff_reaction_section(reaction_section_buffer::Array{String,1})

    # initialize -
    reaction_dictionary = Dict{String,CBMetabolicReaction}()
    reaction_order_array = Array{String,1}()

    # process each line -
    for reaction_line in reaction_section_buffer

        # skip comments and empty lines -
        if (occursin("//", reaction_line) == false && isempty(reaction_line) == false)

            # ok, record: name,{ec},LHS,RHS,R,F

            # initialize -
            reaction_wrapper = CBMetabolicReaction()

            # split around ,
            reaction_record_components_array = split(reaction_line,",")

            # get the name/key -
            reaction_key = reaction_record_components_array[1]

            # cache this (we need to keep the order that appears in the VFF) -
            push!(reaction_order_array, reaction_key)

            # populate -
            reaction_wrapper.record = reaction_line
            reaction_wrapper.reaction_name = reaction_key
            reaction_wrapper.ec_number = reaction_record_components_array[2]
            reaction_wrapper.left_phrase = reaction_record_components_array[3]
            reaction_wrapper.right_phrase = reaction_record_components_array[4]
            reaction_wrapper.reverse = reaction_record_components_array[5]
            reaction_wrapper.forward = reaction_record_components_array[6]

            # cache -
            reaction_dictionary[reaction_key] = reaction_wrapper
        end
    end

    # return -
    return (reaction_order_array, reaction_dictionary)
end

function parse_vff_rules_section(rules_section_buffer::Array{String,1})::Dict{String,CBRule}

    # initialize -
    rules_dictionary = Dict{String,CBRule}()

    # parse rule structure -
    for rule in rules_section_buffer

        # build new wrapper -
        rule_wrapper = CBRule()

        # split around =
        rule_components_array = split(rule,"=")

        # populate the wrapper -
        reaction_name = rule_components_array[1]
        rule_wrapper.reaction_name = reaction_name
        rule_wrapper.rule_text = rule_components_array[2]

        # cache -
        rules_dictionary[reaction_name] = rule_wrapper
    end

    # return -
    return rules_dictionary
end

function parse_vff_gene_order_section(gene_order_section::Array{String,1})::Array{String,1}

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

    # main -
    for gene_symbol in gene_order_section

        # also, grab this symbol for a global list -
        push!(gene_order_array, gene_symbol)
    end

    # return -
    return gene_order_array
end

function parse_vff_reaction_gene_section(reaction_gene_map_section_buffer::Array{String,1})::Dict{String,CBReactionGeneMap}

    # initialize -
    reaction_gene_map = Dict{String,CBReactionGeneMap}()

    # process each of the gene records -
    for gene_record in reaction_gene_map_section_buffer

        # split around =
        gene_record_component_array = split(gene_record,"=")

        # get the reaction tag -
        metabolic_reaction_tag = gene_record_component_array[1]

        # gene data field -
        gene_data_field = gene_record_component_array[2]
        if (isempty(gene_data_field) == false)

            # ok, we have a gene list -

            # Build a mapping object -
            mapping_object = CBReactionGeneMap()
            mapping_object.reaction_name = metabolic_reaction_tag
            gene_symbol_set = Set{String}()

            # split around ::
            local_gene_symbol_list = split(gene_data_field,"::")
            for gene_symbol in local_gene_symbol_list
                if gene_symbol != "[]"

                    # push into local cache -
                    push!(gene_symbol_set, gene_symbol)
                end
            end
            mapping_object.genes = gene_symbol_set

            # insert thos record into the map -
            reaction_gene_map[metabolic_reaction_tag] = mapping_object
        else

            # we don't have a gene record, this is an error.
            # even if we don't have a gene, we should have ::[]
            throw(ErrorException("missing gene record for reaction $(metabolic_reaction_tag)"))
        end
    end

    # return -
    return reaction_gene_map
end

function parse_vff_metabolite_section(metabolite_section_buffer::Array{String,1})::Tuple{Dict{String,CBMetabolite}, Array{String,1}}

    # initialize -
    metabolite_dictionary = Dict{String,CBMetabolite}()
    metabolite_symbol_order_array = Array{String,1}()

    # iterate through each metabolite -
    for metabolite_record in metabolite_section_buffer

        # skip comments and empty lines -
        if (occursin("//", metabolite_record) == false && isempty(metabolite_record) == false)

            # split around =
            metabolite_record_component_array = split(metabolite_record,"=")

            # metabolite symbol is [1]
            metabolite_symbol = metabolite_record_component_array[1]

            # grab -
            push!(metabolite_symbol_order_array, metabolite_symbol)

            # go ..
            metabolite_name_field = metabolite_record_component_array[2]
            if (isempty(metabolite_name_field) == false)

                # split along :: -
                metabolite_name_symbol_array = split(metabolite_name_field,"::")
                full_metabolite_name = metabolite_name_symbol_array[1]
                kegg_id = metabolite_name_symbol_array[2]

                # build wrapper -
                metabolite_wrapper = CBMetabolite()
                metabolite_wrapper.symbol = metabolite_symbol
                metabolite_wrapper.name = full_metabolite_name
                metabolite_wrapper.kegg_id = kegg_id

                # cache -
                metabolite_dictionary[metabolite_symbol] = metabolite_wrapper
            end
        end
    end

    # return -
    return (metabolite_dictionary, metabolite_symbol_order_array)
end

function build_metabolite_symbol_array(reaction_array::Array{CBMetabolicReaction,1})

  # Method variables -
  species_symbol_array = Array{CBMetabolite,1}()

  # Helper function to parse the reaction phrases, split out the symbols
  function _parse_phrase(reaction_phrase::String)

    # Method variables -
    local_species_array = Array{CBMetabolite,1}()

    # Split around + -
    fragment_array = split(reaction_phrase,"+")
    for fragment in fragment_array

      if (contains(fragment,"*"))

          local_fragment_array = split(fragment,"*");
          species_symbol = CBMetabolite();
          species_symbol.symbol = local_fragment_array[2];

      else

        # Build -
        species_symbol = CBMetabolite()
        species_symbol.symbol = fragment
      end

      # check - is this the []?
      if (species_symbol.symbol != "[]")
          # grab -
          push!(local_species_array,species_symbol)
      end
    end

    # return -
    return local_species_array
  end

  function _isequal(species_model_1::CBMetabolite,species_model_2::CBMetabolite)
    if (species_model_1.symbol == species_model_2.symbol)
      return true
    end
    return false
  end

  function _add_symbol!(species_symbol_array,species_symbol)

    contains_species_already = false
    for cached_species_model in species_symbol_array

      if (_isequal(cached_species_model,species_symbol))
        contains_species_already = true
        break
      end
    end

    if (contains_species_already == false)
      push!(species_symbol_array,species_symbol)
    end
  end

  # iterate through and get the symbols -
  for reaction in reaction_array

    tmp_species_array_left = _parse_phrase(reaction.left_phrase)
    tmp_species_array_right = _parse_phrase(reaction.right_phrase)
    append!(tmp_species_array_left,tmp_species_array_right)

    for species_model in tmp_species_array_left
      _add_symbol!(species_symbol_array,species_model)
    end
  end

  # ok, so the species symbol array is *not* sorted)
  # let's sort the species -
  tmp_array = String[]
  for species_symbol::CBMetabolite in species_symbol_array
      push!(tmp_array,species_symbol.symbol)
  end

  # generate permutation array -
  idxa_sorted = sortperm(tmp_array)
  sorted_symbol_array = species_symbol_array[idxa_sorted]
  sorted_tmp_array = tmp_array[idxa_sorted]

  # ok, if a species contains [], then put it at the end -
  c_symbol_array = CBMetabolite[]
  e_symbol_array = CBMetabolite[]
  m_symbol_array = CBMetabolite[]
  r_symbol_array = CBMetabolite[]
  x_symbol_array = CBMetabolite[]
  t_symbol_array = CBMetabolite[]
  no_compartment_array = CBMetabolite[]
  partitioned_symbol_array = CBMetabolite[]

  # find all the species in [e]
  idx_all_e = findall(x->occursin("[e]",x)==true,sorted_tmp_array)
  for index in idx_all_e
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species in [c]
  idx_all_c = findall(x->occursin("[c]",x)==true,sorted_tmp_array)
  for index in idx_all_c
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species in [m]
  idx_all_m = findall(x->occursin("[m]",x)==true,sorted_tmp_array)
  for index in idx_all_m
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species in [r]
  idx_all_r = findall(x->occursin("[r]",x)==true, sorted_tmp_array)
  for index in idx_all_r
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species in [x]
  idx_all_x = findall(x->occursin("[x]",x)==true,sorted_tmp_array)
  for index in idx_all_x
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species in [t]
  idx_all_t = findall(x->occursin("[t]",x)==true,sorted_tmp_array)
  for index in idx_all_t
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all species w/no compartment -
  idx_all_nc = findall(x->occursin(r"\[*\]",x)==false,sorted_tmp_array)
  for index in idx_all_nc
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # return -
  return partitioned_symbol_array
end

function build_stoichiometric_matrix(species_symbol_array::Array{String,1},reaction_array::Array{CBMetabolicReaction,1})

  # Method variables -
  number_of_species = length(species_symbol_array)
  number_of_reactions = length(reaction_array);
  stoichiometric_matrix = zeros(number_of_species,number_of_reactions);

  function _parse_reaction_phrase(lexeme,reaction_phrase)

    # Split around + -
    coefficient = 0.0;
    fragment_array = split(reaction_phrase,"+")
    for fragment in fragment_array

      if (contains(fragment,"*"))
          local_fragment_array = split(fragment,"*");
          test_lexeme = local_fragment_array[2];

          if (lexeme == test_lexeme)
            coefficient = parse(Float64,local_fragment_array[1]);
            break
          end

      else

        # Build -
        test_lexeme = fragment;
        if (lexeme == test_lexeme)
          coefficient = 1.0;
          break
        end
      end
    end

    return coefficient;
  end

  function _find_stoichiometric_coefficient(species_symbol::String,reaction::CBMetabolicReaction)

    # Method variables -
    stoichiometric_coefficient = 0.0

    # Check the left and right phrase -
    stoichiometric_coefficient += -1.0*(_parse_reaction_phrase(species_symbol,reaction.left_phrase))
    stoichiometric_coefficient += _parse_reaction_phrase(species_symbol,reaction.right_phrase)
    return stoichiometric_coefficient;
  end

  # setup counters -
  for (row_index,species_symbol) in enumerate(species_symbol_array)
    for (col_index,reaction) in enumerate(reaction_array)

      # Is this species involved in this reaction?
      stoichiometric_matrix[row_index,col_index] = _find_stoichiometric_coefficient(species_symbol,reaction);

    end
  end

  # return -
  return stoichiometric_matrix
end
# ------------------------------------------------------------------------------ #

# --- PUBLIC METHODS ----------------------------------------------------------- #
function load_vff_model_file(path_to_vff_file::String)::Dict{String,Any}

    # check the file path -
    is_file_path_ok(path_to_vff_file)

    # initialize -
    problem_dictionary = Dict{String,Any}()

    # load the file into an array -
    file_buffer_array = read_file_from_path(path_to_vff_file)

    # parse the reaction section -
    reaction_section = extract_section(file_buffer_array,"#REACTION::START","#REACTION::END")
    (reaction_order_array, reaction_dictionary) = parse_vff_reaction_section(reaction_section)
    problem_dictionary["list_of_reactions"] = reaction_dictionary
    problem_dictionary["reaction_order_array"] = reaction_order_array

    # parse the rules section -
    rule_section = extract_section(file_buffer_array,"#RULES::START","#RULES::END")
    problem_dictionary["list_of_rules"] = parse_vff_rules_section(rule_section)

    # parse the metabolite section -
    metabolite_section = extract_section(file_buffer_array,"#METABOLITES::START","#METABOLITES::END")
    (metabolite_dictionary, metabolite_symbol_order_array) = parse_vff_metabolite_section(metabolite_section)
    problem_dictionary["list_of_metabolites"] = metabolite_dictionary
    problem_dictionary["metabolite_symbol_order_array"] = metabolite_symbol_order_array

    # parse the reaction-gene-mapping section -
    gene_reaction_section = extract_section(file_buffer_array,"#REACTION-GENE-MAP::START","#REACTION-GENE-MAP::END")
    gene_reaction_dictionary = parse_vff_reaction_gene_section(gene_reaction_section)
    problem_dictionary["list_of_genes"] = gene_reaction_dictionary

    # parse the gene-order section -
    gene_order_section = extract_section(file_buffer_array,"#GENE-SYMBOL-ORDER::START","#GENE-SYMBOL-ORDER::END")
    problem_dictionary["gene_order_array"] = parse_vff_gene_order_section(gene_order_section)

    # return -
    return problem_dictionary
end

function generate_vff_from_cobra_dictionary(cobra_dictionary::Dict{String,Any}, path_to_vff_file::String, kegg_organism_code::String; path_to_model_mapping_files::Union{String,Nothing}=nothing)

    # checks -
    is_dir_path_ok(path_to_vff_file)

    # initialize -
    vff_buffer = String[]

    # get the stoichiometric array, and some other stuff
    stoichiometric_matrix = Matrix(cobra_dictionary["S"])
    list_of_reaction_tags = cobra_dictionary["rxns"]
    list_of_genes = cobra_dictionary["genes"]
    list_of_species = cobra_dictionary["mets"]
    list_of_reversible_reactions = cobra_dictionary["rev"]
    list_of_rules = cobra_dictionary["rules"]
    list_of_species_names = cobra_dictionary["metNames"]

    # get the rxn gene mapping matrix -
    RGM = Matrix(cobra_dictionary["rxnGeneMat"])

    # do we have a reaction_to_ec_mapping file?
    reaction_to_ec_mapping_dict = Dict{String,Mapping}()
    if (path_to_model_mapping_files != nothing)

        # create a path to the reaction ec mapping file -
        path_to_reaction_to_ec_mapping_file = joinpath(path_to_model_mapping_files,"ReactionECNumberMap.dat")

        # check path -
        is_file_path_ok(path_to_reaction_to_ec_mapping_file)

        # Ok, build the dictionary -
        reaction_to_ec_mapping_dict = build_mapping_dictionary(path_to_reaction_to_ec_mapping_file)
    end

    # what is the size?
    (number_of_species,number_of_reactions) = size(stoichiometric_matrix)

    # Declare a progress meter for user feedback -
    p = Progress(number_of_reactions,color=:yellow)

    # Add a header -
    filename = splitdir(path_to_vff_file)[2]
    +(vff_buffer,"// ----------------------------------------------------------- //")
    +(vff_buffer,"// $(filename)")
    +(vff_buffer,"// GENERATED BY: CBModelTools")
    +(vff_buffer,"// GENERATED ON: $(Dates.now())")
    +(vff_buffer,"// SOURCE: https://github.com/varnerlab/CBModelTools")
    +(vff_buffer,"//")
    +(vff_buffer,"// VFF (Varner Flat File) format is a simple text based format")
    +(vff_buffer,"// for encoding metabolic (or other) types of biological models.")
    +(vff_buffer,"// Metabolic VFFs have three sections:")
    +(vff_buffer,"// RULES:")
    +(vff_buffer,"// METABOLITES:")
    +(vff_buffer,"// REACTIONS:")
    +(vff_buffer,"// ----------------------------------------------------------- //")
    +(vff_buffer,"")

    # gene-reaction-map section
    +(vff_buffer,"#REACTION-GENE-MAP::START ------------------------------------ //")

    # main loop -
    for reaction_index = 1:number_of_reactions

        # what is the tag for this reaction?
        reaction_tag = list_of_reaction_tags[reaction_index]

        # initialize -
        record = ""
        record *= "$(reaction_tag)="

        # does this reaction have associated genes?
        idx_genes = findall(x->x==1,RGM[reaction_index,:])
        local_gene_set = Set{String}()
        if (!isempty(idx_genes))

            # ok, we have genes, grab them -
            gene_id_array = list_of_genes[idx_genes]

            # process each gene -
            for gene_id in gene_id_array

                # make the KEGG gene id -
                kegg_gene_id = "$(kegg_organism_code):$(gene_id)"

                # cache -
                push!(local_gene_set,kegg_gene_id)
            end
        end

        if (isempty(local_gene_set) == false)

            # make an ec record -
            gene_record = ""
            for gene_number in local_gene_set
                gene_record*="$(gene_number)::"
            end

            # remove trailing ,
            gene_record*="[]"

            # add the record -
            record*="$(gene_record)"

            # push into buffer and go around again -
            +(vff_buffer,record)
        end
    end

    +(vff_buffer,"#REACTION-GENE-MAP::END -------------------------------------- //")
    +(vff_buffer,"")

    # rules section -
    +(vff_buffer,"#RULES::START ------------------------------------------------ //")
    for (index,rule) in enumerate(list_of_rules)

        # init empty line -
        line = ""

        # get tag and reaction text -
        reaction_tag = list_of_reaction_tags[index]
        rule_text = list_of_rules[index]
        if (isempty(rule_text) == true)
            # build the line -
            line *= "$(reaction_tag)=[]"
        else
            # build the line -
            line *= "$(reaction_tag)=$(rule_text)"
        end

        # add rule to the buffer -
        +(vff_buffer,line)
    end

    +(vff_buffer,"#RULES::END -------------------------------------------------- //")
    +(vff_buffer,"")

    # Write the metabolites section -
    +(vff_buffer,"#METABOLITES::START ------------------------------------------ //")

    # build a kegg-id metabolite name map -
    mmt = build_metabolite_matching_table()
    for (index,metabolite) in enumerate(list_of_species)

        # init empty line -
        line = ""

        # get the metabolite name -
        metabolite_name = list_of_species_names[index]

        # lookup KEGG ID -
        kegg_id = "[]"
        if (haskey(mmt,metabolite_name) == true)
            kegg_id = mmt[metabolite_name]
        end

        # write -
        line*="$(metabolite)=$(metabolite_name)::$(kegg_id)"

        # add metabolite to the buffer -
        +(vff_buffer,line)
    end

    +(vff_buffer,"#METABOLITES::END -------------------------------------------- //")
    +(vff_buffer,"")

    # write reactions -
    +(vff_buffer,"#REACTION::START --------------------------------------------- //")
    for reaction_index = 1:number_of_reactions

        # initialize empty buffer -
        line = ""

        # get the reaction tag -
        reaction_tag_string = list_of_reaction_tags[reaction_index]

        # user message -
        msg = "Completed $(reaction_tag_string) ($(reaction_index) of $(number_of_reactions))"

        # update the progress bar -
        ProgressMeter.next!(p; showvalues = [(:status,msg)])

        # add the tag to the buffer -
        line *= "$(reaction_tag_string),"

        # ok, so I need to check, do we have is reaction key in my ec mapping?
        if (haskey(reaction_to_ec_mapping_dict,reaction_tag_string) == true)

            # we have a reaction tag = ec number record
            mapping_record = reaction_to_ec_mapping_dict[reaction_tag_string]

            # get the set of values -
            value_set = mapping_record.value
            for ec_number in value_set
                line*="$(ec_number)::"
            end

            # add [], -
            line*="[],"
        else
            # We have no ecnumber information -
            line *= "[],"
        end

        # find the reactants -
        idx_reactants = findall(stoichiometric_matrix[:,reaction_index].<0.0)
        if (isempty(idx_reactants) == true)
            line *= "[],"
        else

            # how many species do we have?
            number_of_species = length(idx_reactants)
            counter = 1
            for index in idx_reactants

                # get the metabolite -
                metabolite_string = list_of_species[index]
                stcoeff = stoichiometric_matrix[index,reaction_index]

                if (stcoeff != -1.0)
                    # add data to the buffer -
                    line *= "$(abs(stcoeff))*$(metabolite_string)"
                else
                    # add data to the buffer -
                    line *= "$(metabolite_string)"
                end

                # do we have more?
                if (counter < number_of_species)
                    line *= "+"
                else
                    line *= ","
                end

                counter = counter + 1
            end
        end

        # find the products -
        idx_products = findall(stoichiometric_matrix[:,reaction_index].>0.0)
        if (isempty(idx_products) == true)
            line *= "[],"
        else

            # how many species do we have?
            number_of_species = length(idx_products)
            counter = 1
            for index in idx_products

                # get the metabolite -
                metabolite_string = list_of_species[index]
                stcoeff = stoichiometric_matrix[index,reaction_index]

                if (stcoeff != 1.0)
                    # add data to the buffer -
                    line *= "$(stcoeff)*$(metabolite_string)"
                else
                    # add data to the buffer -
                    line *= "$(metabolite_string)"
                end

                # do we have more?
                if (counter < number_of_species)
                    line *= "+"
                else
                    line *= ","
                end

                counter = counter + 1
            end
        end

        # is this reaction reversible?
        rev_value = list_of_reversible_reactions[reaction_index]
        if (rev_value == 1.0)
            line *= "-inf,inf"
        else
            line *= "0,inf"
        end

        # add buffer to list of strings -
        +(vff_buffer,line)
    end
    +(vff_buffer,"#REACTION::END ----------------------------------------------- //")
    +(vff_buffer,"")

    # genes section -
    +(vff_buffer,"#GENE-SYMBOL-ORDER::START ------------------------------------ //")
    for gene_symbol in list_of_genes

        # make a record -
        record = ""
        record *= "$(kegg_organism_code):$(gene_symbol)"

        # write the record -
        +(vff_buffer,record)
    end

    +(vff_buffer,"#GENE-SYMBOL-ORDER::END -------------------------------------- //")

    # Write out the vff file -
    open("$(path_to_vff_file)", "w") do f

        for line_item in vff_buffer
            write(f,"$(line_item)\n")
        end
    end
end
# ------------------------------------------------------------------------------ #

max_stars_repo_name: varnerlab/CBModelTools

max_issues_repo_name: varnerlab/CBModelTools

max_forks_repo_name: varnerlab/CBModelTools


text:
@testset "Mixed Potential Unit Tests" begin
    atoms = [
        :Ar => (@SVector [0.0, 0.0, 0.0])u"bohr",
        :Ar => (@SVector [0.5, 0.0, 0.0])u"bohr",
        :H => (@SVector [0.0, 0.5, 0.0])u"bohr",
        :H => (@SVector [0.5, 0.5, 0.0])u"bohr"
    ]
    box = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]u"bohr"
    system = periodic_system(atoms, box)

    p1 = MockEmpiricalPotential(1.0, 1.0, (:Ar,))
    p2 = MockAbstractPotential(2.0)
    mixed_potentials = [
        +p1,
        -p1,
        p1 + p2,
        p1 - p2,
        2.0 * p1,
        p2 / 2.0,
        2.0 * p1 - 1.0 * p2,
        3.0 * (p1 + p2)
    ]

    true_rcutoffs = [
        1.0,
        1.0,
        Inf,
        Inf,
        1.0,
        Inf,
        Inf,
        Inf
    ]
    true_energies = [
        potential_energy(system, p1),
        -potential_energy(system, p1),
        potential_energy(system, p1) + potential_energy(system, p2),
        potential_energy(system, p1) - potential_energy(system, p2),
        2 * potential_energy(system, p1),
        potential_energy(system, p2) / 2,
        2 * potential_energy(system, p1) - potential_energy(system, p2),
        3 * (potential_energy(system, p1) + potential_energy(system, p2))
    ]
    true_forces = [
        force(system, p1),
        -force(system, p1),
        force(system, p1) + force(system, p2),
        force(system, p1) - force(system, p2),
        2 * force(system, p1),
        force(system, p2) / 2,
        2 * force(system, p1) - force(system, p2),
        3 * (force(system, p1) + force(system, p2))
    ]
    true_virials = [
        virial(system, p1),
        -virial(system, p1),
        virial(system, p1) + virial(system, p2),
        virial(system, p1) - virial(system, p2),
        2 * virial(system, p1),
        virial(system, p2) / 2,
        2 * virial(system, p1) - virial(system, p2),
        3 * (virial(system, p1) + virial(system, p2))
    ]
    true_virial_stresses = [
        virial_stress(system, p1),
        -virial_stress(system, p1),
        virial_stress(system, p1) + virial_stress(system, p2),
        virial_stress(system, p1) - virial_stress(system, p2),
        2 * virial_stress(system, p1),
        virial_stress(system, p2) / 2,
        2 * virial_stress(system, p1) - virial_stress(system, p2),
        3 * (virial_stress(system, p1) + virial_stress(system, p2))
    ]

    for (p, trc, te, tf, tv, tvs) in zip(mixed_potentials, true_rcutoffs, true_energies, true_forces, true_virials, true_virial_stresses)
        @test p isa MixedPotential

        @test get_rcutoff(p) == trc
        @test ismissing(get_species(p))

        @test energy_and_force(system, p).e == te
        @test energy_and_force(system, p).f == tf

        @test potential_energy(system, p) == te
        @test force(system, p) == tf
        @test virial(system, p) == tv
        @test virial_stress(system, p) == tvs
    end
end

ext: jl


text:
notary = Notary()
cypher = CypherSuite(notary)

# Somewhere far far away

mixkey = newsigner(notary)

MIX_PORT = 3000
MIX_ID = id(mixkey)

mix = Mix(MIX_PORT, MIX_ID, cypher)
@async SynchronicBallot.serve(mix, mixkey)

sleep(1.)

# Using the mixer

guardiankey = newsigner(notary)
transactions = Trigger(TransactionVector(), t -> println(typeof(t.document)))
tlog = BraidChain(transactions, Guardian(id(guardiankey)), notary)

registrator = newsigner(notary)
reg = Registrator(id(registrator))
push!(tlog, Transaction(reg, guardiankey))

braider = newsigner(notary)

push!(tlog, Transaction(Braider(id(braider)), guardiankey))

gk = GateKeeper(3001, id(braider), cypher, mix)
sb = BraiderConfig(UInt8(3), UInt8(64), gk)
bo = BraidOfficer(sb, tlog)
botask = @async BraidChains.serve!(bo, braider)
sleep(1.)

### Registration of members ###

member1 = KeyChain(notary)
member2 = KeyChain(notary)
member3 = KeyChain(notary)

# on the server
push!(tlog, Transaction(Member(id(member1)), registrator))
push!(tlog, Transaction(Member(id(member2)), registrator))
push!(tlog, Transaction(Member(id(member3)), registrator))

update!(member1, tlog)
update!(member2, tlog)
update!(member3, tlog)

@show istaskstarted(botask)

### First Braid

@sync begin
    @async braid!(member1, bo)
    @async braid!(member2, bo)
    @async braid!(member3, bo)
end

### After that I shall repeat the procedure

sleep(1.)

update!(member1, tlog)
update!(member2, tlog)
update!(member3, tlog)

@test id(member1[end]) in pseudonyms(tlog)
@test id(member2[end]) in pseudonyms(tlog)
@test id(member3[end]) in pseudonyms(tlog)

### Second Braid

@sync begin
    @async braid!(member1, bo)
    @async braid!(member2, bo)
    @async braid!(member3, bo)
end

sleep(1.)

update!(member1, tlog)
update!(member2, tlog)
update!(member3, tlog)

@test id(member1[end]) in pseudonyms(tlog)
@test id(member2[end]) in pseudonyms(tlog)
@test id(member3[end]) in pseudonyms(tlog)


text:
const ALLOWED_FORMAT = ["gz", "zlib", "zz", "xz", "zst", "bz2"]

"""
    makestream(filepath)

Open a file with an appropriate codec determined by file suffix.
"""
function makestream(filepath, args...; kwargs...)
    io = open(filepath, args...; kwargs...)
    if endswith(filepath, ".gz")
        codec = isreadonly(io) ? GzipDecompressor() : GzipCompressor()
    elseif endswith(filepath, ".zlib")
        codec = isreadonly(io) ? ZlibDecompressor() : ZlibCompressor()
    elseif endswith(filepath, ".zz")
        codec = isreadonly(io) ? DeflateDecompressor() : DeflateCompressor()
    elseif endswith(filepath, ".xz")
        codec = isreadonly(io) ? XzDecompressor() : XzCompressor()
    elseif endswith(filepath, ".zst")
        codec = isreadonly(io) ? ZstdDecompressor() : ZstdCompressor()
    elseif endswith(filepath, ".bz2")
        codec = isreadonly(io) ? Bzip2Decompressor() : Bzip2Compressor()
    else
        return io
    end
    TranscodingStream(codec, io)
end

function makestream(f::Function, args...)
    io = makestream(args...)
    try
        f(io)
    finally
        close(io)
    end
end

"""
    checkplinkfilename(filename, suffix)

Check a filename is valid plink file name. Suffix should be "bed", "fam" or "bim".
"""
function checkplinkfilename(filename::AbstractString, suffix::AbstractString)
    suffix ∈ ["bed", "fam", "bim"] ||
    throw(ArgumentError("suffix should be bed, fam or bim"))
    endswith(filename, "." * suffix) ||
    any(endswith.(filename, ".$suffix." .* ALLOWED_FORMAT)) ||
    throw(ArgumentError("compressed format should be one of $ALLOWED_FORMAT"))
end

"""
    compresse_plink(plink_file, format[="gz"], outfile[=plink_file])

Compress a set of Plink files to `gz` (default), `zlib` or `zz` format.
"""
function compress_plink(
    plinkfile::AbstractString,
    format::AbstractString="gz",
    outfile::AbstractString=plinkfile
    )
    format ∈ ALLOWED_FORMAT || throw(ArgumentError("compressed format should be one of $ALLOWED_FORMAT"))
    for suffix in [".bed", ".fam", ".bim"]
        fname = plinkfile * suffix
        if isfile(fname)
            open(fname) do input
                makestream(fname * "." * format, "w") do output
                    write(output, input)
                end
            end
        end
    end
end

"""
    decompresse_plink(plink_compressed, outfile[=plink_compressed])

Decompress a set of compressed Plink files to plain Plink format.
"""
function decompress_plink(
    plink_compressed::AbstractString,
    format::AbstractString="gz",
    outfile::AbstractString=plink_compressed
    )
    format ∈ ALLOWED_FORMAT || throw(ArgumentError("compressed format should be one of $ALLOWED_FORMAT"))
    for suffix in [".bed", ".fam", ".bim"]
        iname = plink_compressed * suffix * "." * format
        oname = outfile * suffix
        if isfile(iname)
            open(oname, "w") do output
                makestream(iname) do input
                    write(output, input)
                end
            end
        end
    end
end


text:
using BenchmarkTools
using FixedPolynomials
import DynamicPolynomials: @polyvar

@polyvar x y
f1 = x + y + 1
f2 = (x + 1)^4 + y
f3 = 50x^3+83x^2*y+24x*y^2+y^3+392x^2+414x*y+50y^2-28x+59y-100
f4 = x+x^2+y+x*y^3+x^4*y^2+3x^3*y+10*x*y+10x^2*y+10x*y^2+15x^2*y^2+10x^3*y^2

w = rand(2)

function make(f, w)
    p = Polynomial{Float64}(f)
    result = @benchmark evaluate($p, $w)
    println(STDOUT, f)
    show(STDOUT, MIME"text/plain"(), result)
    println("\n")
end

println("FixedPolynomials", "\n")
make(f1, w)
make(f2, w)
make(f3, w)
make(f4, w)


text:
using Test
using Random
using LinearAlgebra
using SparseArrays
using MatrixDepot
using JLD2, FileIO
using MatrixMarket
using TraceEstimation

@testset "DiagonalApprox" begin
    Random.seed!(123543)
    @testset "sprandmatrices" begin
        println("Executing Test 01: sprand matrix")
        A = sprand(3000,3000,0.004)
        A = A + A' + 30I
        M = Matrix(A)
        act = tr(inv(M))
        obv = tr_inv(A, 4, 20)
        @test isapprox(act, obv, rtol = 0.2)
    end
    @testset "poissonmatrix" begin
        println("Executing Test 02: poisson matrix")
        A = matrixdepot("poisson", 50)
        M = Matrix(A)
        act = tr(inv(M))
        obv = tr_inv(A, 8, 40)
        @test isapprox(act, obv, rtol = 0.2)
    end
    @testset "wathenmatrix" begin
        println("Executing Test 03: wathen matrix")
        A = matrixdepot("wathen", 35)
        M = Matrix(A)
        act = tr(inv(M))
        obv = tr_inv(A, 8, 40, pc = "ilu")
        @test isapprox(act, obv, rtol = 0.2)
    end
    #=
    @testset "nasa2146" begin
        println("Executing Test 04: nasa2146 matrix")
        A = MatrixMarket.mmread("nasa2146.mtx")
        M = Matrix(A)
        act = tr(inv(M))
        obv = tr_inv(A, 8, 40)
        @test isapprox(act, obv, rtol = 0.2)
    end
    @testset "kuu" begin
        println("Executing Test 05: kuu matrix")
        A = MatrixMarket.mmread("Kuu.mtx")
        M = Matrix(A)
        act = tr(inv(M))
        obv = tr_inv(A, 8, 80)
        @test isapprox(act, obv, rtol = 0.2)
    end
    =#
    @testset "KMatrix" begin
        println("Executing Test 06: K matrix (high condition)")
        @load "topopt902.jld2" K
        K = SparseMatrixCSC(K)
        M = Matrix(K)
        act = tr(inv(M))
        obv = tr_inv(K, 8, 150)
        @test isapprox(act, obv, rtol = 0.2)
    end
end


text:
"""
sequential function, good for debuging
"""
function setindex_sequential!( ba::BigArray{D,T}, buf::Array{T,N},
                             idxes::Union{UnitRange, Int, Colon} ... ) where {D,T,N}
    idxes = colon2unit_range(buf, idxes)
    C = get_encoding(ba)
    mipLevelName = get_mip_level_name(ba)
    baIter = ChunkIterator(idxes, get_chunk_size(ba); offset=get_offset(ba))
    for (blockID, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) in baIter
        chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = adjust_volume_boundary(ba,
                                    chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer)
        @inbounds chk = buf[rangeInBuffer]
        ba.kvStore[ joinpath(mipLevelName, cartesian_range2string(chunkGlobalRange)) ] =
                                                                            encode( chk, C)
    end
end

"""
    getindex_sequential(ba::BigArray, idxes::Union{UnitRange, Int}...)
sequential implementation for debuging
"""
function getindex_sequential(ba::BigArray{D,T},
                             idxes::Union{UnitRange, Int}...) where {D,T,N}
    t1 = time()
    sz = map(length, idxes)
    buf = zeros(eltype(ba), sz)
    baRange = CartesianIndices(ba)
    C = get_encoding(ba)
    mipLevelName = get_mip_level_name(ba)
    baIter = ChunkIterator(idxes, get_chunk_size(ba); offset=get_offset(ba))
    for (blockId, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) in baIter
        if any(map((x,y)->x>y, first(globalRange).I, last(baRange).I)) ||
            any(map((x,y)->x<y, last(globalRange).I, first(baRange).I))
            @warn("out of volume range, keep it as zeros")
            continue
        end
        chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = adjust_volume_boundary(ba,
                                            chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer)
        chunkSize = (last(chunkGlobalRange) - first(chunkGlobalRange) + one(last(chunkGlobalRange))).I
        try
            #println("global range of chunk: $(cartesian_range2string(chunkGlobalRange))")
            key = cartesian_range2string(chunkGlobalRange)
            keyPath = joinpath(mipLevelName, cartesian_range2string(chunkGlobalRange))
            v = ba.kvStore[keyPath]
            v = Codings.decode(v, C)
            chk = reinterpret(T, v) |> Vector
            chk = reshape(chk, chunkSize)

            @inbounds buf[rangeInBuffer] = chk[rangeInChunk]
        catch err
            if isa(err, KeyError) && ba.fillMissing
                println("no suck key in kvstore: $(err), will fill this block as zeros")
                break
            else
                println("catch an error while getindex in BigArray: $err with type of $(typeof(err))")
                rethrow()
            end
        end
    end

    elapsed = time() - t1 # seconds
    println("cutout speed: $(sizeof(buf)/1024/1024/elapsed) MB/s")
    OffsetArray(buf, idxes...)
end


text:
"""
$(SIGNATURES)

Initialize Sod shock tube
"""
function ib_sod(
    set::AbstractSetup,
    ps::AbstractPhysicalSpace,
    vs::Union{AbstractVelocitySpace,Nothing},
    gas::AbstractProperty,
)

    if set.nSpecies == 1

        primL = [1.0, 0.0, 0.5]
        primR = [0.125, 0.0, 0.625]

        primL = ifelse(set.space[5:6] == "3v", [primL[1:2]; zeros(2); primL[end]], primL)
        primR = ifelse(set.space[5:6] == "3v", [primR[1:2]; zeros(2); primR[end]], primR)

        wL = prim_conserve(primL, gas.γ)
        wR = prim_conserve(primR, gas.γ)

        p = (
            x0 = ps.x0,
            x1 = ps.x1,
            wL = wL,
            wR = wR,
            primL = primL,
            primR = primR,
            γ = gas.γ,
        )

        fw = function (x, p)
            if x <= (p.x0 + p.x1) / 2
                return p.wL
            else
                return p.wR
            end
        end

        bc = function (x, p)
            if x <= (p.x0 + p.x1) / 2
                return p.primL
            else
                return p.primR
            end
        end

        if set.space[1:4] == "1d0f"
            return fw, bc, p
        elseif set.space == "1d1f1v"
            p = (p..., u = vs.u)
            ff = function (x, p)
                w = ifelse(x <= (p.x0 + p.x1) / 2, p.wL, p.wR)
                prim = conserve_prim(w, p.γ)
                h = maxwellian(p.u, prim)
                return h
            end

            return fw, ff, bc, p
        elseif set.space == "1d2f1v"
            p = (p..., u = vs.u, K = gas.K)
            ff = function (x, p)
                w = ifelse(x <= (p.x0 + p.x1) / 2, p.wL, p.wR)
                prim = conserve_prim(w, p.γ)
                h = maxwellian(p.u, prim)
                b = h * p.K / 2 / prim[end]
                return h, b
            end

            return fw, ff, bc, p
        elseif set.space == "1d1f3v"
            p = (p, u = vs.u, v = vs.v, w = vs.w)
            ff = function (x, p)
                w = ifelse(x <= (p.x0 + p.x1) / 2, p.wL, p.wR)
                prim = conserve_prim(w, p.γ)
                h = maxwellian(p.u, p.v, p.w, prim)
                return h
            end

            return fw, ff, bc, p
        end

    elseif set.nSpecies == 2

        primL = zeros(3, 2)
        primL[:, 1] .= [gas.mi, 0.0, gas.mi / 2]
        primL[:, 2] .= [gas.me, 0.0, gas.me / 2]

        primR = zeros(3, 2)
        primR[:, 1] .= [gas.mi * 0.125, 0.0, gas.mi / 1.6]
        primR[:, 2] .= [gas.me * 0.125, 0.0, gas.me / 1.6]

        wL = mixture_prim_conserve(primL, gas.γ)
        wR = mixture_prim_conserve(primR, gas.γ)

        p = (
            x0 = ps.x0,
            x1 = ps.x1,
            wL = wL,
            wR = wR,
            primL = primL,
            primR = primR,
            γ = gas.γ,
        )

        fw = function (x, p)
            if x <= (p.x0 + p.x1) / 2
                return p.wL
            else
                return p.wR
            end
        end

        bc = function (x, p)
            if x <= (p.x0 + p.x1) / 2
                return p.primL
            else
                return p.primR
            end
        end

        if set.space[1:4] == "1d0f"
            return fw, bc, p
        elseif set.space == "1d1f1v"
            hL = mixture_maxwellian(vs.u, primL)
            hR = mixture_maxwellian(vs.u, primR)
            p = (p..., hL = hL, hR = hR)
            ff = function (x, p)
                if x <= (p.x0 + p.x1) / 2
                    return p.hL
                else
                    return p.hR
                end
            end

            return fw, ff, bc, p
        elseif set.space == "1d2f1v"
            hL = mixture_maxwellian(vs.u, primL)
            hR = mixture_maxwellian(vs.u, primR)
            bL = similar(hL)
            bR = similar(hR)
            for j = 1:2
                bL[:, j] .= hL[:, j] .* K ./ (2.0 .* primL[end, j])
                bR[:, j] .= hR[:, j] .* K ./ (2.0 .* primR[end, j])
            end

            p = (p, hL = hL, hR = hR, bL = bL, bR = bR)

            ff = function (x, p)
                if x <= (p.x0 + p.x1) / 2
                    return p.hL, p.bL
                else
                    return p.hR, p.bR
                end
            end

            return fw, ff, bc, p
        end

    end

    return nothing

end


text:
module SkeletonCellFields

using Gridap

export SkeletonPair
export SkeletonCellBasis
export SkeletonCellVector
export SkeletonCellMatrix

export jump
export mean
import Gridap: restrict
import Gridap: inner
import Gridap: integrate
import Base: +, -, *
import Gridap: gradient
import Gridap: symmetric_gradient
import Base: div
import Gridap: trace
import Gridap: curl
import Gridap: outer

function restrict(
  cf::IndexCellFieldLike{Z,T,N},
  desc1::BoundaryDescriptor,
  desc2::BoundaryDescriptor) where {Z,T,N}

  cf1 = restrict(cf,desc1)
  cf2 = restrict(cf,desc2)
  SkeletonPair{Z-1,T,N}(cf1,cf2)
end

function restrict(
  cf::IndexCellFieldLike{D},trian::SkeletonTriangulation{Z}) where {D,Z}
  @assert D == Z + 1
  restrict(cf,trian.descriptor1,trian.descriptor2)
end

struct SkeletonPair{Z,T,N}
  cellfield1
  cellfield2
end

for op in (:+,:-,:(gradient),:(symmetric_gradient),:(div),:(trace),:(curl))
  @eval begin
    function ($op)(a::SkeletonPair{Z,T,N}) where {Z,T,N}
      cf1 = a.cellfield1
      cf2 = a.cellfield2
      SkeletonPair{Z,T,N}($op(cf1),$op(cf2))
    end
  end
end

for op in (:+, :-, :*, :outer)
  @eval begin

    function ($op)(a::SkeletonPair{Z,T,N},b::CellField) where {Z,T,N}
      cf1 = a.cellfield1
      cf2 = a.cellfield2
      SkeletonPair{Z,T,N}($op(cf1,b),$op(cf2,b))
    end

    function ($op)(a::CellField,b::SkeletonPair{Z,T,N}) where {Z,T,N}
      cf1 = b.cellfield1
      cf2 = b.cellfield2
      SkeletonPair{Z,T,N}($op(a,cf1),$op(a,cf2))
    end

  end
end

function jump(sp::SkeletonPair{Z,T,1}) where {Z,T}
  cf1 = sp.cellfield1
  cf2 = sp.cellfield2
  cf1 - cf2
end

function jump(sp::SkeletonPair{Z,T,2}) where {Z,T}
  cf1 = sp.cellfield1
  cf2 = sp.cellfield2
  SkeletonCellBasis{Z,T}(cf1, -cf2)
end

function mean(sp::SkeletonPair{Z,T,1}) where {Z,T}
  cf1 = sp.cellfield1
  cf2 = sp.cellfield2
  0.5*(cf1 + cf2)
end

function mean(sp::SkeletonPair{Z,T,2}) where {Z,T}
  cf1 = sp.cellfield1
  cf2 = sp.cellfield2
  SkeletonCellBasis{Z,T}(0.5*cf1, 0.5*cf2)
end

struct SkeletonCellBasis{Z,T}
  cellbasis1
  cellbasis2
end

function inner(a::SkeletonCellBasis{Z},b::CellField{Z}) where Z
  cb1 = a.cellbasis1
  cb2 = a.cellbasis2
  cm1 = varinner(cb1,b)
  cm2 = varinner(cb2,b)
  SkeletonVarinnerVector{Z}(cm1,cm2)
end

function inner(a::SkeletonCellBasis{Z},b::SkeletonCellBasis{Z}) where Z
  a1 = a.cellbasis1
  a2 = a.cellbasis2
  b1 = b.cellbasis1
  b2 = b.cellbasis2
  cm11 = varinner(a1,b1)
  cm12 = varinner(a1,b2)
  cm21 = varinner(a2,b1)
  cm22 = varinner(a2,b2)
  SkeletonVarinnerMatrix{Z}(cm11,cm12,cm21,cm22)
end

struct SkeletonVarinnerVector{D}
  cellmap1
  cellmap2
end

for op in (:+, :-)
  @eval begin

    function ($op)(a::SkeletonVarinnerVector{Z},b::SkeletonVarinnerVector{Z}) where Z
      c1 = $op(a.cellmap1,b.cellmap1)
      c2 = $op(a.cellmap2,b.cellmap2)
      SkeletonVarinnerVector{Z}(c1,c2)
    end

    function ($op)(a::SkeletonVarinnerVector{Z}) where Z
      c1 = $op(a.cellmap1)
      c2 = $op(a.cellmap2)
      SkeletonVarinnerVector{Z}(c1,c2)
    end

  end
end

function (*)(a::Real,b::SkeletonVarinnerVector{Z}) where Z
  c1 = a*b.cellmap1
  c2 = a*b.cellmap2
  SkeletonVarinnerVector{Z}(c1,c2)
end

struct SkeletonCellVector
  cellvector1
  cellvector2
end

struct SkeletonVarinnerMatrix{D}
  cellmap11
  cellmap12
  cellmap21
  cellmap22
end

for op in (:+, :-)
  @eval begin

    function ($op)(a::SkeletonVarinnerMatrix{Z},b::SkeletonVarinnerMatrix{Z}) where Z
      cellmap11 = $op(a.cellmap11, b.cellmap11)
      cellmap12 = $op(a.cellmap12, b.cellmap12)
      cellmap21 = $op(a.cellmap21, b.cellmap21)
      cellmap22 = $op(a.cellmap22, b.cellmap22)
      SkeletonVarinnerMatrix{Z}(
        cellmap11,
        cellmap12,
        cellmap21,
        cellmap22)
    end

    function ($op)(a::SkeletonVarinnerMatrix{Z}) where Z
      cellmap11 = $op(a.cellmap11)
      cellmap12 = $op(a.cellmap12)
      cellmap21 = $op(a.cellmap21)
      cellmap22 = $op(a.cellmap22)
      SkeletonVarinnerMatrix{Z}(
        cellmap11,
        cellmap12,
        cellmap21,
        cellmap22)
    end

  end
end

function (*)(a::Real,b::SkeletonVarinnerMatrix{Z}) where Z
      cellmap11 = a*b.cellmap11
      cellmap12 = a*b.cellmap12
      cellmap21 = a*b.cellmap21
      cellmap22 = a*b.cellmap22
      SkeletonVarinnerMatrix{Z}(
        cellmap11,
        cellmap12,
        cellmap21,
        cellmap22)
end

struct SkeletonCellMatrix
  cellmatrix11
  cellmatrix12
  cellmatrix21
  cellmatrix22
end

function integrate(
  v::SkeletonVarinnerVector{Z},
  trian::Triangulation{Z},
  quad::CellQuadrature{Z}) where Z

  v1 = integrate(v.cellmap1,trian,quad)
  v2 = integrate(v.cellmap2,trian,quad)
  SkeletonCellVector(v1,v2)
end

function integrate(
  v::SkeletonVarinnerMatrix{Z},
  trian::Triangulation{Z},
  quad::CellQuadrature{Z}) where Z

  m11 = integrate(v.cellmap11,trian,quad)
  m12 = integrate(v.cellmap12,trian,quad)
  m21 = integrate(v.cellmap21,trian,quad)
  m22 = integrate(v.cellmap22,trian,quad)
  SkeletonCellMatrix(m11,m12,m21,m22)
end

end # module


text:
# ---
# title: 1253. Reconstruct a 2-Row Binary Matrix
# id: problem1253
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Greedy
# link: <https://leetcode.com/problems/reconstruct-a-2-row-binary-matrix/description/>
# hidden: true
# ---
#
# Given the following details of a matrix with `n` columns and `2` rows :
#
#   * The matrix is a binary matrix, which means each element in the matrix can be `0` or `1`.
#   * The sum of elements of the 0-th(upper) row is given as `upper`.
#   * The sum of elements of the 1-st(lower) row is given as `lower`.
#   * The sum of elements in the i-th column(0-indexed) is `colsum[i]`, where `colsum` is given as an integer array with length `n`.
#
# Your task is to reconstruct the matrix with `upper`, `lower` and `colsum`.
#
# Return it as a 2-D integer array.
#
# If there are more than one valid solution, any of them will be accepted.
#
# If no valid solution exists, return an empty 2-D array.
#
#
#
# **Example 1:**
#
#
#
#     Input: upper = 2, lower = 1, colsum = [1,1,1]
#     Output: [[1,1,0],[0,0,1]]
#     Explanation: [[1,0,1],[0,1,0]], and [[0,1,1],[1,0,0]] are also correct answers.
#
#
# **Example 2:**
#
#
#
#     Input: upper = 2, lower = 3, colsum = [2,2,1,1]
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: upper = 5, lower = 5, colsum = [2,1,2,0,1,0,1,2,0,1]
#     Output: [[1,1,1,0,1,0,0,1,0,0],[1,0,1,0,0,0,1,1,0,1]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= colsum.length <= 10^5`
#   * `0 <= upper, lower <= colsum.length`
#   * `0 <= colsum[i] <= 2`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Carries out a Kepler step for bodies i & j
function keplerij!(m::Array{Float64,1},x::Array{Float64,2},v::Array{Float64,2},i::Int64,j::Int64,h::Float64,jac_ij::Array{Float64,2})
# The state vector has: 1 time; 2-4 position; 5-7 velocity; 8 r0; 9 dr0dt; 10 beta; 11 s; 12 ds
# Initial state:
state0 = zeros(Float64,12)
# Final state (after a step):
state = zeros(Float64,12)
delx = zeros(Float64,NDIM)
delv = zeros(Float64,NDIM)
# jac_ij should be the Jacobian for going from (x_{0,i},v_{0,i},m_i) &  (x_{0,j},v_{0,j},m_j)
# to  (x_i,v_i,m_i) &  (x_j,v_j,m_j), a 14x14 matrix for the 3-dimensional case.
# Fill with zeros for now:
fill!(jac_ij,0.0)
for k=1:NDIM
  state0[1+k     ] = x[k,i] - x[k,j]
  state0[1+k+NDIM] = v[k,i] - v[k,j]
end
gm = GNEWT*(m[i]+m[j])
# The following jacobian is just computed for the Keplerian coordinates (i.e. doesn't include
# center-of-mass motion, or scale to motion of bodies about their common center of mass):
jac_kepler = zeros(Float64,7,7)
kepler_step!(gm, h, state0, state, jac_kepler)
for k=1:NDIM
  delx[k] = state[1+k] - state0[1+k]
  delv[k] = state[1+NDIM+k] - state0[1+NDIM+k]
end
# Compute COM coords:
mijinv =1.0/(m[i] + m[j])
xcm = zeros(Float64,NDIM)
vcm = zeros(Float64,NDIM)
mi = m[i]*mijinv # Normalize the masses
mj = m[j]*mijinv
for k=1:NDIM
  xcm[k] = mi*x[k,i] + mj*x[k,j]
  vcm[k] = mi*v[k,i] + mj*v[k,j]
end
# Compute the Jacobian:
jac_ij[ 7, 7] = 1.0  # the masses don't change with time!
jac_ij[14,14] = 1.0
for k=1:NDIM
   jac_ij[   k,   k] +=   mi
   jac_ij[   k, 3+k] += h*mi
   jac_ij[   k, 7+k] +=   mj
   jac_ij[   k,10+k] += h*mj
   jac_ij[ 3+k, 3+k] +=   mi
   jac_ij[ 3+k,10+k] +=   mj
   jac_ij[ 7+k,   k] +=   mi
   jac_ij[ 7+k, 3+k] += h*mi
   jac_ij[ 7+k, 7+k] +=   mj
   jac_ij[ 7+k,10+k] += h*mj
   jac_ij[10+k, 3+k] +=   mi
   jac_ij[10+k,10+k] +=   mj
   for l=1:NDIM
# Compute derivatives of \delta x_i with respect to initial conditions:
     jac_ij[   k,   l] += mj*jac_kepler[  k,  l]
     jac_ij[   k, 3+l] += mj*jac_kepler[  k,3+l]
     jac_ij[   k, 7+l] -= mj*jac_kepler[  k,  l]
     jac_ij[   k,10+l] -= mj*jac_kepler[  k,3+l]
# Compute derivatives of \delta v_i with respect to initial conditions:
     jac_ij[ 3+k,   l] += mj*jac_kepler[3+k,  l]
     jac_ij[ 3+k, 3+l] += mj*jac_kepler[3+k,3+l]
     jac_ij[ 3+k, 7+l] -= mj*jac_kepler[3+k,  l]
     jac_ij[ 3+k,10+l] -= mj*jac_kepler[3+k,3+l]
# Compute derivatives of \delta x_j with respect to initial conditions:
     jac_ij[ 7+k,   l] -= mi*jac_kepler[  k,  l]
     jac_ij[ 7+k, 3+l] -= mi*jac_kepler[  k,3+l]
     jac_ij[ 7+k, 7+l] += mi*jac_kepler[  k,  l]
     jac_ij[ 7+k,10+l] += mi*jac_kepler[  k,3+l]
# Compute derivatives of \delta v_j with respect to initial conditions:
     jac_ij[10+k,   l] -= mi*jac_kepler[3+k,  l]
     jac_ij[10+k, 3+l] -= mi*jac_kepler[3+k,3+l]
     jac_ij[10+k, 7+l] += mi*jac_kepler[3+k,  l]
     jac_ij[10+k,10+l] += mi*jac_kepler[3+k,3+l]
   end
# Compute derivatives of \delta x_i with respect to the masses:
   jac_ij[   k, 7] += (x[k,i]+h*v[k,i]-xcm[k]-mj*state[1+k])*mijinv + GNEWT*mj*jac_kepler[  k,7]
   jac_ij[   k,14] += (x[k,j]+h*v[k,j]-xcm[k]+mi*state[1+k])*mijinv + GNEWT*mj*jac_kepler[  k,7]
# Compute derivatives of \delta v_i with respect to the masses:
   jac_ij[ 3+k, 7] += (v[k,i]-vcm[k]-mj*state[4+k])*mijinv + GNEWT*mj*jac_kepler[3+k,7]
   jac_ij[ 3+k,14] += (v[k,j]-vcm[k]+mi*state[4+k])*mijinv + GNEWT*mj*jac_kepler[3+k,7]
# Compute derivatives of \delta x_j with respect to the masses:
   jac_ij[ 7+k, 7] += (x[k,i]+h*v[k,i]-xcm[k]-mj*state[1+k])*mijinv - GNEWT*mi*jac_kepler[  k,7]
   jac_ij[ 7+k,14] += (x[k,j]+h*v[k,j]-xcm[k]+mi*state[1+k])*mijinv - GNEWT*mi*jac_kepler[  k,7]
# Compute derivatives of \delta v_j with respect to the masses:
   jac_ij[10+k, 7] += (v[k,i]-vcm[k]-mj*state[4+k])*mijinv - GNEWT*mi*jac_kepler[3+k,7]
   jac_ij[10+k,14] += (v[k,j]-vcm[k]+mi*state[4+k])*mijinv - GNEWT*mi*jac_kepler[3+k,7]
end
# Advance center of mass & individual Keplerian motions:
centerm!(m,mijinv,x,v,vcm,delx,delv,i,j,h)
return
end


text:
@testset "IR transform" begin

function check_ir(ir_in::IR, ir_out::IR)
    for (x, stmt) in ir_in
        e_in, e_out = stmt.expr, ir_out[x].expr
        if !isexpr(e_in, :call) continue end
        if unwrap(e_in.args[1]) == :rand
            if !(isexpr(e_out, :call) &&
                 e_out.args[1] == GlobalRef(Genify, :trace))
                error("$e_in is not traced")
            end
        elseif e_in.args[1] == GlobalRef(Core, :_apply)
            if !(isexpr(e_out, :call) &&
                 e_out.args[1] == GlobalRef(Core, :_apply) &&
                 e_out.args[2] == GlobalRef(Genify, :trace))
                error("$e_in is not traced")
            end
        end
    end
    return true
end

# Test straight line code
function foo(x::Int)
    p = rand(Beta(1, 1))
    probs = [p, 1-p]
    y = rand(Categorical(probs))
    return (x+y)
end

ir_in = IR(typeof(foo), Int)
ir_out = transform!(copy(ir_in))
@test check_ir(ir_in, ir_out)

# Test branching code
function foo(x::Int)
    if x > 0
        y = rand(Uniform(0, 1))
    else
        y = rand(Normal(0, 1))
    end
    return y
end

ir_in = IR(typeof(foo), Int)
ir_out = transform!(copy(ir_in))
@test check_ir(ir_in, ir_out)

# Test rand with varargs
function foo(dims...)
    p = rand(Beta(1, 1))
    probs = [p, 1-p]
    y = rand(Categorical(probs), dims...)
    return y
end

ir_in = IR(typeof(foo), Any)
ir_out = transform!(copy(ir_in))
@test check_ir(ir_in, ir_out)

end

@testset "Generative function building" begin

function foo(alpha::Real, beta::Real)
    p = rand(Beta(alpha, beta))
    coin = rand(Bernoulli(p))
    return coin
end

genfoo = genify(foo, Int, Float64)
julia_fn = genfoo.julia_function
@test Genify.signature(julia_fn) == Tuple{Real, Real}
@test genfoo.arg_types == [Real, Real]

# Test simulate
trace = simulate(genfoo, (1, 1))
@test 0 <= trace[:p] <= 1
@test trace[:coin] in [true, false]

# Test generate
trace, weight = generate(genfoo, (1, 1), choicemap(:p => 0.25, :coin => false))
@test trace[:p] == 0.25
@test trace[:coin] == false

# Test update
trace, weight, _, discard =
    update(trace, (1, 1), (NoChange(), NoChange()), choicemap(:coin => true))
@test trace[:coin] == true
@test discard[:coin] == false
@test weight == log(0.25) - log(0.75)

# Test gradient propagation
function bar(mu::Real, std::Real)
    z = rand(Normal(mu, std))
    return 2*z
end

genbar = genify(bar, Float64, Float64, gradients=:strict)
@test genbar.has_argument_grads == [true, true]
@test genbar.accepts_output_grad == true

trace, _ = generate(genbar, (0.0, 1.0), choicemap(:z => 0.0))
arg_grads, _, choice_grads = choice_gradients(trace, select(:z), 1.0)
@test arg_grads == (0.0, -1.0)
@test choice_grads[:z] == 2.0

# Test handling of parametric methods
function foo(x::Array{T,N}, start::T, tup::Tuple{Int,T}) where {T, N}
    y = rand(T, size(x))
    x[:] = y .+ start
    return x .* tup[1] .- tup[2]
end

genfoo = genify(foo, Vector{Float64}, Float64, Tuple{Int,Float64})
@test genfoo.arg_types == [Array, Float64, Tuple{Int, T} where T]
@test Genify.signature(genfoo.julia_function) == Tuple{genfoo.arg_types...}
@test all(0 .<= genfoo(zeros(5), 5., (2, 10.)) .<= 2)

end

@testset "Automatic addressing" begin

# Test straight line code
function foo(x::Int)
    p = rand(Beta(1, 1))
    probs = [p, 1-p]
    y = rand(Categorical(probs))
    return (x+y)
end

genfoo = genify(foo, Int; useslots=true)
choices, _, _ = propose(genfoo, (0,))
@test has_value(choices, :p) && has_value(choices, :y)

# Test branching code
function foo(x::Bool)
    if x
        y = rand(Uniform(0, 1))
    else
        y = rand(Normal(0, 1))
    end
end

genfoo = genify(foo, Bool; useslots=true)
choices, _, _ = propose(genfoo, (true,))
@test has_value(choices, :y_1) && !has_value(choices, :y_2)
choices, _, _ = propose(genfoo, (false,))
@test !has_value(choices, :y_1) && has_value(choices, :y_2)

# Test for loops
function foo(n::Int)
    for i in 1:n
        for j in 1:i
            x = rand(Normal(0.0, 1.0))
        end
    end
end

genfoo = genify(foo, Int; useslots=true)
choices, _, _ = propose(genfoo, (5, ))
@test all(choices[:x => i => j] isa Float64 for i in 1:5 for j in 1:i)
choices, _, _ = propose(genfoo, (10, ))
@test all(choices[:x => i => j] isa Float64 for i in 1:10 for j in 1:i)

# Test while loops
function collatz(n::Int)
    seq = Int[]
    while n != 1
        n = mod(n, 2) == 0 ? n ÷ 2 : 3*n + 1
        x = rand(1:n)
        push!(seq, n)
    end
    return seq
end

gencollatz = genify(collatz, Int; useslots=true)
choices, _, seq = propose(gencollatz, (13, ))
@test all(choices[:x => i] isa Int for i in 1:length(seq))
@test all(choices[:x => i] <= n for (i,n) in enumerate(seq))

# Test address generation from arguments
function foo(dist::UnivariateDistribution, callee)
    return callee() + rand(dist)
end

function bar() return rand() end

genfoo = genify(foo, UnivariateDistribution, Any; useslots=false)
choices, _, _ = propose(genfoo, (Normal(0, 1), bar))
@test (:callee, :rand_dist) ⊆ keys(nested_view(choices))

end

@testset "Manual addressing" begin

function model(mu::Real)
    θ = rand(here, submodel, mu) # Test splicing of address namespaces
    x = rand(:obs, observe, θ) # Test nesting of address namespaces
end

function submodel(mu::Real)
    return rand(:theta, LogNormal(mu, 1))
end

function observe(theta::Real)
    x = rand(:m1, Normal(theta, 1))
    y = rand(:m2, Normal(theta, 2))
    z = rand(:m3, Normal(theta, 3))
    return (x, y, z)
end

genmodel = genify(:manual, model, Real)

# Test simulate
trace = simulate(genmodel, (1,))
@test trace[:theta] >= 0
@test trace[:obs] isa Tuple{Float64, Float64, Float64}
@test trace[:obs => :m1] isa Float64
@test trace[:obs => :m2] isa Float64
@test trace[:obs => :m3] isa Float64

# Test generate
observations = choicemap((:obs => :m1, 3), (:obs => :m2, 3), (:obs => :m3, 3))
trace, weight = generate(genmodel, (1,), observations)
@test trace[:obs => :m1] == 3
@test trace[:obs => :m2] == 3
@test trace[:obs => :m3] == 3
@test weight ≈ sum(Gen.logpdf(normal, 3, trace[:theta], i) for i in 1:3)

# Test regenerate
prev_trace = trace
trace, weight, _ = regenerate(trace, (1,), (NoChange(),), select(:theta))
@test trace[:obs] == prev_trace[:obs]

prev_trace = trace
trace, weight, _ = regenerate(trace, (1,), (NoChange(),), select(:obs))
@test trace[:theta] == prev_trace[:theta]

end

@testset "Nested generative functions" begin

function model(mu::Real)
    theta = rand(LogNormal(mu, 1))
    obs = observe(theta)
end

function observe(theta::Real)
    m1 = rand(Normal(theta, 1))
    m2 = rand(Normal(theta, 2))
    m3 = rand(Normal(theta, 3))
    return (m1, m2, m3)
end

genmodel = genify(model, Real)

# Test simulate
trace = simulate(genmodel, (1,))
@test trace[:theta] >= 0
@test trace[:obs] isa Tuple{Float64, Float64, Float64}
@test trace[:obs => :m1] isa Float64
@test trace[:obs => :m2] isa Float64
@test trace[:obs => :m3] isa Float64

# Test generate
observations = choicemap((:obs => :m1, 3), (:obs => :m2, 3), (:obs => :m3, 3))
trace, weight = generate(genmodel, (1,), observations)
@test trace[:obs => :m1] == 3
@test trace[:obs => :m2] == 3
@test trace[:obs => :m3] == 3
@test weight ≈ sum(Gen.logpdf(normal, 3, trace[:theta], i) for i in 1:3)

# Test regenerate
prev_trace = trace
trace, weight, _ = regenerate(trace, (1,), (NoChange(),), select(:theta))
@test trace[:obs] == prev_trace[:obs]

prev_trace = trace
trace, weight, _ = regenerate(trace, (1,), (NoChange(),), select(:obs))
@test trace[:theta] == prev_trace[:theta]

# Test non-recursion
genmodel = genify(model, Real; recurse=false)
choices, _, _ = propose(genmodel, (1,))
@test :obs ∉ keys(nested_view(choices))

end

@testset "Memoized transformation" begin

function foo(alpha::Real, beta::Real)
    p = rand(Beta(alpha, beta))
    coin = rand(Bernoulli(p))
    return coin
end

genfoo1 = genified(foo, Real, Real)
genfoo2 = genify(foo, Real, Real)
genfoo3 = genified(foo, Real, Real)
genfoo4 = genify(foo, Real, Real)
genfoo5 = genified(:minimal, foo, Real, Real)
genfoo6 = genified(foo, Real, Real, recurse=false, useslots=false)

@test genfoo1 !== genfoo2
@test genfoo1 === genfoo3
@test genfoo2 !== genfoo4
@test genfoo1 !== genfoo5
@test genfoo5 === genfoo6

end


text:
# This sets up our own REPL backend as well as some frontend hooks
function supplant_repl_backend(;attempt=0)
    if !isdefined(Base, :active_repl_backend)
        if attempt > 0
            if attempt > 15
                warn("Unable to auto-supplant REPL backend; please try manually by calling supplant_repl_backend()")
                return
            end
            @schedule (sleep(0.1); supplant_repl_backend(attempt=attempt + 1))
            return
        else
            error("Base.active_repl_backend not defined!  This happens when running as a script; don't do that!")
        end
    end

    if all([!isdefined(mode, :on_ready) for mode in Base.active_repl.interface.modes])
        msg = """
        This Julia's REPL doesn't have on_ready hooks!  See README.md for more:
        https://github.com/staticfloat/MakePkgUpdatePrecompileInTheBackground.jl
        """
        error(msg)
    end

    # First, create our own backend
    repl_channel = Channel(1)
    response_channel = Channel(1)
    backend = REPL.REPLBackend(repl_channel, response_channel, false, nothing)
    backend.backend_task = @schedule begin
        global running_mutex
        unlock(running_mutex)
        while true
            tls = task_local_storage()
            tls[:SOURCE_PATH] = nothing
            ast, show_value = take!(backend.repl_channel)
            if show_value == -1
                # exit flag
                break
            end
            # Lock the mutex so stuff can't happen while the REPL is busy, as
            # that would be incredibly messy.  We unlock in on_ready()
            lock(running_mutex)
            REPL.eval_user_input(ast, backend)
        end
    end

    # Now swap it in for the old backend
    old_backend = Base.active_repl_backend
    eval(Base, :(global active_repl_backend; active_repl_backend = $backend))
    eval(Base, :(global active_repl; active_repl.backendref = REPL.REPLBackendRef($repl_channel,$response_channel)))

    # Hook into LineEdit mode `on_enter` and add a CTRL-X keymapping for debugging of cursor position
    failmodes = 0
    for mode in Base.active_repl.interface.modes
        if !isdefined(mode, :on_ready)
            failmodes += 1
            continue
        end

        const old_on_done = mode.on_done
        mode.on_done = (state, buffer, ok) -> begin
            global last_enter

            # If there were lines written before, let's blank them!
            if last_draw_lines > 0
                lock(running_mutex)
                # Skip over the input line(s)
                total_lines = prompt_lines(state)
                skip_lines = cursor_position(state)[1] + 1
                write(STDOUT,"$(CSI)$(skip_lines)A$(CSI)1G")
                for idx in 1:(last_draw_lines)
                    write(STDOUT,"$(CSI)1A$(CSI)0K")
                end
                write(STDOUT, "$(CSI)$(last_draw_lines+1)B$(CSI)1G")
                write(STDOUT, "\n"^(total_lines-1))
                unlock(running_mutex)
            end
            last_enter = time()
            old_on_done(state, buffer, ok)
        end

        const old_on_transition_start = mode.on_transition_start
        mode.on_transition_start = (state) -> begin
            global running_mutex
            lock(running_mutex)
            return old_on_transition_start(state)
        end

        const old_on_transition_stop = mode.on_transition_stop
        mode.on_transition_start = (state) -> begin
            global running_mutex
            unlock(running_mutex)
            return old_on_transition_stop(state)
        end

        # Use a hook specially added for this package to know when the prompt is
        # ready for more input and no longer computing the previous result
        const old_on_ready = mode.on_ready
        mode.on_ready = (state) -> begin
            global prompt_state, running_mutex

            # If we haven't snarfed the prompt_state yet, do so now
            if prompt_state == nothing
                prompt_state = state
            end

            # We're done running whatever we were running before, allow drawing
            unlock(running_mutex)
            return old_on_ready(state)
        end
    end

    # Ask the old backend to die, nicely
    put!(old_backend.repl_channel, (nothing, -1))

    return nothing
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/LibCURL_jll.jl

module LibCURL_jll

using Libdl
using nghttp2_jll

const PATH_list = String[]
const LIBPATH_list = String[]

export libcurl

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libcurl_handle = C_NULL
libcurl_path = ""

if Sys.iswindows()
    const libcurl = "libcurl-4.dll"
elseif Sys.isapple()
    const libcurl = "@rpath/libcurl.4.dylib"
else
    const libcurl = "libcurl.so"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libcurl_handle = dlopen(libcurl)
    global libcurl_path = dlpath(libcurl_handle)
end

is_available() = true

end  # module LibCURL_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
__precompile__()

module heuristicopt

export GA, GATool, SA

module GA
include("ga/ga.jl")
export GAFunction, GAParameter
end

module GATool
include("ga/crossover.jl")
include("ga/elite.jl")
include("ga/mutate.jl")
include("ga/select.jl")
end

module SA
include("sa/simulatedannealing.jl")
end

end

max_stars_repo_name: ariririri/heuristicopt

max_issues_repo_name: ariririri/heuristicopt


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

##########################
# Cholesky Factorization #
##########################

# The dispatch structure in the cholesky, and cholesky! methods is a bit
# complicated and some explanation is therefore provided in the following
#
# In the methods below, LAPACK is called when possible, i.e. StridedMatrices with Float32,
# Float64, Complex{Float32}, and Complex{Float64} element types. For other element or
# matrix types, the unblocked Julia implementation in _chol! is used. For cholesky
# and cholesky! pivoting is supported through a Val(Bool) argument. A type argument is
# necessary for type stability since the output of cholesky and cholesky! is either
# Cholesky or CholeskyPivoted. The latter is only
# supported for the four LAPACK element types. For other types, e.g. BigFloats Val(true) will
# give an error. It is required that the input is Hermitian (including real symmetric) either
# through the Hermitian and Symmetric views or exact symmetric or Hermitian elements which
# is checked for and an error is thrown if the check fails.

# The internal structure is as follows
# - _chol! returns the factor and info without checking positive definiteness
# - cholesky/cholesky! returns Cholesky without checking positive definiteness

# FixMe? The dispatch below seems overly complicated. One simplification could be to
# merge the two Cholesky types into one. It would remove the need for Val completely but
# the cost would be extra unnecessary/unused fields for the unpivoted Cholesky and runtime
# checks of those fields before calls to LAPACK to check which version of the Cholesky
# factorization the type represents.
"""
    Cholesky <: Factorization

Matrix factorization type of the Cholesky factorization of a dense symmetric/Hermitian
positive definite matrix `A`. This is the return type of [`cholesky`](@ref),
the corresponding matrix factorization function.

The triangular Cholesky factor can be obtained from the factorization `F::Cholesky`
via `F.L` and `F.U`.

# Examples
```jldoctest
julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
3×3 Array{Float64,2}:
   4.0   12.0  -16.0
  12.0   37.0  -43.0
 -16.0  -43.0   98.0

julia> C = cholesky(A)
Cholesky{Float64,Array{Float64,2}}
U factor:
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
  ⋅   1.0   5.0
  ⋅    ⋅    3.0

julia> C.U
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
  ⋅   1.0   5.0
  ⋅    ⋅    3.0

julia> C.L
3×3 LowerTriangular{Float64,Array{Float64,2}}:
  2.0   ⋅    ⋅
  6.0  1.0   ⋅
 -8.0  5.0  3.0

julia> C.L * C.U == A
true
```
"""
struct Cholesky{T,S<:AbstractMatrix} <: Factorization{T}
    factors::S
    uplo::Char
    info::BlasInt

    function Cholesky{T,S}(factors, uplo, info) where {T,S<:AbstractMatrix}
        require_one_based_indexing(factors)
        new(factors, uplo, info)
    end
end
Cholesky(A::AbstractMatrix{T}, uplo::Symbol, info::Integer) where {T} =
    Cholesky{T,typeof(A)}(A, char_uplo(uplo), info)
Cholesky(A::AbstractMatrix{T}, uplo::AbstractChar, info::Integer) where {T} =
    Cholesky{T,typeof(A)}(A, uplo, info)

"""
    CholeskyPivoted

Matrix factorization type of the pivoted Cholesky factorization of a dense symmetric/Hermitian
positive semi-definite matrix `A`. This is the return type of [`cholesky(_, Val(true))`](@ref),
the corresponding matrix factorization function.

The triangular Cholesky factor can be obtained from the factorization `F::CholeskyPivoted`
via `F.L` and `F.U`.

# Examples
```jldoctest
julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
3×3 Array{Float64,2}:
   4.0   12.0  -16.0
  12.0   37.0  -43.0
 -16.0  -43.0   98.0

julia> C = cholesky(A, Val(true))
CholeskyPivoted{Float64,Array{Float64,2}}
U factor with rank 3:
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 9.89949  -4.34366  -1.61624
  ⋅        4.25825   1.1694
  ⋅         ⋅        0.142334
permutation:
3-element Array{Int64,1}:
 3
 2
 1
```
"""
struct CholeskyPivoted{T,S<:AbstractMatrix} <: Factorization{T}
    factors::S
    uplo::Char
    piv::Vector{BlasInt}
    rank::BlasInt
    tol::Real
    info::BlasInt

    function CholeskyPivoted{T,S}(factors, uplo, piv, rank, tol, info) where {T,S<:AbstractMatrix}
        require_one_based_indexing(factors)
        new(factors, uplo, piv, rank, tol, info)
    end
end
function CholeskyPivoted(A::AbstractMatrix{T}, uplo::AbstractChar, piv::Vector{<:Integer},
                            rank::Integer, tol::Real, info::Integer) where T
    CholeskyPivoted{T,typeof(A)}(A, uplo, piv, rank, tol, info)
end

# make a copy that allow inplace Cholesky factorization
@inline choltype(A) = promote_type(typeof(sqrt(oneunit(eltype(A)))), Float32)
@inline cholcopy(A) = copy_oftype(A, choltype(A))

# _chol!. Internal methods for calling unpivoted Cholesky
## BLAS/LAPACK element types
function _chol!(A::StridedMatrix{<:BlasFloat}, ::Type{UpperTriangular})
    C, info = LAPACK.potrf!('U', A)
    return UpperTriangular(C), info
end
function _chol!(A::StridedMatrix{<:BlasFloat}, ::Type{LowerTriangular})
    C, info = LAPACK.potrf!('L', A)
    return LowerTriangular(C), info
end
function _chol!(A::StridedMatrix)
    if !ishermitian(A) # return with info = -1 if not Hermitian
        return UpperTriangular(A), convert(BlasInt, -1)
    else
        return _chol!(A, UpperTriangular)
    end
end

## Non BLAS/LAPACK element types (generic)
function _chol!(A::AbstractMatrix, ::Type{UpperTriangular})
    require_one_based_indexing(A)
    n = checksquare(A)
    @inbounds begin
        for k = 1:n
            for i = 1:k - 1
                A[k,k] -= A[i,k]'A[i,k]
            end
            Akk, info = _chol!(A[k,k], UpperTriangular)
            if info != 0
                return UpperTriangular(A), info
            end
            A[k,k] = Akk
            AkkInv = inv(copy(Akk'))
            for j = k + 1:n
                for i = 1:k - 1
                    A[k,j] -= A[i,k]'A[i,j]
                end
                A[k,j] = AkkInv*A[k,j]
            end
        end
    end
    return UpperTriangular(A), convert(BlasInt, 0)
end
function _chol!(A::AbstractMatrix, ::Type{LowerTriangular})
    require_one_based_indexing(A)
    n = checksquare(A)
    @inbounds begin
        for k = 1:n
            for i = 1:k - 1
                A[k,k] -= A[k,i]*A[k,i]'
            end
            Akk, info = _chol!(A[k,k], LowerTriangular)
            if info != 0
                return LowerTriangular(A), info
            end
            A[k,k] = Akk
            AkkInv = inv(Akk)
            for j = 1:k - 1
                @simd for i = k + 1:n
                    A[i,k] -= A[i,j]*A[k,j]'
                end
            end
            for i = k + 1:n
                A[i,k] *= AkkInv'
            end
        end
     end
    return LowerTriangular(A), convert(BlasInt, 0)
end

## Numbers
function _chol!(x::Number, uplo)
    rx = real(x)
    rxr = sqrt(abs(rx))
    rval =  convert(promote_type(typeof(x), typeof(rxr)), rxr)
    rx == abs(x) ? (rval, convert(BlasInt, 0)) : (rval, convert(BlasInt, 1))
end

## for StridedMatrices, check that matrix is symmetric/Hermitian

# cholesky!. Destructive methods for computing Cholesky factorization of real symmetric
# or Hermitian matrix
## No pivoting (default)
function cholesky!(A::RealHermSymComplexHerm, ::Val{false}=Val(false); check::Bool = true)
    C, info = _chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
    check && checkpositivedefinite(info)
    return Cholesky(C.data, A.uplo, info)
end

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
    cholesky!(A, Val(false); check = true) -> Cholesky

The same as [`cholesky`](@ref), but saves space by overwriting the input `A`,
instead of creating a copy. An [`InexactError`](@ref) exception is thrown if
the factorization produces a number not representable by the element type of
`A`, e.g. for integer types.

# Examples
```jldoctest
julia> A = [1 2; 2 50]
2×2 Array{Int64,2}:
 1   2
 2  50

julia> cholesky!(A)
ERROR: InexactError: Int64(6.782329983125268)
Stacktrace:
[...]
```
"""
function cholesky!(A::StridedMatrix, ::Val{false}=Val(false); check::Bool = true)
    checksquare(A)
    if !ishermitian(A) # return with info = -1 if not Hermitian
        check && checkpositivedefinite(-1)
        return Cholesky(A, 'U', convert(BlasInt, -1))
    else
        return cholesky!(Hermitian(A), Val(false); check = check)
    end
end

## With pivoting
### BLAS/LAPACK element types
function cholesky!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix},
                   ::Val{true}; tol = 0.0, check::Bool = true)
    AA, piv, rank, info = LAPACK.pstrf!(A.uplo, A.data, tol)
    C = CholeskyPivoted{eltype(AA),typeof(AA)}(AA, A.uplo, piv, rank, tol, info)
    check && chkfullrank(C)
    return C
end

### Non BLAS/LAPACK element types (generic). Since generic fallback for pivoted Cholesky
### is not implemented yet we throw an error
cholesky!(A::RealHermSymComplexHerm{<:Real}, ::Val{true}; tol = 0.0, check::Bool = true) =
    throw(ArgumentError("generic pivoted Cholesky factorization is not implemented yet"))

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
    cholesky!(A, Val(true); tol = 0.0, check = true) -> CholeskyPivoted

The same as [`cholesky`](@ref), but saves space by overwriting the input `A`,
instead of creating a copy. An [`InexactError`](@ref) exception is thrown if the
factorization produces a number not representable by the element type of `A`,
e.g. for integer types.
"""
function cholesky!(A::StridedMatrix, ::Val{true}; tol = 0.0, check::Bool = true)
    checksquare(A)
    if !ishermitian(A)
        C = CholeskyPivoted(A, 'U', Vector{BlasInt}(),convert(BlasInt, 1),
                            tol, convert(BlasInt, -1))
        check && chkfullrank(C)
        return C
    else
        return cholesky!(Hermitian(A), Val(true); tol = tol, check = check)
    end
end

# cholesky. Non-destructive methods for computing Cholesky factorization of real symmetric
# or Hermitian matrix
## No pivoting (default)
"""
    cholesky(A, Val(false); check = true) -> Cholesky

Compute the Cholesky factorization of a dense symmetric positive definite matrix `A`
and return a [`Cholesky`](@ref) factorization. The matrix `A` can either be a [`Symmetric`](@ref) or [`Hermitian`](@ref)
[`StridedMatrix`](@ref) or a *perfectly* symmetric or Hermitian `StridedMatrix`.
The triangular Cholesky factor can be obtained from the factorization `F` with: `F.L` and `F.U`.
The following functions are available for `Cholesky` objects: [`size`](@ref), [`\\`](@ref),
[`inv`](@ref), [`det`](@ref), [`logdet`](@ref) and [`isposdef`](@ref).

When `check = true`, an error is thrown if the decomposition fails.
When `check = false`, responsibility for checking the decomposition's
validity (via [`issuccess`](@ref)) lies with the user.

# Examples
```jldoctest
julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
3×3 Array{Float64,2}:
   4.0   12.0  -16.0
  12.0   37.0  -43.0
 -16.0  -43.0   98.0

julia> C = cholesky(A)
Cholesky{Float64,Array{Float64,2}}
U factor:
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
  ⋅   1.0   5.0
  ⋅    ⋅    3.0

julia> C.U
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
  ⋅   1.0   5.0
  ⋅    ⋅    3.0

julia> C.L
3×3 LowerTriangular{Float64,Array{Float64,2}}:
  2.0   ⋅    ⋅
  6.0  1.0   ⋅
 -8.0  5.0  3.0

julia> C.L * C.U == A
true
```
"""
cholesky(A::Union{StridedMatrix,RealHermSymComplexHerm{<:Real,<:StridedMatrix}},
    ::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)

## With pivoting
"""
    cholesky(A, Val(true); tol = 0.0, check = true) -> CholeskyPivoted

Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrix `A`
and return a [`CholeskyPivoted`](@ref) factorization. The matrix `A` can either be a [`Symmetric`](@ref)
or [`Hermitian`](@ref) [`StridedMatrix`](@ref) or a *perfectly* symmetric or Hermitian `StridedMatrix`.
The triangular Cholesky factor can be obtained from the factorization `F` with: `F.L` and `F.U`.
The following functions are available for `CholeskyPivoted` objects:
[`size`](@ref), [`\\`](@ref), [`inv`](@ref), [`det`](@ref), and [`rank`](@ref).
The argument `tol` determines the tolerance for determining the rank.
For negative values, the tolerance is the machine precision.

When `check = true`, an error is thrown if the decomposition fails.
When `check = false`, responsibility for checking the decomposition's
validity (via [`issuccess`](@ref)) lies with the user.
"""
cholesky(A::Union{StridedMatrix,RealHermSymComplexHerm{<:Real,<:StridedMatrix}},
    ::Val{true}; tol = 0.0, check::Bool = true) =
    cholesky!(cholcopy(A), Val(true); tol = tol, check = check)

## Number
function cholesky(x::Number, uplo::Symbol=:U)
    C, info = _chol!(x, uplo)
    xf = fill(C, 1, 1)
    Cholesky(xf, uplo, info)
end

function Cholesky{T}(C::Cholesky) where T
    Cnew = convert(AbstractMatrix{T}, C.factors)
    Cholesky{T, typeof(Cnew)}(Cnew, C.uplo, C.info)
end
Factorization{T}(C::Cholesky{T}) where {T} = C
Factorization{T}(C::Cholesky) where {T} = Cholesky{T}(C)
CholeskyPivoted{T}(C::CholeskyPivoted{T}) where {T} = C
CholeskyPivoted{T}(C::CholeskyPivoted) where {T} =
    CholeskyPivoted(AbstractMatrix{T}(C.factors),C.uplo,C.piv,C.rank,C.tol,C.info)
Factorization{T}(C::CholeskyPivoted{T}) where {T} = C
Factorization{T}(C::CholeskyPivoted) where {T} = CholeskyPivoted{T}(C)

AbstractMatrix(C::Cholesky) = C.uplo == 'U' ? C.U'C.U : C.L*C.L'
AbstractArray(C::Cholesky) = AbstractMatrix(C)
Matrix(C::Cholesky) = Array(AbstractArray(C))
Array(C::Cholesky) = Matrix(C)

function AbstractMatrix(F::CholeskyPivoted)
    ip = invperm(F.p)
    U = F.U[1:F.rank,ip]
    U'U
end
AbstractArray(F::CholeskyPivoted) = AbstractMatrix(F)
Matrix(F::CholeskyPivoted) = Array(AbstractArray(F))
Array(F::CholeskyPivoted) = Matrix(F)

copy(C::Cholesky) = Cholesky(copy(C.factors), C.uplo, C.info)
copy(C::CholeskyPivoted) = CholeskyPivoted(copy(C.factors), C.uplo, C.piv, C.rank, C.tol, C.info)

size(C::Union{Cholesky, CholeskyPivoted}) = size(C.factors)
size(C::Union{Cholesky, CholeskyPivoted}, d::Integer) = size(C.factors, d)

function getproperty(C::Cholesky, d::Symbol)
    Cfactors = getfield(C, :factors)
    Cuplo    = getfield(C, :uplo)
    if d === :U
        return UpperTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors'))
    elseif d === :L
        return LowerTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors'))
    elseif d === :UL
        return (Cuplo === 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors))
    else
        return getfield(C, d)
    end
end
Base.propertynames(F::Cholesky, private::Bool=false) =
    (:U, :L, :UL, (private ? fieldnames(typeof(F)) : ())...)

function getproperty(C::CholeskyPivoted{T}, d::Symbol) where T<:BlasFloat
    Cfactors = getfield(C, :factors)
    Cuplo    = getfield(C, :uplo)
    if d === :U
        return UpperTriangular(sym_uplo(Cuplo) == d ? Cfactors : copy(Cfactors'))
    elseif d === :L
        return LowerTriangular(sym_uplo(Cuplo) == d ? Cfactors : copy(Cfactors'))
    elseif d === :p
        return getfield(C, :piv)
    elseif d === :P
        n = size(C, 1)
        P = zeros(T, n, n)
        for i = 1:n
            P[getfield(C, :piv)[i], i] = one(T)
        end
        return P
    else
        return getfield(C, d)
    end
end
Base.propertynames(F::CholeskyPivoted, private::Bool=false) =
    (:U, :L, :p, :P, (private ? fieldnames(typeof(F)) : ())...)

issuccess(C::Cholesky) = C.info == 0

function show(io::IO, mime::MIME{Symbol("text/plain")}, C::Cholesky{<:Any,<:AbstractMatrix})
    if issuccess(C)
        summary(io, C); println(io)
        println(io, "$(C.uplo) factor:")
        show(io, mime, C.UL)
    else
        print(io, "Failed factorization of type $(typeof(C))")
    end
end

function show(io::IO, mime::MIME{Symbol("text/plain")}, C::CholeskyPivoted{<:Any,<:AbstractMatrix})
    summary(io, C); println(io)
    println(io, "$(C.uplo) factor with rank $(rank(C)):")
    show(io, mime, C.uplo == 'U' ? C.U : C.L)
    println(io, "\npermutation:")
    show(io, mime, C.p)
end

ldiv!(C::Cholesky{T,<:AbstractMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
    LAPACK.potrs!(C.uplo, C.factors, B)

function ldiv!(C::Cholesky{<:Any,<:AbstractMatrix}, B::StridedVecOrMat)
    if C.uplo == 'L'
        return ldiv!(adjoint(LowerTriangular(C.factors)), ldiv!(LowerTriangular(C.factors), B))
    else
        return ldiv!(UpperTriangular(C.factors), ldiv!(adjoint(UpperTriangular(C.factors)), B))
    end
end

function ldiv!(C::CholeskyPivoted{T}, B::StridedVector{T}) where T<:BlasFloat
    invpermute!(LAPACK.potrs!(C.uplo, C.factors, permute!(B, C.piv)), C.piv)
end
function ldiv!(C::CholeskyPivoted{T}, B::StridedMatrix{T}) where T<:BlasFloat
    n = size(C, 1)
    for i=1:size(B, 2)
        permute!(view(B, 1:n, i), C.piv)
    end
    LAPACK.potrs!(C.uplo, C.factors, B)
    for i=1:size(B, 2)
        invpermute!(view(B, 1:n, i), C.piv)
    end
    B
end

function ldiv!(C::CholeskyPivoted, B::StridedVector)
    if C.uplo == 'L'
        ldiv!(adjoint(LowerTriangular(C.factors)),
            ldiv!(LowerTriangular(C.factors), permute!(B, C.piv)))
    else
        ldiv!(UpperTriangular(C.factors),
            ldiv!(adjoint(UpperTriangular(C.factors)), permute!(B, C.piv)))
    end
    invpermute!(B, C.piv)
end

function ldiv!(C::CholeskyPivoted, B::StridedMatrix)
    n = size(C, 1)
    for i in 1:size(B, 2)
        permute!(view(B, 1:n, i), C.piv)
    end
    if C.uplo == 'L'
        ldiv!(adjoint(LowerTriangular(C.factors)),
            ldiv!(LowerTriangular(C.factors), B))
    else
        ldiv!(UpperTriangular(C.factors),
            ldiv!(adjoint(UpperTriangular(C.factors)), B))
    end
    for i in 1:size(B, 2)
        invpermute!(view(B, 1:n, i), C.piv)
    end
    B
end

function rdiv!(B::StridedMatrix, C::Cholesky{<:Any,<:AbstractMatrix})
    if C.uplo == 'L'
        return rdiv!(rdiv!(B, adjoint(LowerTriangular(C.factors))), LowerTriangular(C.factors))
    else
        return rdiv!(rdiv!(B, UpperTriangular(C.factors)), adjoint(UpperTriangular(C.factors)))
    end
end

function LinearAlgebra.rdiv!(B::StridedMatrix, C::CholeskyPivoted)
    n = size(C, 2)
    for i in 1:size(B, 1)
        permute!(view(B, i, 1:n), C.piv)
    end
    if C.uplo == 'L'
        rdiv!(rdiv!(B, adjoint(LowerTriangular(C.factors))),
            LowerTriangular(C.factors))
    else
        rdiv!(rdiv!(B, UpperTriangular(C.factors)),
            adjoint(UpperTriangular(C.factors)))
    end
    for i in 1:size(B, 1)
        invpermute!(view(B, i, 1:n), C.piv)
    end
    B
end

isposdef(C::Union{Cholesky,CholeskyPivoted}) = C.info == 0

function det(C::Cholesky)
    dd = one(real(eltype(C)))
    @inbounds for i in 1:size(C.factors,1)
        dd *= real(C.factors[i,i])^2
    end
    return dd
end

function logdet(C::Cholesky)
    dd = zero(real(eltype(C)))
    @inbounds for i in 1:size(C.factors,1)
        dd += log(real(C.factors[i,i]))
    end
    dd + dd # instead of 2.0dd which can change the type
end

function det(C::CholeskyPivoted)
    if C.rank < size(C.factors, 1)
        return zero(real(eltype(C)))
    else
        dd = one(real(eltype(C)))
        for i in 1:size(C.factors,1)
            dd *= real(C.factors[i,i])^2
        end
        return dd
    end
end

function logdet(C::CholeskyPivoted)
    if C.rank < size(C.factors, 1)
        return real(eltype(C))(-Inf)
    else
        dd = zero(real(eltype(C)))
        for i in 1:size(C.factors,1)
            dd += log(real(C.factors[i,i]))
        end
        return dd + dd # instead of 2.0dd which can change the type
    end
end

inv!(C::Cholesky{<:BlasFloat,<:StridedMatrix}) =
    copytri!(LAPACK.potri!(C.uplo, C.factors), C.uplo, true)

inv(C::Cholesky{<:BlasFloat,<:StridedMatrix}) = inv!(copy(C))

function inv(C::CholeskyPivoted)
    ipiv = invperm(C.piv)
    copytri!(LAPACK.potri!(C.uplo, copy(C.factors)), C.uplo, true)[ipiv, ipiv]
end

function chkfullrank(C::CholeskyPivoted)
    if C.rank < size(C.factors, 1)
        throw(RankDeficientException(C.info))
    end
end

rank(C::CholeskyPivoted) = C.rank

"""
    lowrankupdate!(C::Cholesky, v::StridedVector) -> CC::Cholesky

Update a Cholesky factorization `C` with the vector `v`. If `A = C.U'C.U` then
`CC = cholesky(C.U'C.U + v*v')` but the computation of `CC` only uses `O(n^2)`
operations. The input factorization `C` is updated in place such that on exit `C == CC`.
The vector `v` is destroyed during the computation.
"""
function lowrankupdate!(C::Cholesky, v::StridedVector)
    A = C.factors
    n = length(v)
    if size(C, 1) != n
        throw(DimensionMismatch("updating vector must fit size of factorization"))
    end
    if C.uplo == 'U'
        conj!(v)
    end

    for i = 1:n

        # Compute Givens rotation
        c, s, r = givensAlgorithm(A[i,i], v[i])

        # Store new diagonal element
        A[i,i] = r

        # Update remaining elements in row/column
        if C.uplo == 'U'
            for j = i + 1:n
                Aij = A[i,j]
                vj  = v[j]
                A[i,j]  =   c*Aij + s*vj
                v[j]    = -s'*Aij + c*vj
            end
        else
            for j = i + 1:n
                Aji = A[j,i]
                vj  = v[j]
                A[j,i]  =   c*Aji + s*vj
                v[j]    = -s'*Aji + c*vj
            end
        end
    end
    return C
end

"""
    lowrankdowndate!(C::Cholesky, v::StridedVector) -> CC::Cholesky

Downdate a Cholesky factorization `C` with the vector `v`. If `A = C.U'C.U` then
`CC = cholesky(C.U'C.U - v*v')` but the computation of `CC` only uses `O(n^2)`
operations. The input factorization `C` is updated in place such that on exit `C == CC`.
The vector `v` is destroyed during the computation.
"""
function lowrankdowndate!(C::Cholesky, v::StridedVector)
    A = C.factors
    n = length(v)
    if size(C, 1) != n
        throw(DimensionMismatch("updating vector must fit size of factorization"))
    end
    if C.uplo == 'U'
        conj!(v)
    end

    for i = 1:n

        Aii = A[i,i]

        # Compute Givens rotation
        s = conj(v[i]/Aii)
        s2 = abs2(s)
        if s2 > 1
            throw(LinearAlgebra.PosDefException(i))
        end
        c = sqrt(1 - abs2(s))

        # Store new diagonal element
        A[i,i] = c*Aii

        # Update remaining elements in row/column
        if C.uplo == 'U'
            for j = i + 1:n
                vj = v[j]
                Aij = (A[i,j] - s*vj)/c
                A[i,j] = Aij
                v[j] = -s'*Aij + c*vj
            end
        else
            for j = i + 1:n
                vj = v[j]
                Aji = (A[j,i] - s*vj)/c
                A[j,i] = Aji
                v[j] = -s'*Aji + c*vj
            end
        end
    end
    return C
end

"""
    lowrankupdate(C::Cholesky, v::StridedVector) -> CC::Cholesky

Update a Cholesky factorization `C` with the vector `v`. If `A = C.U'C.U`
then `CC = cholesky(C.U'C.U + v*v')` but the computation of `CC` only uses
`O(n^2)` operations.
"""
lowrankupdate(C::Cholesky, v::StridedVector) = lowrankupdate!(copy(C), copy(v))

"""
    lowrankdowndate(C::Cholesky, v::StridedVector) -> CC::Cholesky

Downdate a Cholesky factorization `C` with the vector `v`. If `A = C.U'C.U`
then `CC = cholesky(C.U'C.U - v*v')` but the computation of `CC` only uses
`O(n^2)` operations.
"""
lowrankdowndate(C::Cholesky, v::StridedVector) = lowrankdowndate!(copy(C), copy(v))

max_stars_repo_name: naveenarun/julia

max_issues_repo_name: naveenarun/julia

max_forks_repo_name: naveenarun/julia


text:
function optimize(je_pg_model::JointEstimationPolyaGammaModel)
    local parameters = je_pg_model.parameters
    local latents = je_pg_model.latents
    local n_index = je_pg_model.n_index
    local i_index = je_pg_model.i_index
    local responses_per_item = je_pg_model.responses_per_item
    local responses_per_examinee = je_pg_model.responses_per_examinee

    I = size(parameters, 1)
    N = size(latents, 1)

    batch_den = max(1,(nprocs()-1))
    batch_size_N = Int(ceil(N/batch_den))
    batch_size_I = Int(ceil(I/batch_den))
    tuple_i_n = vcat([[[i, n] for n in n_index[i]] for i in 1:I]...)

    stop = false
    start_time = time()
    W = SharedMatrix{Float64}(I,N)
    iter = 1
    while !stop
        for n = 1:N
            local latent_n = latents[n]
            for i in i_index[n]
                W[i, n] = _generate_w(latent_n, parameters[i])
            end
        end
        #W_per_item = [W[i, n_index[i]] for i in 1:I]

        #W_per_examinee = [W[i_index[n], n] for n in 1:N]

        # parameters = pmap((p, idx, r, w) -> _mcmc_iter_pg(p, latents[idx], r, w[idx]; sampling = je_pg_model.item_sampling),
        #     parameters,
        #     n_index,
        #     je_pg_model.responses_per_item,
        #     eachrow(W),
        #     batch_size = batch_size_I,
        # )
        # latents = pmap((l, idx, r, w) ->
        #             _mcmc_iter_pg(l, parameters[idx], r, w[idx]; sampling = je_pg_model.examinee_sampling),
        #     latents,
        #     i_index,
        #     je_pg_model.responses_per_examinee,
        #     eachcol(W),
        #     batch_size = batch_size_N,
        # )
            for i in 1:I
                if !parameters[i].calibrated
                    parameters[i] = _mcmc_iter_pg(parameters[i], latents[n_index[i]], responses_per_item[i], W[i, n_index[i]]; sampling = je_pg_model.item_sampling)
               end
            end
            for n in 1:N
                if !latents[n].assessed
                   latents[n] = _mcmc_iter_pg(latents[n], parameters[i_index[n]], responses_per_examinee[n], W[i_index[n], n]; sampling = je_pg_model.examinee_sampling)
                end
            end
            if (iter % 200) == 0
                if any([
                    check_iter(iter; max_iter = Int(ext_opt_settings[1]), verbosity = Int(je_pg_model.ext_opt_settings[5])),
                    check_time(start_time; max_time = Int(ext_opt_settings[2]), verbosity = Int(je_pg_model.ext_opt_settings[5])),
                    # check_x_tol_rel!(
                    #     items,
                    #     old_pars;
                    #     x_tol_rel = x_tol_rel
                    # )
                    ])
                    stop = true
                end
            end
        iter += 1
    end
    return parameters::Vector{<:AbstractParameters}, latents::Vector{<:AbstractLatent}
end


text:
clear
load test21.mat
maxQr=max(max(qrain_2400))
maxQc=max(max(qc_2400))
total_qrain=sum(sum(qrain_2400>0.00001));
qrain_area = total_qrain*.25^2
total_qc=sum(sum(qc_2400>0.0001));
qc_area = total_qc*.25^2

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
@connector function Flange(;name)
    sts = @variables begin
        phi(t)
        tau(t), [connect=Flow]
    end
    ODESystem(Equation[], t, sts, [], name=name, defaults=Dict(phi=>0.0, tau=>0.0))
end
Base.@doc """
    Flange(;name)

1-dim. rotational flange of a shaft.

# States:
- `phi`: [rad] Absolute rotation angle of flange
- `tau`: [N.m] Cut torque in the flange
""" Flange

@connector function Support(;name)
    sts = @variables begin
        phi(t)
        tau(t), [connect=Flow]
    end
    ODESystem(Equation[], t, sts, [], name=name, defaults=Dict(phi=>0.0, tau=>0.0))
end
Base.@doc """
    Support(;name)

Support/housing of a 1-dim. rotational shaft

# States:
- `phi`: [rad] Absolute rotation angle of the support/housing
- `tau`: [N.m] Cut torque in the support/housing
""" Support

"""
    PartialCompliant(;name, phi_rel_start=0.0, tau_start=0.0)

Partial model for the compliant connection of two rotational 1-dim. shaft flanges.

# Parameters:
- `phi_rel_start`: [rad] Initial relative rotation angle
- `tau_start`: [N.m] Initial torque between flanges

# States:
- `phi_rel`: [rad] Relative rotation angle (= flange_b.phi - flange_a.phi)
- `tau`: [N.m] Torque between flanges (= flange_b.tau)
"""
function PartialCompliant(;name, phi_rel_start=0.0, tau_start=0.0)
    @named flange_a = Flange()
    @named flange_b = Flange()
    sts = @variables begin
        phi_rel(t)=phi_rel_start
        tau(t)=tau_start
    end
    eqs = [
        phi_rel ~ flange_b.phi - flange_a.phi
        flange_b.tau ~ tau
        flange_a.tau ~ -tau
    ]
    return compose(ODESystem(eqs, t, sts, []; name=name), flange_a, flange_b)
end

"""
    PartialCompliantWithRelativeStates(;name, phi_rel_start=0.0, tau_start=0.0)

Partial model for the compliant connection of two rotational 1-dim. shaft flanges where the relative angle and speed are used as preferred states

# Parameters:
- `phi_rel_start`: [rad] Initial relative rotation angle
- `w_rel_start`: [rad/s] Initial relative angular velocity (= der(phi_rel))
- `a_rel_start`: [rad/s²] Initial relative angular acceleration (= der(w_rel))
- `tau_start`: [N.m] Initial torque between flanges

# States:
- `phi_rel`: [rad] Relative rotation angle (= flange_b.phi - flange_a.phi)
- `w_rel`: [rad/s] Relative angular velocity (= der(phi_rel))
- `a_rel`: [rad/s²] Relative angular acceleration (= der(w_rel))
- `tau`: [N.m] Torque between flanges (= flange_b.tau)
"""
function PartialCompliantWithRelativeStates(;name, phi_rel_start=0.0, w_start=0.0, a_start=0.0, tau_start=0.0)
    @named flange_a = Flange()
    @named flange_b = Flange()
    sts = @variables begin
        phi_rel(t)=phi_rel_start
        w_rel(t)=w_start
        a_rel(t)=a_start
        tau(t)=tau_start
    end
    eqs = [
        phi_rel ~ flange_b.phi - flange_a.phi
        D(phi_rel) ~ w_rel
        D(w_rel) ~ a_rel
        flange_b.tau ~ tau
        flange_a.tau ~ -tau
    ]
    return compose(ODESystem(eqs, t, sts, []; name=name), flange_a, flange_b)
end

"""
    PartialElementaryOneFlangeAndSupport2(;name, use_support=false)

Partial model for a component with one rotational 1-dim. shaft flange and a support used for textual modeling, i.e., for elementary models

# Parameters:
- `use_support`: If support flange enabled, otherwise implicitly grounded

# States:
- `phi_support`: [rad] Absolute angle of support flange"
"""
function PartialElementaryOneFlangeAndSupport2(;name, use_support=false)
    @named flange = Flange()
    sys = [flange]
    @variables phi_support(t)
    if use_support
        @named support = Support()
        eqs = [
            support.phi ~ phi_support
            support.tau ~ -flange.tau
        ]
        push!(sys, support)
    else
        eqs = [phi_support ~ 0]
    end
    return compose(ODESystem(eqs, t, [phi_support], []; name=name), sys)
end

"""
    PartialElementaryTwoFlangesAndSupport2(;name, use_support=false)

Partial model for a component with two rotational 1-dim. shaft flanges and a support used for textual modeling, i.e., for elementary models

# Parameters:
- `use_support`: If support flange enabled, otherwise implicitly grounded

# States:
- `phi_support`: [rad] Absolute angle of support flange"
"""
function PartialElementaryTwoFlangesAndSupport2(;name, use_support=false)
    @named flange_a = Flange()
    @named flange_b = Flange()
    sys = [flange_a, flange_b]
    @variables phi_support(t)=0.0
    if use_support
        @named support = Support()
        eqs = [
            support.phi ~ phi_support
            support.tau ~ -flange_a.tau - flange_b.tau
        ]
        push!(sys, support)
    else
        eqs = [phi_support ~ 0]
    end
    return compose(ODESystem(eqs, t, [phi_support], []; name=name), sys)
end


text:
module SearchDiscoveryReplication
__precompile__()

# Exported functions and types
export Model,estimOptions  # Structs
export Weitzman, FullInfo, RandomSearch, SearchDiscovery # Models
export calcZd, calcZs # Reservation value calculation
export genData, estimateModel, calcWelfareDemand # Main functions for each model

export runSimulation, calcWelfareChanges, gatherResults  # Functions specific for replication
export writeLatexTable # To write to latex table

# Support pkgs
using Parameters: @with_kw
using ArgCheck: @argcheck
using Printf: @sprintf
# Math / Stats pkgs
using Distributions, QuadGK, Roots, ForwardDiff , StatsBase, Random, LinearAlgebra
# Optimizers
using NLopt, GalacticOptim
# Data org pkgs
using BSON: @save, @load

include("definitions_struct.jl") # Own types
include("bivariate.jl") # Functions for Bivariate Normal distribution
include("reservation_values.jl") # Functions for reservation value calculation
include("support.jl")
# Different models
include("models/weitzman.jl")
include("models/randomsearch.jl")
include("models/fullinfo.jl")
include("models/discovery_incomplete.jl")
# Specific for replication
include("replication.jl")

end