text:
using StatsBase, Distributions, IterTools #StatsPlots

function get_freq(h::Histogram, val)
    x = searchsortedfirst.(h.edges, val)
    h.weights[x...]
end

function get_bin(h::Histogram, val)
    bin = searchsortedfirst.(h.edges, val)
    bin[1] - 1
end

function get_probability(h::Histogram, val)
    x = searchsortedfirst.(h.edges, val)
    h.weights[x...]/sum(h.weights)
end

function mutual_information(X::Vector{Vector{T}}, Y::Vector{Vector{T}}) where T <: Number
    hx = fit(Histogram, (X...,))
    hy = fit(Histogram, (Y...,))
    XY = vcat(X,Y)
    hxy = fit(Histogram, (XY...,))
    MI = 0.0
    rx = create_matrix_range(size(hx.weights))
    ry = create_matrix_range(size(hy.weights))
    for x in Iterators.product(rx...)
        p_x = get_bin_probability(hx, x...)
        for y in Iterators.product(ry...)
            p_y = get_bin_probability(hy, y...)
            index_xy = (x..., y...)
            p_xy = get_bin_probability(hxy, index_xy...)
            if p_xy != 0
                MI += p_xy*log(p_xy/(p_x*p_y))
            end
        end
    end
    MI
end

function mutual_information(p_x, p_y, p_xy)
    MI = 0.0
    for x = 1:length(p_x)
        for y = 1:length(p_y)
            prob_xy = p_xy[x, y]
            prob_log = prob_xy/(p_x[x]*p_y[y])
            if prob_xy != 0 && prob_log != NaN
                MI += prob_xy*log(prob_log)
            end
        end
    end
    MI
end

function mutual_information_unit(X::Vector{T}, Y::Vector{T}) where T <: Int
    max_X, max_Y = maximum(X), maximum(Y)
    hx = fit(Histogram, (X), 0:max_X, closed=:right)
    hy = fit(Histogram, (Y), 0:max_Y, closed=:right)
    XY = hcat([X],[Y])
    hxy = fit(Histogram, (X, Y), (0:max_X, 0:max_Y), closed=:right)
    MI = 0.0
    rx = create_matrix_range(size(hx.weights))
    ry = create_matrix_range(size(hy.weights))
    for x in Iterators.product(rx...)
        p_x = get_bin_probability(hx, x...)
        for y in Iterators.product(ry...)
            p_y = get_bin_probability(hy, y...)
            index_xy = (x..., y...)
            p_xy = get_bin_probability(hxy, index_xy...)
            if p_xy != 0
                MI += p_xy*log(p_xy/(p_x*p_y))
            end
        end
    end
    MI
end

function mutual_information_unit(X::Vector{T}, Y::Vector{T}, lim_X, lim_Y) where T <: Int
    min_X, max_X = lim_X
    min_Y, max_Y = lim_Y
    range_X = min_X:max_X
    range_Y = min_Y:max_Y
    hx = fit(Histogram, (X), range_X, closed=:right)
    hy = fit(Histogram, (Y), range_Y, closed=:right)
    XY = hcat([X],[Y])
    hxy = fit(Histogram, (X, Y), (range_X, range_Y), closed=:right)
    MI = 0.0
    rx = create_matrix_range(size(hx.weights))
    ry = create_matrix_range(size(hy.weights))
    for x in Iterators.product(rx...)
        p_x = get_bin_probability(hx, x...)
        for y in Iterators.product(ry...)
            p_y = get_bin_probability(hy, y...)
            index_xy = (x..., y...)
            p_xy = get_bin_probability(hxy, index_xy...)
            if p_xy != 0
                MI += p_xy*log(p_xy/(p_x*p_y))
            end
        end
    end
    MI
end

get_bin_probability(h::Histogram, x...) = h.weights[x...]/sum(h.weights)

function create_matrix_range(t::Tuple)
    tuple_length = length(t)
    ranges = []
    for r = 1:tuple_length
        range = 1:t[r]
        push!(ranges, range)
    end
    ranges
end

function calculate_joint_probability(X::Vector{<:Number}, Y::Vector{<:Number}, lim_X, lim_Y)
    min_X, max_X = lim_X
    min_Y, max_Y = lim_Y
    range_X = min_X:max_X
    range_Y = min_Y:max_Y
    p_x = fit(Histogram, (X), range_X, closed=:right)
    p_y = fit(Histogram, (Y), range_Y, closed=:right)
    p_xy = fit(Histogram, (X, Y), (range_X, range_Y), closed=:right)
    p_X = p_x.weights/sum(p_x.weights)
    p_Y = p_y.weights/sum(p_y.weights)
    p_XY = p_xy.weights/sum(p_xy.weights)
    return p_X, p_Y, p_XY
end

function calculate_joint_probability(X::Vector{<:Number}, Y::Vector{<:Number}, nbins)
    p_x = fit(Histogram, (X), nbins=nbins[1])
    p_y = fit(Histogram, (Y), nbins=nbins[2])
    p_xy = fit(Histogram, (X, Y), nbins=nbins)
    p_X = p_x.weights/sum(p_x.weights)
    p_Y = p_y.weights/sum(p_y.weights)
    p_XY = p_xy.weights/sum(p_xy.weights)
    return p_X, p_Y, p_XY
end

include("layer_histogram.jl")

function calculate_joint_probability(X::Vector{<:Number}, Y::Vector{<:Number}, limits_X, limits_Y, nbins)
    p_x = create_bin_histogram(X, nbins[1], limits_X[1], limits_X[2])
    p_y = create_bin_histogram(Y, nbins[2], limits_Y[1], limits_Y[2])
    p_xy = fit(Histogram, (X, Y), (p_x.edges[1], p_y.edges[1]))
    p_X = p_x.weights/sum(p_x.weights)
    p_Y = p_y.weights/sum(p_y.weights)
    p_XY = p_xy.weights/sum(p_xy.weights)
    return p_X, p_Y, p_XY
end

# Test
#=
X = randn(10000)
Y = randn(10000)

println("Mutual Information: I(X; Y) = ", mutual_information([X], [Y]))
println("Mutual Information: I(Y; X) = ", mutual_information([Y], [X]))
println("Mutual Information: I(X; (Y,X)) = ", mutual_information([X], [Y, X]))
println("Mutual Information: I(X; X) = ", mutual_information([X], [X]))
println("Mutual Information: I(Y; Y) = ", mutual_information([Y], [Y]))
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
const DEFAULT_QUANTIZATION_SAMPLING_RATIO = 1.0
const DEFAULT_QUANTIZAION_K = 2
const DEFAULT_QUANTIZATION_M = 1
const DEFAULT_QUANTIZATION_METHOD = :pq
const DEFAULT_QUANTIZATION_DISTANCE = Distances.SqEuclidean()

ext: jl

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DelayEmbeddings
using DelimitedFiles
using Random

include("../../src/pecuzal_method.jl")
include("../../src/data_analysis_functions.jl")

# compute time series
Random.seed!(1234)

u0_1 = randn(3)
AR3 = ar_process3(u0_1, .45, 0.3, .1, 1.0, 40000)
#AR3 = ar_process3(u0_1, .45, 0.1, .4, 1.0, 40000)

u0_2 = randn(5)
AR5 = ar_process5(u0_2, .2, 0.1, 0.05, .3, 0.2, 1.0, 40000)

u0_3 = randn(7)
AR7 = ar_process7(u0_3, .2, 0.04, 0.0, 0.05, .1, 0.2, .3, 1.0, 40000)
#AR7 = ar_process7(u0_3, .1, 0.04, 0.0, 0.05, .3, 0.2, .3, 1.0, 40000)

figure()
subplot(311)
plot(AR3)
subplot(312)
plot(AR5)
subplot(313)
plot(AR7)

ts_lengths = [500, 1000, 3000, 6000, 10000, 15000, 20000, 30000, 40000]

dim_AR3 = zeros(length(ts_lengths))
dim_AR5 = zeros(length(ts_lengths))
dim_AR7 = zeros(length(ts_lengths))
Ls_3 = []
Ls_5 = []
Ls_7 = []
τ_vals_3 = []
τ_vals_5 = []
τ_vals_7 = []

for (idx,ts_length) in enumerate(ts_lengths)
    println("TS-length: $ts_length")
    w3 = estimate_delay(AR3[1:ts_length],"mi_min")
    w5 = estimate_delay(AR5[1:ts_length],"mi_min")
    w7 = estimate_delay(AR7[1:ts_length],"mi_min")

    Y3, τ3, _, Ls_3_, _ = pecuzal_embedding_update(AR3[1:ts_length]; τs = 0:50 , w = w3)
    Y5, τ5, _, Ls_5_, _ = pecuzal_embedding_update(AR5[1:ts_length]; τs = 0:50 , w = w3)
    Y7, τ7, _, Ls_7_, _ = pecuzal_embedding_update(AR7[1:ts_length]; τs = 0:50 , w = w7)

    push!(Ls_3, Ls_3_)
    push!(Ls_5, Ls_5_)
    push!(Ls_7, Ls_7_)
    push!(τ_vals_3, τ3)
    push!(τ_vals_5, τ5)
    push!(τ_vals_7, τ7)
    dim_AR3[idx] = size(Y3,2)
    dim_AR5[idx] = size(Y5,2)
    dim_AR7[idx] = size(Y7,2)
    if idx == 5
        break
    end
end

# writedlm("./scripts/AR Processes/AR3_dims.csv", dim_AR3)
# writedlm("./scripts/AR Processes/AR5_dims.csv", dim_AR5)
# writedlm("./scripts/AR Processes/AR7_dims.csv", dim_AR7)


text:
"""
    SimpleEdgeIter

The function [`edges`](@ref) returns a `SimpleEdgeIter` for `AbstractSimpleGraph`s.
The iterates are in lexicographical order, smallest first. The iterator is valid for
one pass over the edges, and is invalidated by changes to the graph.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = path_graph(3);

julia> es = edges(g)
SimpleEdgeIter 2

julia> e_it = iterate(es)
(Edge 1 => 2, SimpleEdgeIterState [2, 2])

julia> iterate(es, e_it[2])
(Edge 2 => 3, SimpleEdgeIterState [0, 1])
```
"""
struct SimpleEdgeIter{G} <: AbstractEdgeIter
    g::G
end

eltype(::Type{SimpleEdgeIter{SimpleGraph{T}}}) where {T} = SimpleGraphEdge{T}
eltype(::Type{SimpleEdgeIter{SimpleDiGraph{T}}}) where {T} = SimpleDiGraphEdge{T}

@traitfn @inline function iterate(eit::SimpleEdgeIter{G}, state=(one(eltype(eit.g)), 1) ) where {G <: AbstractSimpleGraph; !IsDirected{G}}
    g = eit.g
    fadjlist = fadj(g)
    T = eltype(g)
    n = T(nv(g))
    u, i = state

    @inbounds while u < n
        list_u = fadjlist[u]
        if i > length(list_u)
            u += one(u)
            i = searchsortedfirst(fadjlist[u], u)
            continue
        end
        e = SimpleEdge(u, list_u[i])
        state = (u, i + 1)
        return e, state
    end

    @inbounds (n == 0 || i > length(fadjlist[n])) && return nothing

    e = SimpleEdge(n, n)
    state = (u, i + 1)
    return e, state
end

@traitfn @inline function iterate(eit::SimpleEdgeIter{G}, state=(one(eltype(eit.g)), 1) ) where {G <: AbstractSimpleGraph; IsDirected{G}}
    g = eit.g
    fadjlist = fadj(g)
    T = eltype(g)
    n = T(nv(g))
    u, i = state

    n == 0 && return nothing

    @inbounds while true
        list_u = fadjlist[u]
        if i > length(list_u)
            u == n && return nothing

            u += one(u)
            list_u = fadjlist[u]
            i = 1
            continue
        end
        e = SimpleEdge(u, list_u[i])
        state = (u, i + 1)
        return e, state
    end

    return nothing
end

length(eit::SimpleEdgeIter) = ne(eit.g)

function _isequal(e1::SimpleEdgeIter, e2)
    k = 0
    for e in e2
        has_edge(e1.g, e) || return false
        k += 1
    end
    return k == ne(e1.g)
end
==(e1::SimpleEdgeIter, e2::AbstractVector{SimpleEdge}) = _isequal(e1, e2)
==(e1::AbstractVector{SimpleEdge}, e2::SimpleEdgeIter) = _isequal(e2, e1)
==(e1::SimpleEdgeIter, e2::Set{SimpleEdge}) = _isequal(e1, e2)
==(e1::Set{SimpleEdge}, e2::SimpleEdgeIter) = _isequal(e2, e1)

function ==(e1::SimpleEdgeIter, e2::SimpleEdgeIter)
    g = e1.g
    h = e2.g
    ne(g) == ne(h) || return false
    m = min(nv(g), nv(h))
    for i in 1:m
        fadj(g, i) == fadj(h, i) || return false
    end
    nv(g) == nv(h) && return true
    for i in m+1:nv(g)
        isempty(fadj(g, i)) || return false
    end
    for i in m+1:nv(h)
        isempty(fadj(h, i)) || return false
    end
    return true
end

in(e, es::SimpleEdgeIter) = has_edge(es.g, e)

show(io::IO, eit::SimpleEdgeIter) = write(io, "SimpleEdgeIter $(ne(eit.g))")


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 0782768c-951f-4189-954d-cba8b401314d
begin

	# load some external packages
	using PlutoUI
	using DataFrames
	using BSON
	using GLPK
	using PrettyTables
	using Plots
	using CSV
	using Optim
	using Statistics
	using LinearAlgebra
	using Plots
	using TypedTables

	# setup my paths (where are my files?)
	_PATH_TO_ROOT = pwd()
	_PATH_TO_SRC = joinpath(_PATH_TO_ROOT,"src")
	_PATH_TO_MODEL = joinpath(_PATH_TO_ROOT,"model")
	_PATH_TO_FIGS = joinpath(_PATH_TO_ROOT,"figs")

	# load the ENGRI 1120 project code library -
	include(joinpath(_PATH_TO_SRC,"Include.jl"))

	# load the model -
	MODEL = BSON.load(joinpath(_PATH_TO_MODEL,"model_v2.bson"), @__MODULE__)

	# show -
	nothing
end

# ╔═╡ 4855467c-3670-4e0b-a64e-0e09effa6e0d
md"""
## ENGRI 1120: Design and Analysis of a Sustainable Cell-Free Production Process for Industrially Important Small Molecules
"""

# ╔═╡ 8ca99e1f-afd6-4c1c-8918-db6d9747099c
html"""

Team name: Super Swag Super Heated Steam Engines

Team Members: Lauren de Silva, Dennis Wu, Henry Lin

Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca NY 14850

"""

# ╔═╡ a0ad3474-1844-41bc-bd95-242aa94a5ff1
md"""
### Introduction
"""

# ╔═╡ 27d0013a-3c23-4f92-a334-5e4b55e2447e
md"""PDGN (1,2-propanediol dinitrate or propane glycol dinitrate) is a compound commonly found in propellants and explosives. It is an nitrate ester high explosive that is found both in the liquid propellant and the liquid explosive form [5]. The US navy is a main consumer of Otto Fuel II, a torpedo fuel that comprises approximately 76% PGDN[2] . In this work, we analyze how to manufacture PDGN at 95% purity and stream rate of at least 1.0g/hour in an E. Coli cell-free bioreactor from alpha-d-glucose, a renewable sugar feedstock, in addition to oxygen and potassium nitrate.

Through this glucose based method, which is accessible all around the world,  we aim to propose a more economically feasible method of producing PDGN. Using oxygen from the air and potassium nitrate, we can develop a relatively cost efficient process that is inexpensive in comparison to current production methods.

Although there is no direct route from sucrose to PGDN, the bioreactor allows us to take advantage of the glycolysis of sucrose in the bioreactor. This reaction mechanism follows the follow chemical equations:

* (1) ATP + alpha-D-Glucose ⇌ ADP + alpha-D-Glucose 6-phosphate-

* (2)Alpha-D-Glucose 6-phosphate ⇌ beta-D-Fructose 6-phosphate

* (3)beta-D-Fructose 1,6-bisphosphate ⇌ Glycerone phosphate + D-Glyceraldehyde 3-phosphate

* (4)D-Glyceraldehyde 3-phosphate ⇌ Glycerone phosphate

* (5)Glycerone phosphate ⇌ Methylglyoxal + Orthophosphate

* (6)Methylglyoxal + NADPH + H+ ⇌ (S)-Lactaldehyde + NADP+

* (7)(S)-Lactaldehyde + NADH + H+ ⇌ Propane-1,2-diol + NAD+

* (8)Propane-1,2-diol + 2KNO3 ⇌ PDGN + 2KOH

The only addition to this biological reaction is potassium nitrate in order to produce the PDGN. The overall reaction follows the equation:

$$glucose + 2KNO_3 ⇌ PDGN + 2KOH + CO_2 + excess-materials$$

Our goal is to determine the most economically efficient pathway to synthesize the PGDN specifically looking at the sucrose reaction mechanism above.
"""

# ╔═╡ 40da982c-1cc4-4881-a2ea-fbeef5c46d2d
md"""
### Materials and Methods
"""

# ╔═╡ 446bd227-6579-4c87-8588-0de111c5a529
md""" ##### Flux Balance Analysis Assumptions and Methods
The primary method used to compute the optimal open extent of reaction was a flux balance analysis. The flux balance is a practice of optimizing the genetic processes within the cell to increase the production of what we need. This form of metabolic engineering is modeled in the flux balance code used below. The program uses a linear equation as follows:
max vi=sum c1v1
The equation seeks to take the metabolic flux of the objective coefficient and subject it to two constraints of the metabolic flux to find the optimal flux. To use the flux balance analysis we make the simplifying assumptions that the microfluidic chip used is well-mixed and operates at steady-state. Also that the chip is being run at a constant temperature and pressure and that the liquid phase is ideal.

The program reads that the optimal reaction extent is 2.466 mmol/hr. This is our maximum flow of PGDN we can create from just one chip. After this is established we move onto setting up a series of these reactors to maximize the flow of our stream.
"""

# ╔═╡ 4867b51c-fb6c-42a9-b87d-4131e014b402
md"""
##### Setting up the Flux Balance Analysis (FBA) calculation to estimate of optimal reaction rates
"""

# ╔═╡ 059e3b4e-4e84-4934-b787-32d0f42a0247
begin

	# setup the FBA calculation for the project -

	# === SELECT YOUR PRODUCT HERE ==================================================== #
	# What rate are trying to maximize? (select your product)
	# rn:R08199 = isoprene
	# rn:28235c0c-ec00-4a11-8acb-510b0f2e2687 = PGDN
	# rn:rn:R09799 = Hydrazine
	# rn:R03119 = 3G
	idx_target_rate = find_reaction_index(MODEL,:reaction_number=>"rn:28235c0c-ec00-4a11-8acb-510b0f2e2687")
	# ================================================================================= #

	# First, let's build the stoichiometric matrix from the model object -
	(cia,ria,S) = build_stoichiometric_matrix(MODEL);

	# Next, what is the size of the system? (ℳ = number of metabolites, ℛ = number of reactions)
	(ℳ,ℛ) = size(S)

	# Next, setup a default bounds array => update specific elements
	# We'll correct the directionality below -
	Vₘ = (13.7)*(3600)*(50e-9)*(1000) # units: mmol/hr
	flux_bounds = [-Vₘ*ones(ℛ,1) Vₘ*ones(ℛ,1)]

	# update the flux bounds -> which fluxes can can backwards?
	# do determine this: sgn(v) = -1*sgn(ΔG)
	updated_flux_bounds = update_flux_bounds_directionality(MODEL,flux_bounds)

	# hard code some bounds that we know -
	updated_flux_bounds[44,1] = 0.0  # ATP synthesis can't run backwards

	# What is the default mol flow input array => update specific elements
	# strategy: start with nothing in both streams, add material(s) back
	n_dot_input_stream_1 = zeros(ℳ,1)	# stream 1
	n_dot_input_stream_2 = zeros(ℳ,1)	# stream 2

	# === YOU NEED TO CHANGE BELOW HERE ====================================================== #
	# Let's lookup stuff that we want/need to supply to the chip to get the reactiont to go -
	# what you feed *depends upon your product*
	compounds_that_we_need_to_supply_feed_1 = [
		"alpha-d-glucose", "potassium nitrate"
	]

	# what are the amounts that we need to supply to chip in feed stream 1 (units: mmol/hr)?
	mol_flow_values_feed_1 = [
		13 ; # alpha-d-glucose mmol/hr
		30 ; #potassium nitrate mmol/hr
	]

	# what is coming into feed stream 2?
	compounds_that_we_need_to_supply_feed_2 = [
		"oxygen"
	]

	# let's always add Vₘ into feed stream 2
	mol_flow_values_feed_2 = [
		0.615; # oxygen mmol/hr
	]

	# === YOU NEED TO CHANGE ABOVE HERE ====================================================== #

	# stream 1:
	idx_supply_stream_1 = Array{Int64,1}()
	for compound in compounds_that_we_need_to_supply_feed_1
		idx = find_compound_index(MODEL,:compound_name=>compound)
		push!(idx_supply_stream_1,idx)
	end

	# stream 2:
	idx_supply_stream_2 = Array{Int64,1}()
	for compound in compounds_that_we_need_to_supply_feed_2
		idx = find_compound_index(MODEL,:compound_name=>compound)
		push!(idx_supply_stream_2,idx)
	end

	# supply for stream 1 and stream 2
	n_dot_input_stream_1[idx_supply_stream_1] .= mol_flow_values_feed_1
	n_dot_input_stream_2[idx_supply_stream_2] .= mol_flow_values_feed_2

	# setup the species bounds array -
	species_bounds = [-1.0*(n_dot_input_stream_1.+n_dot_input_stream_2) 1000.0*ones(ℳ,1)]

	# Lastly, let's setup the objective function -
	c = zeros(ℛ)
	c[idx_target_rate] = -1.0

	# show -
	nothing
end

# ╔═╡ 7410ecfc-0280-4ff2-8932-cfb2022e47c4
begin

	# compute the optimal flux -
	result = calculate_optimal_flux_distribution(S, updated_flux_bounds, species_bounds, c);

	# get the open extent vector -
	ϵ_dot = result.calculated_flux_array

	# what is the composition coming out of the first chip?
	n_dot_out_chip_1 = (n_dot_input_stream_1 + n_dot_input_stream_2 + S*ϵ_dot);

	# did this converge?
	with_terminal() do

		# get exit/status information from the solver -
		exit_flag = result.exit_flag
		status_flag = result.status_flag

		# display -
		println("Computed optimal flux distribution w/exit_flag = 0: $(exit_flag==0) and status_flag = 5: $(status_flag == 5)")
	end
end

# ╔═╡ 64daa21a-ac42-4b20-9e6b-ec2d19cd50fc
md"""
###### Table 1: Optimal reaction extent table computed by flux balance analysis.

Each row corresponds to a reaction in the model where $\dot{\epsilon}_{i}$ denotes the optimal open reaction extent computed by flux balance analysis.
"""

# ╔═╡ 637c5a39-63db-4d76-bb4c-ecc819482036
md""" ##### Connecting the Chip Reactions in Series"""

# ╔═╡ ae718670-bd74-48f5-96ed-d45a8d3bccb9
md"""Though one chip produces 2.466 mmol/hr, the output stream is  not solely composed of PGDN. We must arrange the chips in a series to achieve our goal of having a pure output of 1g/hr. From the program outlined below, a output of 14.796mmol/hr is found to produce enough PGDN in the stream to achieve our goal. We need approximately 6 chips to have enough PGDN after purification. Placing the chips in series is overall more economically efficient as it yields a higher output with fewer chips. The assumptions we make are that the microfluidic chip is well-mixed and operates at steady-state, the chip remains at constant temperature, pressure and the liquid phase is ideal.
"""

# ╔═╡ 3d84169a-3a85-4d9f-b7bb-b7264017fbb8
begin

	# setup calculation for chips i = 2,....,N
	N = 6 # number of chips

	# initialize some space to store the mol flow rates -
	series_mol_state_array = zeros(ℳ,N)
	exit_flag_array = Array{Int64,1}()
	status_flag_array = Array{Int64,1}()

	# the initial col of this array is the output of from chip 1
	for species_index = 1:ℳ
		series_mol_state_array[species_index,1] = n_dot_out_chip_1[species_index]
	end

	# assumption: we *always* feed glycerol into port 2 - so we only need to update the input flow into port 1
	for chip_index = 2:N

		# update the input into the chip -
		n_dot_input_port_1 = series_mol_state_array[:,chip_index - 1] 		# the input to chip j comes from j - 1

		# setup the species bounds array -
		species_bounds_next_chip = [-1.0*(n_dot_input_port_1.+n_dot_input_stream_2) 1000.0*ones(ℳ,1)]

		# run the optimal calculation -
		result_next_chip = calculate_optimal_flux_distribution(S, updated_flux_bounds, species_bounds_next_chip, c);

		# grab the status and exit flags ... so we can check all is right with the world ...
		push!(exit_flag_array, result_next_chip.exit_flag)
		push!(status_flag_array, result_next_chip.status_flag)

		# Get the flux from the result object -
		ϵ_dot_next_chip = result_next_chip.calculated_flux_array

		# compute the output from chip j = chip_index
		n_dot_out_next_chip = (n_dot_input_port_1 + n_dot_input_stream_2 + S*ϵ_dot_next_chip);

		# copy this state vector into the state array
		for species_index = 1:ℳ
			series_mol_state_array[species_index,chip_index] = n_dot_out_next_chip[species_index]
		end

		# go around again ...
	end
end

# ╔═╡ d541352c-4956-4bbe-a11a-27f78f3b8afe
exit_flag_array

# ╔═╡ 538f7122-e6ae-4256-8756-cd41d27ca085
status_flag_array

# ╔═╡ e933ddd9-8fd8-416a-8710-a64d3eb36f79
md"""
###### Table 2: State table from a single chip (species mol flow rate mmol/hr at exit).

The mass flow for all species at the exit of the chip is encoded in the `mass_dot_output_array` array.
"""

# ╔═╡ 7d35f315-927c-44b2-948a-c4e3d273a5e1
begin

	# what chip r we looking at?
	n_dot_output_chip = series_mol_state_array[:,end]

	# get the array of MW -
	MW_array = MODEL[:compounds][!,:compound_mw]

	# convert the output mol stream to a mass stream -
	mass_dot_output = (n_dot_output_chip.*MW_array)*(1/1000)

	# what is the total coming out?
	total_mass_out = sum(mass_dot_output)

	# display code makes the table -
	with_terminal() do

		# what are the compound names and code strings? -> we can get these from the MODEL object
		compound_name_strings = MODEL[:compounds][!,:compound_name]
		compound_id_strings = MODEL[:compounds][!,:compound_id]

		# how many molecules are in the state array?
		ℳ_local = length(compound_id_strings)

		# initialize some storage -
		number_of_cols = 3 + N + 2
		state_table = Array{Any,2}(undef,ℳ_local,number_of_cols)

		# get the uptake array from the result -
		uptake_array = result.uptake_array

		# populate the state table -
		for compound_index = 1:ℳ_local
			state_table[compound_index,1] = compound_index
			state_table[compound_index,2] = compound_name_strings[compound_index]
			state_table[compound_index,3] = compound_id_strings[compound_index]

			for chip_index = 1:N
				tmp_value = abs(series_mol_state_array[compound_index, chip_index])
				state_table[compound_index,chip_index + 3] = (tmp_value) <= 1e-6 ? 0.0 :
					series_mol_state_array[compound_index, chip_index]
			end

			# show the mass -
			tmp_value = abs(mass_dot_output[compound_index])
			state_table[compound_index,(N + 3 + 1)] = (tmp_value) <= 1e-6 ? 0.0 : mass_dot_output[compound_index]

			# show the mass fraction -
			# show the mass -
			tmp_value = abs(mass_dot_output[compound_index])
			state_table[compound_index, (N + 3 + 2)] = (tmp_value) <= 1e-6 ? 0.0 :
				(1/total_mass_out)*mass_dot_output[compound_index]
		end

		# build the table header -
		id_header_row = Array{String,1}()
		units_header_row = Array{String,1}()

		# setup id row -
		push!(id_header_row, "i")
		push!(id_header_row, "name")
		push!(id_header_row, "id")
		for chip = 1:N
			push!(id_header_row, "Chip $(chip)")
		end
		push!(id_header_row, "m_dot")
		push!(id_header_row, "ωᵢ_output")

		# setup units header row -
		push!(units_header_row, "")
		push!(units_header_row, "")
		push!(units_header_row, "")
		for chip = 1:N
			push!(units_header_row, "mmol/hr")
		end
		push!(units_header_row, "g/hr")
		push!(units_header_row, "")

		# header row -
		state_table_header_row = (id_header_row, units_header_row)

		# write the table -
		pretty_table(state_table; header=state_table_header_row)
	end
end

# ╔═╡ 7720a5a6-5d7e-4c79-8aba-0a4bb04973af
md"""
##### Compute the downstream separation using Magical Separation Units (MSUs)
"""

# ╔═╡ 8530bb61-ab28-4747-8037-0b9d481685a7
md"""After we obtained the 14.796mmol/hr stream, we then set up magical separator units that are 75% efficient in purifying our stream. The magical separator unit program works by establishing a theta value to which we assign the goal purity of our stream. In this scieneitro it is 0.75.  Then the program computes the composition of the upstream. Each level of separation added, we achieve a high purity. A minimum of 4 levels are needed to achieve  95.6 % purity at 1.08g/hr. To reduce cost, we decided to not filter out the waste stream for the main byproducts, potassium hydroxide and acetate, do not provide a significant return on value. Note: The MSU's are purely hypotheical and make the assumptions that the seperation of each stream occurs at exaclty 75% efficency and that the each of the streams are ideal
"""

# ╔═╡ 87e21b0b-5f6b-4402-baf3-68a150ef0fc2
md"""
Reaction string format: $(@bind rxn_string_ver Select(["first"=>"KEGG", "second"=>"HUMAN"]))
"""

# ╔═╡ 77df0589-a5a0-45cf-a0db-156621634262
with_terminal() do

	# initialize some storage -
	flux_table = Array{Any,2}(undef,ℛ,6)

	# what are the reaction strings? -> we can get these from the MODEL object
	reaction_strings = MODEL[:reactions][!,:reaction_markup]
	reaction_id = MODEL[:reactions][!,:reaction_number]

	# translate the reaction string to HUMAN -
	human_rxn_string_array = translation_reaction_string_to_human(MODEL)

	# populate the state table -
	for reaction_index = 1:ℛ
		flux_table[reaction_index,1] = reaction_index
		flux_table[reaction_index,2] = reaction_id[reaction_index]

			if (rxn_string_ver == "first")
				flux_table[reaction_index,3] = reaction_strings[reaction_index]
			else
				flux_table[reaction_index,3] = human_rxn_string_array[reaction_index]
		end
		flux_table[reaction_index,4] = flux_bounds[reaction_index,1]
		flux_table[reaction_index,5] = flux_bounds[reaction_index,2]

		# clean up the display -
		tmp_value = abs(ϵ_dot[reaction_index])
		flux_table[reaction_index,6] = tmp_value < 1e-6 ? 0.0 : ϵ_dot[reaction_index]
	end
		# header row -
	flux_table_header_row = (["i","RID","R","ϵ₁_dot LB", "ϵ₁_dot UB", "ϵᵢ_dot"],
		["","","", "mmol/hr", "mmol/hr", "mmol/hr"]);

	# write the table -
	pretty_table(flux_table; header=flux_table_header_row)
end

# ╔═╡ e1da943a-0d11-4776-bb67-2d8caad4cb18
# how many levels are we going to have in the separation tree?
number_of_levels = 4

# ╔═╡ 28a9763c-c5e1-41ae-b52c-445bcb839755
begin

	# define the split -
	θ = 0.75

	# most of the "stuff" has a 1 - θ in the up, and a θ in the down
	u = (1-θ)*ones(ℳ,1)
	d = θ*ones(ℳ,1)

	# However: the desired product has the opposite => correct for my compound of interest -> this is compound i = ⋆
	idx_target_compound = find_compound_index(MODEL,:compound_name=>"PGDN")

	# correct defaults -
	u[idx_target_compound] = θ
	d[idx_target_compound] = 1 - θ

	# let's compute the composition of the *always up* stream -

	# initialize some storage -
	species_mass_flow_array_top = zeros(ℳ,number_of_levels)
	species_mass_flow_array_bottom = zeros(ℳ,number_of_levels)

	for species_index = 1:ℳ
		value = mass_dot_output[species_index]
		species_mass_flow_array_top[species_index,1] = value
		species_mass_flow_array_bottom[species_index,1] = value
	end

	for level = 2:number_of_levels

		# compute the mass flows coming out of the top -
		m_dot_top = mass_dot_output.*(u.^(level-1))
		m_dot_bottom = mass_dot_output.*(d.^(level-1))

		# update my storage array -
		for species_index = 1:ℳ
			species_mass_flow_array_top[species_index,level] = m_dot_top[species_index]
			species_mass_flow_array_bottom[species_index,level] = m_dot_bottom[species_index]
		end
	end

	# what is the mass fraction in the top stream -
	species_mass_fraction_array_top = zeros(ℳ,number_of_levels)
	species_mass_fraction_array_bottom = zeros(ℳ,number_of_levels)

	# array to hold the *total* mass flow rate -
	total_mdot_top_array = zeros(number_of_levels)
	total_mdot_bottom_array = zeros(number_of_levels)

	# this is a dumb way to do this ... you're better than that JV come on ...
	T_top = sum(species_mass_flow_array_top,dims=1)
	T_bottom = sum(species_mass_flow_array_bottom,dims=1)
	for level = 1:number_of_levels

		# get the total for this level -
		T_level_top = T_top[level]
		T_level_bottom = T_bottom[level]

		# grab -
		total_mdot_top_array[level] = T_level_top
		total_mdot_bottom_array[level] = T_level_bottom

		for species_index = 1:ℳ
			species_mass_fraction_array_top[species_index,level] = (1/T_level_top)*
				(species_mass_flow_array_top[species_index,level])
			species_mass_fraction_array_bottom[species_index,level] = (1/T_level_bottom)*
				(species_mass_flow_array_bottom[species_index,level])
		end
	end
end

# ╔═╡ 24d220cd-0ead-44f1-9327-9db647b8108b
begin

	stages = (1:number_of_levels) |> collect
	plot(stages,species_mass_fraction_array_top[idx_target_compound,:], linetype=:steppre,lw=2,legend=:bottomright,
		label="Mass fraction i = PDO Tops")
	xlabel!("Stage index l",fontsize=18)
	ylabel!("Tops mass fraction ωᵢ (dimensionless)",fontsize=18)

	# make a 0.95 line target line -
	target_line = 0.95*ones(number_of_levels)
	plot!(stages, target_line, color="red", lw=2,linestyle=:dash, label="Target 95% purity")
end

# ╔═╡ cf367be1-ebde-4f46-974c-e2fdd1fdd903
with_terminal() do

	# initialize some space -
	state_table = Array{Any,2}(undef, number_of_levels, 3)
	for level_index = 1:number_of_levels
		state_table[level_index,1] = level_index
		state_table[level_index,2] = species_mass_fraction_array_top[idx_target_compound, level_index]
		state_table[level_index,3] = total_mdot_top_array[level_index]
	end

	# header -
	state_table_header_row = (["stage","ωᵢ i=⋆ top","mdot"],
			["","","g/hr"]);

	# write the table -
	pretty_table(state_table; header=state_table_header_row)
end

# ╔═╡ ea3305ae-cb8e-4a2b-83e2-ef882223e961
html"""
"""

# ╔═╡ 77107b01-dae9-4900-b9f7-f0c0224a492b
md"""
### Results and Discussion
"""

# ╔═╡ dcb915aa-e721-4080-a3fb-71b15d0b866d
md"""Although originally it was suggested that the sugar feedstocks be sucrose, we would advise the use of alpha-d-glucose instead. This sugar not only leads to a lower materials cost, but also produces a lower concentration of CO2 gas. Through a flux balance analysis, it was found that 1.08429g/hr of PDGN of 95.6% purity could be produced using 30mmol/hr of KNO3 and 13mmol/hr alpha glucose in stream 1, and .615mmol/hr in stream 2. With the reactor running for 24 hours a day, and the PDGN being sold for $1.20/gram, this process provides a $31.23 return each day.
The complete process requires 6 reactor chips, 4 separators, and 1 pump, with total equipment cost of $1680.00. The total cost of materials per day was calculated to be $0.72. Using this information, it was found that the entire process would take approximately 55 days to break even on this investment, and the total profit after the first year (based off of 252 trading days) is $6009.32. After this first year, the total profit per year is $7689.63.

The net present value of these final results and this reaction mechanism set at a 10% discount rate was then compared to an alternative investment with a 10% risk free interest rate that is constant over the payback lifetime, 5 years. The net present value of this process was calculated to be $25,932, a more favorable performance in comparison to investing the original $1,680 equipment cost in the alternative.
"""

# ╔═╡ 218fa01e-5230-49bb-a48a-ebafe81dac4a
md"""
##### *Financial Anylsis Using Net Present Value (NPV) Reuslts*
"""

# ╔═╡ 4030c9b7-b30d-462d-9d28-b8c2ba11a5e9
md"""
 ###### **Table 1** summarizes the cost and inventory needed for the process on Day 1.
"""

# ╔═╡ cf16db05-b530-4693-8607-f1f45aa702cc
with_terminal() do
data = ["Reactors" "100.00" "100.00" "6" "N/A" "N/A" "600.00"
		"Separators" "20.00" "20.00" "4" "N/A" "N/A" "80.00"
		"Pump" "1,000.00" "1,000.00" "1" "N/A" "N/A" "1,000.00"
		"Alpha D-Glucose" "17.90 per 25 grams" "0.72 per gram" "N/A" "0.00882764" "0.21186336" "0.15"
		"Oxygen" "Free" "Free" "N/A" "Free" "Free" "Free"
		"Potassium Nitrate" "57 per 100 grams" "0.57 per gram" "N/A" "0.0412501" "0.9900024" "0.56"
		"PGDN" "1196.5 per 1 kg" "1.20 per gram" "N/A" "1.08429" 		"26.02296" "31.23"];

header=(["Items", "Purchase Cost", "Calculation Cost", "Number of Equipment Purchased", "Mass Flow Rates", "Cost to run for 24 hours ", "Total Costs after 1 day"],["dollars","dollars","dollars","","g/hr","dollars","dollars"])

hl_lastRow = Highlighter(f      = (data, i, j) -> i %7 == 0,
                            crayon = Crayon(background = :light_blue))
hl_lastColumn = Highlighter(f      = (data, i, j) -> j %7 == 0,
                            crayon = Crayon(background = :light_blue))
pretty_table(data; header,highlighters =(hl_lastColumn, hl_lastRow))
end

# ╔═╡ 22178cca-c670-4d4f-9021-22328447f9d3
md"""
###### *Note: The last highlighted row represents the revenue made on a day-to-day basis from selling the final product (1,2-propanediol dinitrate).*

From Table 1 above, the total initial (capital) cost for equipment and to run the process on Day 1 accumulates to a total of $1680.72. Then, the cost to run the process will be an additional cost of $0.72 each day. This simple relationship is summarized in the linear equation below.

$$TotalSpendingOnDayX = 0.76x+1680.72$$

On the revenue side, the product that is sold each day provides a return of $31.23.

Total Revenue on Day

$$TotalSpendingOnDayX= 31.23x$$

As a result, it would take roughly 55 days of running the entire process to break even on the investment (by calculating for x). Additionally, at the end of the first year (following the typical 252 trading days rule), the total profit will sit at $6009.32. The results are summarized in Table 2.

"""

# ╔═╡ d6e3755f-8c25-4627-90a8-e056fc2b640b
md"""
 ###### **Table 2**
"""

# ╔═╡ cbe36492-a695-4514-bbd8-c2df553a120d
with_terminal() do
data = ["Total Spending" "1,680.72" "1681.44" "1682.16" "1,719.43" "1,720.14" "1,860.64"
"Total Revenue" "1.23" "31.23" "31.23" "1,717.65" "1,748.88" "7,869.96"
"Profit" "-1649.49" "-1650.21" "-1650.93" "-1.78" "28.74" "6,009.32"
];

header=(["" ,"Inital Cost", "Cost on Day 1", "Cost on Day 2", "After 55 days", "After 56 Days", "After 1 year(252 days)"], ["","dollars", "dollars", "dollars", "dollars", "dollars", "dollars"])

pretty_table(data; header)
end

# ╔═╡ af08e78f-2910-47d6-a4c9-0e01ec0a7697
md"""  **Table 3**: Total Profit Process Generates Over the Course of 5 years.
"""

# ╔═╡ eef3b8c8-410e-4880-a55e-2a09c065bef2
with_terminal() do
data = [
"Total Spending" "1,860.64" "2,041.29" "2,221.93" "2,403.29" "2,583.94"
"Total Revenues" "7,869.96" "15,739.92" "23,609.88" "31,479.84" "39,349.80"
"Total Profit" "6,009.32" "13,698.63" "21,387.95" "29,076.55" "36,765.86"
"Difference from prior year" "6,009.32" "7689.63" "7689.63" "7689.63" "7689.63"
];

header=(["","After 1 Year","After 2 Years","After 3 Years", "After 4 Years","After 5 Years"],["","(252 days)","(504 days)","(756 days)", "(1008 days)","(1260 days)"],["", "dollars", "dollars", "dollars", "dollars", "dollars"])

hl_lastRow = Highlighter(f = (data, i, j) -> i %3 == 0,
                            crayon = Crayon(background = :light_blue))

pretty_table(data; header,highlighters =hl_lastRow)
end

# ╔═╡ 7850b71d-d834-43ce-9775-77e8c6e7513b
md"""

As noted after the first year profit of $6,009.32, the process produces a constant profit of $7689.63 each year.

These final results were compared to an alternative investment at both a constant 10% discount rate and at a constant 1% discount rate via the Net Present Value method.
* The lifetime of the comparison was set to 5 years.
* The constant discount rate was set to 10% or 0.1 in the parameters, indicating a comparison to a 10% risk-free interest rate that is constant over the payback lifetime.
* Since the cash flow array was in the thousands, the values were imputed as listed: [-1.68, 6.0, 7.689, 7.689, 7.689, 7.689]. The initial -1.68 represents the initial (capital) investment that was made, and each of the following terms represents the total profit made each year following this investment.
* The Present Value of the cash values are as listed: [-1.68, 5.45455, 6.35455, 5.77686, 5.25169, 4.774326]. And then by summing up these values, the Net Present Value is calculated to be 25.9319.
* Since the Net Present Value is positive, the project would be a better investment of the initial amount of $1,680 than an alternative investment at a 10% interest rate.
* The same logic can be applied to an alternative investment at a 1% interest rate.

Using the Net Present Value approach, the process performance is indeed more economically feasible than both the alternative investments at a constant yield of 1 % per year or at a constant yield of 10% per year.

"""

# ╔═╡ ab9084e1-57a2-4c6b-ad32-8fee8c142c43
md"""
### Conclusions
"""

# ╔═╡ 3bea3df5-a395-4831-ab43-1c713deb69f2
md"""
While the sucrose to PGDN pathway is a feasible option, it is more economically efficient to use alpha-d-glucose as the starting material. When transferring from sucrose to alpha-d-glucose there was a reduction in co2 production resulting making this pathway a slightly more environmentally favorable option. The pathway analyzed in this reaction is advisable as an alternative to the direct production feeding in 1-2propaindiol and KNO3. The downsides of this reaction is that in creating the PGDN, the Co2 byproduct is approximately 1.25x the amount of PGDN synthesized. [4]Methods such as absorption, adsorption, chemical looping, membrane gas separation or gas hydration have been seen to be useful in capturing Co2. The Co2 emissions should be readily monitored and to ensure it does not exceed a certain threshold. These measures must be taken into account, and thus are one downside to the process. [1] Co2, if captured, however, can be sold back to the food production industry for beverages and other sources. The scale of production of the PGDN would have to be significantly greater for this to be an economically viable option. In addition to the Co2 byproduct, the reaction also produces significant amounts of potassium hydroxide and acetate. [7] Potassium Hydroxide is relatively cheap thus is not worth filtering and selling. This chemical however must be carefully disposed of by diluting and neutralizing. Potassium hydroxide is potentially detrimental to air, soil, and water quality and regulations must be carefully abided by during the production process. [6] The acetate is most easily disposed of by controlled incineration. The management of the byproduct disposal is a tedious yet vital part of the production process however, the overall yield and turnover rate is great enough to still profit off the process while maintaining regulatory compliance.
"""

# ╔═╡ 836e69f7-9a2d-4674-8c20-51b07d13b7ab
md"""
### References
"""

# ╔═╡ f8e863a7-76df-4a78-b86b-9aacc5490663
md"""
[1] Chameides, E. R. and B. (2008, April 17). CO2: They Should Bottle That Stuff. Time. Retrieved December 10, 2021, from http://content.time.com/time/specials/2007/article/0,28804,1730759_1731383_1731989,00.html.

[2] The Custodian. (2009, July 19). Otto Fuel II. Otto fuel II - Everything2.com. Retrieved December 11, 2021, from https://everything2.com/title/Otto+fuel+II.

[3] D -glucose anhydrous, 96 492-62-6. Sigma Aldrich. (n.d.). Retrieved December 10, 2021, from https://www.sigmaaldrich.com/US/en/product/aldrich/158968.

[4] Environmental Protection Agency. (2021, January). Pollution Prevention and Waste Management. EPA. Retrieved December 10, 2021, from https://www.epa.gov/trinationalanalysis/pollution-prevention-and-waste-management. TRI National Analysis

[5] Hichem FETTAKA,* Michel H. LEFEBVRE. (2016). Propylene glycol dinitrate (PGDN) as an explosive taggant. Central European Journal of Energetic Materials. Retrieved December 11, 2021, from https://ipo.lukasiewicz.gov.pl/wydawnictwa/wp-content/uploads/2021/04/Fettaka-1.pdf.

[6] National Institutes of Health. (n.d.). Ethyl acetate - disposal methods. U.S. National Library of Medicine. Retrieved December 10, 2021, from https://webwiser.nlm.nih.gov/substance?substanceId=446&identifier=Ethyl+Acetate&identifierType=name&menuItemId=59&catId=76.

[7] National Institutes of Health. (n.d.). Potassium hydroxide - disposal methods. U.S. National Library of Medicine. Retrieved December 10, 2021, from https://webwiser.nlm.nih.gov/substance?substanceId=401&identifier=Potassium+Hydroxide&identifierType=name&menuItemId=59&catId=76.

[8] Sucrose for Molecular Biology, = 99.5 GC ... - sigma-aldrich. Sigma Aldrich. (n.d.). Retrieved December 10, 2021, from https://www.sigmaaldrich.com/US/en/product/sigma/s0389.

"""

# ╔═╡ 18b29a1a-4787-11ec-25e3-5f29ebd21430
html"""


text:
module Unreal

export unreal

"""
    unreal(::T) where T <: Real

The unreal part of a real number is an empty array.
"""
function unreal(::T) where T <: Real
    T[]
end

end


text:
#Almost always precompile
__precompile__()

module IntegralTransforms
# Use our Dependencies
#
using ForwardDiff, MacroTools, StaticArrays, StructsOfArrays, SpecialFunctions, Compat

# Define our abstract types
abstract type AbstractIntegralTransform end
abstract type AbstractTransform <: AbstractIntegralTransform end
abstract type AbstractLaplaceTransform <: AbstractIntegralTransform end
abstract type  AbstractForwardLaplaceTransform <: AbstractLaplaceTransform end
abstract type AbstractInverseLaplaceTransform <: AbstractLaplaceTransform end
# We export the public API, not the internal functionality.
# \U+2212
export ℒ, invℒ, func, func2, func3, func4, func5, func6

include("laplace.jl")

end # module


text:
abstract type AbstractSubstance end
abstract type AbstractPhaseState end
abstract type AbstractMixture end
abstract type AbstractThermoBackend end
abstract type AbstractEquationOfState <: AbstractThermoBackend end
abstract type AbstractCubicEOS <: AbstractEquationOfState end
abstract type AbstractMixingRules <: AbstractThermoBackend end
abstract type AbstractActivityCoefficientModel <: AbstractThermoBackend end
abstract type AbstractProperty end

@enum ThermoProp pressure=1 volume=2 temperature=3 entropy=4 internalenergy=5 enthalpy=6 helmholtzenergy=7  gibbsenergy=8 heatcapacity=9 vaporpressure=10 chemicalpotential=11 fugacity=12 fugacitycoefficient=13 activity=14 activitycoefficient=15
@enum MixingRules none=0 wongsandler=1

function thermocalc(st::AbstractPhaseState, eos::AbstractEquationOfState, p::AbstractProperty)
println("We're sorry, the property $show(p) isn't implemented for $show(st) , using the $show(eos) equation of state.") end
function thermocalc!(st::AbstractPhaseState, eos::AbstractEquationOfState, p::AbstractProperty)
println("We're sorry, the property $show(p) isn't implemented for $show(st) , using the $show(eos) equation of state.") end


text:
using AdFem
using PyPlot

function g_func(x, y)
    -2*x*(1-x) - 2*y*(1-y) + 2*(1-2x)*y*(1-y) + 3*(1-2y)*x*(1-x)
end

m = 30
n = 30
h = 1/n

xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
u = @. (x*(1-x)*y*(1-y))

A = constant(compute_fem_laplace_matrix1(m, n, h))
u0 = 2ones(2(m+1)*(n+1))
v0 = 3ones(2(m+1)*(n+1))
u0 = fem_to_gauss_points(u0, m, n, h)
v0 = fem_to_gauss_points(v0, m, n, h)
B = compute_fem_advection_matrix1(constant(u0), constant(v0), m, n, h)
L = -A+B
g = eval_f_on_gauss_pts(g_func, m, n, h)
rhs = compute_fem_source_term1(g, m, n, h)
bd = bcnode("all", m, n, h)
L, _ = fem_impose_Dirichlet_boundary_condition1(L, bd, m, n, h)
rhs[bd] .= 0.0
sol = L\rhs

sess = Session(); init(sess)
S = run(sess, sol)
figure(figsize=(10,4))
subplot(121)
visualize_scalar_on_fem_points(S, m, n, h)
title("computed solution")
subplot(122)
visualize_scalar_on_fem_points(u, m, n, h)
title("exact solution")

# subplot(133)
# visualize_scalar_on_fem_points(S-u, m, n, h)
# title("difference between computed and exact")


text:
using IJuliaTimeMachine
using Documenter

cp(normpath(@__FILE__, "../../README.md"), normpath(@__FILE__, "../src/overview.md"); force=true)
cp(normpath(@__FILE__, "../../examples/Demo.pdf"), normpath(@__FILE__, "../src/Demo.pdf"); force=true)

makedocs(;
    modules=[IJuliaTimeMachine],
    authors="Daniel Spielman <daniel.spielman@yale.edu> and contributors",
    repo="https://github.com/danspielman/IJuliaTimeMachine.jl/blob/{commit}{path}#L{line}",
    sitename="IJuliaTimeMachine.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://danspielman.github.io/IJuliaTimeMachine.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Docstrings" => "docstrings.md"
    ],
)

deploydocs(;
    repo="github.com/danspielman/IJuliaTimeMachine.jl",
)

rm(normpath(@__FILE__, "../src/overview.md"))
rm(normpath(@__FILE__, "../src/Demo.pdf"))

max_stars_repo_name: danspielman/IJuliaTimeMachine

max_issues_repo_name: danspielman/IJuliaTimeMachine

max_forks_repo_name: danspielman/IJuliaTimeMachine


text:
import ChainRulesCore:
    rrule, @non_differentiable, NoTangent, Tangent, backing, AbstractTangent, ZeroTangent

function create_circuit_tangent(circuit, params)
    gc = dispatch(circuit, params)
    res = recursive_create_tangent(gc)
    return res
end

# fallback
function recursive_create_tangent(c::AbstractBlock)
    if nparameters(c) == 0
        return NoTangent()
    else
        error("`Tangent` for type $(typeof(c)) is not defined!")
    end
end
# primitive blocks
unsafe_primitive_tangent(::Any) = NoTangent()
unsafe_primitive_tangent(x::Number) = x
for GT in [:RotationGate, :ShiftGate, :TimeEvolution, :PhaseGate]
    @eval function recursive_create_tangent(c::$GT)
        lst = map(fieldnames(typeof(c))) do fn
            fn => unsafe_primitive_tangent(getfield(c, fn))
        end
        nt = NamedTuple(lst)
        Tangent{typeof(c),typeof(nt)}(nt)
    end
end
# composite blocks
unsafe_composite_tangent(::Any) = NoTangent()
unsafe_composite_tangent(c::AbstractVector{<:AbstractBlock}) = recursive_create_tangent.(c)
unsafe_composite_tangent(c::AbstractBlock) = recursive_create_tangent(c)
for GT in [
    :ChainBlock,
    :Add,
    :KronBlock,
    :RepeatedBlock,
    :PutBlock,
    :Subroutine,
    :CachedBlock,
    :Daggered,
    :Scale,
]
    @eval function recursive_create_tangent(c::$GT)
        lst = map(fieldnames(typeof(c))) do fn
            fn => unsafe_composite_tangent(getfield(c, fn))
        end
        nt = NamedTuple(lst)
        Tangent{typeof(c),typeof(nt)}(nt)
    end
end

extract_circuit_gradients!(c::Number, output) = push!(output, c)
extract_circuit_gradients!(::Nothing, output) = output
extract_circuit_gradients!(::NoTangent, output) = output
extract_circuit_gradients!(::ZeroTangent, output) = output
function extract_circuit_gradients!(c::AbstractVector, output)
    for ci in c
        extract_circuit_gradients!(ci, output)
    end
    return output
end
function extract_circuit_gradients!(c::Tangent, output)
    for fn in getfield(c, :backing)
        extract_circuit_gradients!(fn, output)
    end
    return output
end
function extract_circuit_gradients!(c::NamedTuple, output)
    for fn in c
        extract_circuit_gradients!(fn, output)
    end
    return output
end

function rrule(::typeof(apply), reg::ArrayReg, block::AbstractBlock)
    out = apply(reg, block)
    out, function (outδ)
        (in, inδ), paramsδ = apply_back((copy(out), outδ), block)
        return (NoTangent(), inδ, create_circuit_tangent(block, paramsδ))
    end
end
function rrule(::typeof(apply), reg::ArrayReg, block::Add)
    out = apply(reg, block)
    out, function (outδ)
        (in, inδ), paramsδ = apply_back((copy(out), outδ), block; in = reg)
        return (NoTangent(), inδ, create_circuit_tangent(block, paramsδ))
    end
end

function rrule(::typeof(dispatch), block::AbstractBlock, params)
    out = dispatch(block, params)
    out, function (outδ::AbstractTangent)
        g = extract_circuit_gradients!(outδ, empty(params))
        res = (NoTangent(), NoTangent(), g)
        return res
    end
end

function rrule(::typeof(expect), op::AbstractBlock, reg::AbstractRegister{B}) where {B}
    out = expect(op, reg)
    out, function (outδ)
        greg = expect_g(op, reg)
        for b in 1:B
            viewbatch(greg, b).state .*= 2 * outδ[b]
        end
        return (NoTangent(), NoTangent(), greg)
    end
end

function rrule(
    ::typeof(expect),
    op::AbstractBlock,
    reg_and_circuit::Pair{<:ArrayReg{B},<:AbstractBlock},
) where {B}
    out = expect(op, reg_and_circuit)
    out,
    function (outδ)
        greg, gcircuit = expect_g(op, reg_and_circuit)
        for b in 1:B
            viewbatch(greg, b).state .*= 2 * outδ[b]
        end
        return (
            NoTangent(),
            NoTangent(),
            Tangent{typeof(reg_and_circuit)}(;
                first = greg,
                second = create_circuit_tangent(reg_and_circuit.second, gcircuit),
            ),
        )
    end
end

function rrule(::Type{T}, block::AbstractBlock) where {T<:Matrix}
    out = T(block)
    out, function (outδ)
        paramsδ = mat_back(block, outδ)
        return (NoTangent(), create_circuit_tangent(block, paramsδ))
    end
end

function rrule(::typeof(mat), ::Type{T}, block::AbstractBlock) where {T}
    out = mat(T, block)
    out, function (outδ)
        paramsδ = mat_back(block, outδ)
        return (NoTangent(), NoTangent(), create_circuit_tangent(block, paramsδ))
    end
end

function rrule(::Type{ArrayReg{B}}, raw::AbstractArray) where {B}
    ArrayReg{B}(raw), adjy -> (NoTangent(), reshape(adjy.state, size(raw)))
end

function rrule(::Type{ArrayReg}, raw::AbstractArray)
    ArrayReg(raw), adjy -> (NoTangent(), reshape(adjy.state, size(raw)))
end

function rrule(::typeof(copy), reg::ArrayReg) where {B}
    copy(reg), adjy -> (NoTangent(), adjy)
end

for (BT, BLOCKS) in [(:Add, :(outδ.list)) (:ChainBlock, :(outδ.blocks))]
    for ST in [:AbstractVector, :Tuple]
        @eval function rrule(::Type{BT}, source::$ST) where {N,BT<:$BT}
            out = BT(source)
            out, function (outδ)
                return (NoTangent(), $ST($BLOCKS))
            end
        end
    end
    @eval function rrule(::Type{BT}, args::AbstractBlock...) where {N,BT<:$BT}
        out = BT(args...)
        out, function (outδ)
            return (NoTangent(), $BLOCKS...)
        end
    end
end

_totype(::Type{T}, x::AbstractArray{T}) where {T} = x
_totype(::Type{T}, x::AbstractArray{T2}) where {T,T2} = convert.(T, x)
rrule(::typeof(state), reg::ArrayReg{B,T}) where {B,T} =
    state(reg), adjy -> (NoTangent(), ArrayReg(_totype(T, adjy)))
rrule(::typeof(statevec), reg::ArrayReg{B,T}) where {B,T} =
    statevec(reg), adjy -> (NoTangent(), ArrayReg(_totype(T, adjy)))
rrule(::typeof(state), reg::AdjointArrayReg{B,T}) where {B,T} =
    state(reg), adjy -> (NoTangent(), ArrayReg(_totype(T, adjy)')')
rrule(::typeof(statevec), reg::AdjointArrayReg{B,T}) where {B,T} =
    statevec(reg), adjy -> (NoTangent(), ArrayReg(_totype(T, adjy)')')
rrule(::typeof(parent), reg::AdjointArrayReg) = parent(reg), adjy -> (NoTangent(), adjy')
rrule(::typeof(Base.adjoint), reg::ArrayReg) = Base.adjoint(reg), adjy -> (NoTangent(), parent(adjy))
@non_differentiable nparameters(::Any)
@non_differentiable zero_state(args...)
@non_differentiable rand_state(args...)
@non_differentiable uniform_state(args...)
@non_differentiable product_state(args...)

max_stars_repo_name: QuantumBFS/YaoStdCircuit

max_issues_repo_name: QuantumBFS/YaoStdCircuit


text:
module MetadataWindow

using Gtk, Gtk.ShortNames, Gtk.GConstants
using ..Ahorn, Maple

metadataWindow = nothing

const metaDropdownOptions = Dict{String, Any}(
    "IntroType" => Maple.intro_types,
    "ColorGrade" => Maple.color_grades,
    "CoreMode" => Maple.core_modes,
    "CassetteSong" => sort(collect(values(Maple.CassetteSongs.songs))),
    "Wipe" => Maple.wipe_names
)

const modeDropdownOptions = Dict{String, Any}(
    "Inventory" => Maple.inventories,
    "StartLevel" => String[]
)

const audiostateDropdownOptions = Dict{String, Any}(
    "Music" => sort(collect(values(Maple.Songs.songs))),
    "Ambience" => sort(collect(values(Maple.AmbientSounds.sounds)))
)

const metaFieldOrder = String[
    "Name", "SID", "Icon",
    "CompleteScreenName", "CassetteCheckpointIndex",
    "CassetteNoteColor", "CassetteSong",
    "TitleBaseColor", "TitleAccentColor",
    "TitleTextColor", "IntroType",
    "ColorGrade", "Wipe",
    "DarknessAlpha", "BloomBase", "BloomStrength",
    "Jumpthru", "CoreMode",
    "ForegroundTiles", "BackgroundTiles", "AnimatedTiles",
    "Portraits", "Sprites"
]

const modeFieldOrder = String[]
const audiostateFieldOrder = String[]

function getOptions(data::Dict{String, Any}, dropdownOptions::Dict{String, Any}, langdata::Ahorn.LangData, defaults::Dict{String, Any})
    res = Ahorn.Form.Option[]

    names = get(langdata, :names)
    tooltips = get(langdata, :tooltips)

    for (attr, value) in data
        # Due to the way Metadata is structured we don't want to generate an option for Dicts
        if isa(value, Dict)
            continue
        end

        sattr = Symbol(attr)
        keyOptions = get(dropdownOptions, attr, nothing)
        displayName = haskey(names, sattr) ? names[sattr] : Ahorn.humanizeVariableName(attr)
        tooltip = Ahorn.expandTooltipText(get(tooltips, sattr, ""))
        preferredType = typeof(get(defaults, attr, value))

        push!(res, Ahorn.Form.suggestOption(displayName, value, tooltip=tooltip, dataName=attr, choices=keyOptions, editable=true, preferredType=preferredType))
    end

    return res
end

function writeIfExists!(data::Dict{K, V}, key::K, value::V) where {K, V}
    if value !== nothing && !isempty(value)
        data[key] = value

    else
        # Delete the key to make sure any old values are cleared
        delete!(data, key)
    end
end

function addKeyIfAbsent!(data::Dict{K, V}, key::K) where {K, V}
    if !haskey(data, key)
        data[key] = Dict{K, V}()
    end
end

function deleteKeyIfEmpty!(data::Dict{K, V}, key::K) where {K, V}
    if haskey(data, key) && isempty(data[key])
        delete!(data, key)
    end
end

function createWindow()
    targetSide = Ahorn.loadedState.side

    displayCustomValues = get(Ahorn.debug.config, "DISPLAY_ALL_META_ATTRIBUTES", false)

    # Fill in possible start room values manually
    modeDropdownOptions["StartLevel"] = String[room.name for room in targetSide.map.rooms]

    langdata = get(Ahorn.langdata, :meta_window)

    metaLangdata = get(langdata, :meta)
    metaData = merge(Maple.default_meta, get(Dict{String, Any}, targetSide.data, "meta"))

    modeLangdata = get(langdata, :mode)
    modeData = merge(Maple.default_mode, get(Dict{String, Any}, metaData, "mode"))

    audiostateLangdata = get(langdata, :audiostate)
    audiostateData = merge(Maple.default_audiostate, get(Dict{String, Any}, modeData, "audiostate"))

    if !displayCustomValues
        filter!(p -> haskey(Maple.default_meta, p.first), metaData)
        filter!(p -> haskey(Maple.default_mode, p.first), modeData)
        filter!(p -> haskey(Maple.default_audiostate, p.first), audiostateData)
    end

    modeOptions = getOptions(modeData, modeDropdownOptions, modeLangdata, Maple.default_mode)
    modeSection = Ahorn.Form.Section("Mode", modeOptions, dataName="mode", fieldOrder=modeFieldOrder)

    metaOptions = getOptions(metaData, metaDropdownOptions, metaLangdata, Maple.default_meta)
    metaSection = Ahorn.Form.Section("Meta", metaOptions, dataName="meta", fieldOrder=metaFieldOrder)

    audiostateOptions = getOptions(audiostateData, audiostateDropdownOptions, audiostateLangdata, Maple.default_audiostate)
    audiostateSection = Ahorn.Form.Section("Audio", audiostateOptions, dataName="audiostate", fieldOrder=audiostateFieldOrder)

    sections = Ahorn.Form.Section[
        metaSection, modeSection, audiostateSection
    ]

    function callback(data::Dict{String, Dict{String, Any}})
        addKeyIfAbsent!(targetSide.data, "meta")
        addKeyIfAbsent!(targetSide.data["meta"], "mode")
        addKeyIfAbsent!(targetSide.data["meta"]["mode"], "audiostate")

        for (attr, value) in data["meta"]
            writeIfExists!(targetSide.data["meta"], attr, value)
        end

        for (attr, value) in data["mode"]
            writeIfExists!(targetSide.data["meta"]["mode"], attr, value)
        end

        for (attr, value) in data["audiostate"]
            writeIfExists!(targetSide.data["meta"]["mode"]["audiostate"], attr, value)
        end

        deleteKeyIfEmpty!(targetSide.data["meta"]["mode"], "audiostate")
        deleteKeyIfEmpty!(targetSide.data["meta"], "mode")
        deleteKeyIfEmpty!(targetSide.data, "meta")

        Ahorn.loadXMLMeta()

        name = get(data["meta"], "Name", targetSide.map.package)
        GAccessor.title(Ahorn.window, "$(Ahorn.baseTitle) - $name")
    end

    metadataWindow = Ahorn.Form.createFormWindow("$(Ahorn.baseTitle) - Configure Metadata", sections, callback=callback)

    return metadataWindow
end

function configureMetadata(widget::Union{Ahorn.MenuItemsTypes, Nothing}=nothing)
    if Ahorn.loadedState.side !== nothing
        if metadataWindow !== nothing
            Gtk.destroy(metadataWindow)
        end

        @Ahorn.catchall global metadataWindow = createWindow()

        showall(metadataWindow)
    end
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
#= ============================================
Generate 𝐹 and ∇ₓ𝐹 from user input
============================================ =#

# TODO replace this with DiffEqOperators when possible
struct LinearOperators{T<:Tuple}
    ops::T
end
LinearOperators(ops...) = LinearOperators((ops...,))
SparseArrays.blockdiag(As::LinearOperators...) = blockdiag([sum(A.ops) for A in As]...)
function LinearAlgebra.mul!(du, A::LinearOperators, u, α, β)
    for (i, op) in enumerate(A.ops)
        (i==1) ? mul!(du, op, u, α, β) : mul!(du, op, u, α, 1)
    end
    du
end
LinearAlgebra.factorize(A::LinearOperators) = factorize(sum(A.ops))
Base.:\(A::LinearOperators, u::AbstractArray) = factorize(A) \ u
Base.:+(A::LinearOperators, B::AbstractSparseArray) = LinearOperators(A.ops..., B)
Base.:+(B::AbstractSparseArray, A::LinearOperators) = A + B
function Base.:*(A::LinearOperators, u::AbstractArray) # has to be same type!
    du = similar(u)
    mul!(du, A, u, 1, 0)
end
export LinearOperators

# AIBECSFunction creates an ODEFunction from the Ts and Gs
AIBECSFunction(T::AbstractSparseArray, G::Function) = AIBECSFunction(p -> T, G, T.n)
AIBECSFunction(T::Function, G::Function, nb::Int) = AIBECSFunction(T, (G,), nb)
function AIBECSFunction(T::Function, Gs::Tuple, nb::Int)
    nt = length(Gs) # if a single T is given, then nt is given by the number of Gs
    Tidx = ones(Int64, nt) # and all the tracers share the same T
    AIBECSFunction((T,), Gs, nb, nt, Tidx)
end
function AIBECSFunction(Ts::Tuple, Gs::Tuple, nb::Int, nt::Int=length(Ts), Tidx::AbstractVector=1:nt)
    if all(isinplace(G, nt + 2) for G in Gs) # +2 to account for dx and p
        iipAIBECSFunction(Ts, Gs, nb, nt, Tidx)
    else
        oopAIBECSFunction(Ts, Gs, nb, nt, Tidx)
    end
end

function iipAIBECSFunction(Ts, Gs, nb, nt=length(Ts), Tidx=1:nt)
    tracers(u) = state_to_tracers(u, nb, nt)
    tracer(u, i) = state_to_tracer(u, nb, nt, i)
    function G(du, u, p)
        for j in 1:nt
            Gs[j](tracer(du, j), tracers(u)..., p)
        end
        du
    end
    function f(du, u, p, t=0)
        G(du, u, p)
        for jT in eachindex(Ts)
            op = Ts[jT](p)
            for j in findall(Tidx .== jT)
                mul!(tracer(du, j), op, tracer(u, j), -1, 1)
            end
        end
        du
    end
    # Jacobian
    ∇ₓG(u, p) = inplace_local_jacobian(Gs, u, p, nt, nb)
    function T(p)
        uniqueTs = [Tⱼ(p) for Tⱼ in Ts]
        blockdiag([uniqueTs[Tidx[j]] for j in 1:nt]...)
    end
    jac(u, p, t=0) = ∇ₓG(u, p) - T(p)
    return ODEFunction{true}(f, jac=jac)
end
function oopAIBECSFunction(Ts, Gs, nb, nt=length(Ts), Tidx=1:nt)
    tracers(u) = state_to_tracers(u, nb, nt)
    tracer(u, i) = state_to_tracer(u, nb, nt, i)
    G(u, p) = reduce(vcat, Gⱼ(tracers(u)..., p) for Gⱼ in Gs)
    function f(u, p, t=0)
        du = copy(G(u, p))
        for jT in eachindex(Ts)
            op = Ts[jT](p)
            for j in findall(Tidx .== jT)
                mul!(tracer(du, j), op, tracer(u, j), -1, 1)
            end
        end
        du
    end
    # Jacobian
    ∇ₓG(u, p) = local_jacobian(Gs, u, p, nt, nb)
    function T(p)
        uniqueTs = [Tⱼ(p) for Tⱼ in Ts]
        blockdiag([uniqueTs[Tidx[j]] for j in 1:nt]...)
    end
    jac(u, p, t=0) = ∇ₓG(u, p) - T(p)
    return ODEFunction(f, jac=jac)
end
# AIBECSFunction calls itself to allow for both λ::Vector and p::APar
function AIBECSFunction(Ts, Gs, nb::Int, ::Type{P}) where {P <: APar}
    AIBECSFunction(AIBECSFunction(Ts, Gs, nb), P)
end
function AIBECSFunction(fun::ODEFunction{false}, ::Type{P}) where {P <: APar}
    jac(u, p::P, t=0) = fun.jac(u, p, t)
    jac(u, λ::Vector, t=0) = fun.jac(u, λ2p(P, λ), t)
    f(u, p::P, t=0) = fun.f(u, p, t)
    f(u, λ::Vector, t=0) = fun.f(u, λ2p(P, λ), t)
    return ODEFunction{false}(f, jac=jac)
end
function AIBECSFunction(fun::ODEFunction{true}, ::Type{P}) where {P <: APar}
    jac(u, p::P, t=0) = fun.jac(u, p, t)
    jac(u, λ::Vector, t=0) = fun.jac(u, λ2p(P, λ), t)
    f(du, u, p::P, t=0) = fun.f(du, u, p, t)
    f(du, u, λ::Vector, t=0) = fun.f(du, u, λ2p(P, λ), t)
    return ODEFunction{true}(f, jac=jac)
end

export AIBECSFunction

"""
    F, ∇ₓF = F_and_∇ₓF(Ts, Gs, nb)

Returns the state function `F` and its Jacobian, `∇ₓF`.

    F, ∇ₓF = F_and_∇ₓF(T, Gs, nb)

Returns the state function `F` and its Jacobian, `∇ₓF` (with all tracers transported by single `T`).

This function is deprecated. Use

    F = AIBECSFunction(Ts, Gs, nb)

instead.
You can then call `F(x,p)` for the tendencies, and `F(Val{:jac},x,p)` for the Jacobian.
"""
function F_and_∇ₓF(fun::ODEFunction)
    Base.depwarn("""Deprecation:
        F, ∇ₓF = F_and_∇ₓF(args...)
    is deprecated. Use
        F = AIBECSFunction(args...)
    instead! (And then you can still call `F(x,p)`.)
    """, :F_and_∇ₓF, force=true)
    fun.f, fun.jac
end
F_and_∇ₓF(args...) = F_and_∇ₓF(AIBECSFunction(args...))
export F_and_∇ₓF

"""
    localderivative(G, x, p)
    localderivative(Gᵢ, xs, i, p)
    localderivative(Gᵢ, dx, xs, i, p)

Returns the "local" derivative of `G` (or `Gᵢ`), i.e., equivalent to the vector

```
∇ₓG(x,p) * ones(size(x))
```

but using ForwardDiff's Jacobian instead.
"""
function localderivative(G, x, p) # for single tracer
    return ForwardDiff.derivative(λ -> G(x .+ λ, p), 0.0)
end
function localderivative(Gᵢ, xs, j, p) # for multiple tracers
    return ForwardDiff.derivative(λ -> Gᵢ(perturb_tracer(xs, j, λ)..., p), 0.0)
end
function localderivative(Gᵢ!, dx, xs, j, p) # if Gᵢ are in-place
    return ForwardDiff.derivative((dx, λ) -> Gᵢ!(dx, perturb_tracer(xs, j, λ)..., p), dx, 0.0)
end
perturb_tracer(xs, j, λ) = (xs[1:j - 1]..., xs[j] .+ λ, xs[j + 1:end]...)

"""
    F, L, NL, ∇ₓF, ∇ₓL, ∇ₓNL, T = split_state_function_and_Jacobian(Ts, Ls, NLs, nb)

Returns the state function `F` and its jacobian, `∇ₓF`, as well as a collection of split operators. This is experimental. Use at your own risk!
"""
function split_state_function_and_Jacobian(Ts::Tuple, Ls::Tuple, NLs::Tuple, nb)
    nt = length(Ts)
    tracers(x) = state_to_tracers(x, nb, nt)
    T(p) = blockdiag([Tⱼ(p) for Tⱼ in Ts]...) # Big T (linear part)
    NL(x, p) = reduce(vcat, NLⱼ(tracers(x)..., p) for NLⱼ in NLs) # nonlinear part
    L(x, p) = reduce(vcat, Lⱼ(tracers(x)..., p) for Lⱼ in Ls) # nonlinear part
    F(x, p) = NL(x, p) + L(x, p) - T(p) * x                     # full 𝐹(𝑥) = -T 𝑥 + 𝐺(𝑥)
    ∇ₓNL(x, p) = local_jacobian(NLs, x, p, nt, nb)     # Jacobian of nonlinear part
    ∇ₓL(p) = local_jacobian(Ls, zeros(nt * nb), p, nt, nb)     # Jacobian of nonlinear part
    ∇ₓF(x, p) = ∇ₓNL(x, p) + ∇ₓL(p) - T(p)       # full Jacobian ∇ₓ𝐹(𝑥) = -T + ∇ₓ𝐺(𝑥)
    return F, L, NL, ∇ₓF, ∇ₓL, ∇ₓNL, T
end
function split_state_function_and_Jacobian(T, L, NL, nb)
    F(x, p) = NL(x, p) + L(x, p) - T(p) * x                     # full 𝐹(𝑥)
    ∇ₓNL(x, p) = sparse(Diagonal(localderivative(NL, x, p))) # Jacobian of nonlinear part
    ∇ₓL(p) = sparse(Diagonal(localderivative(L, zeros(nb), p)))     # Jacobian of nonlinear part
    ∇ₓF(x, p) = ∇ₓNL(x, p) + ∇ₓL(p) - T(p)       # full Jacobian ∇ₓ𝐹(𝑥) = -T + ∇ₓ𝐺(𝑥)
    return F, L, NL, ∇ₓF, ∇ₓL, ∇ₓNL, T
end
export split_state_function_and_Jacobian

function local_jacobian(Gs, x, p, nt, nb)
    return reduce(vcat, local_jacobian_row(Gⱼ, x, p, nt, nb) for Gⱼ in Gs)
end
function inplace_local_jacobian(Gs, x, p, nt, nb)
    return reduce(vcat, inplace_local_jacobian_row(Gⱼ!, x, p, nt, nb) for Gⱼ! in Gs)
end

function local_jacobian_row(Gᵢ, x, p, nt, nb)
    tracers(x) = state_to_tracers(x, nb, nt)
    return reduce(hcat, sparse(Diagonal(localderivative(Gᵢ, tracers(x), j, p))) for j in 1:nt)
end
function inplace_local_jacobian_row(Gᵢ!, x, p, nt, nb)
    tracers(x) = state_to_tracers(x, nb, nt)
    dx = Vector{Float64}(undef, nb)
    return reduce(hcat, sparse(Diagonal(localderivative(Gᵢ!, dx, tracers(x), j, p))) for j in 1:nt)
end

#= ============================================
Generate 𝑓 and derivatives from user input
============================================ =#

function generate_f(ωs, μx, σ²x, v, ωp, ::Type{T}) where {T <: APar}
    nt, nb = length(ωs), length(v)
    tracers(x) = state_to_tracers(x, nb, nt)
    f(x, λorp) = ωp * mismatch(T, λorp) +
        sum([ωⱼ * mismatch(xⱼ, μⱼ, σⱼ², v) for (ωⱼ, xⱼ, μⱼ, σⱼ²) in zip(ωs, tracers(x), μx, σ²x)])
    return f
end
function generate_f(ωs, ωp, grd, obs, ::Type{T}; kwargs...) where {T <: APar}
    nt, nb = length(ωs), count(iswet(grd))
    tracers(x) = state_to_tracers(x, nb, nt)
    Ms = [interpolationmatrix(grd, obsⱼ) for obsⱼ in obs]
    cs = get(kwargs, :cs, (collect(identity for i in 1:nt)...,))
    f(x, λorp) = ωp * mismatch(T, λorp) +
        sum([ωⱼ * mismatch(xⱼ, grd, obsⱼ, M=Mⱼ, c=cⱼ) for (ωⱼ, xⱼ, obsⱼ, Mⱼ, cⱼ) in zip(ωs, tracers(x), obs, Ms, cs)])
    return f
end
function generate_f(ωs, ωp, grd, modify::Function, obs, ::Type{T}) where {T <: APar}
    nt, nb = length(ωs), count(iswet(grd))
    Ms = [interpolationmatrix(grd, obsⱼ) for obsⱼ in obs]
    iwets = [iswet(grd, obsⱼ) for obsⱼ in obs]
    function f(x, λorp)
        xs = unpack_tracers(x, grd)
        return ωp * mismatch(T, λorp) + sum([ωᵢ * indirectmismatch(xs, grd, modify, obs, i, Mᵢ, iwetᵢ) for (i, (ωᵢ, Mᵢ, iwetᵢ)) in enumerate(zip(ωs, Ms, iwets))])
    end
    return f
end

function generate_∇ₓf(ωs, μx, σ²x, v)
    nt, nb = length(ωs), length(v)
    tracers(x) = state_to_tracers(x, nb, nt)
    ∇ₓf(x) = reduce(hcat, ωⱼ * ∇mismatch(xⱼ, μⱼ, σⱼ², v) for (ωⱼ, xⱼ, μⱼ, σⱼ²) in zip(ωs, tracers(x), μx, σ²x))
    ∇ₓf(x, p) = ∇ₓf(x)
    return ∇ₓf
end
function generate_∇ₓf(ωs, grd, obs; kwargs...)
    nt, nb = length(ωs), count(iswet(grd))
    tracers(x) = state_to_tracers(x, nb, nt)
    Ms = [interpolationmatrix(grd, obsⱼ) for obsⱼ in obs]
    cs = get(kwargs, :cs, (collect(identity for i in 1:nt)...,))
    ∇ₓf(x) = reduce(hcat, ωⱼ * ∇mismatch(xⱼ, grd, obsⱼ, M=Mⱼ, c=cⱼ) for (ωⱼ, xⱼ, obsⱼ, Mⱼ, cⱼ) in zip(ωs, tracers(x), obs, Ms, cs))
    ∇ₓf(x, p) = ∇ₓf(x)
    return ∇ₓf
end
function generate_∇ₓf(ωs, grd, modify::Function, obs)
    nt, nb = length(ωs), count(iswet(grd))
    Ms = [interpolationmatrix(grd, obsⱼ) for obsⱼ in obs]
    iwets = [iswet(grd, obsⱼ) for obsⱼ in obs]
    function ∇ₓf(x)
        xs = unpack_tracers(x, grd)
        sum([ωᵢ * ∇indirectmismatch(unpack_tracers(x, grd), grd, modify, obs, i, Mᵢ, iwetᵢ) for (i, (ωᵢ, Mᵢ, iwetᵢ)) in enumerate(zip(ωs, Ms, iwets))])
    end
    ∇ₓf(x, p) = ∇ₓf(x)
    return ∇ₓf
end

function f_and_∇ₓf(ωs, μx, σ²x, v, ωp, ::Type{T}) where {T <: APar}
    generate_f(ωs, μx, σ²x, v, ωp, T), generate_∇ₓf(ωs, μx, σ²x, v)
end
function f_and_∇ₓf(ωs, ωp, grd, obs, ::Type{T}; kwargs...) where {T <: APar}
    generate_f(ωs, ωp, grd, obs, T; kwargs...), generate_∇ₓf(ωs, grd, obs; kwargs...)
end
function f_and_∇ₓf(ωs, ωp, grd, modify::Function, obs, ::Type{T}) where {T <: APar}
    generate_f(ωs, ωp, grd, modify, obs, T), generate_∇ₓf(ωs, grd, modify, obs)
end
export f_and_∇ₓf

"""
    mismatch(x, xobs, σ²xobs, v)

Volume-weighted mismatch of modelled tracer `x` against observed mean, `xobs`, given observed variance, `σ²xobs`, and volumes `v`.
"""
function mismatch(x, xobs, σ²xobs, v)
    δx = x - xobs
    W = Diagonal(v ./ σ²xobs)
    return 0.5 * transpose(δx) * W * δx / (transpose(xobs) * W * xobs)
end
mismatch(x, ::Missing, args...) = 0

"""
    ∇mismatch(x, xobs, σ²xobs, v)

Adjoint of the gradient of `mismatch(x, xobs, σ²xobs, v)`.
"""
function ∇mismatch(x, xobs, σ²xobs, v)
    δx = x - xobs
    W = Diagonal(v ./ σ²xobs)
    return transpose(W * δx) / (transpose(xobs) * W * xobs)
end
∇mismatch(x, ::Missing, args...) = transpose(zeros(length(x)))

## new functions for more generic obs packages
# TODO Add an optional function argument to transform the data before computingn the mismatch
# Example if for isotope tracers X where one ususally wants to minimize the mismatch in δ or ε.
function mismatch(x, grd::OceanGrid, obs; c=identity, W=I, M=interpolationmatrix(grd, obs), iwet=iswet(grd, obs))
    o = view(obs, iwet)
    δx = M * c(x) - o
    return 0.5 * transpose(δx) * W * δx / (transpose(o) * W * o)
end
mismatch(x, grd::OceanGrid, ::Missing; kwargs...) = 0
function ∇mismatch(x, grd::OceanGrid, obs; c=identity, W=I, M=interpolationmatrix(grd, obs), iwet=iswet(grd, obs))
    ∇c = Diagonal(ForwardDiff.derivative(λ -> c(x .+ λ), 0.0))
    o = view(obs, iwet)
    δx = M * c(x) - o
    return transpose(W * δx) * M * ∇c / (transpose(o) * W * o)
end
∇mismatch(x, grd::OceanGrid, ::Missing; kwargs...) = transpose(zeros(length(x)))

# In case the mismatch is not based on the tracer but on some function of it
# TODO Add option for correlation matrix for Bayesian inferenece
function indirectmismatch(xs::Tuple, grd::OceanGrid, modify::Function, obs, i, M=interpolationmatrix(grd, obs[i]), iwet=iswet(grd, obs[i]))
    x2 = modify(xs...)
    out = 0.0
    o = obs[i][iwet, :value]
    δx = M * x2[i] - o
    return 0.5 * transpose(δx) * δx / (transpose(o) * o)
end
function ∇indirectmismatch(xs::Tuple, grd::OceanGrid, modify::Function, obs, i, M=interpolationmatrix(grd, obs[i]), iwet=iswet(grd, obs[i]))
    x2 = modify(xs...)
    o = obs[i][iwet, :value]
    δx = M * x2[i] - o
    ∇modᵢ = ∇modify(modify, xs, i)
    return transpose(δx) * M * ∇modᵢ / (transpose(o) * o)
end
# TODO think of more efficient way to avoid recomputing ∇modify whole for each i
function ∇modify(modify, xs, i, j)
    return sparse(Diagonal(ForwardDiff.derivative(λ -> modify(perturb_tracer(xs, j, λ)...)[i], 0.0)))
end
∇modify(modify, xs, i) = reduce(hcat, ∇modify(modify, xs, i, j) for j in 1:length(xs))

#= ============================================
multi-tracer norm
============================================ =#

function volumeweighted_norm(nt, v)
    w = repeat(v, nt)
    return nrm(x) = transpose(x) * Diagonal(w) * x
end

#= ============================================
unpacking of multi-tracers
============================================ =#

state_to_tracers(x, nb, nt) = ntuple(i -> state_to_tracer(x, nb, nt, i), nt)
state_to_tracer(x, nb, nt, i) = view(x, tracer_indices(nb, nt, i))
function state_to_tracers(x, grd)
    nb = number_of_wet_boxes(grd)
    nt = Int(round(length(x) / nb))
    return state_to_tracers(x, nb, nt)
end
tracer_indices(nb, nt, i) = (i - 1) * nb + 1:i * nb
tracers_to_state(xs) = reduce(vcat, xs)
export state_to_tracers, state_to_tracer, tracers_to_state, tracer_indices
# Alias for better name
unpack_tracers = state_to_tracers
export unpack_tracers


text:
type GraphLearningResult
    fileroot     :: AbstractString
    target_lat   :: AbstractFeature
    target_lon   :: AbstractFeature
    parents_lat  :: Vector{AbstractFeature}
    parents_lon  :: Vector{AbstractFeature}
    features     :: Vector{AbstractFeature}

    adj          :: BitMatrix               # this is of the size of the resulting network (ie, |f_inds|)
    stats        :: Vector{Matrix{Float64}} # NOTE(tim): this does not include prior counts
    bayescore    :: Float64

    function GraphLearningResult(
        basefolder     :: AbstractString,
        features       :: Vector{AbstractFeature},
        ind_lat        :: Int,
        ind_lon        :: Int,
        parentinds_lat :: Vector{Int}, # indices are given in terms of input features vector
        parentinds_lon :: Vector{Int},
        bayescore      :: Float64,
        r              :: AbstractVector{Int},
        d              :: AbstractMatrix{Int}
        )

        # 1 - build the set of features used in the net, ordered by [target_lat, target_lon, indicators]
        # 2 - pull the indeces for the parents of target_lat
        # 3 - pull the indeces for the parents of target_lon

        net_feature_indeces = _feature_indeces_in_net(ind_lat, ind_lon, parentinds_lat, parentinds_lon)
        net_features = features[net_feature_indeces]
        net_ind_lat = 1
        net_ind_lon = 2
        net_parent_indices_lat = _find_index_mapping(parentinds_lat, net_feature_indeces)
        net_parent_indices_lon = _find_index_mapping(parentinds_lon, net_feature_indeces)
        num_net_features = length(net_feature_indeces)

        # println("net_feature_indeces:    ", net_feature_indeces)
        # println("net_features:           ", net_features)
        # println("net_parent_indices_lat: ", net_parent_indices_lat)
        # println("net_parent_indices_lon: ", net_parent_indices_lon)
        # println("num_net_features:       ", num_net_features)

        adj   = construct_model_adjacency_matrix(net_ind_lat, net_ind_lon, net_parent_indices_lat, net_parent_indices_lon, num_net_features)
        stats = convert(Vector{Matrix{Float64}}, SmileExtra.statistics(adj, r[net_feature_indeces], d[net_feature_indeces,:]))

        # println("adj: ")
        # println(adj)

        target_lat = features[ind_lat]
        target_lon = features[ind_lon]

        # need to permute these appropriately

        parents_lat = features[net_parent_indices_lat]
        parents_lon = features[net_parent_indices_lon]

        new(basefolder, target_lat, target_lon, parents_lat, parents_lon,
            net_features, adj, stats, bayescore)
    end
end

function get_emstats(res::GraphLearningResult, binmapdict::Dict{Symbol, AbstractDiscretizer})
    emstats = Dict{AbstractString, Any}()
    emstats["binmaps"] = binmapdict
    emstats["targets"] = [res.target_lat, res.target_lon]
    emstats["indicators"] = res.features[3:end]
    emstats["statistics"] = res.stats
    emstats["adjacency"] = res.adj
    emstats
end
function get_emstats{D<:AbstractDiscretizer, F<:AbstractFeature}(res::GraphLearningResult, binmaps::Vector{D}, features::Vector{F})
    binmapdict = Dict{Symbol,AbstractDiscretizer}()
    for (b,f) in zip(binmaps, features)
        binmapdict[symbol(f)] = b
    end
    get_emstats(res, binmapdict)
end


text:
@enum RequestType TypeFollowers TypeFollowing TypeUserinfo

const USERAGENT = "Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_4) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/80.0.3987.149 Safari/537.36"
const HOST = "api.github.com"
const TIMEZONE = "Asia/Shanghai"

max_stars_repo_name: tosone/julia-spider

max_issues_repo_name: tosone/julia-spider

max_forks_repo_name: tosone/julia-spider


text:
using DiffEqBase, OrdinaryDiffEq
using StaticArrays: SVector, SizedArray
using RecursiveArrayTools: ArrayPartition

import DiffEqBase: ODEProblem, solve, ODE_DEFAULT_NORM, @.., addsteps!

import OrdinaryDiffEq: OrdinaryDiffEqAdaptiveAlgorithm,
OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache,
alg_order, alg_cache, initialize!, perform_step!, @unpack,
@cache, stepsize_controller!, step_accept_controller!

# TODO: check which keywords work fine
const warnkeywords = (:save_idxs, :d_discontinuities, :unstable_check, :save_everystep,
    :save_end, :initialize_save, :adaptive, :dt, :reltol, :dtmax,
    :dtmin, :force_dtmin, :internalnorm, :gamma, :beta1, :beta2,
    :qmax, :qmin, :qsteady_min, :qsteady_max, :qoldinit, :failfactor,
    :isoutofdomain, :unstable_check,
    :calck, :progress, :timeseries_steps, :dense)

global warnlist = Set(warnkeywords)

abstract type TaylorAlgorithm <: OrdinaryDiffEqAdaptiveAlgorithm end
struct TaylorMethod <: TaylorAlgorithm
    order::Int
    parse_eqs::Bool
end

TaylorMethod(order; parse_eqs=true) = TaylorMethod(order, parse_eqs) # set `parse_eqs` to `true` by default

alg_order(alg::TaylorMethod) = alg.order

TaylorMethod() = error("Maximum order must be specified for the Taylor method")

export TaylorMethod

# overload DiffEqBase.ODE_DEFAULT_NORM for Taylor1 arrays
ODE_DEFAULT_NORM(x::AbstractArray{Taylor1{T}, N},y) where {T<:Number, N} = norm(x, Inf)

### cache stuff
struct TaylorMethodCache{uType, rateType, tTType, uTType} <: OrdinaryDiffEqMutableCache
    u::uType
    uprev::uType
    tmp::uType
    k::rateType
    fsalfirst::rateType
    tT::tTType
    uT::uTType
    duT::uTType
    uauxT::uTType
    parse_eqs::Ref{Bool}
    tmpTaylor::uTType
    arrTaylor::Vector{uTType}
end

full_cache(c::TaylorMethodCache) = begin
    tuple(c.u, c.uprev, c.tmp, c.k, c.fsalfirst, c.tT, c.uT, c.duT, c.uauxT, c.parse_eqs,
        c.tmpTaylor, c.arrTaylor)
end

struct TaylorMethodConstantCache{uTType} <: OrdinaryDiffEqConstantCache
    uT::uTType
    parse_eqs::Ref{Bool}
    tmpTaylor::Vector{uTType}
    arrTaylor::Vector{Vector{uTType}}
end

function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p,
        calck,::Val{true})
    order = alg.order
    tT = Taylor1(typeof(t), order)
    tT[0] = t
    uT = Taylor1.(u, order)
    duT = zero.(Taylor1.(u, order))
    uauxT = similar(uT)
    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(alg.parse_eqs, f, tT, uT, duT, p)
    TaylorMethodCache(
        u,
        uprev,
        similar(u),
        zero(rate_prototype),
        zero(rate_prototype),
        tT,
        uT,
        duT,
        uauxT,
        Ref(parse_eqs),
        tmpTaylor,
        arrTaylor
        )
end

# This method is used for DynamicalODEFunction's (`parse_eqs=false`): tmpT1 and arrT1
# must have the proper type to build `TaylorMethodCache`
function alg_cache(alg::TaylorMethod, u::ArrayPartition, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p,
        calck, ::Val{true})
    order = alg.order
    tT = Taylor1(typeof(t), order)
    tT[0] = t
    uT = Taylor1.(u, order)
    duT = zero.(Taylor1.(u, order))
    uauxT = similar(uT)
    parse_eqs, _, _ = _determine_parsing!(alg.parse_eqs, f, tT, uT, duT, p)
    ## `tmpT1` must have the same type than `uT` and `arrT1` is a `[tmpT1]`
    tmpT1 = similar(uT)
    arrT1 = similar([tmpT1])
    TaylorMethodCache(
        u,
        uprev,
        similar(u),
        zero(rate_prototype),
        zero(rate_prototype),
        tT,
        uT,
        duT,
        uauxT,
        Ref(parse_eqs),
        tmpT1,
        arrT1
        )
end

function alg_cache(alg::TaylorMethod, u, rate_prototype, uEltypeNoUnits,
        uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, dt, reltol, p, calck,
        ::Val{false})
    order = alg.order
    tT = Taylor1(typeof(t), order)
    tT[0] = t
    uT = Taylor1(u, order)
    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(alg.parse_eqs, f, tT, uT, p)
    TaylorMethodConstantCache(Taylor1(u, alg.order), Ref(parse_eqs), tmpTaylor, arrTaylor)
end

function initialize!(integrator, c::TaylorMethodConstantCache)
    @unpack u, t, f, p = integrator
    tT = Taylor1(typeof(t), integrator.alg.order)
    tT[0] = t
    c.uT .= Taylor1(u, tT.order)
    __jetcoeffs!(Val(c.parse_eqs.x), f, tT, c.uT, p, c.tmpTaylor, c.arrTaylor)
    # FSAL stuff
    integrator.kshortsize = 2
    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
    integrator.destats.nf += 1
    # Avoid undefined entries if k is an array of arrays
    integrator.fsallast = zero(integrator.fsalfirst)
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator,cache::TaylorMethodConstantCache)
    @unpack u, t, dt, f, p = integrator
    tT = Taylor1(typeof(t), integrator.alg.order)
    tT[0] = t+dt
    u = evaluate(cache.uT, dt)
    cache.uT[0] = u
    __jetcoeffs!(Val(cache.parse_eqs.x), f, tT, cache.uT, p, cache.tmpTaylor, cache.arrTaylor)
    k = f(u, p, t+dt) # For the interpolation, needs k at the updated point
    integrator.destats.nf += 1
    integrator.fsallast = k
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.u = u
end

function initialize!(integrator, cache::TaylorMethodCache)
    @unpack u, t, f, p = integrator
    @unpack k, fsalfirst, tT, uT, duT, uauxT, parse_eqs = cache
    __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, cache.tmpTaylor, cache.arrTaylor)
    # FSAL for interpolation
    integrator.fsalfirst = fsalfirst
    integrator.fsallast = k
    integrator.kshortsize = 1
    resize!(integrator.k, integrator.kshortsize)
    integrator.k[1] = integrator.fsalfirst
    # integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t)
    integrator.fsalfirst = constant_term.(duT)
    integrator.destats.nf += 1
end

function perform_step!(integrator, cache::TaylorMethodCache)
    @unpack t, dt, u, f, p = integrator
    @unpack k, tT, uT, duT, uauxT, parse_eqs, tmpTaylor, arrTaylor = cache
    evaluate!(uT, dt, u)
    tT[0] = t+dt
    for i in eachindex(u)
        @inbounds uT[i][0] = u[i]
        duT[i].coeffs .= zero(duT[i][0])
    end
    __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, tmpTaylor, arrTaylor)
    k = constant_term.(duT) # For the interpolation, needs k at the updated point
    integrator.destats.nf += 1
end

stepsize_controller!(integrator,alg::TaylorMethod) =
    stepsize(integrator.cache.uT, integrator.opts.abstol)
step_accept_controller!(integrator, alg::TaylorMethod, q) = q

function DiffEqBase.solve(
        prob::DiffEqBase.AbstractODEProblem{uType, tupType, isinplace},
        alg::TaylorMethod, args...;
        verbose=true,
        kwargs...) where
        {uType, tupType, isinplace}

    if verbose
        warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
        warned && warn_compat()
    end

    f = prob.f
    parse_eqs = haskey(kwargs, :parse_eqs) ? kwargs[:parse_eqs] : true # `true` is the default
    if !isinplace && typeof(prob.u0) <: AbstractArray
        ### TODO: allow `parse_eqs=true` for DynamicalODEFunction
        if prob.f isa DynamicalODEFunction
            f1! = (dv, v, u, p, t) -> (dv .= prob.f.f1(v, u, p, t); 0)
            f2! = (du, v, u, p, t) -> (du .= prob.f.f2(v, u, p, t); 0)
            _alg = TaylorMethod(alg.order, parse_eqs = false)
            ### workaround use of `SVector` with oop `DynamicalODEProblem`
            ### TODO: add proper support for oop problems with arrays
            if eltype(prob.u0.x) <: SVector
                _u0 = ArrayPartition(SizedArray{Tuple{length(prob.u0.x[1])}}.(prob.u0.x))
            else
                _u0 = prob.u0
            end
            _prob = DynamicalODEProblem(f1!, f2!, _u0.x[1], _u0.x[2], prob.tspan, prob.p; prob.kwargs...)
        else
            f! = (du, u, p, t) -> (du .= f(u, p, t); 0)
            _alg = TaylorMethod(alg.order, parse_eqs = parse_eqs)
            _prob = ODEProblem(f!, prob.u0, prob.tspan, prob.p; prob.kwargs...)
        end
    else
        _alg = TaylorMethod(alg.order, parse_eqs = parse_eqs)
        _prob = prob
    end

    # DiffEqBase.solve(prob, _alg, args...; kwargs...)
    integrator = DiffEqBase.__init(_prob, _alg, args...; kwargs...)
    integrator.dt = stepsize(integrator.cache.uT, integrator.opts.abstol) # override handle_dt! setting of initial dt
    DiffEqBase.solve!(integrator)
    integrator.sol
end

# used in continuous callbacks and related methods to update Taylor expansions cache
function update_jetcoeffs_cache!(u,f,p,cache::TaylorMethodCache)
    @unpack tT, uT, duT, uauxT, parse_eqs, tmpTaylor, arrTaylor = cache
    for i in eachindex(u)
        @inbounds uT[i][0] = u[i]
        duT[i].coeffs .= zero(duT[i][0])
    end
    __jetcoeffs!(Val(parse_eqs.x), f, tT, uT, duT, uauxT, p, tmpTaylor, arrTaylor)
    return nothing
end

# This function was modified from OrdinaryDiffEq.jl; MIT-licensed
# DiffEqBase.addsteps! overload for ::TaylorMethodCache to handle continuous
# and vector callbacks with TaylorIntegration.jl via the common interface
function DiffEqBase.addsteps!(k, t, uprev, u, dt, f, p, cache::TaylorMethodCache,
        always_calc_begin = false, allow_calc_end = true,force_calc_end = false)
    if length(k)<2 || always_calc_begin
        if typeof(cache) <: OrdinaryDiffEqMutableCache
            rtmp = similar(u, eltype(eltype(k)))
            f(rtmp,uprev,p,t)
            copyat_or_push!(k,1,rtmp)
            f(rtmp,u,p,t+dt)
            copyat_or_push!(k,2,rtmp)
        else
            copyat_or_push!(k,1,f(uprev,p,t))
            copyat_or_push!(k,2,f(u,p,t+dt))
        end
    end
    update_jetcoeffs_cache!(u,f,p,cache)
    nothing
end

@inline __jetcoeffs!(::Val{false}, f::ODEFunction, t, x, params, tmpTaylor, arrTaylor) =
    __jetcoeffs!(Val(false), f.f, t, x, params, tmpTaylor, arrTaylor)
@inline __jetcoeffs!(::Val{true},  f::ODEFunction, t, x, params, tmpTaylor, arrTaylor) =
    __jetcoeffs!(Val(true), f.f, t, x, params, tmpTaylor, arrTaylor)
@inline __jetcoeffs!(::Val{false}, f::ODEFunction, t, x, dx, xaux, params, tmpTaylor, arrTaylor) =
    __jetcoeffs!(Val(false), f.f, t, x, dx, xaux, params, tmpTaylor, arrTaylor)
@inline __jetcoeffs!(::Val{true},  f::ODEFunction, t, x, dx, xaux, params, tmpTaylor, arrTaylor) =
    __jetcoeffs!(Val(true), f.f, t, x, dx, xaux, params, tmpTaylor, arrTaylor)

_determine_parsing!(parse_eqs::Bool, f::ODEFunction, t, x, params) =
    _determine_parsing!(parse_eqs, f.f, t, x, params)
_determine_parsing!(parse_eqs::Bool, f::ODEFunction, t, x, dx, params) =
    _determine_parsing!(parse_eqs, f.f, t, x, dx, params)


text:
using JuMP;
using Ipopt;
using DataFrames;

my_table0 = readtable( "C:/Users/mdzhd/Documents/partial_linear/lin_12_07_2016/newdata.csv");

my_table = convert( DataFrame, my_table0 );

( N, K1 ) = size( my_table );

m = Model(solver=IpoptSolver());

@defVar( m, 0.0 <= b[1:(K1-1)] <= 100.0 )

@NLobjective( m, Min, sum{ ( my_table[i,1] - sum{ my_table[i,j] * b[ j - 1] , j in 2:K1} )^2 , i in 1:N} )

statusM = solve( m )

df = DataFrame( beta = getvalue( b ) )

writetable( "C:/Users/mdzhd/Documents/partial_linear/lin_12_07_2016/solution_lin.csv", df)

# ----------------------- LOGISTIC --------------------------------

using JuMP, Ipopt, DataFrames
my_table = readtable( "/home/iriadmin/g/JStack/src/tst/newdata.csv");
( N, K1 ) = size( my_table );
m = Model(solver=IpoptSolver());
@defVar( m, 0.0 <= b[1:(K1-1)] <= 100.0 )
@NLobjective( m, Min, sum{ ( my_table[i,1] - sum{ my_table[i,j] * b[ j - 1] , j in 2:K1} )^2 , i in 1:N} )
statusM = solve( m )
df = DataFrame( beta = getvalue( b ) )

dfd = readtable( "/home/iriadmin/g/JStack/src/tst/logreg.csv")
dfd=dfd[setdiff(names(dfd),[:x0])]
mm = ModelMatrix(ModelFrame(y ~ 1 + x1 + x2 + x3 + x4 + x6 + x7, dfd))
g = glm(y ~ 1 + x1 + x2 + x3 + x4 + x6 + x7, dfd,  Bernoulli())   #m = glm(y ~ 1 + x1 + x2 + x3 + x4 + x6 + x7, dfd, Poisson())
a = convert(Array,dfd) #array(dfd[:, 1:end])
B = zeros(length(mm.m[1,2:end]))
x = a[:,2:end]
y = a[:,1]
# sum( -1 * y .* (x*beta) + log(1+exp((x*beta))) )

m = Model(solver=IpoptSolver());
@variable(m, B[i=1:6], start=0)

@defVar( m, 0.0 <= b[1:(K1-1)] <= 100.0 )
@NLobjective( m, Min, sum{ ( my_table[i,1] - sum{ my_table[i,j] * b[ j - 1] , j in 2:K1} )^2 , i in 1:N} )
statusM = solve( m )
df = DataFrame( beta = getvalue( b ) )

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using GeneralizedLeastSquares: coef, LinRegQR, LinRegCholesky, LinRegBunchKaufman, LinRegSweep
using LinearAlgebra
using Test

n, p = 10^6, 50

x = randn(n, p)

y = x * (1:50) + randn(n)

w = Diagonal(rand(n))

β̂ = x \ y

β̂2 = (x' * w * x) \ (x' * w * y)

@testset "OLS" begin
    @testset "LinRegQR" begin
        model = LinRegQR(x, y)
        @test coef(model) ≈ β̂
    end
    @testset "LinRegCholesky" begin
        model = LinRegCholesky(x, y)
        @test coef(model) ≈ β̂
    end
    @testset "LinRegBunchKaufman" begin
        model = LinRegBunchKaufman(x, y)
        @test coef(model) ≈ β̂
    end
    @testset "LinRegSweep" begin
        model = LinRegSweep(x, y)
        @test coef(model) ≈ β̂
    end
end
@testset "GLS" begin
    @testset "LinRegQR" begin
        wsqrt = sqrt(w)
        x2 = wsqrt * x
        y2 = wsqrt * y
        model = LinRegQR(x2, y2)
        @test coef(model) ≈ β̂2
    end
    @testset "LinRegCholesky" begin
        model = LinRegCholesky(x, y, w)
        @test coef(model) ≈ β̂2
    end
    @testset "LinRegBunchKaufman" begin
        model = LinRegBunchKaufman(x, y, w)
        @test coef(model) ≈ β̂2
    end
    @testset "LinRegSweep" begin
        model = LinRegSweep(x, y, w)
        @test coef(model) ≈ β̂2
    end
end


text:
using IncrementalInference
using Test

##

@testset "test default userdata::FactorMetadata..." begin

##

fgt = initfg()

addVariable!(fgt, :x1, ContinuousScalar)
addFactor!(fgt, [:x1], Prior(Normal()))

addVariable!(fgt, :x2, ContinuousScalar)
addFactor!(fgt, [:x1;:x2], LinearRelative(Normal(10,1)))

fc = DFG.getSolverData(getFactor(fgt, :x1x2f1))

@test length(IIF._getCCW(fc).cpt[1].factormetadata.fullvariables) == 2
@test IIF._getCCW(fc).cpt[1].factormetadata.solvefor == :x2

##

end

#