text:
function CR3BPDifferentialCorrection(arg0::PVCoordinates, arg1::CR3BPSystem, arg2::jdouble)
    return CR3BPDifferentialCorrection((PVCoordinates, CR3BPSystem, jdouble), arg0, arg1, arg2)
end

function CR3BPDifferentialCorrection(arg0::PVCoordinates, arg1::CR3BPSystem, arg2::jdouble, arg3::AttitudeProvider, arg4::TimeScale)
    return CR3BPDifferentialCorrection((PVCoordinates, CR3BPSystem, jdouble, AttitudeProvider, TimeScale), arg0, arg1, arg2, arg3, arg4)
end

function compute(obj::CR3BPDifferentialCorrection, arg0::LibrationOrbitType)
    return jcall(obj, "compute", PVCoordinates, (LibrationOrbitType,), arg0)
end

function compute_lyapunov(obj::CR3BPDifferentialCorrection, arg0::STMEquations)
    return jcall(obj, "computeLyapunov", PVCoordinates, (STMEquations,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_orbital_period(obj::CR3BPDifferentialCorrection)
    return jcall(obj, "getOrbitalPeriod", jdouble, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
# Compute tracks from entity edge data.

using JLD, PyPlot
include("findtracks.jl")
include("findtrackgraph.jl")

# Load edge incidence matrix.
# Load the data file
file_dir = joinpath(Base.source_dir(),"../1EntityAnalysis/Entity.jld")
E = load(file_dir)["E"]
#E = loadassoc(file_dir)
E = logical(E)

# Set prefixes
p = StartsWith("PERSON/,")
t = StartsWith("TIME/,")
x = StartsWith("LOCATION/,")

# Build entity tracks with routine.
A = findtracks(E,t,p,x)

# Track Graph
G = findtrackgraph(A)
print(G > 5)
figure()
spy(G)

# Track graph pattern
o = "ORGANIZATION/international monetary fund,"
p = StartsWith("PERSON/,")
Go = findtrackgraph(A[:,getcol(E[getrow(E[:,o]),p])]);

print((Go > 2) & ((Go ./ G) > 0.2))

ext: jl


text:
using ChiralGraphs
using Test

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


text:
const ybard = [1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1,
               3.5e-1, 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1,
               1.34e0, 2.10e0, 4.39e0]

const bard = let res_init=zeros(15), jac_init=zeros(15,3), x_init=[1.0, 1.0, 1.0]

    function res(x, r)
        for i = 1:15
            u = Float64(i)
            v = 1.6e1 - u
            w = min(u, v)
            r[i] = ybard[i] - (x[1] + u / (x[2] * v + x[3] * w))
        end
        return r
    end

    function jac(x, j)
        for i = 1:15
            j[i, 1] = -1.e0
            u = Float64(i)
            v = 1.6e1 - u
            w = min(u, v)
            t = u/(x[2]*v + x[3]*w)^2
            j[i, 2] = v*t
            j[i, 3] = w*t
        end
        return j
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("bard", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end


text:
using Genie.Router
using Genie.Requests
using Genie.Renderer.Json
using MainPageController
using DataFrames

const DROPDOWNITEMS = Dict("1" => "Entry 1", "2" => "Entry 2", "3" => "Entry 3")
const TABLEDATA = DataFrame(ID=[1, 2, 3, 4], Name=["AA", "BB", "CC", "DD"],
                Country=["DE", "GB", "DK", "US"])

"""
    dict_to_dropdown(data)

Transform the `Dict` `data` to a dictionary with the correct data format for
the dropdowns.
"""
function dict_to_dropdown(dict::Dict)
    out = []
    for key in keys(dict)
        push!(out, Dict("id" => key,
            "value" => dict[key]))
    end
    return out
end

"""
    df_to_table(data, columns, id)

Prepare a dictionary for the web frontend to show tabular data from the data source `data`,
with column names from the vector `columns`.
"""
function df_to_table(df::DataFrame, columns, id)
    out = []
    for i = 1:size(df, 1)
        dict = Dict("key" => df[i,1],
        "ID" => string(df[i,1]))
        for (j, name) in enumerate(columns)
            dict[name] = df[i,j + 1] * " " * string(id)
        end
        push!(out, dict)
    end
    return out
end

"""
    to_table_columns(colnames)

Convert the vector of column titles `colnames` to the dictionary used by the tables in the
web frontend.
"""
function to_table_columns(colnames)
    out = []
    for name in colnames
        push!(out,Dict("title" => name,
                    "dataIndex" => name,
                    "key" => name))
    end
    return out
end

route("/", MainPageController.mainpage)

route("/upload/fileSubmit", method=POST) do
    payload = filespayload()
    filename = payload["Upload,file"].name
    @info filename
    write(payload["Upload,file"]; filename="./assets/uploads/" * filename)
end

route("/api/getData") do
    data = Dict("DDData" => dict_to_dropdown(DROPDOWNITEMS))
    return json(data)
end

route("/api/formSubmit", method=POST) do
    data = jsonpayload()
    id = parse(Int, data["DataRequest"])
    columns = ["Name", "Country"]
    Data = Dict("dataSource" => df_to_table(TABLEDATA, columns, id),
            "columns" => to_table_columns(columns))
    return json(Data)
end


text:
"""
    find_blocked_reactions(model, [tolerance = 1e-10])

Return a list of reactions that exhibit zero flux
"""
function find_blocked_reactions(model::Model; tolerance::Number = 1e-10)
	return find_blocked_reactions(model, tolerance)
end

function find_blocked_reactions(model::Model, tolerance::Number = 1e-10)

	minF,maxF = fva(model,0);
	index = Int64[];

	# Find all min/max fluxes that approach 0
	for i = 1:length(minF)
		if (abs(minF[i]) < tolerance) && (abs(maxF[i]) < tolerance)
			index = [index; model.rxns[i]];
		end
	end

	return index;
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    find_deadend_metabolites(model, lbfilter = true)

Locate all metabolites that are only either produced or consumed but
not both

### Optional arguements:
    lbfilter: true / false
Filter the result for metabolites who participate in reactions that have `positive` lower bounds
"""
function find_deadend_metabolites(model::Model; lbfilter::Bool = true)
	return find_deadend_metabolites(model, lbfilter)
end

function find_deadend_metabolites(model::Model, lbfilter::Bool = true)
	dead_end = zeros(model.S.m)
	reactions = get_sparse_columns(model.S)

	# determine pos/neg numbers
	negative = (x -> model.S.nzval[x] < 0.0)
	positive = (x -> model.S.nzval[x] > 0.0)

	# reactions with a POSITIVE lower bound
	lb_filter = (x -> model.lb[reactions[x]] >= 0.0)

	for i = 1:length(model.b)

		# Need to find for each met, the reactions it participates in
		# Then filter out every met thats either always a producer or always
		# a consumer
		# Optionally I filter out those that have a positive lower bound >= 0.0
		met_s_location = findin(model.S.rowval, i)

		if lbfilter
			met_s_location = filter(lb_filter, met_s_location)
		end

		if length(met_s_location) < 2
			continue
		end

		# arrays of neg/pos occurrences
		consume = map(negative, met_s_location)
		produce = map(positive, met_s_location)

		# any(a) returns true if any element of a is true
		if any(consume)
			dead_end[i] -= 1
		end
		if any(produce)
			dead_end[i] += 1
		end
	end

	return find(dead_end)
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    find_essential_genes(model, min_growth)

Find those genes which, if knocked out, would result in a biomass-production below that
specified by **min_growth**

### Examples
    julia> find_essential_genes(model,0.01)
    2-element Array{String,1}:
     "b2416"
     "b2415"

Note that **min_growth** represents the actual growth, not a percentage. For min growth of
0.1%, do

    julia> find_essential_genes(model,0.001 * fba(model).obj)
    2-element Array{String,1}:
     "b2416"
     "b2415"
"""
function find_essential_genes(model::Model, min_growth::Number)
    # find all single gene deletions
    ess_genes = gene_deletion(model,1).g_f

    # find those single gene deletions that bring growth below "min_growth"
    for (key,val) in ess_genes
        if val > min_growth
            delete!(ess_genes,key)
        end
    end
    ess_genes = collect(keys(ess_genes))
    ess_genes = vcat(ess_genes...)
    return ess_genes
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    find_exchange_reactions(model[, output = 'index')

Find all exhange reactions in the model

### Optional arguements:
    output:
* 'index': return an array of indices of the exchange reactions
* 'name' : return an array of the names of the exchange reactions

"""
function find_exchange_reactions(model::Model; output::String = "index")
	return find_exchange_reactions(model, output)
end

function find_exchange_reactions(model::Model, output::String = "index")
	# Need to find all reactions that only appear ONCE in the reaction indices
	# as these are those that only map to A SINGLE metabolite
	count_single_occurences = (x -> length(nzrange(model.S,x)) == 1)

	index = map(count_single_occurences, collect(1:length(model.c)))
	index = find(index)
	if contains(lowercase(output),"name")
		return model.rxns[index]
	end

	return index
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    find_reactions_from_genes(model, gene [, output = 'index'])

Find the reactions in which a specific ***gene*** participates

### Optional arguements:
    output:
* 'index': return an array of indices of the reactions
* 'name' : return an array of the names of the reactions

"""
function find_reactions_from_gene(model::Model, gene::String; output::String = "index")
	return find_reactions_from_gene(model, gene, output)
end

function find_reactions_from_gene(model::Model, gene::Int64; output::String = "index")
	return find_reactions_from_gene(model, gene, output)
end

function find_reactions_from_gene(model::Model, gene::String, output::String = "index")
	index = find_gene_effects(model, gene)

	if in(lowercase(output),["n", "name"])
		return model.rxns[index]
	end

	return index;
end

function find_reactions_from_gene(model::Model, gene::Int64, output::String = "index")
	gene = model.genes[gene]
	return find_reactions_from_gene(model, gene, output)
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    find_reactions_from_metabolites(model, metabolite [, output = 'index'])

Find the reactions in which a specific ***metabolite*** participates

### Optional arguements:
    output:
* 'index': return an array of indices of the reactions
* 'name' : return an array of the names of the reactions

"""
function find_reactions_from_metabolite(model::Model, metabolite::String; output::String = "index")
	return find_reactions_from_metabolite(model, metabolite, output)
end

function find_reactions_from_metabolite(model::Model, metabolite::String, output::String = "index")
	index = []

	if metabolite_exists(model,metabolite)
		push!(index,findin(model.S.rowval, find_metabolite(model,metabolite)))
	end

	index = vcat(index...)

	reactions = get_sparse_columns(model.S)

	index = reactions[index]

	if lowercase(output) == "formula" || lowercase(output[1]) == 'f'
		for i in index
			print_reaction_formula(model,model.rxns[i])
		end

	elseif lowercase(output) == "name" || lowercase(output[1]) == 'n'
		return model.rxns[index]

	else
		return index
	end
end

function find_reactions_from_metabolite(model::Model, metabolite::Int64; output::String = "index")
	return find_reactions_from_metabolite(model, metabolite, output)
end

function find_reactions_from_metabolite(model::Model, metabolite::Int64, output::String = "index")
	metabolite = model.mets[metabolite]
	return find_reactions_from_metabolite(model, metabolite, output)
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    gene_deletion(model, num_genes)

Perform fba for all gene combinations knock-outs from 1:n

### Examples
    knockouts = gene_deletion(model,2)

Return a GeneDeletion object which has the following fields:
*r_g - Reaction/Gene which displays the reactions disabled by gene combinations
*r_f - Reaction/Flow which displays the flow after disabling reactions
*g_f - Gene/Flow which displays the flow for each gene combination knockout
"""
function gene_deletion(model::Model, n::Number)
    gene_list, rule_list, saved_names, _ = abstract_gene_names(model)

    lp = setup_lp(model)
    wt = fba(model).obj

    # create a dictionary containing the reactions that are affected
    # by each geme, so that geneX => [rxnX1,...rxnXn]
    rxn_affected_by_gene = find_gene_effects(rule_list, gene_list)
    rxn_affected_by_gene = Dict(zip(gene_list,rxn_affected_by_gene))

    # empty dicts, yet to be filled
    r_g_dict = Dict()
    r_f_dict = Dict()
    truth_table = Dict()

    for i in 1:n

        # every gene combination
        combos = combinations(gene_list,i)
        num_combos = length(combos)

        # 1: Finding knockout effects
        # every gene combination affects some combination
        # of reactions, rxns_affected contains the indices of
        # these reactions
        rxns_affected = Array(Array{Any, 1}, num_combos)
        for (index,combo) in enumerate(combos)
            rxns = map(x -> rxn_affected_by_gene[x], combo)
            rxns = sort(unique(vcat(rxns...)))

            rxns_affected[index] = rxns
        end

        # 28% of looptime
        # 2: Convert rules to true/false statements
        # convert rules to "safe mode", where every gene is
        # simply named gene_xxx, for string matching safety
        converted_rules = Array(Array{Any,1}, num_combos)
        for (index,combo) in enumerate(combos)
            rules = rule_list[rxns_affected[index]]
            convert_rule_to_boolean(combo, rules)

            # TEST speedup big rules
            rules = map(x -> length(x) > 1000 ? compressor(x) : x, rules)

            converted_rules[index] = rules
        end
        conv_rule_list = deepcopy(converted_rules)
        converted_rules = vcat(converted_rules...)

        # 19% of looptime
        # 3: Evaluating every true/false statement
        # les út allar convertaðar reglur og eval(parse)
        new_truth_table = Dict(zip(converted_rules, fill(true, length(converted_rules))))
        rules_to_test = setdiff(collect(keys(new_truth_table)), collect(keys(truth_table)))
        truth_table = merge(new_truth_table, truth_table)

        for key in rules_to_test
            try
            truth_table[key] = eval(parse(key))
            catch
                println(key)
            end
        end

        # 31% of looptime
        # 4: Finding disabled reactions and connecting to gene knockouts
        # map reactions to gene combinations and insert into a dictionary
        dead_reactions = Array(Array{Any,1}, num_combos)
        for (index, combo) in enumerate(combos)
            rules = conv_rule_list[index]

            dead = find(map(x -> !truth_table[x], rules))
            dead_reactions[index] = rxns_affected[index][dead]
        end

        # 5: not sure
        # create a dictionary that comtains the reactions that
        # get disabled for every particular gene combination
        unique_rxns = unique(dead_reactions)
        combo_arrays = Array(Array{Any, 1}, length(unique_rxns))
        (y -> map(x -> y[x] = [], 1:length(y)))(combo_arrays)
        rxn_combo_dict = Dict(zip(unique_rxns, combo_arrays))
        for (index,combo) in enumerate(combos)
            push!(rxn_combo_dict[dead_reactions[index]], combo)
        end

        # 10% of total
        # 6: all the solver stuff
        # calculate fba for every combination of disabled reactions
        if length(procs()) != 1
            println("Multiple cores active. Using parallel...")
            rxn_flow_dict = PCBM.pgd(model, unique_rxns)
            rxn_flow_dict = merge(rxn_flow_dict...)
        else
            rxn_flow_dict = analyze_list_of_blocked_reaction_lp(lp, unique_rxns)
        end
        # 7: Unioning dictionaries
        # find the union of the rxn_combo_dict and former loops
        intersects = intersect(collect(keys(rxn_combo_dict)), collect(keys(r_g_dict)))
        for clash in intersects
            r_g_dict[clash] = vcat(r_g_dict[clash], rxn_combo_dict[clash])
        end
        r_g_dict = merge(rxn_combo_dict, r_g_dict)

        # 8: unioning more dictionaries
        # find the union of rxns_flow from other loops and this loop
        intersects = intersect(collect(keys(rxn_flow_dict)), collect(r_f_dict))
        for clash in intersects
            r_f_dict[clash] = vcat(r_f_dict[clash], rxn_flow_dict[clash])
        end
        r_f_dict = merge(rxn_flow_dict, r_f_dict)

    end

    for (key,val) in r_g_dict, (i,v) in enumerate(val), (j,k) in enumerate(v)
        v[j] = saved_names[k]
    end

    c_f_dict = Dict()
    for (key,val) in r_g_dict, gene in val
        c_f_dict[convert(Array{UTF8String},gene)] = r_f_dict[key]
    end

    return GeneDeletion(r_g_dict, r_f_dict, c_f_dict)
end

function compressor(rule)
    rule = replace(rule, "( true )", "true")
    rule = replace(rule, "( false )", "false")

    b = length(rule)

    while b != 0
        rule = replace(rule, "true && true", "true")
        a = length(rule)
        b = a == b ? 0 : a
    end

    b = length(rule)
    while b != 0
        rule = replace(rule, "true || true || true", "true || true")
        a = length(rule)
        b = a == b ? 0 : a
    end

    b = length(rule)
    while b != 0
        rule = replace(rule, "false && true", "false")
        rule = replace(rule, "true && false", "false")
        a = length(rule)
        b = a == b ? 0 : a
    end

    b = length(rule)
    while b != 0
        rule = replace(rule, "|| false || false ||", "|| false ||")
        rule = replace(rule, "true && false", "false")
        a = length(rule)
        b = a == b ? 0 : a
    end

    return rule
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    knockout_genes(model, genes)

Perform fba after a knocking out selected genes genes

Return a GeneKnockout type with fields ``growth``, ``flux``, ``disabled`` and ``affected``

### Example

    julia> knockout_genes(model, [model.genes[4], model.genes[16], model.genes[99]])

    Type: GeneKnockout
                            Fields
                            growth :                      0.374230
                              flux :              Array{Float64,1}
                          disabled :                       [12,71]
                          affected ::
                                                   gene:      reactions:
                                              gene_016 :            [12]
                                              gene_004 :             [3]
                                              gene_099 :            [71]

So knocking out genes 4, 16 and 99, would **touch** reactions 12, 3 and 71, respectively, but only reactions 12 and 71 would be disabled

"""
function knockout_genes(model::Model, gene::String)
    return knockout_genes(model, [gene])
end

function knockout_genes(model::Model, gene_index::Number)
    if (gene_index <= 0) | (gene_index > length(model.genes))
        warn("Gene index is invalid")
    end

    if !isinteger(gene_index)
        warn("Gene index must be integer")
    end

    return knockout_genes(model, [gene_index])
end

function knockout_genes{T <: Number}(model::Model, gene_indices::Array{T})
    gene_indices = unique(gene_indices)
    if !isempty(find(x -> (x <= 0) | (x > length(model.genes)), gene_indices))
        warn("One or more indices are invalid")
        return
    end

    return knockout_genes(model, model.genes[gene_indices])
end

function knockout_genes{T <: String}(model::Model, genes::Array{T})
    if any(x -> !in(x,model.genes), genes)
        warn("One or more gene does not appear in the model")
        return
    end

    # find wild type growth
    wt = fba(model).obj

    lp = setup_lp(model)

    # create a safe collection of rules
    gene_list, rule_list, saved_names, _ = abstract_gene_names(model)
    inverse_saved_names = invert_dict(saved_names)

    safe_genes = deepcopy(genes)
    convert_from_abstract(inverse_saved_names, safe_genes)
    rxns_affected = find_gene_effects(rule_list, safe_genes)
    gene_effect = Dict(zip(map(x -> saved_names[x], safe_genes), rxns_affected))
    rxns_affected = sort(unique(vcat(rxns_affected...)))

    rules_to_check = rule_list[rxns_affected]
    genes_to_kill = map(x -> inverse_saved_names[x], genes)
    convert_rule_to_boolean(genes_to_kill, rules_to_check)

    rxns_to_kill = []
    for (index, value) in enumerate(rules_to_check)
        push!(rxns_to_kill, !eval(parse(value)))
    end

    killed_reactions = rxns_affected[find(rxns_to_kill)]

    map(x -> set_col_bounds(lp, x, 0.0, 0.0), killed_reactions)
    solve(lp)

    return GeneKnockout(get_objective_value(lp), get_solution(lp), killed_reactions, gene_effect)
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    knockout_reactions(model, reactions)

Attempt to find the simplest combination of genes that need to be knocked out to disable the selected reactions
"""
function knockout_reactions(model::Model, reaction::String)
    rxn_idx = find(model.rxns .== reaction)[1]

    if isempty(rxn_idx)
        warn("Reaction " * reaction * " not found")
        return
    end

    return knockout_reactions(model, [rxn_idx])
end

function knockout_reactions{T <: String}(model::Model, reactions::Array{T})
    rxn_idx = find(model.rxns .== reactions)

    if length(rxn_idx) != length(reactions)
        warn("Not all reactions found in model")
    end

    return knockout_reactions(model, rxn_idx)
end

function knockout_reactions(model::Model, reactions::Number)
    return knockout_reactions(model, [reactions])
end

function knockout_reactions{T <: Number}(model::Model, reactions::Array{T})
    gene_list, rule_list, saved_names, _ = abstract_gene_names(model)
    gene_list = Dict(zip(gene_list, collect(1:length(gene_list))))

    # get the gene reaction rules for the reactions I want to disable
    knockout_rules = rule_list[reactions]

    # check if any of the rules is "", i.e empty i.e not "knockoutable"
    sanity = find(map(x -> !contains(x, "gene_"), knockout_rules))
    if length(sanity) > 0
        for index in sanity
            warn("Reaction " * string(reactions[index]), " cannot be knocked out")
        end
        deleteat!(knockout_rules, sanity)
    end

    # find the genes that appear in the reactions, a1,a2,a3...b1,b2,b3....n1,n2...
    # put that list of genes into an array of arrays, RG = [[a1,a2,a3...], [b1,b2,b3...],......[n1,n2...]]
    r_g_array = Array[]
    for (index, rule) in enumerate(knockout_rules)
        push!(r_g_array, genes_in_reaction(rule))
    end

    # locate the smallest gene knockout every reaction needs
    Combo_Array = Array[]
    for (index, rule) in enumerate(knockout_rules)#, n in 1:length(r_g_array[index])
        push!(Combo_Array, [])
        knockout_found = false

        n = 0
        while (n <= length(r_g_array[index])) & (!knockout_found)
            n += 1
            for combo in combinations(r_g_array[index], n)
                splt = split(rule)
                for gene in combo, (i,v) in enumerate(splt)
                    splt[i] = replace(splt[i], gene, "false")
                end

                for (i,v) in enumerate(splt)
                    if startswith(v,"gene_")
                        splt[i] = "true"
                    end
                end

                if !eval(parse(join(splt, " ")))
                    knockout_found = true
                    push!(Combo_Array[index],combo)
                end

            end
        end
    end

    Reaction_Matrices = Array[]
    for (i,v) in enumerate(Combo_Array)
        push!(Reaction_Matrices, zeros(Int64, length(v[1]), length(v)))
        for (j,k) in enumerate(v)
            for(o,p) in enumerate(k)
                Reaction_Matrices[i][o,j] = gene_list[p]
            end
        end
    end
    return Reaction_Matrices
    # With help from Tasos Papastylianou
    # http://stackoverflow.com/a/39087691/4551168
    # Find the smalest combinations of columns between matrices
    # where every matrix represents reactions and columns represent the
    # smallest set of genes that can disable the reaction
    #########################################################################
    toJagged(x) = [x[:,i] for i in 1:size(x,2)];
    JaggedMatrices = [toJagged(x) for x in Reaction_Matrices];

    Combined = [unique(i) for i in JaggedMatrices[1]];
    for n in 2:length(JaggedMatrices)
        Combined = [unique([i;j]) for i in Combined, j in JaggedMatrices[n]];
    end

    # Find the minimum
    Lengths         = [length(s) for s in Combined];
    Minima          = findin(Lengths, min(Lengths...));
    SubscriptsArray = ind2sub(size(Lengths), Minima);
    #########################################################################

    genes = Combined[SubscriptsArray[1]]

    gene_list = invert_dict(gene_list)

    Genes = Array(AbstractString, length(genes[1]), length(genes))

    for (i,v) in enumerate(genes), (j,gene) in enumerate(v)
        Genes[j,i] = saved_names[gene_list[gene]]
    end

    return Genes
end

"""
    print_reaction_formula(model,reaction::String)
    print_reaction_formula(model,index::Number)

Print out the formula for a reaction
"""
function print_medium(model)
	# List medium contents
	idx = find_exchange_reactions(model)
	for i in idx
	  if model.lb[i] < 0
	    @printf("%s\t%1.4f", model.rxns[i], model.lb[i])
	  end
	end
end

function print_reaction_formula(model::Model,reaction::String)
	id = find_reaction(model, reaction)
	id = vcat(id...)

	for i in id
		print_reaction_formula(model, i)
	end

end

function print_reaction_formula(model::Model, index::Number)

		sparse_indices = collect(nzrange(model.S,index))
		coefficients = model.S.nzval[sparse_indices]
		metabolites = model.S.rowval[sparse_indices]
		key = deepcopy(model.mets[index])

		right = ""
		left = ""

		for (i,v) in enumerate(coefficients)

			if v > 0
				right = string(right,v," ", model.mets[metabolites[i]], " + ")
			end
			if v < 0
				#left = string(left,metabolites[i]," + ")
				left = string(left,abs(v)," ", model.mets[metabolites[i]], " + ")
			end
		end

		right = right[1:end-2]
		left = left[1:end-2]
		println(left, " -> ", right)
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    solution = robustness_analysis(model, reactions; objective[, pts = [], direction = "max"])

Perform a robustness analysis for any number of fixed reactions, specified by indices or reaction names

**Example**
To perform a robustness analysis on reaction 13 againts reactions 5,8 and 11, where the point resolution
for reactions 5,8 and 11 is 4,2 and 10, respectively

    robustness_analysis(model, [5,8,11], 13, [4,2,10], "max")
    Robustness Analysis
                          result :   (2,3,4) Array
                          ranges :
                                  reaction :          range:
                                         5 :     (-0.0,20.0)
                                         8 :      (0.0,20.0)
                                        11 :    (8.39,175.0)

The method returns a **Robustness** object with fields

 * **solution.result**

    a multidimensional matrix which contains the flux value of the objective for each point

 * **solution.reactions**

    names of the reactions

 * **solution.ranges**

    the flux values of the reactions at specific points

        solution.ranges[1]
        4-element Array{Float64,1}:
         -1.3884e-28
          6.66667
         13.3333
         20.0

    So reaction 1 (reaction 5, ACONTb) is fixed at 6.6666 at point 2

To see the value of the objective, reaction 13, when reaction 5
is fixed at point 1, reaction 8 at point 2 and reaction 11 at
point 10

    solution[1,2,10]
    6.365820260756199e-16

To see every point where the objective is between 0.5 and 0.9

    range(solution, 0.5, 0.9)
    4-element Array{Tuple{Int64,Int64,Int64},1}:
     (1,2,3)
     (2,2,3)
     (3,2,3)
     (4,2,3)

To get the entire vector for reaction 5 when reaction 8 is fixed at
point 2 and reaction 11 is fixed at point 6

    solution[:,2,6]
    4-element Array{Float64,1}:
     0.425313
     0.314254
     0.203194
     0.0921351

Return a matrix with the objective reactions optimal value at each
value of the control reaction

**Plotting**

It is reccomended to use PyPlot to plot in Julia.

using PyPlot, to plot a 2D image of the effect of reaction 5
while reaction 8 is fixed at point 2 and reaction 11 at point 6

    PyPlot.plot(a[:,2,6])

to plot a 3D image in PyPlot for reactions 5 and 8 while
reaction 11 is fixed at point 6

    PyPlot.surf(a[:,:,6])
"""
function robustness_analysis(model::Model, reactions::Array{String}; objective::String = "", pts = [], direction::String = "max")
    return robustness_analysis(model, reacitons, objective, pts, direction)
end

function robustness_analysis(model::Model, rxns::Array{String}, objective::String = "", pts = [], direction::String = "max")
    rxns_idx = map(x -> find(x .== model.rxns), rxns)
    rxns_idx = filter(x -> length(x) == 1, rxns_idx)
    rxns_idx = vcat(rxns_idx...)

    if length(rxns_idx) != length(rxns)
        warn("One or more reactions not found in model.rxns")
        return
    end

    if objective == ""
        objective = model.rxns[find(model.c)[1]]
    end

    obj_idx = find(model.rxns .== objective)

    if isempty(obj_idx)
        warn("Objective reaction not found in model.rxns")
    end

    obj_idx = obj_idx[1]

    if pts == []
        pts = fill(20, length(rxns_idx))
    end

    if length(pts) != length(rxns)
        warn("Please supply as many points as you supply reactions")
        return
    end

    return robustness_analysis(model, rxns_idx, obj_idx, pts, direction)
end

function robustness_analysis{T <: Number}(model::Model, rxns::Array{T}; objective::Number = 0, pts = [], direction::String = "max")
    robustness_analysis(model, rxns; objective, pts, direction)
end

function robustness_analysis{T <: Number}(model::Model, rxns::Array{T}, objective::Number = 0, pts = [], direction::String = "max")
	if objective == 0
		objective = find(model.c)[1]
	end
    if pts == []
        pts = fill(20, length(rxns))
    end

    # just to relieve the code of endless length(model.c)
    num_rxns = length(rxns)

    # create a collection of indices for every single "fixation" of the reactions
    # so that if I mean to fix reaction 5 in 4 different points and reaction 8 in 2 points
    # combos will initially be [1,1], [1,2], [1,3], [1,4], ... [4,3], [4,4]
    # but idx will filter out those indices that I wont want, such as [1,4], point 4 is
    # out of range for reaction 8, as I only want 1 and 2 for 8
    # so in the end, I will have
    # [1,1],  [1,2],  [2,1], ... , [4,1], [4,2]
    combos = collect(combinations(repmat(1:maximum(pts), 1, num_rxns), num_rxns))
    combos = unique(combos)
    combos = sort(combos, lt=lexless)

    idx = map(y -> map(x -> combos[x][y] > pts[y], 1:length(combos)), 1:length(pts))
    idx = map(x -> find(x), idx)
    idx = union(idx...)
    idx = sort(idx, rev=true)
    map(x -> splice!(combos, x), idx)
    num_combos = length(combos)

    linspaces = map(x -> Array(Float64, x), pts)

    results = Array(Float64, num_combos)

    lp = setup_lp(model, direction)

    for (i,v) in enumerate(rxns)
        minF, maxF = min_max_lp(lp, v)
        linspaces[i] = collect(linspace(minF, maxF, pts[i]))
    end

    change_objective_coef(lp, objective)

    last_combo = zeros(Int64, length(rxns))
    for (i, combo) in enumerate(combos)
        combos_that_change = find(last_combo .!= combo)
        map((x,y) -> set_col_bounds(lp, rxns[x], linspaces[x][y]), combos_that_change, combo)

        results[i] = answer_lp(lp)
    end

    # format the solution
    result_matrix = reshape(results, tuple(combos[end]...))

    return Robustness(result_matrix, rxns, linspaces)
end

# -------------------------------------------------------------------
# -------------------------------------------------------------------

"""
    synthetic_lethal_genes(model, cutoff = 0.1, num_runs = 100)

Check the model for essential genes, conditionally essential genes and non-essential genes

Returns a ``SLG`` type, which contains the fields ``ess``, ``cond_ess`` and ``non_ess``,::

* ``cutoff`` represents the minimum biomass flux as a fraction of the wild-type flux.

* ``num_runs`` indicates how many times the algorithm runs, higher number gives better results, but takes longer.

**Example**

To find essentiality with biomass fixed at ``0.1`` ::

    julia> slg = synthetic_lethal_genes(model, 0.1, 200)
    run 200 of 200 Number of conditionally essential genes found: 115

    Type: SLG
         Field                  Contains    Size
           ess                Essential:       7
      cond_ess    Conditional Essential:     115
       non_ess            Non-Essential:      15
"""
function synthetic_lethal_genes(model, cutoff = 0.1, nruns = 100)
    # just a counter
    cond_count = 0

    # calculate the wt_growth and gene_flow first.
    # see description of prep_genes() for information. The purpose is to
    # format the data for quicker and easier lookup
    wt_growth = fba(model).obj * cutoff
    gene_flow = gene_deletion(model,1).g_f
    rules, new_genes, gene_lookup, effects, rxn_genes = prep_genes(model)

    ngenes = length(model.genes)
    nrxns = length(model.rxns)

    # quick_gene is a dict with keys 1 => model.genes[1] ....
    # for faster lookup when replacing them in strings
    quick_gene = Dict(zip(1:ngenes, new_genes))
    quick_essential = Dict()
    quick_cond = Dict()

    # will be filled with evaluated booleans
    # rxn_rules are converted to true/false statements
    # that need to be evaluated, and there is tremendous benefit
    # in only evaluating once instdead of repeadidly evaluating fx
    # the expression "true | false"
    truth_table = Dict()

    # will be filled with shuffled up indices of genes
    perms = zeros(Int64, ngenes, nruns)

    # fill up the essential  array of essential genes
    # and sort them. In addition, create a dictionary quickgene
    # for quick lookup
    for (key,val) in gene_flow
        if val < wt_growth
            quick_essential[find(key .== model.genes)[1]] = val
        end
    end

    # fill up the permutation matrix "perm" of gene indices and
    # turn all "essentials" inside to 0 so that we dont
    # repeatedly check genes that we know will kill
    for i = 1:nruns
        perms[:,i] = randperm(ngenes)
    end
    perms = map(x -> haskey(quick_essential, x) ? 0 : x, perms)

    lp = setup_lp(model)
    lb, ub = get_variable_bounds(lp)
    for i in 1:nruns
        print(string("\n\u1b[1Frun ", i, " of ", nruns, " Number of conditionally essential genes found: ", cond_count))

        map(x -> set_col_bounds(lp, x, lb[x], ub[x]), 1:nrxns)

        tmp_rules = deepcopy(rules)

        # must do this every round, restard the lp and get its bounds

        disabled = []

        for j in 1:ngenes
            gene_idx = perms[j,i]
            if gene_idx == 0
                continue
            end

            affected_idx = effects[gene_idx]
            affected_rules = deepcopy(tmp_rules[affected_idx])

            for (k,rule) in enumerate(affected_rules)
                affected_rules[k] = replace(affected_rules[k], quick_gene[gene_idx], " false ")
                affected_rules[k] = replace(affected_rules[k], r"x\[[0-9]*\]", " true ")
            end

            for expression in affected_rules
                if !haskey(truth_table, expression)
                    truth_table[expression] = eval(parse(expression))
                end
            end

            # find the indices of the reactions that a gene
            # disables, but only the indices of those reactions
            # that haven yet been disabled
            dead = []
            for (k,r) in enumerate(affected_rules)
                if !truth_table[r]
                    push!(dead, affected_idx[k])
                end
            end
            dead = setdiff(dead, disabled)
            disabled = vcat(disabled, dead)

            map(x -> set_col_bounds(lp, x, 0.0, 0.0), dead)

            if answer_lp(lp) < wt_growth
                if !haskey(quick_cond, gene_idx)
                    quick_cond[gene_idx] = []
                    cond_count += 1
                end
                push!(quick_cond[gene_idx], disabled)
                map(x -> set_col_bounds(lp, x, lb[x], ub[x]), dead)
            else
                for (k,rule) in enumerate(tmp_rules)
                    tmp_rules[k] = replace(tmp_rules[k], quick_gene[gene_idx], " false ")
                end
            end
        end
    end
    println()
    cond_ess = sort(collect(keys(quick_cond)))
    essential = sort(collect(keys(quick_essential)))

    non_essential = setdiff(collect(1:ngenes), cond_ess)
    non_essential = setdiff(non_essential, essential)

    return SLG(quick_essential, quick_cond, non_essential)
end

max_stars_repo_name: isebarn/CBM

max_issues_repo_name: isebarn/CBM

max_forks_repo_name: isebarn/CBM


text:
# ------------------------------------------------------------------------------------------
# # Bucles
#
# Temas:
# 1. bucles `while`
# 2. bucles `for`
#

#
# ## bucles while
#
# La sintaxis para un `while` es
#
# ```julia
# while *condición*
#     *cuerpo de bucle*
# end
# ```
#
# Por ejemplo, usaremos un `while` para contar o iterar sobre un arreglo.
# ------------------------------------------------------------------------------------------

n = 0
while n < 10
    n += 1
    println(n)
end

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

i = 1
while i <= length(myfriends)
    friend = myfriends[i]
    println("Hola $friend, ¡gusto en verte!")
    i += 1
end

# ------------------------------------------------------------------------------------------
# ## bucles for
#
# La sintaxis para un bucle `for` es
#
# ```julia
# for *var* in *iterable de bucle*
#     *cuerpo de bucle*
# end
# ```
#
# Podemos usar un bucle for para generar los datos anteriores
# ------------------------------------------------------------------------------------------

for n in 1:10
    println(n)
end

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

for friend in myfriends
    println("Hi $friend, it's great to see you!")
end

# ------------------------------------------------------------------------------------------
# Nota: se puede reemplazar `in` con `=` ó `∈`.
# ------------------------------------------------------------------------------------------

for n = 1:10
    println(n)
end

for n ∈ 1:10
    println(n)
end

# ------------------------------------------------------------------------------------------
# Ahora usemos bucles `for` para crear tablas de sumas donde el valor de cada entrada es la
# suma de los índices del renglón y la columna.

#
# Primero, inicializamos el arreglo con puros 0s.
# ------------------------------------------------------------------------------------------

m, n = 5, 5
A = fill(0, (m, n))

for i in 1:m
    for j in 1:n
        A[i, j] = i + j
    end
end
A

# ------------------------------------------------------------------------------------------
# Aquí va un poquito de azúcar sintáctica para el bucle `for`
# ------------------------------------------------------------------------------------------

B = fill(0, (m, n))

for i in 1:m, j in 1:n
    B[i, j] = i + j
end
B

# ------------------------------------------------------------------------------------------
# La manera mas "Juliana" de crear esta tabla es por medio de un *array comprehension /
# comprehensión de arreglo*.
# ------------------------------------------------------------------------------------------

C = [i + j for i in 1:m, j in 1:n]

# ------------------------------------------------------------------------------------------
# En el próximo ejemplo, embebimos un arreglo de comprehensión en un bucle `for`, generando
# tablas de adición de tamaño creciente.
# ------------------------------------------------------------------------------------------

for n in 1:10
   A = [i + j for i in 1:n, j in 1:n]
   display(A)
end

# ------------------------------------------------------------------------------------------
# ### Ejercicios
#
# 4.1 Crea un diccionario `squares`, que tiene llaves de valores de 1 a 100. El valor
# asociado a cada llave es el cuadrado de la llave. Guarda los valores asociados a las
# llaves pares como enteros y las impares como cadenas. Por ejemplo,
#
# ```julia
# squares[10] == 100
# squares[11] == "121"
# ```
#
# (¡No necesitas condicionales para esto!)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 4.2 Usa `fill` para crea una matriz de `10x10` de solo `0`'s. Pobla las primeras 10
# entradas con el índice de esa entrada. ¿Julia usa el orden de primero columna o primero
# renglón? (O sea, ¿el "segundo" elemento es el de la primera columna en el primer renglón,
# ó es el de el primer renglón en la segunda columna?)
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
export mutate

function mutate_add(parent::GRN)
    ids = [copy(parent.ids); rand()]
    enh = [copy(parent.enh); rand()]
    inh = [copy(parent.inh); rand()]
    nreg = parent.nreg + 1
    GRN(parent.nin, parent.nout, nreg, ids, enh, inh, parent.beta, parent.delta)
end

function mutate_delete(parent::GRN)
    to_remove = rand((parent.nin+parent.nout+1):length(parent.ids))
    ids = deleteat!(copy(parent.ids), to_remove)
    enh = deleteat!(copy(parent.enh), to_remove)
    inh = deleteat!(copy(parent.inh), to_remove)
    nreg = parent.nreg - 1
    GRN(parent.nin, parent.nout, nreg, ids, enh, inh, parent.beta, parent.delta)
end

function mutate_modify(parent::GRN, config::Config)
    target = rand(1:(length(parent.ids)+2))
    ids = copy(parent.ids); enh = copy(parent.enh); inh = copy(parent.inh)
    beta = parent.beta; delta = parent.delta
    if target <= length(parent.ids)
        tag = rand(1:3)
        if tag == 1
            ids[target] = rand()
        elseif tag == 2
            enh[target] = rand()
        else
            inh[target] = rand()
        end
    elseif target == (length(parent.ids) + 1)
        beta = rand() * (config.beta_max - config.beta_min) + config.beta_min
    else
        delta = rand() * (config.delta_max - config.delta_min) + config.delta_min
    end
    GRN(parent.nin, parent.nout, parent.nreg, ids, enh, inh, beta, delta)
end

function mutate(parent::GRN, config::Config)
    mutation = rand()
    if mutation < config.mutate_add_rate
        child = mutate_add(parent)
    elseif (mutation < (config.mutate_add_rate + config.mutate_delete_rate) &&
            length(parent.ids) > (parent.nin + parent.nout))
        child = mutate_delete(parent)
    else
        child = mutate_modify(parent, config)
    end
    child
end


text:
@testset "Sorts" begin
    sorts = [
        bubble_sort!
        bucket_sort!
        counting_sort!
        exchange_sort!
        insertion_sort!
        merge_sort!
        quick_sort!
        selection_sort!
    ]

    for f in sorts
        x = [3,5,1,4,2]
        y = [12,4,3,2,1]
        empty = []
        one = [1]

        f(x)
        @test x == [1,2,3,4,5]

        f(y)
        @test y == [1,2,3,4,12]

        f(empty)
        @test empty == []

        f(one)
        @test one == [1]
    end
end

max_stars_repo_name: gohpeijin/Julia

max_issues_repo_name: gohpeijin/Julia

max_forks_repo_name: gohpeijin/Julia


text:
s = open("input.txt") do file
    read(file, String)
end

passports = split(s, "\n\n");

function ValidPassport(passport)
    requiredFields = [
        r"byr:(19[2-9]\d|200[0-2])(\s|$)",
        r"iyr:20(1\d|20)(\s|$)",
        r"eyr:20(2\d|30)(\s|$)",
        r"hgt:(1([5-8]\d|9[0-3])cm|(59|6\d|7[0-6])in)(\s|$)",
        r"hcl:#[0-9a-f]{6}(\s|$)",
        r"ecl:(amb|blu|brn|gry|grn|hzl|oth)(\s|$)",
        r"pid:\d{9}(\s|$)"
    ];

    for field in requiredFields
        if !occursin(field, passport)
            return 0;
        end
    end

    return 1;
end

sum(map(ValidPassport, passports))
ValidPassport("eyr:1972 cid:100\nhcl:#18171d ecl:amb hgt:170 pid:186cm iyr:2018 byr:1926")

ValidPassport("iyr:2019 hcl:#602927 eyr:1967 hgt:170cm ecl:grn pid:012533040 byr:1946")

ValidPassport("hcl:dab227 iyr:2012
ecl:brn hgt:182cm pid:021572410 eyr:2020 byr:1992 cid:277")

ValidPassport("hgt:59cm ecl:zzz
eyr:2038 hcl:74454a iyr:2023
pid:3556412378 byr:2007")

ValidPassport("pid:087499704 hgt:74in ecl:grn iyr:2012 eyr:2030 byr:1980
hcl:#623a2f")

ValidPassport("eyr:2029 ecl:blu cid:129 byr:1989
iyr:2014 pid:896056539 hcl:#a97842 hgt:165cm")

ValidPassport("hcl:#888785
hgt:164cm byr:2001 iyr:2015 cid:88
pid:545766238 ecl:hzl
eyr:2022")

ValidPassport("iyr:2010 hgt:158cm hcl:#b6652a ecl:blu byr:1944 eyr:2021 pid:093154719")

    # byr (Birth Year) - four digits; at least 1920 and at most 2002.
    # iyr (Issue Year) - four digits; at least 2010 and at most 2020.
    # eyr (Expiration Year) - four digits; at least 2020 and at most 2030.
    # hgt (Height) - a number followed by either cm or in:
    #     If cm, the number must be at least 150 and at most 193.
    #     If in, the number must be at least 59 and at most 76.
    # hcl (Hair Color) - a # followed by exactly six characters 0-9 or a-f.
    # ecl (Eye Color) - exactly one of: amb blu brn gry grn hzl oth.
    # pid (Passport ID) - a nine-digit number, including leading zeroes.
    # cid (Country ID) - ignored, missing or not.


text:
using PlutoStaticHTML
using Pluto:
    Cell,
    Notebook,
    ServerSession
using Test:
    @testset,
    @test

const PKGDIR = string(pkgdir(PlutoStaticHTML))::String
const NOTEBOOK_DIR = joinpath(pkgdir(PlutoStaticHTML), "docs", "src")
const NOTEBOOK_PATH = joinpath(NOTEBOOK_DIR, "notebook.jl")

function pluto_notebook_content(code)
    return """
        ### A Pluto.jl notebook ###
        # v0.17.1

        using Markdown
        using InteractiveUtils

        # ╔═╡ a6dda572-3f2c-11ec-0eeb-69e2323a92de
        $(code)

        # ╔═╡ Cell order:
        # ╠═a6dda572-3f2c-11ec-0eeb-69e2323a92de
        """
end

function notebook2html!(notebook::Notebook; append_cells=Cell[], kwargs...)
    session = ServerSession()
    PlutoStaticHTML._append_cell!(notebook, append_cells)
    run_notebook!(notebook, session)
    return notebook2html(notebook; kwargs...)
end

include("context.jl")
include("html.jl")
include("build.jl")


text:
type ActiveSet
  indexes::Vector{Int64}
  numActive::Int64
end

function ActiveSet{T<:AbstractFloat}(
    x::StridedVector{T};
    options::LassoOptions = LassoOptions()
    )
  zero_thr = options.zerothr
  numElem = length(x)
  activeset = [1:numElem;]
  numActive = 0
  @inbounds for j = 1:numElem
    if abs(x[j]) >= zero_thr
      numActive += 1
      activeset[numActive], activeset[j] = activeset[j], activeset[numActive]
    end
  end
  ActiveSet(activeset, numActive)
end

function _add_violator!{T<:AbstractFloat}(
    activeset::ActiveSet,
    x::StridedVector{T},
    A::StridedMatrix{T},
    b::StridedVector{T},
    λ::StridedVector{T};
    options::LassoOptions   =   LassoOptions()
    )
  changed = false

  zerothr = options.zerothr
  numElem = length(x)
  numActive = activeset.numActive
  indexes = activeset.indexes
  # check for things to be removed from the active set
  i = 0
  @inbounds while i < numActive
    i = i + 1
    t = indexes[i]
    if abs(x[t]) < zerothr
      x[t] = zero(T)
      changed = true
      indexes[numActive], indexes[i] = indexes[i], indexes[numActive]
      numActive -= 1
      i = i - 1
    end
  end
  activeset.numActive = numActive

  I = 0
  V = zero(T)
  @inbounds for i=numActive+1:numElem
    t = indexes[i]
    tmp = _Axk(A, x, t, activeset) + b[t]
    nV = abs(tmp) - λ[t]
    if V < nV
      I = i
      V = nV
    end
  end
  if I > 0 && V > options.gradtol
    changed = true
    numActive += 1
    indexes[numActive], indexes[I] = indexes[I], indexes[numActive]
  end
  activeset.numActive = numActive
  changed
end

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

shrink{T<:AbstractFloat}(v::T, c::T) = v > c ? v - c : (v < -c ? v + c : zero(T))

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

# computes x'A[:, j]
function _Axk{T<:AbstractFloat}(A::StridedMatrix{T}, x::StridedVector{T}, j::Int64, activeset::ActiveSet)
  s = zero(T)
  indexes = activeset.indexes
  @inbounds for i=1:activeset.numActive
    ind = indexes[i]
    s += x[ind] * A[ind, j]
  end
  s
end

function _lasso!{T<:AbstractFloat}(
    x::StridedVector{T},
    A::StridedMatrix{T},
    b::StridedVector{T},
    λ::StridedVector{T},
    activeset::ActiveSet;
    options::LassoOptions   =   LassoOptions()
    )

  indexes = activeset.indexes
  iter = 0
  while true
    iter += 1
    maxUpdate = zero(T)
    @inbounds for i=1:activeset.numActive
      j = indexes[i]
      a = A[j, j]
      S0 = _Axk(A, x, j, activeset) - a * x[j] + b[j]
      newValue = shrink(-S0 / a, abs(λ[j] / a))
      if abs(x[j] - newValue) > maxUpdate
        maxUpdate = abs(x[j] - newValue)
      end
      x[j] = newValue
    end
    if iter > options.max_inner_iter || maxUpdate < options.xtol
      break
    end
  end
  x
end

lasso!{T<:AbstractFloat}(
    x::StridedVector{T},
    A::StridedMatrix{T},
    b::StridedVector{T},
    λ::T;
    options::LassoOptions   =   LassoOptions()
    ) = lasso!(x, A, b, fill(λ, length(x)); options=options)

function lasso!{T<:AbstractFloat}(
    x::StridedVector{T},
    A::StridedMatrix{T},
    b::StridedVector{T},
    λ::StridedVector{T};
    options::LassoOptions   =   LassoOptions()
    )

  maxoutiter = options.max_outer_iter
  activeset = ActiveSet(x; options=options)
  if activeset.numActive == 0
    _add_violator!(activeset, x, A, b, λ; options=options)
  end

  for outiter=1:maxoutiter
    _lasso!(x, A, b, λ, activeset; options=options)
    if ~_add_violator!(activeset, x, A, b, λ; options=options)
      break
    end
  end
  x, activeset
end

max_stars_repo_name: mlakolar/CDLasso

max_issues_repo_name: mlakolar/CDLasso

max_forks_repo_name: mlakolar/CDLasso


text:
"""
    AdvancedFileLogger(
        dir::AbstractString,
        file_name_pattern::AbstractString;
        log_format_function::Function=print_standard_format,
        min_level=Logging.Info
    )

Return an file logger with extra functionality (in contrast to a `FileLogger` which is only a interface to `SimpleLogger`).

### Parameters
- `dir` -- log file directory
- `file_name_pattern` -- file name pattern for rotating log files (e.g. `raw"\\a\\c\\c\\e\\s\\s-YYYY-mm-dd-HH-MM.\\l\\o\\g")`)
- `log_format_function` -- optional, default=`print_standard_format`
- `min_level` -- optional, default=`Logging.Info`

### Features
- Define a `DateFormat`/`String` file pattern for rotating log messages to a new file → argument `file_name_pattern`
- Provide a custom formatting function for writting to log files → argument `log_format_function`

### Returns
A `MinLevelLogger` which includes a `DatetimeRotatingFileLogger`.
"""
function AdvancedFileLogger(
    dir_path::AbstractString,
    file_name_pattern::AbstractString;
    log_format_function::Function=print_standard_format,
    min_level::Logging.LogLevel=Logging.Info,
)
    file_logger = DatetimeRotatingFileLogger(log_format_function, dir_path, file_name_pattern)
    return MinLevelLogger(file_logger, min_level)
end

# Adjust
const DISPLAY_SIZE = (100, 100)
"""
    print_standard_format(io::IO, log_args::NamedTuple)

Print log message and keyword arguments in a format similar to the one used by the [`ConsoleLogger`](@ref).
"""
function print_standard_format(io::IO, log_args::NamedTuple)
    # TODO: replace by eachsplit for Julia > 1.6
    levelstr = log_args.level == Logging.Warn ? "Warning" : string(log_args.level)
    msglines = split(chomp(string(log_args.message)::String), '\n')

    println(io, "[", levelstr, "] ", msglines[1]) #┌

    for i in 2:length(msglines)
        println(io, "│ ", msglines[i])
    end
    n_kwargs = length(log_args.kwargs)
    for (idx, (key, val)) in enumerate(log_args.kwargs)
        kwarg_prefix = (idx == n_kwargs) ? "└ " : "├ "

        if key == :back_trace
            msglines = split(chomp(string(val)::String), '\n')
        else
            # If the kwarg is not an Exception back_trace, use the ConsoleLogger logic to display the variable
            # NOTE: for the future, it might be reasonable, to also exclude kwargs of type String from this branch
            iob = IOBuffer()
            # DisplaySize defines "how much" of the val will be printed to file
            show(IOContext(iob, :limit => true, :displaysize => DISPLAY_SIZE), "text/plain", val)
            val_str = String(take!(iob))
            msglines = split(chomp(val_str), '\n')
        end
        if length(msglines) == 1
            println(io, kwarg_prefix, key, " = ",  msglines[1])
        else
            println(io, "├ ", key, " = ", msglines[1])
            for i in 2:(length(msglines) - 1)
                println(io, "│  ", msglines[i])
            end
            if idx == n_kwargs
                println(io, "└  ", msglines[end])
            else
                println(io, "│  ", msglines[end])
            end
        end
    end
    return nothing
end


text:
# Note: Functions surrounded by a comment blocks are there because `Vararg` is still allocating.
# When Vararg is fast enough, they can simply be removed.

#######################
# UndefBlocksInitializer #
#######################

"""
    UndefBlocksInitializer

Singleton type used in block array initialization, indicating the
array-constructor-caller would like an uninitialized block array. See also
undef_blocks (@ref), an alias for UndefBlocksInitializer().

Examples

≡≡≡≡≡≡≡≡≡≡
```julia
julia> BlockArray(undef_blocks, Matrix{Float32}, [1,2], [3,2])
2×2-blocked 3×5 BlockArrays.BlockArray{Float32,2,Array{Float32,2}}:
 #undef  #undef  #undef  │  #undef  #undef
 ------------------------┼----------------
 #undef  #undef  #undef  │  #undef  #undef
 #undef  #undef  #undef  │  #undef  #undef
 ```
"""
struct UndefBlocksInitializer end

"""
    undef_blocks

Alias for UndefBlocksInitializer(), which constructs an instance of the singleton
type UndefBlocksInitializer (@ref), used in block array initialization to indicate the
array-constructor-caller would like an uninitialized block array.

Examples

≡≡≡≡≡≡≡≡≡≡
```julia
julia> BlockArray(undef_blocks, Matrix{Float32}, [1,2], [3,2])
2×2-blocked 3×5 BlockArrays.BlockArray{Float32,2,Array{Float32,2}}:
 #undef  #undef  #undef  │  #undef  #undef
 ------------------------┼----------------
 #undef  #undef  #undef  │  #undef  #undef
 #undef  #undef  #undef  │  #undef  #undef
 ```
"""
const undef_blocks = UndefBlocksInitializer()

##############
# BlockArray #
##############

function _BlockArray end

"""
    BlockArray{T, N, R <: AbstractArray{T, N}} <: AbstractBlockArray{T, N}

A `BlockArray` is an array where each block is stored contiguously. This means that insertions and retrieval of blocks
can be very fast and non allocating since no copying of data is needed.

In the type definition, `R` defines the array type that each block has, for example `Matrix{Float64}`.
"""
struct BlockArray{T, N, R <: AbstractArray{T, N}} <: AbstractBlockArray{T, N}
    blocks::Array{R, N}
    block_sizes::BlockSizes{N}

    global function _BlockArray(blocks::Array{R, N}, block_sizes::BlockSizes{N}) where {T, N, R <: AbstractArray{T, N}}
        new{T, N, R}(blocks, block_sizes)
    end
end

# Auxilary outer constructors
function _BlockArray(blocks::Array{R, N}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R <: AbstractArray{T, N}}
    return _BlockArray(blocks, BlockSizes(block_sizes...))
end

const BlockMatrix{T, R <: AbstractMatrix{T}} = BlockArray{T, 2, R}
const BlockVector{T, R <: AbstractVector{T}} = BlockArray{T, 1, R}
const BlockVecOrMat{T, R} = Union{BlockMatrix{T, R}, BlockVector{T, R}}

################
# Constructors #
################

@inline function _BlockArray(::Type{R}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R <: AbstractArray{T, N}}
    _BlockArray(R, BlockSizes(block_sizes...))
end

function _BlockArray(::Type{R}, block_sizes::BlockSizes{N}) where {T, N, R <: AbstractArray{T, N}}
    n_blocks = nblocks(block_sizes)
    blocks = Array{R, N}(undef, n_blocks)
    _BlockArray(blocks, block_sizes)
end

@inline function undef_blocks_BlockArray(::Type{R}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R <: AbstractArray{T, N}}
    _BlockArray(R, block_sizes...)
end

"""
Constructs a `BlockArray` with uninitialized blocks from a block type `R` with sizes defind by `block_sizes`.

```jldoctest
julia> BlockArray(undef_blocks, Matrix{Float64}, [1,3], [2,2])
2×2-blocked 4×4 BlockArrays.BlockArray{Float64,2,Array{Float64,2}}:
 #undef  │  #undef  #undef  #undef  │
 --------┼--------------------------┼
 #undef  │  #undef  #undef  #undef  │
 #undef  │  #undef  #undef  #undef  │
 --------┼--------------------------┼
 #undef  │  #undef  #undef  #undef  │
```
"""
@inline function BlockArray(::UndefBlocksInitializer, ::Type{R}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R <: AbstractArray{T, N}}
    undef_blocks_BlockArray(R, block_sizes...)
end

@inline function BlockArray{T}(::UndefBlocksInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N}
    BlockArray(undef_blocks, Array{T,N}, block_sizes...)
end

@inline function BlockArray{T,N}(::UndefBlocksInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N}
    BlockArray(undef_blocks, Array{T,N}, block_sizes...)
end

@inline function BlockArray{T,N,R}(::UndefBlocksInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R <: AbstractArray{T, N}}
    BlockArray(undef_blocks, R, block_sizes...)
end

@generated function initialized_blocks_BlockArray(::Type{R}, block_sizes::BlockSizes{N}) where {T, N, R <: AbstractArray{T, N}}
    return quote
        block_arr = _BlockArray(R, block_sizes)
        @nloops $N i i->(1:nblocks(block_sizes, i)) begin
            block_index = @ntuple $N i
            setblock!(block_arr, similar(R, blocksize(block_sizes, block_index)), block_index...)
        end

        return block_arr
    end
end

function initialized_blocks_BlockArray(::Type{R}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R <: AbstractArray{T, N}}
    initialized_blocks_BlockArray(R, BlockSizes(block_sizes...))
end

@inline function BlockArray{T}(::UndefInitializer, block_sizes::BlockSizes{N}) where {T, N}
    initialized_blocks_BlockArray(Array{T, N}, block_sizes)
end

@inline function BlockArray{T, N}(::UndefInitializer, block_sizes::BlockSizes{N}) where {T, N}
    initialized_blocks_BlockArray(Array{T, N}, block_sizes)
end

@inline function BlockArray{T, N, R}(::UndefInitializer, block_sizes::BlockSizes{N}) where {T, N, R <: AbstractArray{T, N}}
    initialized_blocks_BlockArray(R, block_sizes)
end

@inline function BlockArray{T}(::UndefInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N}
    initialized_blocks_BlockArray(Array{T, N}, block_sizes...)
end

@inline function BlockArray{T, N}(::UndefInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N}
    initialized_blocks_BlockArray(Array{T, N}, block_sizes...)
end

@inline function BlockArray{T, N, R}(::UndefInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R <: AbstractArray{T, N}}
    initialized_blocks_BlockArray(R, block_sizes...)
end

function BlockArray(arr::AbstractArray{T, N}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T,N}
    for i in 1:N
        if sum(block_sizes[i]) != size(arr, i)
            throw(DimensionMismatch("block size for dimension $i: $(block_sizes[i]) does not sum to the array size: $(size(arr, i))"))
        end
    end
    BlockArray(arr, BlockSizes(block_sizes...))
end

@generated function BlockArray(arr::AbstractArray{T, N}, block_sizes::BlockSizes{N}) where {T,N}
    return quote
        block_arr = _BlockArray(typeof(arr), block_sizes)
        @nloops $N i i->(1:nblocks(block_sizes, i)) begin
            block_index = @ntuple $N i
            indices = globalrange(block_sizes, block_index)
            setblock!(block_arr, arr[indices...], block_index...)
        end

        return block_arr
    end
end

# Convert AbstractArrays that conform to block array interface
convert(::Type{BlockArray{T,N,R}}, A::BlockArray{T,N,R}) where {T,N,R} = A
convert(::Type{BlockArray{T,N}}, A::BlockArray{T,N}) where {T,N} = A
convert(::Type{BlockArray{T}}, A::BlockArray{T}) where {T} = A
convert(::Type{BlockArray}, A::BlockArray) = A

BlockArray{T, N}(A::AbstractArray{T2, N}) where {T,T2,N} =
    BlockArray(Array{T, N}(A), blocksizes(A))
BlockArray{T1}(A::AbstractArray{T2, N}) where {T1,T2,N} = BlockArray{T1, N}(A)
BlockArray(A::AbstractArray{T, N}) where {T,N} = BlockArray{T, N}(A)

convert(::Type{BlockArray{T, N}}, A::AbstractArray{T2, N}) where {T,T2,N} =
    BlockArray(convert(Array{T, N}, A), blocksizes(A))
convert(::Type{BlockArray{T1}}, A::AbstractArray{T2, N}) where {T1,T2,N} =
    convert(BlockArray{T1, N}, A)
convert(::Type{BlockArray}, A::AbstractArray{T, N}) where {T,N} =
    convert(BlockArray{T, N}, A)

copy(A::BlockArray) = _BlockArray(copy.(A.blocks), copy(A.block_sizes))

################################
# AbstractBlockArray Interface #
################################
@inline blocksizes(block_array::BlockArray) = block_array.block_sizes

@inline function getblock(block_arr::BlockArray{T,N}, block::Vararg{Int, N}) where {T,N}
    @boundscheck blockcheckbounds(block_arr, block...)
    block_arr.blocks[block...]
end

@inline function Base.getindex(block_arr::BlockArray{T,N}, blockindex::BlockIndex{N}) where {T,N}
    @boundscheck checkbounds(block_arr.blocks, blockindex.I...)
    @inbounds block = block_arr.blocks[blockindex.I...]
    @boundscheck checkbounds(block, blockindex.α...)
    @inbounds v = block[blockindex.α...]
    return v
end

###########################
# AbstractArray Interface #
###########################

@inline function Base.similar(block_array::BlockArray{T,N}, ::Type{T2}) where {T,N,T2}
    _BlockArray(similar(block_array.blocks, Array{T2, N}), copy(blocksizes(block_array)))
end

@inline function Base.getindex(block_arr::BlockArray{T, N}, i::Vararg{Int, N}) where {T,N}
    @boundscheck checkbounds(block_arr, i...)
    @inbounds v = block_arr[global2blockindex(blocksizes(block_arr), i)]
    return v
end

@inline function Base.setindex!(block_arr::BlockArray{T, N}, v, i::Vararg{Int, N}) where {T,N}
    @boundscheck checkbounds(block_arr, i...)
    @inbounds block_arr[global2blockindex(blocksizes(block_arr), i)] = v
    return block_arr
end

############
# Indexing #
############

function _check_setblock!(block_arr::BlockArray{T, N}, v, block::NTuple{N, Int}) where {T,N}
    for i in 1:N
        if size(v, i) != blocksize(block_arr, i, block[i])
            throw(DimensionMismatch(string("tried to assign $(size(v)) array to ", blocksize(block_arr, block), " block")))
        end
    end
end

@inline function setblock!(block_arr::BlockArray{T, N}, v, block::Vararg{Int, N}) where {T,N}
    @boundscheck blockcheckbounds(block_arr, block...)
    @boundscheck _check_setblock!(block_arr, v, block)
    @inbounds block_arr.blocks[block...] = v
    return block_arr
end

@propagate_inbounds function Base.setindex!(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N}) where {T,N}
    getblock(block_array, block_index.I...)[block_index.α...] = v
end

########
# Misc #
########

@generated function Base.Array(block_array::BlockArray{T, N, R}) where {T,N,R}
    # TODO: This will fail for empty block array
    return quote
        block_sizes = blocksizes(block_array)
        arr = similar(block_array.blocks[1], size(block_array)...)
        @nloops $N i i->(1:nblocks(block_sizes, i)) begin
            block_index = @ntuple $N i
            indices = globalrange(block_sizes, block_index)
            arr[indices...] = getblock(block_array, block_index...)
        end

        return arr
    end
end

@generated function copyto!(block_array::BlockArray{T, N, R}, arr::R) where {T,N,R <: AbstractArray}
    return quote
        block_sizes = blocksizes(block_array)

        @nloops $N i i->(1:nblocks(block_sizes, i)) begin
            block_index = @ntuple $N i
            indices = globalrange(block_sizes, block_index)
            copyto!(getblock(block_array, block_index...), arr[indices...])
        end

        return block_array
    end
end

function Base.fill!(block_array::BlockArray, v)
    for block in block_array.blocks
        fill!(block, v)
    end
    block_array
end

function lmul!(α::Number, block_array::BlockArray)
    for block in block_array.blocks
        lmul!(α, block)
    end
    block_array
end

function rmul!(block_array::BlockArray, α::Number)
    for block in block_array.blocks
        rmul!(block, α)
    end
    block_array
end


text:
#=
Discrete Decision Processes

@author : Daisuke Oyama, Spencer Lyon, Matthew McKay

@date: 24/Sep/2015

References
----------

http://quant-econ.net/jl/ddp.html

Notes
-----
1. This currently implements:
    a. Value Iteration,
    b. Policy Iteration, and
    c. Modified Policy Iteration
   For:
    a. Dense Matrices
    b. State-Action Pair Formulation

=#

import Base: *

#------------------------#
#-Types and Constructors-#
#------------------------#

"""
DiscreteDP type for specifying paramters for discrete dynamic programming model

##### Parameters

- `R::Array{T,NR}` : Reward Array
- `Q::Array{T,NQ}` : Transition Probability Array
- `beta::Float64`  : Discount Factor
- `a_indices::Nullable{Vector{Tind}}`: Action Indices. Null unless using
  SA formulation
- `a_indptr::Nullable{Vector{Tind}}`: Action Index Pointers. Null unless using
  SA formulation

##### Returns

- `ddp::DiscreteDP` : DiscreteDP object

"""
type DiscreteDP{T<:Real,NQ,NR,Tbeta<:Real,Tind}
    R::Array{T,NR}                     # Reward Array
    Q::Array{T,NQ}                     # Transition Probability Array
    beta::Tbeta                        # Discount Factor
    a_indices::Nullable{Vector{Tind}}  # Action Indices
    a_indptr::Nullable{Vector{Tind}}   # Action Index Pointers

    function DiscreteDP(R::Array, Q::Array, beta::Real)
        # verify input integrity 1
        if NQ != 3
            msg = "Q must be 3-dimensional without s-a formulation"
            throw(ArgumentError(msg))
        end
        if NR != 2
            msg = "R must be 2-dimensional without s-a formulation"
            throw(ArgumentError(msg))
        end
        (beta < 0 || beta >= 1) &&  throw(ArgumentError("beta must be [0, 1)"))

        # verify input integrity 2
        num_states, num_actions = size(R)
        if size(Q) != (num_states, num_actions, num_states)
            throw(ArgumentError("shapes of R and Q must be (n,m) and (n,m,n)"))
        end

        # check feasibility
        R_max = s_wise_max(R)
        if any(R_max .== -Inf)
            # First state index such that all actions yield -Inf
            s = find(R_max .== -Inf) #-Only Gives True
            throw(ArgumentError("for every state the reward must be finite for
                some action: violated for state $s"))
        end

        # here the indices and indptr are null.
        _a_indices = Nullable{Vector{Int}}()
        a_indptr = Nullable{Vector{Int}}()

        new(R, Q, beta, _a_indices, a_indptr)
    end

    # Note: We left R, Q as type Array to produce more helpful error message with regards to shape.
    # R and Q are dense Arrays
    function DiscreteDP(R::AbstractArray, Q::AbstractArray, beta::Real,
                        s_indices::Vector, a_indices::Vector)
        # verify input integrity 1
        if NQ != 2
            throw(ArgumentError("Q must be 2-dimensional with s-a formulation"))
        end
        if NR != 1
            throw(ArgumentError("R must be 1-dimensional with s-a formulation"))
        end
        (beta < 0 || beta >= 1) && throw(ArgumentError("beta must be [0, 1)"))

        # verify input integrity (same length)
        num_sa_pairs, num_states = size(Q)
        if length(R) != num_sa_pairs
            throw(ArgumentError("shapes of R and Q must be (L,) and (L,n)"))
        end
        if length(s_indices) != num_sa_pairs
            msg = "length of s_indices must be equal to the number of s-a pairs"
            throw(ArgumentError(msg))
        end
        if length(a_indices) != num_sa_pairs
            msg = "length of a_indices must be equal to the number of s-a pairs"
            throw(ArgumentError(msg))
        end

        if _has_sorted_sa_indices(s_indices, a_indices)
            a_indptr = Array(Int64, num_states+1)
            _a_indices = copy(a_indices)
            _generate_a_indptr!(num_states, s_indices, a_indptr)
        else
            # transpose matrix to use Julia's CSC; now rows are actions and
            # columns are states (this is why it's called as_ptr not sa_ptr)
            m = maximum(a_indices)
            n = maximum(s_indices)
            msg = "Duplicate s-a pair found"
            as_ptr = sparse(a_indices, s_indices, 1:num_sa_pairs, m, n,
                            (x,y)->throw(ArgumentError(msg)))
            _a_indices = as_ptr.rowval
            a_indptr = as_ptr.colptr

            R = R[as_ptr.nzval]
            Q = Q[as_ptr.nzval, :]
        end

        # check feasibility
        aptr_diff = diff(a_indptr)
        if any(aptr_diff .== 0.0)
            # First state index such that no action is available
            s = find(aptr_diff .== 0.0)  # Only Gives True
            throw(ArgumentError("for every state at least one action
                must be available: violated for state $s"))
        end

        # package into nullables before constructing type
        _a_indices = Nullable{Vector{Tind}}(_a_indices)
        a_indptr = Nullable{Vector{Tind}}(a_indptr)

        new(R, full(Q), beta, _a_indices, a_indptr)
    end
end

"""
DiscreteDP type for specifying parameters for discrete dynamic programming
model Dense Matrix Formulation

##### Parameters

- `R::Array{T,NR}` : Reward Array
- `Q::Array{T,NQ}` : Transition Probability Array
- `beta::Float64`  : Discount Factor

##### Returns

- `ddp::DiscreteDP` : Constructor for DiscreteDP object
"""
DiscreteDP{T,NQ,NR,Tbeta}(R::Array{T,NR}, Q::Array{T,NQ}, beta::Tbeta) =
    DiscreteDP{T,NQ,NR,Tbeta,Int}(R, Q, beta)

"""
DiscreteDP type for specifying parameters for discrete dynamic programming
model State-Action Pair Formulation

##### Parameters

- `R::Array{T,NR}` : Reward Array
- `Q::Array{T,NQ}` : Transition Probability Array
- `beta::Float64`  : Discount Factor
- `s_indices::Nullable{Vector{Tind}}`: State Indices. Null unless using
  SA formulation
- `a_indices::Nullable{Vector{Tind}}`: Action Indices. Null unless using
  SA formulation
- `a_indptr::Nullable{Vector{Tind}}`: Action Index Pointers. Null unless using
  SA formulation

##### Returns

- `ddp::DiscreteDP` : Constructor for DiscreteDP object

"""
function DiscreteDP{T,NQ,NR,Tbeta,Tind}(R::AbstractArray{T,NR},
                                        Q::AbstractArray{T,NQ},
                                        beta::Tbeta, s_indices::Vector{Tind},
                                        a_indices::Vector{Tind})
    DiscreteDP{T,NQ,NR,Tbeta,Tind}(R, Q, beta, s_indices, a_indices)
end

#--------------#
#-Type Aliases-#
#--------------#

typealias DDP{T,Tbeta,Tind} DiscreteDP{T,3,2,Tbeta,Tind}
typealias DDPsa{T,Tbeta,Tind} DiscreteDP{T,2,1,Tbeta,Tind}

#--------------------#
#-Property Functions-#
#--------------------#

num_states(ddp::DDP) = size(ddp.R, 1)
num_states(ddp::DDPsa) = size(ddp.Q, 2)
num_actions(ddp::DiscreteDP) = size(ddp.R, 2)

abstract DDPAlgorithm
"""
This refers to the Value Iteration solution algorithm.

References
----------

http://quant-econ.net/jl/ddp.html

"""
immutable VFI <: DDPAlgorithm end

"""
This refers to the Policy Iteration solution algorithm.

References
----------

http://quant-econ.net/jl/ddp.html

"""
immutable PFI <: DDPAlgorithm end

"""
This refers to the Modified Policy Iteration solution algorithm.

References
----------

http://quant-econ.net/jl/ddp.html

"""
immutable MPFI <: DDPAlgorithm end

"""
DPSolveResult is an object for retaining results and associated metadata after
solving the model

##### Parameters

- `ddp::DiscreteDP` : DiscreteDP object

##### Returns

- `ddpr::DPSolveResult` : DiscreteDP Results object

"""
type DPSolveResult{Algo<:DDPAlgorithm,Tval<:Real}
    v::Vector{Tval}
    Tv::Array{Tval}
    num_iter::Int
    sigma::Array{Int,1}
    mc::MarkovChain

    function DPSolveResult(ddp::DiscreteDP)
        v = s_wise_max(ddp, ddp.R) # Initialise v with v_init
        ddpr = new(v, similar(v), 0, similar(v, Int))

        # fill in sigma with proper values
        compute_greedy!(ddp, ddpr)
        ddpr
    end

    # method to pass initial value function (skip the s-wise max)
    function DPSolveResult(ddp::DiscreteDP, v::Vector)
        ddpr = new(v, similar(v), 0, similar(v, Int))

        # fill in sigma with proper values
        compute_greedy!(ddp, ddpr)
        ddpr
    end
end

# ------------------------ #
# Bellman operator methods #
# ------------------------ #

"""
The Bellman operator, which computes and returns the updated value function Tv
for a value function v.

##### Parameters

- `ddp::DiscreteDP` : Object that contains the model parameters
- `v::Vector{T<:AbstractFloat}`: The current guess of the value function
- `Tv::Vector{T<:AbstractFloat}`: A buffer array to hold the updated value
  function. Initial value not used and will be overwritten
- `sigma::Vector`: A buffer array to hold the policy function. Initial
  values not used and will be overwritten

##### Returns

- `Tv::Vector` : Updated value function vector
- `sigma::Vector` : Updated policiy function vector
"""
function bellman_operator!(ddp::DiscreteDP, v::Vector, Tv::Vector, sigma::Vector)
    vals = ddp.R + ddp.beta * (ddp.Q * v)
    s_wise_max!(ddp, vals, Tv, sigma)
    Tv, sigma
end

"""
Apply the Bellman operator using `v=ddpr.v`, `Tv=ddpr.Tv`, and `sigma=ddpr.sigma`

##### Notes

Updates `ddpr.Tv` and `ddpr.sigma` inplace

"""
bellman_operator!(ddp::DiscreteDP, ddpr::DPSolveResult) =
    bellman_operator!(ddp, ddpr.v, ddpr.Tv, ddpr.sigma)

"""
The Bellman operator, which computes and returns the updated value function Tv
for a given value function v.

This function will fill the input `v` with `Tv` and the input `sigma` with the
corresponding policy rule

##### Parameters

- `ddp::DiscreteDP`: The ddp model
- `v::Vector{T<:AbstractFloat}`: The current guess of the value function. This
  array will be overwritten
- `sigma::Vector`: A buffer array to hold the policy function. Initial
  values not used and will be overwritten

##### Returns

- `Tv::Vector`: Updated value function vector
- `sigma::Vector{T<:Integer}`: Policy rule
"""
bellman_operator!{T<:AbstractFloat}(ddp::DiscreteDP, v::Vector{T}, sigma::Vector) =
    bellman_operator!(ddp, v, v, sigma)

# method to allow dispatch on rationals
# TODO from albep: not sure how to update this to the state-action pair formulation
function bellman_operator!{T1<:Rational,T2<:Rational,NR,NQ,T3<:Rational}(ddp::DiscreteDP{T1,NR,NQ,T2},
                                                                         v::Vector{T3},
                                                                         sigma::Vector)
    bellman_operator!(ddp, v, v, sigma)
end

"""
The Bellman operator, which computes and returns the updated value function Tv
for a given value function v.

##### Parameters

- `ddp::DiscreteDP`: The ddp model
- `v::Vector`: The current guess of the value function

##### Returns

- `Tv::Vector` : Updated value function vector
"""
bellman_operator(ddp::DiscreteDP, v::Vector) =
    s_wise_max(ddp, ddp.R + ddp.beta * (ddp.Q * v))

# ---------------------- #
# Compute greedy methods #
# ---------------------- #

"""
Compute the v-greedy policy

##### Parameters

- `ddp::DiscreteDP` : Object that contains the model parameters
- `ddpr::DPSolveResult` : Object that contains result variables

##### Returns

- `sigma::Vector{Int}` : Array containing `v`-greedy policy rule

##### Notes

modifies ddpr.sigma and ddpr.Tv in place

"""
compute_greedy!(ddp::DiscreteDP, ddpr::DPSolveResult) =
    (bellman_operator!(ddp, ddpr); ddpr.sigma)

"""
Compute the v-greedy policy.

#### Arguments

- `v::Vector` Value function vector of length `n`
- `ddp::DiscreteDP` Object that contains the model parameters

#### Returns

- `sigma:: v-greedy policy vector, of length `n`

"""
function compute_greedy{TV<:Real}(ddp::DiscreteDP, v::Vector{TV})
    Tv = similar(v)
    sigma = ones(Int, length(v))
    bellman_operator!(ddp, v, Tv, sigma)
    sigma
end

# ----------------------- #
# Evaluate policy methods #
# ----------------------- #

"""
Method of `evaluate_policy` that extracts sigma from a `DPSolveResult`

See other docstring for details
"""
evaluate_policy(ddp::DiscreteDP, ddpr::DPSolveResult) =
    evaluate_policy(ddp, ddpr.sigma)

"""
Compute the value of a policy.

##### Parameters

- `ddp::DiscreteDP` : Object that contains the model parameters
- `sigma::Vector{T<:Integer}` : Policy rule vector

##### Returns

- `v_sigma::Array{Float64}` : Value vector of `sigma`, of length n.

"""
function evaluate_policy{T<:Integer}(ddp::DiscreteDP, sigma::Vector{T})
    R_sigma, Q_sigma = RQ_sigma(ddp, sigma)
    b = R_sigma
    A = I - ddp.beta * Q_sigma
    return A \ b
end

# ------------- #
# Solve methods #
# ------------- #

"""
Solve the dynamic programming problem.

##### Parameters

- `ddp::DiscreteDP` : Object that contains the Model Parameters
- `method::Type{T<Algo}(VFI)`: Type name specifying solution method. Acceptable
arguments are `VFI` for value function iteration or `PFI` for policy function
iteration or `MPFI` for modified policy function iteration
- `;max_iter::Int(250)` : Maximum number of iterations
- `;epsilon::Float64(1e-3)` : Value for epsilon-optimality. Only used if
`method` is `VFI`
- `;k::Int(20)` : Number of iterations for partial policy evaluation in modified
policy iteration (irrelevant for other methods).

##### Returns

- `ddpr::DPSolveResult{Algo}` : Optimization result represented as a
DPSolveResult. See `DPSolveResult` for details.
"""
function solve{Algo<:DDPAlgorithm,T}(ddp::DiscreteDP{T}, method::Type{Algo}=VFI;
                                     max_iter::Integer=250, epsilon::Real=1e-3,
                                     k::Integer=20)
    ddpr = DPSolveResult{Algo,T}(ddp)
    _solve!(ddp, ddpr, max_iter, epsilon, k)
    ddpr.mc = MarkovChain(ddp, ddpr)
    ddpr
end

function solve{Algo<:DDPAlgorithm,T}(ddp::DiscreteDP{T}, v_init::Vector{T},
                                   method::Type{Algo}=VFI; max_iter::Integer=250,
                                   epsilon::Real=1e-3, k::Integer=20)
    ddpr = DPSolveResult{Algo,T}(ddp, v_init)
    _solve!(ddp, ddpr, max_iter, epsilon, k)
    ddpr.mc = MarkovChain(ddp, ddpr)
    ddpr
end

# --------- #
# Other API #
# --------- #

"""
Returns the controlled Markov chain for a given policy `sigma`.

##### Parameters

- `ddp::DiscreteDP` : Object that contains the model parameters
- `ddpr::DPSolveResult` : Object that contains result variables

##### Returns

mc : MarkovChain
     Controlled Markov chain.
"""
QuantEcon.MarkovChain(ddp::DiscreteDP, ddpr::DPSolveResult) =
    MarkovChain(RQ_sigma(ddp, ddpr)[2])

"""
Method of `RQ_sigma` that extracts sigma from a `DPSolveResult`

See other docstring for details
"""
RQ_sigma(ddp::DiscreteDP, ddpr::DPSolveResult) = RQ_sigma(ddp, ddpr.sigma)

"""
Given a policy `sigma`, return the reward vector `R_sigma` and
the transition probability matrix `Q_sigma`.

##### Parameters

- `ddp::DiscreteDP` : Object that contains the model parameters
- `sigma::Vector{Int}`: policy rule vector

##### Returns

- `R_sigma::Array{Float64}`: Reward vector for `sigma`, of length n.

- `Q_sigma::Array{Float64}`: Transition probability matrix for `sigma`,
  of shape (n, n).

"""
function RQ_sigma{T<:Integer}(ddp::DDP, sigma::Array{T})
    R_sigma = [ddp.R[i, sigma[i]] for i in 1:length(sigma)]
    Q_sigma = hcat([getindex(ddp.Q, i, sigma[i], Colon())[:] for i=1:num_states(ddp)]...)
    return R_sigma, Q_sigma'
end

# TODO: express it in a similar way as above to exploit Julia's column major order
function RQ_sigma{T<:Integer}(ddp::DDPsa, sigma::Array{T})
    sigma_indices = Array(T, num_states(ddp))
    _find_indices!(get(ddp.a_indices), get(ddp.a_indptr), sigma, sigma_indices)
    R_sigma = ddp.R[sigma_indices]
    Q_sigma = ddp.Q[sigma_indices, :]
    return R_sigma, Q_sigma
end

# ---------------- #
# Internal methods #
# ---------------- #

## s_wise_max for DDP

s_wise_max(ddp::DiscreteDP, vals::AbstractMatrix) = s_wise_max(vals)

s_wise_max!(ddp::DiscreteDP, vals::AbstractMatrix, out::Vector, out_argmax::Vector) =
    s_wise_max!(vals, out, out_argmax)

"""
Return the `Vector` `max_a vals(s, a)`,  where `vals` is represented as a
`AbstractMatrix` of size `(num_states, num_actions)`.
"""
s_wise_max(vals::AbstractMatrix) = vec(maximum(vals, 2))

"""
Populate `out` with  `max_a vals(s, a)`,  where `vals` is represented as a
`AbstractMatrix` of size `(num_states, num_actions)`.
"""
s_wise_max!(vals::AbstractMatrix, out::Vector) = (println("calling this one! "); maximum!(out, vals))

"""
Populate `out` with  `max_a vals(s, a)`,  where `vals` is represented as a
`AbstractMatrix` of size `(num_states, num_actions)`.

Also fills `out_argmax` with the column number associated with the indmax in
each row
"""
function s_wise_max!(vals::AbstractMatrix, out::Vector, out_argmax::Vector)
    # naive implementation where I just iterate over the rows
    nr, nc = size(vals)
    for i_r in 1:nr
        # reset temporaries
        cur_max = -Inf
        out_argmax[i_r] = 1

        for i_c in 1:nc
            @inbounds v_rc = vals[i_r, i_c]
            if v_rc > cur_max
                out[i_r] = v_rc
                out_argmax[i_r] = i_c
                cur_max = v_rc
            end
        end

    end
    out, out_argmax
end

## s_wise_max for DDPsa

function s_wise_max(ddp::DDPsa, vals::Vector)
    s_wise_max!(get(ddp.a_indices), get(ddp.a_indptr), vals,
                 Array(Float64, num_states(ddp)))
end

s_wise_max!(ddp::DDPsa, vals::Vector, out::Vector, out_argmax::Vector) =
    s_wise_max!(get(ddp.a_indices), get(ddp.a_indptr), vals, out, out_argmax)

"""
Populate `out` with  `max_a vals(s, a)`,  where `vals` is represented as a
`Vector` of size `(num_sa_pairs,)`.
"""
function s_wise_max!(a_indices::Vector, a_indptr::Vector, vals::Vector,
                      out::Vector)
    n = length(out)
    for i in 1:n
        if a_indptr[i] != a_indptr[i+1]
            m = a_indptr[i]
            for j in a_indptr[i]+1:(a_indptr[i+1]-1)
                if vals[j] > vals[m]
                    m = j
                end
            end
            out[i] = vals[m]
        end
    end
    return out
end

"""
Populate `out` with  `max_a vals(s, a)`,  where `vals` is represented as a
`Vector` of size `(num_sa_pairs,)`.

Also fills `out_argmax` with the cartesiean index associated with the indmax in
each row
"""
function s_wise_max!(a_indices::Vector, a_indptr::Vector, vals::Vector,
                      out::Vector, out_argmax::Vector)
    n = length(out)
    for i in 1:n
        if a_indptr[i] != a_indptr[i+1]
            m = a_indptr[i]
            for j in a_indptr[i]+1:(a_indptr[i+1]-1)
                if vals[j] > vals[m]
                    m = j
                end
            end
            out[i] = vals[m]
            out_argmax[i] = a_indices[m]
        end
    end
    out, out_argmax
end

"""
Check whether `s_indices` and `a_indices` are sorted in lexicographic order.

Parameters
----------
s_indices, a_indices : Vectors

Returns
-------
bool: Whether `s_indices` and `a_indices` are sorted.
"""
function _has_sorted_sa_indices(s_indices::Vector, a_indices::Vector)
    L = length(s_indices)
    for i in 1:L-1
        if s_indices[i] > s_indices[i+1]
            return false
        end
        if s_indices[i] == s_indices[i+1]
            if a_indices[i] >= a_indices[i+1]
                return false
            end
        end
    end
    return true
end

"""
Generate `a_indptr`; stored in `out`. `s_indices` is assumed to be
in sorted order.

Parameters
----------
num_states : Int

s_indices : Vector{Int}

out : Vector{Int} with length = num_states+1
"""
function _generate_a_indptr!(num_states::Int, s_indices::Vector, out::Vector)
    idx = 1
    out[1] = 1
    for s in 1:num_states-1
        while(s_indices[idx] == s)
            idx += 1
        end
        out[s+1] = idx
    end
    # need this +1 to be consistent with Julia's sparse pointers:
    # colptr[i]:(colptr[i+1]-1)
    out[num_states+1] = length(s_indices)+1
    out
end

function _find_indices!(a_indices::Vector, a_indptr::Vector, sigma::Vector,
                        out::Vector)
    n = length(sigma)
    for i in 1:n, j in a_indptr[i]:(a_indptr[i+1]-1)
        if sigma[i] == a_indices[j]
            out[i] = j
        end
    end
end

"""
Define Matrix Multiplication between 3-dimensional matrix and a vector

Matrix multiplication over the last dimension of A

"""
function *{T}(A::Array{T,3}, v::Vector)
    shape = size(A)
    size(v, 1) == shape[end] || error("wrong dimensions")

    B = reshape(A, (prod(shape[1:end-1]), shape[end]))
    out = B * v

    return reshape(out, shape[1:end-1])
end

"""
Impliments Value Iteration
NOTE: See `solve` for further details
"""
function _solve!(ddp::DiscreteDP, ddpr::DPSolveResult{VFI}, max_iter::Integer,
                 epsilon::Real, k::Integer)
    if ddp.beta == 0.0
        tol = Inf
    else
        tol = epsilon * (1-ddp.beta) / (2*ddp.beta)
    end

    for i in 1:max_iter
        # updates Tv in place
        bellman_operator!(ddp, ddpr)

        # compute error and update the v inside ddpr
        err = maxabs(ddpr.Tv .- ddpr.v)
        copy!(ddpr.v, ddpr.Tv)
        ddpr.num_iter += 1

        if err < tol
            break
        end
    end

    ddpr
end

"""
Policy Function Iteration

NOTE: The epsilon is ignored in this method. It is only here so dispatch can
      go from `solve(::DiscreteDP, ::Type{Algo})` to any of the algorithms.
      See `solve` for further details
"""
function _solve!(ddp::DiscreteDP, ddpr::DPSolveResult{PFI}, max_iter::Integer,
                 epsilon::Real, k::Integer)
    old_sigma = copy(ddpr.sigma)

    for i in 1:max_iter
       ddpr.v = evaluate_policy(ddp, ddpr)
       compute_greedy!(ddp, ddpr)

       ddpr.num_iter += 1
       if all(old_sigma .== ddpr.sigma)
           break
       end
       copy!(old_sigma, ddpr.sigma)

    end

    ddpr
end

span(x::Vector) = maximum(x) - minimum(x)
midrange(x::Vector) = mean(extrema(x))

"""
Modified Policy Function Iteration
"""
function _solve!(ddp::DiscreteDP, ddpr::DPSolveResult{MPFI}, max_iter::Integer,
                 epsilon::Real, k::Integer)
    beta = ddp.beta
    fill!(ddpr.v, minimum(ddp.R[ddp.R .> -Inf]) / (1.0 - beta))
    old_sigma = copy(ddpr.sigma)

    tol = beta > 0 ? epsilon * (1-beta) / beta : Inf

    for i in 1:max_iter
        bellman_operator!(ddp, ddpr)  # updates Tv, sigma inplace
        dif = ddpr.Tv - ddpr.v

        ddpr.num_iter += 1

        # check convergence
        if span(dif) < tol
            ddpr.v = ddpr.Tv + midrange(dif) * beta / (1-beta)
            break
        end

        # now update v to use the output of the bellman step when entering
        # policy loop
        copy!(ddpr.v, ddpr.Tv)

        # now do k iterations of policy iteration
        R_sigma, Q_sigma = RQ_sigma(ddp, ddpr)
        for i in 1:k
            ddpr.Tv = R_sigma + beta * Q_sigma * ddpr.v
            copy!(ddpr.v, ddpr.Tv)
        end

    end

    ddpr
end


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--generate_data"
            arg_type = Bool
            default = false
        "--version"
            arg_type = String
            default = "0000"
        "--verbose"
            arg_type = Bool
            default = false
    end

    return parse_args(s)
end

args = parse_commandline()
if !isdir("./$(args["version"])")
    mkdir("./$(args["version"])")
end

using PyTensorFlow
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)
include("ops_imseq.jl")
include("../Ops/FWI/fwi_util.jl")
include("fwi_util_op.jl")
np = pyimport("numpy")
# NOTE  Parameters
# const ALPHA = 0.006323996017182
# const SRC_CONST = 5.6146
# const GRAV_CONST = 1.0/144.0
const ALPHA = 1.0
const SRC_CONST = 86400.0
const GRAV_CONST = 1.0

# NOTE Hyperparameter for flow simulation
m = 15
n = 30
h = 30.0 # meter

NT  = 50
dt_survey = 5
Δt = 20.0 # day
z = (1:m)*h|>collect
x = (1:n)*h|>collect
X, Z = np.meshgrid(x, z)

# ρw = 996.9571
# ρo = 640.7385
# μw = 1.0
# μo = 3.0
ρw = 501.9
ρo = 1053.0
μw = 0.1
μo = 1.0

K_init = 20.0 .* ones(m,n)

g = 9.8*GRAV_CONST
ϕ = 0.25 .* ones(m,n)
qw = zeros(NT, m, n)
qw[:,9,3] .= 0.005 * (1/h^2)/10.0 * SRC_CONST
qo = zeros(NT, m, n)
qo[:,9,28] .= -0.005 * (1/h^2)/10.0 * SRC_CONST
sw0 = zeros(m, n)
survey_indices = collect(1:dt_survey:NT+1) # 10 stages
n_survey = length(survey_indices)

# NOTE Hyperparameter for fwi_op

# argsparse.jl
# ENV["CUDA_VISIBLE_DEVICES"] = 1
# ENV["PARAMDIR"] = "Src/params/"
# config = tf.ConfigProto(device_count = Dict("GPU"=>0))

dz = 3 # meters
dx = 3
nz = Int64(round((m * h) / dz)) + 1
nx = Int64(round((n * h) / dx)) + 1
nPml = 64
nSteps = 3001
dt = 0.00025
f0 = 50.0
nPad = 32 - mod((nz+2*nPml), 32)
nz_pad = nz + 2*nPml + nPad
nx_pad = nx + 2*nPml

# reflection
# x_src = collect(5:20:nx-5)
# z_src = 5ones(Int64, size(x_src))
# x_rec = collect(5:1:nx-5)
# z_rec = 5 .* ones(Int64, size(x_rec))

# xwell
# # z_src = collect(5:10:nz-5) #14->11srcs 10->15srcs
# # z_src = collect(5:10:nz-5)
z_src = collect(5:10:nz-5)
x_src = 5ones(Int64, size(z_src))
z_rec = collect(5:1:nz-5)
x_rec = (nx-5) .* ones(Int64, size(z_rec))

# para_fname = "./$(args["version"])/para_file.json"
# survey_fname = "./$(args["version"])/survey_file.json"
# data_dir_name = "./$(args["version"])/Data"
# paraGen(nz, nx, dz, dx, nSteps, dt, f0, nPml, nPad, filter_para, isAc, para_fname, survey_fname, data_dir_name)
# surveyGen(z_src, x_src, z_rec, x_rec, survey_fname)

cp_nopad = 3500.0 .* ones(nz, nx) # initial cp
cs = cp_nopad ./ sqrt(3.0)
den = 2200.0 .* ones(nz, nx)
cp_pad = 3500.0 .* ones(nz_pad, nx_pad) # initial cp
cs_pad = cp_pad ./ sqrt(3.0)
den_pad = 2200.0 .* ones(nz_pad, nx_pad)
cp_pad_value = 3500.0

# tf_cp = constant(cp)
tf_cs = constant(cs_pad)
tf_den = constant(den_pad)

# src = Matrix{Float64}(undef, 1, 2001)
# # src[1,:] = Float64.(reinterpret(Float32, read("../Ops/FWI/Src/params/ricker_10Hz.bin")))
# src[1,:] = Float64.(reinterpret(Float32, read("../Ops/FWI/Src/params/Mar_source_2001.bin")))
src = sourceGene(f0, nSteps, dt)
tf_stf = constant(repeat(src, outer=length(z_src)))
# tf_para_fname = tf.strings.join([para_fname])
tf_gpu_id0 = constant(0, dtype=Int32)
tf_gpu_id1 = constant(1, dtype=Int32)
nGpus = 1
# tf_gpu_id_array = constant(collect(0:nGpus-1), dtype=Int32)
tf_gpu_id_array = constant([0], dtype=Int32)
tf_shot_ids0 = constant(collect(Int32, 0:length(x_src)-1), dtype=Int32)
tf_shot_ids1 = constant(collect(Int32, 13:25), dtype=Int32)

# NOTE Hyperparameter for rock physics
tf_bulk_fl1 = constant(2.735e9)
tf_bulk_fl2 = constant(0.125e9) # to displace fl1
tf_bulk_sat1 = constant(den .* (cp_nopad.^2 .- 4.0/3.0 .* cp_nopad.^2 ./3.0)) # vp/vs ratio as sqrt(3)
tf_bulk_min = constant(36.6e9)
tf_shear_sat1 = constant(den .* cp_nopad.^2 ./3.0)
tf_ϕ_pad = tf.image.resize_bilinear(tf.reshape(constant(ϕ), (1, m, n, 1)), (nz, nx)) # upsample the porosity
tf_ϕ_pad = cast(tf_ϕ_pad, Float64)
tf_ϕ_pad = squeeze(tf_ϕ_pad)
tf_shear_pad = tf.pad(tf_shear_sat1, [nPml (nPml+nPad); nPml nPml],
    constant_values=den[1,1] * cp_nopad[1,1]^2 /3.0) / 1e6

function Gassman(sw)
    tf_bulk_fl_mix = 1.0/( (1-sw)/tf_bulk_fl1 + sw/tf_bulk_fl2 )
    temp = tf_bulk_sat1/(tf_bulk_min - tf_bulk_sat1) - tf_bulk_fl1/tf_ϕ_pad /(tf_bulk_min - tf_bulk_fl1) + tf_bulk_fl_mix/tf_ϕ_pad /(tf_bulk_min - tf_bulk_fl_mix)

    tf_bulk_new = tf_bulk_min / (1.0/temp + 1.0)
    # tf_den_new = constant(den) + tf_ϕ_pad .* sw * (ρw - ρo) *16.018463373960138;
    tf_den_new = constant(den) + tf_ϕ_pad .* sw * (ρw - ρo)
    # tf_cp_new = sqrt((tf_bulk_new + 4.0/3.0 * tf_shear_sat1)/tf_den_new)
    tf_lambda_new = tf_bulk_new - 2.0/3.0 * tf_shear_sat1
    return tf_lambda_new, tf_den_new
end

tf_brie_coef = Variable(2.0*30.0)
# tf_brie_coef = constant(3.0)
# tf_brie_coef = constant(2.0)
function Brie(sw)
    tf_bulk_fl_mix = (tf_bulk_fl1-tf_bulk_fl2)*(1-sw)^(tf_brie_coef/30.0) + tf_bulk_fl2
    temp = tf_bulk_sat1/(tf_bulk_min - tf_bulk_sat1) - tf_bulk_fl1/tf_ϕ_pad /(tf_bulk_min - tf_bulk_fl1) + tf_bulk_fl_mix/tf_ϕ_pad /(tf_bulk_min - tf_bulk_fl_mix)

    tf_bulk_new = tf_bulk_min / (1.0/temp + 1.0)
    # tf_den_new = constant(den) + tf_ϕ_pad .* sw * (ρw - ρo) *16.018463373960138;
    tf_den_new = constant(den) + tf_ϕ_pad .* sw * (ρw - ρo)
    # tf_cp_new = sqrt((tf_bulk_new + 4.0/3.0 * tf_shear_sat1)/tf_den_new)
    tf_lambda_new = tf_bulk_new - 2.0/3.0 * tf_shear_sat1
    return tf_lambda_new, tf_den_new
end

function RockLinear(sw)
    # tf_lambda_new = constant(7500.0*1e6 .* ones(nz,nx)) + (17400.0-7500.0)*1e6 * sw
    tf_lambda_new = constant(7500.0*1e6 .* ones(nz,nx)) + (9200.0-7500.0)*1e6 * sw
    tf_den_new = constant(den) + tf_ϕ_pad .* sw * (ρw - ρo)
    return tf_lambda_new, tf_den_new
end

tf_patch_temp = tf_bulk_sat1/(tf_bulk_min - tf_bulk_sat1) -
    tf_bulk_fl1/tf_ϕ_pad /(tf_bulk_min - tf_bulk_fl1) +
    tf_bulk_fl2/tf_ϕ_pad /(tf_bulk_min - tf_bulk_fl2)
tf_bulk_sat2 = tf_bulk_min/(1.0/tf_patch_temp + 1.0)
function Patchy(sw)
    tf_bulk_new = 1/( (1-sw)/(tf_bulk_sat1+4.0/3.0*tf_shear_sat1)
            + sw/(tf_bulk_sat2+4.0/3.0*tf_shear_sat1) ) - 4.0/3.0*tf_shear_sat1
    tf_lambda_new = tf_bulk_new - 2.0/3.0 * tf_shear_sat1
    tf_den_new = constant(den) + tf_ϕ_pad .* sw * (ρw - ρo)
    return tf_lambda_new, tf_den_new
end


text:
module MMA

include("math.jl")

end # module


text:
"write 1 byte as `UInt8` to `io`"
write_uint1(io::IO, arg) = write(io, UInt8(arg))

"write 2 bytes as `UInt16` to `io`"
write_uint2(io::IO, arg) = write(io, UInt16(arg))

"write 4 bytes as `Int32` to `io`"
write_int4(io::IO, arg) = write(io, Int32(arg))

"write 4 bytes as `UInt32` to `io`"
write_uint4(io::IO, arg) = write(io, UInt32(arg))

"write 8 bytes as `UInt64` to `io`"
write_uint8(io::IO, arg) = write(io, UInt64(arg))

@doc raw"""
write unicode string as plain ascii string to `io`

# Examples
```jldoctest
julia> sprint(Pickle.write_plain_str, "abcα😀")
"abc\\u03b1\\U0001f600"
```
"""
function write_plain_str(io::IO, arg)
  len = 0
  for c in arg
    if isascii(c)
      if c == '\\'
        len += write(io, "\\u005c")
      elseif c == '\0'
        len += write(io, "\\u0000")
      elseif c == '\n'
        len += write(io, "\\u000a")
      elseif c == '\r'
        len += write(io, "\\u000d")
      elseif c == '\u1a'
        len += write(io, "\\u001a")
      else
        len += write(io, c)
      end
    else
      if codepoint(c) >> 16 |> !iszero # large unicode
        len += write(io, "\\U")
        len += write(io, lpad(string(codepoint(c); base=16), 8, '0'))
      else
        len += write(io, "\\u")
        len += write(io, lpad(string(codepoint(c); base=16), 4, '0'))
      end
    end
  end
  len
end

"""
convert a Integer into 2's complement byte array (Vector{UInt8}).

# Examples
```jldoctest
julia> Pickle.int_to_bytes(0) |> isempty # 0-element Vector{UInt8}
true

julia> Pickle.int_to_bytes(255)
2-element Vector{UInt8}:
 0xff
 0x00

julia> Pickle.int_to_bytes(32767)
2-element Vector{UInt8}:
 0xff
 0x7f

julia> Pickle.int_to_bytes(-256)
2-element Vector{UInt8}:
 0x00
 0xff

julia> Pickle.int_to_bytes(-32768)
2-element Vector{UInt8}:
 0x00
 0x80

julia> Pickle.int_to_bytes(-128)
1-element Vector{UInt8}:
 0x80

julia> Pickle.int_to_bytes(127)
1-element Vector{UInt8}:
 0x7f
```
"""
function int_to_bytes(x)
  iszero(x) && return Vector{UInt8}()
  islittle = Base.ENDIAN_BOM == 0x04030201
  isneg = x < 0

  if isneg
    x += 1
  end

  hexstr = x > 0 ? string(x; base=16) : string(-x; base=16)
  if length(hexstr) |> isodd
    hexstr = lpad(hexstr, length(hexstr)+1, '0')
  else
    if !(first(hexstr) < '8')
      hexstr = lpad(hexstr, length(hexstr)+2, '0')
    end
  end

  bytes = hex2bytes(hexstr)
  islittle && (bytes = reverse(bytes))

  if isneg
    bytes .⊻= 0xff
  end

  return bytes
end


text:
@testset "Core Tests" begin
  include("core-test.jl")
  @test_broken "Interface Tests"

  @test_skip "GenericNeurons"
  @test_skip "GenericSynapses"
  @test_skip "ConvNeurons"
  @test_skip "ConvSynapses"
end


text:
using Test

@testset "Conjugate" begin

    using ProximalOperators
    using ProximalOperators: Zero
    using ProximalAlgorithms

    x = [1.0, -2.0, 3.0, -4.0, 5.0, -6.0]

    f = Conjugate(IndZero()) # = IndFree
    grad_f_x, f_x = gradient(f, x)
    @test iszero(grad_f_x)
    @test iszero(f_x)

    g = Conjugate(Zero()) # = IndZero
    prox_g_x, g_y = prox(g, x)
    @test iszero(prox_g_x)
    @test iszero(g_y)

    l = Conjugate(SqrNormL2())
    grad_l_x, l_x = gradient(l, x)
    @test isequal(grad_l_x, x)

    lam = 2
    l = Conjugate(SqrNormL2(lam))
    grad_l_x, l_x = gradient(l, x)
    @test isequal(grad_l_x, x / lam)

end