text:
import Libdl

const depsfile = joinpath(@__DIR__, "deps.jl")

function find_matlab_root()
    # Determine MATLAB library path and provide facilities to load libraries with this path
    matlab_root = get(ENV, "MATLAB_ROOT",
                        get(ENV, "MATLAB_HOME", nothing))
    if isnothing(matlab_root)
        matlab_exe = Sys.which("matlab")
        if !isnothing(matlab_exe)
            matlab_root = dirname(dirname(matlab_exe))
        else
            if Sys.isapple()
                default_dir = "/Applications"
                if isdir(default_dir)
                    dirs = readdir(default_dir)
                    filter!(app -> occursin(r"^MATLAB_R[0-9]+[ab]\.app$", app), dirs)
                    if !isempty(dirs)
                        matlab_root = joinpath(default_dir, maximum(dirs))
                    end
                end
            elseif Sys.iswindows()
                default_dir = Sys.WORD_SIZE == 32 ? "C:\\Program Files (x86)\\MATLAB" : "C:\\Program Files\\MATLAB"
                if isdir(default_dir)
                    dirs = readdir(default_dir)
                    filter!(dir -> occursin(r"^R[0-9]+[ab]$", dir), dirs)
                    if !isempty(dirs)
                        matlab_root = joinpath(default_dir, maximum(dirs))
                    end
                end
            end
        end
    end
    !isnothing(matlab_root) && isdir(matlab_root) && @info("Detected MATLAB root folder at \"$matlab_root\"")
    return matlab_root
end

function find_matlab_libpath(matlab_root)
    # get path to MATLAB libraries
    matlab_libdir = if Sys.islinux()
        Sys.WORD_SIZE == 32 ? "glnx86" : "glnxa64"
    elseif Sys.isapple()
        Sys.WORD_SIZE == 32 ? "maci" : "maci64"
    elseif Sys.iswindows()
        Sys.WORD_SIZE == 32 ? "win32" : "win64"
    end
    matlab_libpath = joinpath(matlab_root, "bin", matlab_libdir)
    isdir(matlab_libpath) && @info("Detected MATLAB library path at \"$matlab_libpath\"")
    return matlab_libpath
end

function find_matlab_cmd(matlab_root)
    if Sys.iswindows()
        matlab_cmd = joinpath(matlab_root, "bin", (Sys.WORD_SIZE == 32 ? "win32" : "win64"), "matlab.exe")
        isfile(matlab_cmd) && @info("Detected MATLAB executable at \"$matlab_cmd\"")
    else
        matlab_exe = joinpath(matlab_root, "bin", "matlab")
        isfile(matlab_exe) && @info("Detected MATLAB executable at \"$matlab_exe\"")
        matlab_cmd = "exec $(Base.shell_escape(matlab_exe))"
    end
    return matlab_cmd
end

matlab_root = find_matlab_root()

if !isnothing(matlab_root)
    matlab_libpath = find_matlab_libpath(matlab_root)
    matlab_cmd = find_matlab_cmd(matlab_root)
    libmx_size = filesize(Libdl.dlpath(joinpath(matlab_libpath, "libmx")))
    open(depsfile, "w") do io
        println(io,
            """
            # This file is automatically generated, do not edit.

            function check_deps()
                if libmx_size != filesize(Libdl.dlpath(joinpath(matlab_libpath, "libmx")))
                    error("MATLAB library has changed, re-run Pkg.build(\\\"MATLAB\\\")")
                end
            end
            """
            )
            println(io, "const matlab_libpath = \"$(escape_string(matlab_libpath))\"")
            println(io, "const matlab_cmd = \"$(escape_string(matlab_cmd))\"")
            println(io, "const libmx_size = $libmx_size")
    end
elseif get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", nothing) == "true" || get(ENV, "CI", nothing) == "true"
    # We need to be able to install and load this package without error for
    # Julia's registry AutoMerge to work, so we just use dummy values.
    # Similarly we want to also be able to install and load this package for CI.
    matlab_libpath = ""
    matlab_cmd = ""
    libmx_size = 0

    open(depsfile, "w") do io
        println(io,
            """
            # This file is automatically generated, do not edit.

            check_deps() = nothing
            """
            )
            println(io, "const matlab_libpath = \"$(escape_string(matlab_libpath))\"")
            println(io, "const matlab_cmd = \"$(escape_string(matlab_cmd))\"")
            println(io, "const libmx_size = $libmx_size")
    end
else
    error("MATLAB cannot be found. Set the \"MATLAB_ROOT\" environment variable to the MATLAB root directory and re-run Pkg.build(\"MATLAB\").")
end


text:
"""Note: Does not (currently) escape all of latex."""
function escape_latex(s::AbstractString)
    s #not actually escaping anything  infact
end

latex_brace(x::AbstractString) = "{" * x * "}"

to_latex(s::MathSymbolInner) = to_latex("$s") #Must define first so will be over writtern
#TODO: Special case Floats, to handle e
to_latex(s::AbstractString) = escape_latex(s)

to_latex(x::MathSymbol) = to_latex(x.sym)

function to_latex(x::MathSuperscript)
    to_latex(x.base) * "^" * latex_brace(to_latex(x.expo))
end

function to_latex(x::MathSubscript)
    to_latex(x.base) * "_" * latex_brace(to_latex(x.sub))
end

function to_latex(x::MathFrac)
    "\\frac" * latex_brace(to_latex(x.numer)) * latex_brace(to_latex(x.denom))
end

function to_latex(x::MathFunctionCall)
    name = x.fname |> string |> escape_latex
    if length(name) > 1
        name = "\\mathrm"*latex_brace(name) #Multicharacter names are typeset upright
    end
    name*to_latex(x.arg)
end

function to_latex(x::MathMatrix)
    desc = "c"^size(x.mat,1)
    rows = mapslices(y->join(y," & "), to_latex.(x.mat),2) |> vec
    latexmat = join(rows, "\\\\")

    """\\left[\\begin{array}{$desc}
    $latexmat
    \\end{array}\\right]"""
end

function to_latex(x::MathSequence)
    inner = join(to_latex.(x.fragments), " ")
    if length(x) == 1 && typeof(first(x.fragments))==MathSequence
       # skip outer-most brackets on Sequence that will be breackets internally anyway
       x.fragments |> first |> to_latex
    else
        "(" * join(to_latex.(x.fragments)," ") * ")"
    end
end

ext: jl


text:
using CPIMC, DataStructures
import CPIMC: orbs, adjacent_kinks_affecting_orbs, kinks_affecting_orbs, τ_borders, isunaffected, time_ordered_orbs, occupations_at, longest_type_1_chain_length, right_type_1_count
using CPIMC.PlaneWaves

a = PlaneWave((-2,0,0))
b = PlaneWave((3,0,0))
c = PlaneWave((0,0,0))
d = PlaneWave((1,0,0))
e = PlaneWave((5,9,9))

sd = Kinks( ImgTime(0.2) => T4(a,b,c,d),
           ImgTime(0.5) => T4(c,d,a,b),
           ImgTime(0.6) => T4(b,a,d,c),
           ImgTime(0.8) => T4(d,c,b,a) )

conf = Configuration(sphere(PlaneWave((0,0,0),Up),dk=1),sd)

@testset "flip(o::PlaneWave)" begin
    @test flip(PlaneWave((1,0,0),Down)) == PlaneWave((1,0,0),Up)
    x = PlaneWave((0,),Down)
    @test flip(flip(x)) == x
end

@testset "dimension(o::PlaneWave{D}) where {D}" begin
    @test dimension(a) == 3
    @test dimension(PlaneWave((1,))) == 1
    @test dimension(conf.occupations) == 3
end

@testset "fractional_spin_polarization(occ::Set{PlaneWave)" begin
    # equal number of Spin Up and Spin Down orbitals
    occ = Set([ PlaneWave((1,0,0),Down), PlaneWave((2,0,0),Up),
                PlaneWave((3,0,0),Down), PlaneWave((4,0,0),Up) ])
    @test fractional_spin_polarization(occ) == 0
    # 1 Spin Up and 3 Spin Down orbitals
    occ = Set([ PlaneWave((1,0,0),Down), PlaneWave((2,0,0),Down),
                PlaneWave((3,0,0),Down), PlaneWave((4,0,0),Up) ])
    @test fractional_spin_polarization(occ) == 0.5
    @test fractional_spin_polarization(occ) == fractional_spin_polarization(Set([flip(o) for o in occ]))
end

@testset "fractional_spin_polarization(c::Configuration{PlaneWave})" begin
    @test fractional_spin_polarization(conf) == 0
end

@testset "find_fourth_orb_for_kink" for _ in (1:5)
    r = (1:100)
    orb1 = PlaneWave((rand(r),rand(r),rand(r)),rand([Up,Down]))
    orb2 = orb1
    while orb1 == orb2
        orb2 = PlaneWave((rand(r),rand(r),rand(r)),rand([Up,Down]))
    end
    orb3 = orb1
    while (orb3 == orb2) || (orb3 == orb1)
        orb3 = PlaneWave((rand(r),rand(r),rand(r)),rand([orb1.spin,orb2.spin]))
    end
    @test in(find_fourth_orb_for_kink(orb3, orb1, orb2).spin, [orb1.spin, orb2.spin])
    @test (find_fourth_orb_for_kink(orb3, orb1, orb2).spin == orb3.spin) == (orb1.spin == orb2.spin)
    @test find_fourth_orb_for_kink(orb3, orb1, orb2).vec + orb3.vec == orb1.vec + orb2.vec
    orb4 =  find_fourth_orb_for_kink(orb3, orb1, orb2)
    @test find_fourth_orb_for_kink(orb3, orb1, orb2) == find_fourth_orb_for_kink(orb3, orb3, orb4)
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Test tight bound on GINI impurity index with remaining
"""
module GINIOptimTest2

export testmain

using Distributions #mode
using Reexport
@reexport using DataFrames
using RLESUtils, MathUtils, BitCounters

function testmain(ntrials=50; d=3)
  results = DataFrame([Int64,Int64,Int64,Int64,Int64,Int64,Int64,Int64,Int64,
                       Int64,Int64,Int64,Int64,Int64,Int64,Float64,
                       Int64,Int64,Int64,Int64,Int64,Int64,Float64,
                       Int64,Int64,Int64,Int64,Int64,Int64,Float64],
                      [:Ca, :Cb, :Cc, :n1a, :n1b, :n1c, :n2a, :n2b,:n2c,
                       :samp_u1a, :samp_u1b, :samp_u1c, :samp_u2a, :samp_u2b, :samp_u2c, :samp_score,
                       :allc_u1a, :allc_u1b, :allc_u1c, :allc_u2a, :allc_u2b, :allc_u2c, :allc_score,
                       :heur_u1a, :heur_u1b, :heur_u1c, :heur_u2a, :heur_u2b, :heur_u2c, :heur_score], 0)
  for i = 1:ntrials
    C = rand(1:10, d)
    n1 = rand(1:10, d)
    n2 = rand(1:10, d)

    samp_u = samplingbest(C, n1, n2)
    samp_score = evalscore(n1, n2, samp_u...)

    allc_u = allcombos(C, n1, n2)
    allc_score = evalscore(n1, n2, allc_u...)

    heur_u = heuristic1(C, n1, n2)
    heur_score = evalscore(n1, n2, heur_u...)

    trial = [C..., n1..., n2...,
             samp_u..., samp_score,
             allc_u..., allc_score,
             heur_u..., heur_score]
    push!(results, trial)
  end

  results
end

function evalscore(n1, n2, u1, u2)
  gini_from_counts(n1 + u1, n2 + u2)
end

function samplingbest(C, n1, n2; nsamples=10000)
  u1 = u1_ = similar(C)
  u2 = u2_ = similar(C)
  score = 1.0

  #draw samples
  for i = 1:nsamples
    for j = 1:length(C)
      u1_[j] = rand(0:C[j])
      u2_[j] = C[j] - u1_[j]
    end
    score_ = evalscore(n1, n2, u1_, u2_)

    if score_ < score
      score = score_
      u1 = u1_
      u2 = u2_
    end
  end

  u1, u2
end

function allcombos(C, n1, n2)
  ncases = 2^(length(C))
  case_u1 = Array(Vector, ncases)
  case_u2 = Array(Vector, ncases)
  case_score = Array(Float64, ncases)
  for (i, v) in enumerate(BitCounter(length(C)))
    case_u1[i] = v .* C
    case_u2[i] = !v .* C
    case_score[i] = evalscore(n1, n2, case_u1[i], case_u2[i])
  end

  i = indmin(case_score)
  return case_u1[i], case_u2[i]
end

#apply "remaining" bound greedily based on descending order on C[i]
function heuristic1(C, n1, n2)
  u1 = zeros(Int64, length(C))
  u2 = zeros(Int64, length(C))
  n1_ = deepcopy(n1)
  n2_ = deepcopy(n2)
  C_ = deepcopy(C)

  for i in sortperm(C_)
    j = evalcriteria(n1_, n2_, i, C_[i])
    if j == 1
      n1_[i] += C_[i]
      u1[i] += C_[i]
    else #j==2
      n2_[i] += C_[i]
      u2[i] += C_[i]
    end
    C_[i] = 0
  end
  u1, u2
end

function evalcriteria(c1, c2, i, N)
  n1 = sum(c1)
  n2 = sum(c2)

  c1_ = deepcopy(c1)
  c1_[i] = 0 #don't count the i of interest
  c2_ = deepcopy(c2)
  c2_[i] = 0 #don't count the i of interest

  b1 = (sum(c1_)^2 + sumabs2(c1_)) / (n1 * (n1 + N))
  b2 = (sum(c2_)^2 + sumabs2(c2_)) / (n2 * (n2 + N))

  if b1 < b2
    return 1
  else
    return 2
  end
end

end #module


text:
# # Loading Data into Comrade

# The VLBI field does not have a standarized data format, and the EHT uses a
# particular uvfits format that is similar to the optical interferometry oifits format.
# As a result, we reuse the excellent `eht-imaging` package to load data into `Comrade`.

# Once the data is loaded we then load the data into specific data formats that `Comrade`
# expects. Note that in the future this may change to a Julia package as the Julia radio
# astronomy group grows.

# To get started we will load `Comrade` and `Plots` to enable visualizations of the data

using Comrade
using Plots

# Now we load ehtim. This assumes you have a working installation of `eht-imaging`.
# To install eht-imaging see the [ehtim](https://github.com/achael/ehtim) github repo.

load_ehtim()

# Now we load the data. We will use the 2017 public M87 data which can be downloaded from
# [cyverse](https://datacommons.cyverse.org/browse/iplant/home/shared/commons_repo/curated/EHTC_FirstM87Results_Apr2019)

obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "../assets/SR1_M87_2017_096_lo_hops_netcal_StokesI.uvfits"))
# Add scan and coherently average over them. The eht data has been phase calibrated so that
# this is fine to do.
obs.add_scans()
obs = obs.avg_coherent(0.0, scan_avg=true)

# We can now extract data products that `Comrade` can use
vis = extract_vis(obs) #complex visibilites
amp = extract_amp(obs) # visibility amplitudes
cphase = extract_cphase(obs; count="min") # extract minimal set of closure phases
lcamp = extract_lcamp(obs; count="min") # extract minimal set of log-closure amplitudes

# We can also recover the array used in the observation using
ac = arrayconfig(vis)
plot(ac) # Plot the baseline coverage

# To plot the data we just call

l = @layout [ a b; c d]
pv = plot(vis)
pa = plot(amp)
pcp = plot(cphase)
plc = plot(lcamp)

plot(pv, pa, pcp, plc; layout=l)


text:
"""
    raytracing(fragments;args)

Returns a list of fragments and a list of rays.
The rays represent the fragments trajectories and can be used to intersect targets.
Optional argument is the v_threshold. Which is a threshold velocity used to discard fragments
for which the velocity is lower. The default value is 0.0.
"""
function raytracing(fragments::Vector;v_threshold=0.0)
    #v_threshold = 500.0
    fr =[]
    ray =[]
    for fragment in fragments
    #number=number + length(fragment)
        for f in fragment
            if norm(f.velocity) > v_threshold
                push!(fr,f)
                #r = Ray(Point(f.position),Point(f.velocity)-Point(0.0,0.0,0.0))
                r =OpticSim.Ray(OpticSim.Vec3(f.position),OpticSim.Vec3(f.velocity))
                push!(ray,r)

            end

        end
    end
    return fr,ray
end


text:
using HTTP: request
using LazyJSON: value
using SQLite: DB, Query, drop!
using Gumbo: parsehtml

function parse_license(license)
    if occursin("Apache", license)
        output = "Apache"
    elseif occursin("MIT", license)
        output = "MIT"
    elseif occursin("BSD", license)
        output = "BSD"
    elseif occursin(r"(CPL|EPL|Common Public)", license)
        output = "EPL"
    elseif occursin("Artistic", license)
        output = "Artistic"
    elseif occursin("CeCILL", license)
        output = "CECILL"
    elseif occursin(r"(Mozilla|MPL)", license)
        output = "MPL"
    elseif occursin(r"((?<!A)GPL|General Public)", license)
        output = "GPL"
    elseif occursin(r"(AGPL|Affero)", license)
        output = "AGPL"
    elseif occursin("EUPL", license)
        output = "EUPL"
    elseif occursin("Lucent", license)
        output = "LPL"
    elseif occursin("Unlimited", license)
        output = "FSFUL"
    elseif occursin("BSL", license)
        output = "BSL"
    elseif occursin(r"(CC|Creative)", license)
        output = "CC"
    elseif occursin("ACM", license)
        output = "ACM"
    elseif occursin(r"<a href=\"LICEN[CS]E\">", license)
        output = "Custom"
    elseif isempty(license)
        output = "GPL"
    else
        output = license
    end
    output
end
osi(license) = license ∈ ["Apache", "MIT", "BSD", "EPL", "Artistic", "CECILL",
                          "MPL", "GPL", "AGPL", "EUPL", "LPL", "BSL"]

function CRAN!()
    # Connect to the database
    db = DB("data/OSS.db")
    drop!(db, "cran_pkg", ifexists = true)
    drop!(db, "cran_contributors", ifexists = true)
    drop!(db, "cran_deps", ifexists = true)
    Query(db, "create table cran_pkg (
               pkg varchar not null,
               version varchar not null,
               license varchar not null,
               osi bool not null,
               slug varchar,
               primary key (pkg)
               );")
    Query(db, "create table cran_contributors (
               pkg varchar not null,
               name varchar not null,
               author bool not null,
               compiler bool not null,
               contributor bool not null,
               copyright_holder bool not null,
               creator bool not null,
               thesis_advisor bool not null,
               translator bool not null,
               foreign key (pkg) references cran_pkg(pkg)
               );")
    Query(db, "create table cran_deps (
               pkg varchar,
               deps varchar,
               foreign key (pkg) references cran_pkg(pkg)
               );")
    # Timestamp: 2012-12-17
    # This code is to query available packages in CRAN
    # response = request("GET", "https://cran.r-project.org/web/packages/available_packages_by_name.html")
    # response.status == 200 || throw(ArgumentError("Could not connect to CRAN"))
    # body = String(response.body)
    # pkgs = eachmatch(r"(?<=\.{2}/\.{2}/web/packages/).*(?=/index\.html)",
    #                  body) |>
    #        (x -> getproperty.(x, :match))
    # write("data/CRAN/_index.txt", join(pkgs, '\n'))
    # Read available packages in CRAN
    pkgs = "data/CRAN/_index.txt" |> read |> String |> (x -> split(x, '\n'))
    # Downloads and save the package information
    for pkg ∈ pkgs
        file = "data/CRAN/$pkg.html"
        if ~isfile(file) # Change to true when updating
            response = request("GET", "https://cran.r-project.org/web/packages/$pkg/index.html")
            String(response.body) |> (x -> write(file, x))
        end
    end
    for pkg ∈ pkgs
        html = "data/CRAN/$pkg.html" |> read |> String |> parsehtml
        vals = string.(getindex.(getindex.(html.root[2][3][1][:], 1), 1))
        version = findfirst(isequal("Version:"), vals) |>
            (x -> html.root[2][3][1][x][2]) |>
            string |>
            (x -> match(r"\p{N}+(\.\p{N}+){0,2}", x).match) |>
            VersionNumber
        license = findfirst(isequal("License:"), vals) |>
            (x -> html.root[2][3][1][x][2]) |>
            string |>
            (x -> eachmatch(r"(?<=>).+?(?=</)", x)) |>
            (x -> getproperty.(x, :match)) |>
            (x -> join(x, ',')) |>
            parse_license
        osi_approved = osi(license)
        dependency = Vector{String}()
        deps = findfirst(isequal("Depends:"), vals)
        imps = findfirst(isequal("Imports:"), vals)
        if ~isa(deps, Nothing)
            html.root[2][3][1][deps][2] |>
            string |>
            (x -> replace(x, r"(<.*?>|\(.*?\))" => "")) |>
            (x -> split(x, ',')) |>
            (x -> strip.(x)) |>
            (x -> append!(dependency, x))
        end
        if ~isa(imps, Nothing)
            html.root[2][3][1][imps][2] |>
            string |>
            (x -> replace(x, r"(<.*?>|\n)" => "")) |>
            (x -> replace(x, r"\(.*?\)" => "")) |>
            (x -> split(x, ',')) |>
            (x -> strip.(x)) |>
            (x -> append!(dependency, x))
        end
        filter!(!isequal("R"), dependency) |> sort!
        author = findfirst(isequal("Author:"), vals)
        authors = html.root[2][3][1][author][2][1] |>
            string |>
            (x -> replace(x, r"<.*?>" => "")) |>
            (x -> split(x, r",(?!\s*(aut|com|ctb|cph|cre|ths|trl)[,\]])\s*"))
        maintainer = findfirst(isequal("Maintainer:"), vals)
        maintainers = html.root[2][3][1][maintainer][2] |>
            string |>
            (x -> replace(x, r"</*td>" => "")) |>
            (x -> replace(x, " at " => "@"))
        email = match(r"(?<=<).*(?=>)", maintainers) |>
            (x -> isa(x, Nothing) ? "null" : x.match)
        slug = "null"
        bugs = findfirst(isequal("BugReports:"), vals)
        if ~isa(bugs, Nothing)
            link = html.root[2][3][1][bugs][2] |>
                string |>
                (x -> replace(x, r"/+" => "/")) |>
                (x -> match(r"(?<=github\.com/).*?/.*?(?=(/|\"))", x))
            if isa(link, RegexMatch)
                slug = "'$(link.match)'"
            end
        end
        URL = findfirst(isequal("URL:"), vals)
        if ~isa(URL, Nothing)
            link = html.root[2][3][1][URL][2] |>
                string |>
                (x -> replace(x, r"/+" => "/")) |>
                (x -> match(r"(?<=github\.com/).*?/.*?(?=(/|\"))", x))
            if isa(link, RegexMatch)
                slug = "'$(link.match)'"
            end
        end
        Query(db, "insert into cran_pkg values (" *
                  "'$pkg', '$version', '$license', $osi_approved, $slug)")
        for author ∈ authors
            # author = authors[1]
            name = match(r".*?(?=(\[|$))", author).match |>
                (x -> strip(x, ['\'',' '])) |>
                (x -> replace(x, "'" => "''"))
            role = match(r"(?<=\[).*(?=\])", author)
            if isa(role, Nothing)
                author = compiler = contributor = copyright_holder = creator =
                thesis_advisor = translator = false
            else
                role = role.match
                author = occursin("aut", role)
                compiler = occursin("com", role)
                contributor = occursin("ctb", role)
                copyright_holder = occursin("copyright_holder", role)
                creator = occursin("cre", role)
                thesis_advisor = occursin("ths", role)
                translator = occursin("trl", role)
            end
            Query(db, "insert into cran_contributors values (
                       '$pkg', '$name', $author, $compiler,
                       $contributor, $copyright_holder, $creator,
                       $thesis_advisor, $translator)")
        end
        for deps ∈ dependency
            Query(db, "insert into cran_deps values (" *
                      "'$pkg', '$deps')")
        end
    end
end
CRAN!()


text:
"""
`Maker.Targets`

An Associative type for storing `AbstractTargets`. A simple wrapper
of a Dict to allow more suitable display of the set of AbstractTargets.
"""
immutable Targets <: Associative
    dict::Dict{UTF8String, AbstractTarget}
end

Base.length(x::Targets) = length(x.dict)
Base.keys(x::Targets) = keys(x.dict)
Base.getindex(x::Targets, key) = getindex(x.dict, key)
Base.get(x::Targets, key, default) = get(x.dict, key, default)
Base.setindex!(x::Targets, key, val) = setindex!(x.dict, key, val)
Base.empty!(x::Targets) = empty!(x.dict)
Base.start(x::Targets) = start(x.dict)
Base.done(x::Targets, i) = done(x.dict, i)
Base.next(x::Targets, i) = next(x.dict, i)

"""
```julia
Maker.TARGETS
```

The global `Targets` dictionary for registered AbstractTargets.
"""
const TARGETS = Targets(Dict{UTF8String, AbstractTarget}())
NEXTDOC = UTF8String[""]

"""
```julia
@desc
```

Define a description for the next Maker task.
"""
macro desc(arg)
    :(Maker.NEXTDOC[1] = $arg)
end

"""
```julia
tasks()
tasks(name)
```

Return the target `name` or return all targets. Often used to show summaries of
all registered targets or a specific target.
"""
tasks() = TARGETS

tasks(name) = TARGETS.dict[name]

"""
```julia
tasks!(name, target)
```

Add the `target` task with name `name` to the global database of targets.
"""
function tasks!{T <: AbstractTarget}(name, val::T)
    TARGETS[name] = val
end

function funloc(f::Function)
    if isgeneric(f)
        a = functionloc(f, has1arg(f) ? (AbstractTarget,) : ())
        return string(a[1], " line ", a[2])
    else
        return funloc(Base.uncompressed_ast(f.code))
    end
end

function funloc(e::Expr)
    for a in e.args[3].args
         if isa(a, LineNumberNode)
             return string(a.file, " line ", a.line)
         end
   end
   return ""
end

function Base.show{T <: AbstractTarget}(io::IO, t::T)
    print(io, "\"", t.name, "\"")
    println(io, " => ", replace(string(t.dependencies), "UTF8String", ""))
    println(io, T)
    !isempty(t.description) && println(io, t.description)
    println(io, "Timestamp: ", t.timestamp, " UTC")
    println(io, "Action defined at: ", funloc(t.action))
    t.isstale && println(io, "STALE task")
end

function markdown(X::Targets)
    ks = collect(keys(X))
    md = """
         task | dependencies | description
         :----|:-------------|:-----------
         """
    for k in sort(ks)
        t = X[k]
        md *= "\"" * t.name * "\" | " *
              replace(string(t.dependencies), "UTF8String", "") *
              " | " * t.description * "\n"
    end
    Base.Markdown.parse(md)
end

function Base.show(io::IO, X::Targets)
    writemime(io, MIME"text/plain"(), markdown(X))
end

function Base.writemime(io::IO, T::MIME"text/plain", X::Targets)
    writemime(io, T, markdown(X))
end

function Base.display(d::Base.REPL.REPLDisplay, X::Targets)
    Base.Markdown.term(Base.REPL.outstream(d.repl), markdown(X))
end


text:
using Base.Dates
using Base.Test

using MarketData

using TimeSeries

@testset "apply" begin

@testset "time series methods" begin
    @testset "lag takes previous day and timestamps it to next day" begin
        @test isapprox(lag(cl).values[1]   , 111.94, atol=.01)
        @test lag(cl).timestamp[1] == Date(2000,1,4)
    end

    @testset "lag accepts other offset values" begin
        @test lag(cl, 9).timestamp[1] == Date(2000,1,14)
    end

    @testset "lag operates on 2d arrays" begin
        @test lag(ohlc, 9).timestamp[1] == Date(2000,1,14)
    end

    @testset "lag returns 1d from 1d time arrays" begin
        @test ndims(lag(cl).values) == 1
    end

    @testset "lag returns 2d from 2d time arrays" begin
        @test ndims(lag(ohlc).values) == 2
    end

    @testset "lead takes next day and timestamps it to current day" begin
        @test isapprox(lead(cl).values[1]   , 102.5, atol=.1)
        @test lead(cl).timestamp[1] == Date(2000,1,3)
    end

    @testset "lead accepts other offset values" begin
        @test lead(cl, 9).values[1]    == 100.44
        @test lead(cl, 9).timestamp[1] == Date(2000,1,3)
    end

    @testset "lead operates on 2d arrays" begin
        @test lead(ohlc, 9).timestamp[1] == Date(2000,1,3)
    end

    @testset "lead returns 1d from 1d time arrays" begin
        @test ndims(lead(cl).values) == 1
    end

    @testset "lead returns 2d from 2d time arrays" begin
        @test ndims(lead(ohlc).values) == 2
    end

    @testset "diff calculates 1st-order differences" begin
        @test diff(op).timestamp                              == diff(op, padding=false).timestamp
        @test diff(op).values                                 == diff(op, padding=false).values
        @test diff(op, padding=false).values[1]               == op[2].values[1] .- op[1].values[1]
        @test isequal(diff(op, padding=true).values[1], NaN)  == true
        @test diff(op, padding=true).values[2]                == diff(op).values[1]
        @test diff(op, padding=true).values[2]                == op[2].values[1] .- op[1].values[1]
    end

    @testset "simple return value" begin
        @test percentchange(cl, :simple).values == percentchange(cl).values
        @test percentchange(cl).values          == percentchange(cl, padding=false).values
        @test isapprox(percentchange(cl).values[1]                   , (102.5-111.94)/111.94, atol=.01)
        @test isapprox(percentchange(ohlc).values[1, :]              , (ohlc.values[2,:] - ohlc.values[1,:]) ./ ohlc.values[1,:])
        @test isnan.(percentchange(cl, padding=true).values[1])
        @test isapprox(percentchange(cl, padding=true).values[2]     , (102.5-111.94)/111.94, atol=.01)
        @test isapprox(percentchange(ohlc, padding=true).values[2, :], (ohlc.values[2,:] - ohlc.values[1,:]) ./ ohlc.values[1,:])
    end

    @testset "log return value" begin
        @test percentchange(cl, :log).values == percentchange(cl, :log, padding=false).values
        @test isapprox(percentchange(cl, :log).values[1]                   , log(102.5) - log(111.94), atol=.01)
        @test isapprox(percentchange(ohlc, :log).values[1, :]              , log.(ohlc.values[2,:]) .- log.(ohlc.values[1,:]), atol=.01)
        @test isnan.(percentchange(cl, :log, padding=true).values[1])
        @test isapprox(percentchange(cl, :log, padding=true).values[2]     , log(102.5) - log(111.94))
        @test isapprox(percentchange(ohlc, :log, padding=true).values[2, :], log.(ohlc.values[2,:]) .- log.(ohlc.values[1,:]))
    end

    @testset "moving supplies correct window length" begin
        @test moving(cl, mean, 10).values                                == moving(cl, mean, 10, padding=false).values
        @test moving(cl, mean, 10).timestamp[1]                          == Date(2000,1,14)
        @test isapprox(moving(cl, mean, 10).values[1], mean(cl.values[1:10]))
        @test moving(cl, mean, 10, padding=true).timestamp[1]            == Date(2000,1, 3)
        @test moving(cl, mean, 10, padding=true).timestamp[10]           == Date(2000,1,14)
        @test isequal(moving(cl, mean, 10, padding=true).values[1], NaN) == true
        @test moving(cl, mean, 10, padding=true).values[10]              == moving(cl, mean, 10).values[1]
        @test moving(ohlc, mean, 10).values                              == moving(ohlc, mean, 10, padding=false).values
        @test isapprox(moving(ohlc, mean, 10).values[1, :]', mean(ohlc.values[1:10, :], 1))
        @test isequal(moving(ohlc, mean, 10, padding=true).values[1, :], [NaN, NaN, NaN, NaN]) == true
        @test moving(ohlc, mean, 10, padding=true).values[10, :]         == moving(ohlc, mean, 10).values[1, :]
    end

    @testset "upto method accumulates" begin
        @test isapprox(upto(cl, sum).values[10]       , sum(cl.values[1:10]))
        @test isapprox(upto(cl, mean).values[10]      , mean(cl.values[1:10]))
        @test upto(cl, sum).timestamp[10] == Date(2000,1,14)
        # transpose the upto value output from column to row vector but values are identical
        @test isapprox(upto(ohlc, sum).values[10, :]' , sum(ohlc.values[1:10, :], 1))
        @test isapprox(upto(ohlc, mean).values[10, :]', mean(ohlc.values[1:10, :], 1))
    end
end

@testset "base element-wise operators on TimeArray values" begin
    @testset "only values on intersecting Dates computed" begin
        @test isapprox((cl[1:2] ./ op[2:3]).values[1], 0.94688222)
        @test isapprox((cl[1:4] .+ op[4:7]).values[1], 201.12, atol=.01)
        @test length(cl[1:2] ./ op[2:3]) == 1
        @test length(cl[1:4] .+ op[4:7]) == 1
    end

    @testset "unary operation on TimeArray values" begin
        @test (+cl).values[1]                  == cl.values[1]
        @test (-cl).values[1]                  == -cl.values[1]
        @test (.!(cl .== op)).values[1]        == true
        @test log.(cl).values[1]               == log.(cl.values[1])
        @test sqrt.(cl).values[1]              == sqrt.(cl.values[1])
        @test (+ohlc).values[1,:]              == ohlc.values[1,:]
        @test (-ohlc).values[1,:]              == -(ohlc.values[1,:])
        @test (.!(ohlc .== ohlc)).values[1, 1] == false
        @test log.(ohlc).values[1, :]          == log.(ohlc.values[1, :])
        @test sqrt.(ohlc).values[1, :]         == sqrt.(ohlc.values[1, :])
        @test_throws DomainError sqrt.(-ohlc)
    end

    @testset "dot operation between TimeArray values and Int/Float and viceversa" begin
        @test isapprox((cl .- 100).values[1]    , 11.94, atol=.01)
        @test isapprox((cl .+ 100).values[1]    , 211.94, atol=.01)
        @test isapprox((cl .* 100).values[1]    , 11194, atol=1)
        @test isapprox((cl ./ 100).values[1]    , 1.1194, atol=.0001)
        @test isapprox((cl .^ 2).values[1]      , 12530.5636, atol=.0001)
        @test (cl .% 2).values[1] == cl.values[1] % 2
        @test isapprox((100 .- cl).values[1]    , -11.94, atol=.01)
        @test isapprox((100 .+ cl).values[1]    , 211.94, atol=.01)
        @test isapprox((100 .* cl).values[1]    , 11194, atol=.01)
        @test isapprox((100 ./ cl).values[1]    , 0.8933357155619082)
        @test (2 .^ cl).values[1]       == 4980784073277740581384811358191616
        @test (2 .% cl).values[1]       == 2
        @test (ohlc .- 100).values[1,:] == ohlc.values[1,:] .- 100
        @test (ohlc .+ 100).values[1,:] == ohlc.values[1,:] .+ 100
        @test (ohlc .* 100).values[1,:] == ohlc.values[1,:] .* 100
        @test (ohlc ./ 100).values[1,:] == ohlc.values[1,:] ./ 100
        @test (ohlc .^ 2).values[1,:]   == ohlc.values[1,:] .^ 2
        @test (ohlc .% 2).values[1,:]   == ohlc.values[1,:] .% 2
        @test (100 .- ohlc).values[1,:] == 100 .- ohlc.values[1,:]
        @test (100 .+ ohlc).values[1,:] == 100 .+ ohlc.values[1,:]
        @test (100 .* ohlc).values[1,:] == 100 .* ohlc.values[1,:]
        @test (100 ./ ohlc).values[1,:] == 100 ./ ohlc.values[1,:]
        @test (2 .^ ohlc).values[1,:]   == 2 .^ ohlc.values[1,:]
        @test (2 .% ohlc).values[1,:]   == 2 .% ohlc.values[1,:]
    end

    @testset "mathematical operations between two same-column-count TimeArrays" begin
        @test isapprox((cl .+ op).values[1], 216.82, atol=.01)
        @test isapprox((cl .- op).values[1], 7.06, atol=.01)
        @test isapprox((cl .* op).values[1], 11740.2672, atol=0.0001)
        @test isapprox((cl ./ op).values[1], 1.067315, atol=0.001)
        @test (cl .% op).values                                 == cl.values .% op.values
        @test (cl .^ op).values                                 == cl.values .^ op.values
        @test (cl .* (cl.> 200)).values                         == cl.values .* (cl.values .> 200)
        @test (round.(cl) .* cl).values                         == round.(Int, cl.values) .* cl.values
        @test (ohlc .+ ohlc).values                             == ohlc.values .+ ohlc.values
        @test (ohlc .- ohlc).values                             == ohlc.values .- ohlc.values
        @test (ohlc .* ohlc).values                             == ohlc.values .* ohlc.values
        @test (ohlc ./ ohlc).values                             == ohlc.values ./ ohlc.values
        @test (ohlc .% ohlc).values                             == ohlc.values .% ohlc.values
        @test (ohlc .^ ohlc).values                             == ohlc.values .^ ohlc.values
        @test (ohlc .* (ohlc .> 200)).values                    == ohlc.values .* (ohlc.values .> 200)
        @test (round.(ohlc) .* ohlc).values                     == round.(Int, ohlc.values) .* ohlc.values
    end

    @testset "broadcasted mathematical operations between different-column-count TimeArrays" begin
        @test (ohlc .+ cl).values                             == ohlc.values .+ cl.values
        @test (ohlc .- cl).values                             == ohlc.values .- cl.values
        @test (ohlc .* cl).values                             == ohlc.values .* cl.values
        @test (ohlc ./ cl).values                             == ohlc.values ./ cl.values
        @test (ohlc .% cl).values                             == ohlc.values .% cl.values
        @test (ohlc .^ cl).values                             == ohlc.values .^ cl.values
        @test (ohlc .* (cl.> 200)).values                     == ohlc.values .* (cl.values .> 200)
        @test (round.(ohlc) .* cl).values                     == round.(Int, ohlc.values) .* cl.values
        @test (cl .+ ohlc).values                             == cl.values .+ ohlc.values
        @test (cl .- ohlc).values                             == cl.values .- ohlc.values
        @test (cl .* ohlc).values                             == cl.values .* ohlc.values
        @test (cl ./ ohlc).values                             == cl.values ./ ohlc.values
        @test (cl .% ohlc).values                             == cl.values .% ohlc.values
        @test (cl .^ ohlc).values                             == cl.values .^ ohlc.values
        @test (cl .* (ohlc .> 200)).values                    == cl.values .* (ohlc.values .> 200)
        @test (round.(cl) .* ohlc).values                     == round.(Int, cl.values) .* ohlc.values
        # One array must have a single column
        @test_throws ErrorException (ohlc["Open", "Close"] .+ ohlc)
    end

    @testset "comparison operations between TimeArray values and Int/Float (and viceversa)" begin
        @test (cl .> 111.94).values[1]      == false
        @test (cl .< 111.94).values[1]      == false
        @test (cl .>= 111.94).values[1]     == true
        @test (cl .<= 111.94).values[1]     == true
        @test (cl .== 111.94).values[1]     == true
        @test (cl .!= 111.94).values[1]     == false
        @test (111.94 .> cl).values[1]      == false
        @test (111.94 .< cl).values[1]      == false
        @test (111.94 .>= cl).values[1]     == true
        @test (111.94 .<= cl).values[1]     == true
        @test (111.94 .== cl).values[1]     == true
        @test (111.94 .!= cl).values[1]     == false
        @test (ohlc .> 111.94).values[1,:]  == [false, true, false, false]
        @test (ohlc .< 111.94).values[1,:]  == [true, false, true, false]
        @test (ohlc .>= 111.94).values[1,:] == [false, true, false, true]
        @test (ohlc .<= 111.94).values[1,:] == [true, false, true, true]
        @test (ohlc .== 111.94).values[1,:] == [false, false, false, true]
        @test (ohlc .!= 111.94).values[1,:] == [true, true, true, false]
        @test (111.94 .> ohlc).values[1,:]  == [true, false, true, false]
        @test (111.94 .< ohlc).values[1,:]  == [false, true, false, false]
        @test (111.94 .>= ohlc).values[1,:] == [true, false, true, true]
        @test (111.94 .<= ohlc).values[1,:] == [false, true, false, true]
        @test (111.94 .== ohlc).values[1,:] == [false, false, false, true]
        @test (111.94 .!= ohlc).values[1,:] == [true, true, true, false]
    end

    @testset "comparison operations between TimeArray values and Bool (and viceversa)" begin
        @test ((cl .> 111.94) .== true).values[1]     == false
        @test ((cl .> 111.94) .!= true).values[1]     == true
        @test (true .== (cl .> 111.94)).values[1]     == false
        @test (true .!= (cl .> 111.94)).values[1]     == true
        @test ((ohlc .> 111.94).== true).values[1,:]  == [false, true, false, false]
        @test ((ohlc .> 111.94).!= true).values[1,:]  == [true, false, true, true]
        @test (true .== (ohlc .> 111.94)).values[1,:] == [false, true, false, false]
        @test (true .!= (ohlc .> 111.94)).values[1,:] == [true, false, true, true]
    end

    @testset "comparison operations between same-column-count TimeArrays" begin
        @test (cl .> op).values[1]        == true
        @test (cl .< op).values[1]        == false
        @test (cl .<= op).values[1]       == false
        @test (cl .>= op).values[1]       == true
        @test (cl .== op).values[1]       == false
        @test (cl .!= op).values[1]       == true
        @test (ohlc .> ohlc).values[1,:]  == [false, false, false, false]
        @test (ohlc .< ohlc).values[1,:]  == [false, false, false, false]
        @test (ohlc .<= ohlc).values[1,:] == [true, true, true, true]
        @test (ohlc .>= ohlc).values[1,:] == [true, true, true, true]
        @test (ohlc .== ohlc).values[1,:] == [true, true, true, true]
        @test (ohlc .!= ohlc).values[1,:] == [false, false, false, false]
    end

    @testset "comparison operations between different-column-count TimeArrays" begin
        @test (ohlc .> cl).values  == (ohlc.values .> cl.values)
        @test (ohlc .< cl).values  == (ohlc.values .< cl.values)
        @test (ohlc .>= cl).values == (ohlc.values .>= cl.values)
        @test (ohlc .<= cl).values == (ohlc.values .<= cl.values)
        @test (ohlc .== cl).values == (ohlc.values .== cl.values)
        @test (ohlc .!= cl).values == (ohlc.values .!= cl.values)
        @test (cl .> ohlc).values  == (cl.values .> ohlc.values)
        @test (cl .< ohlc).values  == (cl.values .< ohlc.values)
        @test (cl .>= ohlc).values == (cl.values .>= ohlc.values)
        @test (cl .<= ohlc).values == (cl.values .<= ohlc.values)
        @test (cl .== ohlc).values == (cl.values .== ohlc.values)
        @test (cl .!= ohlc).values == (cl.values .!= ohlc.values)
        # One array must have a single column
        @test_throws ErrorException (ohlc["Open", "Close"] .== ohlc)
    end

    @testset "bitwise elementwise operations between bool and TimeArrays' values" begin
        @test ((cl .> 100) .& true).values[1]      == true
        @test ((cl .> 100) .| true).values[1]      == true
        @test ((cl .> 100) .⊻ true).values[1]      == false
        @test (false .& (cl .> 100)).values[1]     == false
        @test (false .| (cl .> 100)).values[1]     == true
        @test (false .⊻ (cl .> 100)).values[1]     == true
        @test ((ohlc .> 100) .& true).values[4,:]  == [true, true, false, false]
        @test ((ohlc .> 100) .| true).values[4,:]  == [true, true, true, true]
        @test ((ohlc .> 100) .⊻ true).values[4,:]  == [false, false, true, true]
        @test (false .& (ohlc .> 100)).values[4,:] == [false, false, false, false]
        @test (false .| (ohlc .> 100)).values[4,:] == [true, true, false, false]
        @test (false .⊻ (ohlc .> 100)).values[4,:] == [true, true, false, false]
      end

    @testset "bitwise elementwise operations between same-column-count TimeArrays' boolean values" begin
        @test ((cl .> 100) .& (cl .< 120)).values[1]       == true
        @test ((cl .> 100) .| (cl .< 120)).values[1]       == true
        @test ((cl .> 100) .⊻ (cl .< 120)).values[1]       == false
        @test ((ohlc .> 100) .& (ohlc .< 120)).values[4,:] == [true, true, false, false]
        @test ((ohlc .> 100) .| (ohlc .< 120)).values[4,:] == [true, true, true, true]
        @test ((ohlc .> 100) .⊻ (ohlc .< 120)).values[4,:] == [false, false, true, true]
        @test ((ohlc .> 100) .⊻ (cl .< 120)).values[4,:]   == [false, false, true, true]
    end
end

#TODO: add test cases of non-standard function on dot-call

@testset "basecall works with Base methods" begin
    @testset "cumsum works" begin
        @test basecall(cl, cumsum).values[2] == cl.values[1] + cl.values[2]
    end
end

@testset "adding/removing missing rows works" begin
    uohlc = uniformspace(ohlc[1:8])

    @testset "uniform spacing detection works" begin
        @test uniformspaced(datetime1)
        @test uniformspaced(uohlc)
        @test !uniformspaced(cl)
        @test !uniformspaced(ohlc)
    end

    @testset "forcing uniform spacing works" begin
        @test length(uohlc)                               == 10
        @test uohlc[5].values                             == ohlc[5].values
        @test isequal(uohlc[6].values, [NaN NaN NaN NaN]) == true
        @test isequal(uohlc[7].values, [NaN NaN NaN NaN]) == true
        @test uohlc[8].values                             == ohlc[6].values
    end

    @testset "dropnan works" begin
        nohlc = TimeArray(ohlc.timestamp, copy(ohlc.values), ohlc.colnames, ohlc.meta)
        nohlc.values[7:12, 2] = NaN

        @test dropnan(uohlc).timestamp       == dropnan(uohlc, :all).timestamp
        @test dropnan(uohlc).values          == dropnan(uohlc, :all).values

        @test dropnan(ohlc, :all).values     == ohlc.values
        @test dropnan(nohlc, :all).timestamp == ohlc.timestamp
        @test dropnan(uohlc, :all).timestamp == ohlc[1:8].timestamp
        @test dropnan(uohlc, :all).values    == ohlc[1:8].values

        @test dropnan(ohlc, :any).values     == ohlc.values
        @test dropnan(nohlc, :any).timestamp == ohlc.timestamp[[1:6;13:end]]
        @test dropnan(nohlc, :any).values    == ohlc.values[[1:6;13:end], :]
        @test dropnan(uohlc, :any).timestamp == ohlc[1:8].timestamp
        @test dropnan(uohlc, :any).values    == ohlc[1:8].values
    end
end

end  # @testset "apply"


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

using Markdown
using InteractiveUtils

# ╔═╡ 8336bf6e-0fa8-11eb-1f48-37cfe47a6d0f
using Pkg

# ╔═╡ 8ae04a5a-0fa8-11eb-2821-8dc46ee73269
Pkg.activate("../")

# ╔═╡ 66923a50-0fa8-11eb-1547-cf8e289818b0
using Revise

# ╔═╡ 3f9bdb02-7028-449b-93fc-01c4d3632db4
using PolaronMobility

# ╔═╡ 2024f235-aca8-427b-8579-2c869a0ba901
using QuadGK

# ╔═╡ f0529460-0fa8-11eb-1bec-1feb30645d42
using Plots

# ╔═╡ ffff1d5b-e9db-4da4-b01c-8c8075023997
using Gnuplot

# ╔═╡ d263806b-0851-4854-ba06-df1071b3d38b
const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s

# ╔═╡ c6606f73-3c31-4df8-8585-81698966008a
const Boltzmann = const kB =  1.3806504e-23;                  # kg m2 / K s2

# ╔═╡ 8b357378-449b-45cb-aa8e-fa2a5bf5feed
T=10

# ╔═╡ c3fd3779-7c5c-4370-8378-66c4b0d963af
MAPIe=polaronmobility(T, 4.5, 24.1, 2.25E12, 0.12)

# ╔═╡ 386fd0a4-0faa-11eb-14e3-a3cdf3f3cb5d
#α=5
α=MAPIe.α[1]

# ╔═╡ a86c04ae-0faa-11eb-3565-7bd3108906b2
#v,w=feynmanvw(α) # Variational solution, athermal (original) action

# ╔═╡ 4306f2f4-0faa-11eb-1d32-dd8ade1d930a
f=2.25E12

# ╔═╡ b9a25625-38bd-43d7-a4b7-5a0da5a28772
ω=2π*f

# ╔═╡ 947530ce-e248-47af-ae49-34d698312d12
βred=ħ*ω/(kB*T)

# ╔═╡ b3e10abe-0faa-11eb-308e-91c4c0886b8b
v,w=feynmanvw(α, βred) # Variational solution, thermal (Osaka) action

# ╔═╡ 4b26ec8c-0faa-11eb-2af1-e3f5e39fb358
meff=0.12

# ╔═╡ 34a7c84e-0faa-11eb-08b6-83f4ab86b689
nurange=0.0:0.1:3.5 #I know this says nu, but due to the rescaling below, this is actually THz

# ╔═╡ ad10a50c-0fa8-11eb-30a9-1b1bb4fc26b7
oldX=PolaronMobility.ImX(nurange,v,w,βred,α,ω,meff*PolaronMobility.MassElectron)

# ╔═╡ 9ca217bf-985f-48fc-8a35-913c18fd79fc
plot(oldX.nu, oldX.ImX, label="ImX")

# ╔═╡ e45514ef-e24f-46ed-9844-b525e1042a80
βred1=ħ*2π*1E12/(kB*T)

# ╔═╡ 6ffe0fba-4146-43d4-90bd-dc325da74482
feynmanvw(α,βred1)

# ╔═╡ 2f09ce0a-8243-4e3f-9a3a-48ddf0a5b12a
βred2=ħ*2π*2.25E12/(kB*T)

# ╔═╡ 172d0553-d260-491a-9b94-79b865e3ee9a
feynmanvw(α,βred2)

# ╔═╡ 28040f04-9d08-449d-bc6c-079fae287657
X=[PolaronMobility.χ(nu, βred, α, v,w) for nu in nurange]

# ╔═╡ d721c0a3-234d-4824-9619-8a8159be9220
begin
	plot(nurange.*f/1E12, imag.(X), label="ImX")
	plot!(nurange.*f/1E12, real.(X), label="ReX")
end

# ╔═╡ 73ac4039-4675-4dce-8eab-b6635481406e
begin
	# re-scaling to frequency directly within the selection of nurange
	#  i.e. now nurange is in THz
	X1=[PolaronMobility.χ(nu, βred1, α,feynmanvw(α,βred1)...)/nu for nu in nurange./(1.0)]
	X2=[PolaronMobility.χ(nu, βred2, α,feynmanvw(α,βred2)...)/nu for nu in nurange./(2.25)]
end

# ╔═╡ 399bb77e-13a9-41dc-b432-a9bc605af25a
begin
	plot(nurange, imag.(X1), label="ImX 1")
	plot!(nurange, real.(X1), label="ReX 1")
	plot!(nurange, imag.(X2), label="ImX 2")
	plot!(nurange, real.(X2), label="ReX 2")
end

# ╔═╡ d33d2c9d-ac56-47ac-a0ad-b769c2fcf245
begin
	plot(nurange, imag.(X1+1.5*X2), label="ImX")
	plot!(nurange, real.(X1+1.5*X2), label="ReX")
	xlabel!("Frequeny (THz)")
	ylabel!("Χ(ν)/ν")
end

# ╔═╡ 3474df75-bd91-4803-acf3-bd6f697a6163
begin
	@gp nurange imag.(0.4*X1+0.6*X2) "w l t 'Im'"
	@gp :- nurange real.(0.4*X1+0.6*X2) "w l t 'Re'"
	@gp :- "set xlabel 'Frequency (THz)' "
	@gp :- "set ylabel 'Χ(ν)/ν' "
end

# ╔═╡ 04a89c80-7ecc-49fe-90f1-e1e081824d87
Gnuplot.save(term="pngcairo size 1024,768", output="~/tmp_plot.png")

# ╔═╡ 7e796fb8-1149-44fa-9076-9d7c2641be03
#Notable Mode f= 2.4380741812443247e12    α_partial=0.33509445526579384
#Notable Mode f= 2.249091763771941e12    α_partial=0.46413279655805495
#Notable Mode f= 2.0796321906344238e12    α_partial=0.5034016572517709
#Notable Mode f= 2.0336707697261187e12    α_partial=0.14188851068346853
#Notable Mode f= 1.5673011873879714e12    α_partial=0.16044621957165453
#Notable Mode f= 1.0188379384951798e12    α_partial=0.16189515169321733
#Notable Mode f= 9.970130778462072e11    α_partial=0.14036635422200763
#Notable Mode f= 9.20178190638621e11    α_partial=0.18115470952505333
#Notable Mode f= 5.738689505255511e11    α_partial=0.3479224374216951
#Sum alpha: 2.658367465100708
0.33509445526579384+0.46413279655805495+0.5034016572517709+0.14188851068346853

# ╔═╡ cafcff56-6e56-4b92-8ef0-87c5df309c7a
0.16044621957165453+0.16189515169321733+0.14036635422200763+0.18115470952505333+0.3479224374216951

# ╔═╡ 8530958f-7390-470a-9b61-ce88046bc084
1.5/2.5

# ╔═╡ 41398cba-bb88-4d4a-a817-1de8f551cc9c
function σ(χ)
	ϵr=χ+1
	n=√ϵr # assuming relative permeability of MAPI is 1...
	σ=im*(1-n)
	σ
end

# ╔═╡ c21ad7eb-9d9b-4202-b4ad-b1206dc94b3c
σ(1+0im)

# ╔═╡ 775c2be5-5aa4-41e2-9956-8e4131bff73b
σ(0+1im)

# ╔═╡ dab2fc71-d8d8-4e6d-9c43-d85688ca7558
σ1=σ.([PolaronMobility.χ(nu, βred1, α,feynmanvw(α,βred1)...)/nu for nu in nurange./(1.0)])

# ╔═╡ 9b8b2293-2395-4285-86fe-1b7fb7895297
σ2=σ.([PolaronMobility.χ(nu, βred2, α,feynmanvw(α,βred2)...)/nu for nu in nurange./(2.25)])

# ╔═╡ 53303534-308e-448f-9051-2bc3d859dc14
Gnuplot.options.mime

# ╔═╡ cdba35db-b258-433c-9387-42b505e2efb4
begin
	@gp nurange imag.(0.4*σ1+0.6*σ2) "w lp t 'Im'"
	@gp :- nurange real.(0.4*σ1+0.6*σ2) "w lp t 'Re'"
	@gp :- "set xlabel 'Frequency (THz)' "
	@gp :- "set ylabel 'σ' "

	Gnuplot.save("~/complex_conductivity_MAPI.gp")
	Gnuplot.save(term="pngcairo size 1024,768 fontscale 2 pointscale 2 linewidth 2", output="~/complex_conductivity_MAPI.png")

end

# ╔═╡ e0c76116-8c8a-4053-86c7-2f7c48f808bf
begin
	@gp nurange imag.(σ1) "w lp t '1THz phonons Im'"
	@gp :- nurange real.(σ1) "w lp t '1THz phonons Re'"
	@gp :- nurange imag.(σ2) "w lp t '2.25THz phonons Im'"
	@gp :- nurange real.(σ2) "w lp t '2.25THz phonons Re'"
	@gp :- nurange imag.(0.4*σ1+0.6*σ2) "w lp t 'Combined Im'"
	@gp :- nurange real.(0.4*σ1+0.6*σ2) "w lp t 'Combined Re'"

	@gp :- "set xlabel 'Frequency (THz)' "
	@gp :- "set ylabel 'σ' "
end

# ╔═╡ 1880d08f-f808-4690-8ad1-af1616cda0bf
# Multiple branch free energy optimisation

# ╔═╡ e39b7ce9-5396-4065-ba6a-eac83e626794
MAPI= [
# 96.20813558773261 0.4996300522819191
# 93.13630357703363 1.7139631746083817
# 92.87834578121567 0.60108592692181
# 92.4847918585963 0.0058228799414729
# 92.26701437594754 0.100590086574602
# 89.43972834606603 0.006278895133832249
# 46.89209141511332 0.2460894564364346
# 46.420949316788 0.14174282581124137
# 44.0380222871706 0.1987196948553428
# 42.89702947649343 0.011159939465770681
# 42.67180170168193 0.02557751102757614
# 41.46971205834201 0.012555230726601503
# 37.08982543385215 0.00107488277468418
# 36.53555265689563 0.02126940080871224
# 30.20608114002676 0.009019481779712388
# 27.374810898415028 0.03994453721421388
# 26.363055017011728 0.05011922682554448
# 9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
# 0.022939578929507105 8.355742795827834e-06   # Acoustic modes!
# 0.04882611767873102 8.309858592685e-06
# 0.07575149723846182 2.78248540373041e-05
]

# ╔═╡ 3f2969ce-53c4-44ce-b442-c6c1e4c8e63c
v_params, w_params = PolaronMobility.multi_variation(0.01, 4.5, 0.12, (6.29e-10)^3, MAPI; N = 1)

# ╔═╡ 9f8e375a-e45d-430e-89f0-e62573f2f45b
F_j = [PolaronMobility.multi_free_energy(v_params[j, :], w_params[j, :], 10, 4.5, 0.12, (6.29e-10)^3, MAPI, j) for j in 1:length(MAPI[:, 1])]

# ╔═╡ b0fc3a52-b4da-4ac4-a36e-e6229cd1f20e
F_total = sum(F_j)

# ╔═╡ d36a215e-9411-4332-82cf-caa87a56c7f4
# multiple branch susceptibility and conductivity

# ╔═╡ c2e007f4-dac0-4be5-ab55-10186fb34d61
phonon_mode_freqs = MAPI[:, 1]

# ╔═╡ ecd96809-9808-46fa-9405-884aeb8c7e6e
ir_activities = MAPI[:, 2]

# ╔═╡ f0d30402-ee5f-4cf9-b4bd-af35a102f71d
ϵ_ionic = [PolaronMobility.ϵ_ionic_mode(f, r, (6.29e-10)^3) for (f, r) in zip(phonon_mode_freqs, ir_activities)]

# ╔═╡ e2faa878-0ad1-42f7-99e8-53699c1c3468
ϵ_total = sum(ϵ_ionic)

# ╔═╡ 2d161b89-5594-4601-bb8e-c1748ade6269
α_j = [PolaronMobility.frohlich_α_j(4.5, ϵ_i, ϵ_total, f, 0.12) for (ϵ_i, f) in zip(ϵ_ionic, phonon_mode_freqs)]

# ╔═╡ 928e7f3b-6eab-4081-94bf-a4fa3a953f35
α_eff = sum(α_j)

# ╔═╡ 1dd3fcda-6818-4601-9fa4-7b75a07ee8d8
βr = [ħ * 2π * 1e12 * f / kB / 40 for f in phonon_mode_freqs]

# ╔═╡ 9b410a25-1c09-4bcb-b06d-63c928ab5bbc
Ω_range = 0.01:0.01:2.5

# ╔═╡ d20fd376-b5a6-48ce-9380-fccf420bdddd
Plots.PyPlotBackend()

# ╔═╡ c69faa18-759e-465c-b720-d8c6d429780d
χ = [PolaronMobility.multi_impedence(Ω, βr, α_j, v_params, w_params, phonon_mode_freqs, 0.12) for Ω in Ω_range]

# ╔═╡ e052454b-f885-4d76-8eb1-baedeff24e1e
begin
	χ_plot = plot(Ω_range, (imag.(χ)), xlabel = "Frequency (THz)", ylabel = "z", label = "Imz", minorgrid = true, markershape = :circle)
	plot!(χ_plot, Ω_range, (real.(χ)), label = "Rez", markershape = :circle)
end

# ╔═╡ dbe241ba-5b5c-45d0-9d63-3ccefc82766f
σ_new = PolaronMobility.multi_conductivity(Ω_range, χ)

# ╔═╡ ec793073-968b-447b-a445-9182f1cfdf0f
begin
	σ_plot = plot(Ω_range, (real.(σ_new)), xlabel = "Frequency (THz)", ylabel = "σ", label = "Reσ", legend = :bottomright, minorgrid = true, markershape = :circle)
	plot!(σ_plot, Ω_range, (imag.(σ_new)), label = "Imσ", markershape = :circle)
end

# ╔═╡ 36f09479-cf29-4a0f-9f97-81169abbcd9b
σ_abs_plot = plot(Ω_range, (sqrt.(abs.(σ_new))), xlabel = "Frequency (THz)", ylabel = "abs(σ)", legend = false, minorgrid = true, markershape = :circle)

# ╔═╡ Cell order:
# ╠═d263806b-0851-4854-ba06-df1071b3d38b
# ╠═c6606f73-3c31-4df8-8585-81698966008a
# ╠═66923a50-0fa8-11eb-1547-cf8e289818b0
# ╠═8336bf6e-0fa8-11eb-1f48-37cfe47a6d0f
# ╠═8ae04a5a-0fa8-11eb-2821-8dc46ee73269
# ╠═3f9bdb02-7028-449b-93fc-01c4d3632db4
# ╠═2024f235-aca8-427b-8579-2c869a0ba901
# ╠═8b357378-449b-45cb-aa8e-fa2a5bf5feed
# ╠═c3fd3779-7c5c-4370-8378-66c4b0d963af
# ╠═386fd0a4-0faa-11eb-14e3-a3cdf3f3cb5d
# ╠═a86c04ae-0faa-11eb-3565-7bd3108906b2
# ╠═947530ce-e248-47af-ae49-34d698312d12
# ╠═b3e10abe-0faa-11eb-308e-91c4c0886b8b
# ╠═4306f2f4-0faa-11eb-1d32-dd8ade1d930a
# ╠═b9a25625-38bd-43d7-a4b7-5a0da5a28772
# ╠═4b26ec8c-0faa-11eb-2af1-e3f5e39fb358
# ╠═34a7c84e-0faa-11eb-08b6-83f4ab86b689
# ╠═ad10a50c-0fa8-11eb-30a9-1b1bb4fc26b7
# ╠═f0529460-0fa8-11eb-1bec-1feb30645d42
# ╠═9ca217bf-985f-48fc-8a35-913c18fd79fc
# ╠═e45514ef-e24f-46ed-9844-b525e1042a80
# ╠═6ffe0fba-4146-43d4-90bd-dc325da74482
# ╠═2f09ce0a-8243-4e3f-9a3a-48ddf0a5b12a
# ╠═172d0553-d260-491a-9b94-79b865e3ee9a
# ╠═28040f04-9d08-449d-bc6c-079fae287657
# ╠═d721c0a3-234d-4824-9619-8a8159be9220
# ╠═73ac4039-4675-4dce-8eab-b6635481406e
# ╠═399bb77e-13a9-41dc-b432-a9bc605af25a
# ╠═d33d2c9d-ac56-47ac-a0ad-b769c2fcf245
# ╠═ffff1d5b-e9db-4da4-b01c-8c8075023997
# ╠═3474df75-bd91-4803-acf3-bd6f697a6163
# ╠═04a89c80-7ecc-49fe-90f1-e1e081824d87
# ╠═7e796fb8-1149-44fa-9076-9d7c2641be03
# ╠═cafcff56-6e56-4b92-8ef0-87c5df309c7a
# ╠═8530958f-7390-470a-9b61-ce88046bc084
# ╠═41398cba-bb88-4d4a-a817-1de8f551cc9c
# ╠═c21ad7eb-9d9b-4202-b4ad-b1206dc94b3c
# ╠═775c2be5-5aa4-41e2-9956-8e4131bff73b
# ╠═dab2fc71-d8d8-4e6d-9c43-d85688ca7558
# ╠═9b8b2293-2395-4285-86fe-1b7fb7895297
# ╠═53303534-308e-448f-9051-2bc3d859dc14
# ╠═cdba35db-b258-433c-9387-42b505e2efb4
# ╠═e0c76116-8c8a-4053-86c7-2f7c48f808bf
# ╠═1880d08f-f808-4690-8ad1-af1616cda0bf
# ╠═e39b7ce9-5396-4065-ba6a-eac83e626794
# ╠═3f2969ce-53c4-44ce-b442-c6c1e4c8e63c
# ╠═9f8e375a-e45d-430e-89f0-e62573f2f45b
# ╠═b0fc3a52-b4da-4ac4-a36e-e6229cd1f20e
# ╠═d36a215e-9411-4332-82cf-caa87a56c7f4
# ╠═c2e007f4-dac0-4be5-ab55-10186fb34d61
# ╠═ecd96809-9808-46fa-9405-884aeb8c7e6e
# ╠═f0d30402-ee5f-4cf9-b4bd-af35a102f71d
# ╠═e2faa878-0ad1-42f7-99e8-53699c1c3468
# ╠═2d161b89-5594-4601-bb8e-c1748ade6269
# ╠═928e7f3b-6eab-4081-94bf-a4fa3a953f35
# ╠═1dd3fcda-6818-4601-9fa4-7b75a07ee8d8
# ╠═9b410a25-1c09-4bcb-b06d-63c928ab5bbc
# ╠═d20fd376-b5a6-48ce-9380-fccf420bdddd
# ╠═c69faa18-759e-465c-b720-d8c6d429780d
# ╠═e052454b-f885-4d76-8eb1-baedeff24e1e
# ╠═dbe241ba-5b5c-45d0-9d63-3ccefc82766f
# ╠═ec793073-968b-447b-a445-9182f1cfdf0f
# ╠═36f09479-cf29-4a0f-9f97-81169abbcd9b


text:
"""
    _pretty_substances(ss::Vector{String})::String

Nicely format a substance list.
"""
function _pretty_substances(ss::Vector{String})::String
    if isempty(ss)
        "∅"
    elseif length(ss) > 5
        join([ss[1], ss[2], "...", ss[end-1], ss[end]], " + ")
    else
        join(ss, " + ")
    end
end

function Base.show(io::IO, ::MIME"text/plain", r::Reaction)
    if r.ub > 0.0 && r.lb < 0.0
        arrow = " ⟷  "
    elseif r.ub <= 0.0 && r.lb < 0.0
        arrow = " ⟵  "
    elseif r.ub > 0.0 && r.lb >= 0.0
        arrow = " ⟶  "
    else
        arrow = " →∣←  " # blocked reaction
    end
    substrates =
        ["$(-v) $k" for (k, v) in Iterators.filter(((_, v)::Pair -> v < 0), r.metabolites)]
    products =
        ["$v $k" for (k, v) in Iterators.filter(((_, v)::Pair -> v >= 0), r.metabolites)]

    for fname in fieldnames(Reaction)
        if fname == :metabolites
            _pretty_print_keyvals(
                io,
                "Reaction.$(string(fname)): ",
                _pretty_substances(substrates) * arrow * _pretty_substances(products),
            )
        elseif fname == :grr
            _pretty_print_keyvals(
                io,
                "Reaction.$(string(fname)): ",
                _maybemap(x -> _unparse_grr(String, x), r.grr),
            )
        elseif fname in (:lb, :ub, :objective_coefficient)
            _pretty_print_keyvals(
                io,
                "Reaction.$(string(fname)): ",
                string(getfield(r, fname)),
            )
        else
            _pretty_print_keyvals(io, "Reaction.$(string(fname)): ", getfield(r, fname))
        end
    end
end


text:
# Imports
using OpenMDAO: AbstractExplicitComp, VarData, PartialsData, get_rows_cols
import OpenMDAO: setup, compute!, compute_partials!
using Interpolations: LinearInterpolation
using PCHIPInterpolation
using Statistics: mean
using ForwardDiff: jacobian!, JacobianConfig
using ComponentArrays
using ConcreteStructs

# Define source struct
@concrete struct Source{} <: AbstractExplicitComp
    settings
    data
    ac
    shield
    n_t
    idx_src
    source!
    X
    Y
    J
end

function Source(settings, data, ac, shield, n_t::Int, idx_src)

    include(settings.pyNA_directory * "/src/noise_src_jl/fan.jl")

    include(settings.pyNA_directory * "/src/noise_src_jl/core.jl")

    include(settings.pyNA_directory * "/src/noise_src_jl/jet.jl")

    include(settings.pyNA_directory * "/src/noise_src_jl/airframe.jl")

    function source!(settings, data, ac, shield, n_t, idx_src, y, x)

        # Create input_src
        input_src = x[:TS]
        input_src = vcat(input_src, x[:M_0])
        input_src = vcat(input_src, x[:c_0])
        input_src = vcat(input_src, x[:rho_0])
        input_src = vcat(input_src, x[:mu_0])
        input_src = vcat(input_src, x[:T_0])
        input_src = vcat(input_src, x[:theta])
        input_src = vcat(input_src, x[:phi])
        if settings.jet_mixing == true && settings.jet_shock == false
            input_src = vcat(input_src, x[:V_j_star])
            input_src = vcat(input_src, x[:rho_j_star])
            input_src = vcat(input_src, x[:A_j_star])
            input_src = vcat(input_src, x[:Tt_j_star])
        elseif settings.jet_shock == true && settings.jet_mixing == false
            input_src = vcat(input_src, x[:V_j_star])
            input_src = vcat(input_src, x[:M_j])
            input_src = vcat(input_src, x[:A_j_star])
            input_src = vcat(input_src, x[:Tt_j_star])
        elseif settings.jet_shock ==true && settings.jet_mixing == true
            input_src = vcat(input_src, x[:V_j_star])
            input_src = vcat(input_src, x[:rho_j_star])
            input_src = vcat(input_src, x[:A_j_star])
            input_src = vcat(input_src, x[:Tt_j_star])
            input_src = vcat(input_src, x[:M_j])
        end
        if settings.core
            if settings.method_core_turb == "GE"
                input_src = vcat(input_src, x[:mdoti_c_star])
                input_src = vcat(input_src, x[:Tti_c_star])
                input_src = vcat(input_src, x[:Ttj_c_star])
                input_src = vcat(input_src, x[:Pti_c_star])
                input_src = vcat(input_src, x[:DTt_des_c_star])
            elseif settings.method_core_turb == "PW"
                input_src = vcat(input_src, x[:mdoti_c_star])
                input_src = vcat(input_src, x[:Tti_c_star])
                input_src = vcat(input_src, x[:Ttj_c_star])
                input_src = vcat(input_src, x[:Pti_c_star])
                input_src = vcat(input_src, x[:rho_te_c_star])
                input_src = vcat(input_src, x[:c_te_c_star])
                input_src = vcat(input_src, x[:rho_ti_c_star])
                input_src = vcat(input_src, x[:c_ti_c_star])
            end
        end
        if settings.airframe
            input_src = vcat(input_src, x[:theta_flaps])
            input_src = vcat(input_src, x[:I_landing_gear])
        end
        if settings.fan_inlet==true || settings.fan_discharge==true
            input_src = vcat(input_src, x[:DTt_f_star])
            input_src = vcat(input_src, x[:mdot_f_star])
            input_src = vcat(input_src, x[:N_f_star])
            input_src = vcat(input_src, x[:A_f_star])
            input_src = vcat(input_src, x[:d_f_star])
        end

        # Extract TS
        TS = x[:TS]

        # Get type of input vector
        T = eltype(input_src)

        # Initialize source mean-square acoustic pressure
        msap_source = zeros(T, (n_t, settings.N_f))

        if settings.fan_inlet
            msap_fan_inlet = fan(settings, data, ac, n_t, shield, idx_src, input_src, "fan_inlet")
            msap_source = msap_source .+ msap_fan_inlet
        end
        if settings.fan_discharge
            msap_fan_discharge = fan(settings, data, ac, n_t, shield, idx_src, input_src, "fan_discharge")
            msap_source = msap_source .+ msap_fan_discharge
        end

        if settings.core
            msap_core = core(settings, data, ac, n_t, idx_src, input_src)
            if settings.suppression && settings.case_name in ["NASA STCA Standard", "stca_enginedesign_standard"]
                msap_core[findall(TS.*ones(1, settings.N_f).>0.8)] = (10. ^(-2.3 / 10.) * msap_core)[findall(TS.*ones(1, settings.N_f).>0.8)]
            end
            msap_source = msap_source .+ msap_core
        end

        if settings.jet_mixing
            msap_jet_mixing = jet_mixing(settings, data, ac, n_t, idx_src, input_src)
            if settings.suppression && settings.case_name in ["NASA STCA Standard", "stca_enginedesign_standard"]
                msap_jet_mixing[findall(TS.*ones(1, settings.N_f).>0.8)] = (10. ^(-2.3 / 10.) * msap_jet_mixing)[findall(TS.*ones(1, settings.N_f).>0.8)]
            end
            msap_source = msap_source .+ msap_jet_mixing
        end
        if settings.jet_shock
            msap_jet_shock = jet_shock(settings, data, ac, n_t, idx_src, input_src)
            if settings.suppression && settings.case_name in ["NASA STCA Standard", "stca_enginedesign_standard"]
                msap_jet_shock[findall(TS.*ones(1, settings.N_f).>0.8)] = (10. ^(-2.3 / 10.) * msap_jet_shock)[findall(TS.*ones(1, settings.N_f).>0.8)]
            end
            msap_source = msap_source .+ msap_jet_shock
        end

        if settings.airframe
            msap_af = airframe(settings, data, ac, n_t, idx_src, input_src)
            msap_source = msap_source .+ msap_af
        end

         y[:msap_source] = 10*log10.(msap_source)
        #y[:msap_source] = msap_source

    end

    X = ComponentArray(TS=1.,
                       M_0=1.,
                       c_0=1.,
                       rho_0=1.,
                       mu_0=1.,
                       T_0=1.,
                       theta=1.,
                       phi=1.)
    if settings.jet_mixing == true && settings.jet_shock == false
        X = vcat(X, ComponentArray(V_j_star=1.,
                                   rho_j_star=1.,
                                   A_j_star=1.,
                                   Tt_j_star=1.))
    elseif settings.jet_shock == true && settings.jet_mixing == false
        X = vcat(X, ComponentArray(V_j_star=1.,
                                   M_j=1.,
                                   A_j_star=1.,
                                   Tt_j_star=1.))
    elseif settings.jet_shock ==true && settings.jet_mixing == true
        X = vcat(X, ComponentArray(V_j_star=1.,
                                   rho_j_star=1.,
                                   A_j_star=1.,
                                   Tt_j_star=1.,
                                   M_j=1.))
    end
    if settings.core
        if settings.method_core_turb == "GE"
            X = vcat(X, ComponentArray(mdoti_c_star=1.,
                                       Tti_c_star=1.,
                                       Ttj_c_star=1.,
                                       Pti_c_star=1.,
                                       DTt_des_c_star=1.))
        elseif settings.method_core_turb == "PW"
            X = vcat(X, ComponentArray(mdoti_c_star=1.,
                                       Tti_c_star=1.,
                                       Ttj_c_star=1.,
                                       Pti_c_star=1.,
                                       rho_te_c_star=1.,
                                       c_te_c_star=1.,
                                       rho_ti_c_star=1.,
                                       c_ti_c_star=1.))
        end
    end
    if settings.airframe
        X = vcat(X, ComponentArray(theta_flaps=1.,
                                   I_landing_gear=1.))
    end
    if settings.fan_inlet==true || settings.fan_discharge==true
        X = vcat(X, ComponentArray(DTt_f_star=1.,
                                   mdot_f_star=1.,
                                   N_f_star=1.,
                                   A_f_star=1.,
                                   d_f_star=1.))
    end

    Y = ComponentArray(msap_source=zeros(Float64, settings.N_f))
    J = Y.*X'
    #'

    return Source(settings, data, ac, shield, n_t, idx_src, source!, X, Y, J)
end

function setup(self::Source)
    # Load options
    n_t = self.n_t
    settings = self.settings

    # Create inputs
    inputs = Vector{VarData}()
    if settings.jet_mixing == true && settings.jet_shock == false
        push!(inputs, VarData("V_j_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("rho_j_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("A_j_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("Tt_j_star", shape=(n_t,), val=ones(n_t)))
    elseif settings.jet_shock == true && settings.jet_mixing == false
        push!(inputs, VarData("V_j_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("M_j", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("A_j_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("Tt_j_star", shape=(n_t,), val=ones(n_t)))
    elseif settings.jet_shock ==true && settings.jet_mixing == true
        push!(inputs, VarData("V_j_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("rho_j_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("A_j_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("Tt_j_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("M_j", shape=(n_t,), val=ones(n_t)))
    end
    if settings.core
        if settings.method_core_turb == "GE"
            push!(inputs, VarData("mdoti_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("Tti_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("Ttj_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("Pti_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("DTt_des_c_star", shape=(n_t,), val=ones(n_t)))
        elseif settings.method_core_turb == "PW"
            push!(inputs, VarData("mdoti_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("Tti_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("Ttj_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("Pti_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("rho_te_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("c_te_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("rho_ti_c_star", shape=(n_t,), val=ones(n_t)))
            push!(inputs, VarData("c_ti_c_star", shape=(n_t,), val=ones(n_t)))
        end
    end
    if settings.airframe
        push!(inputs, VarData("theta_flaps", shape=(n_t,), val=ones(n_t), units="deg"))
        push!(inputs, VarData("I_landing_gear", shape=(n_t,), val=ones(n_t)))
    end
    if settings.fan_inlet==true || settings.fan_discharge==true
        push!(inputs, VarData("DTt_f_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("mdot_f_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("N_f_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("A_f_star", shape=(n_t,), val=ones(n_t)))
        push!(inputs, VarData("d_f_star", shape=(n_t,), val=ones(n_t)))
    end

    # Power setting input
    push!(inputs, VarData("TS", shape=(n_t,), val=ones(n_t)))
    push!(inputs, VarData("M_0", shape=(n_t,), val=ones(n_t)))
    push!(inputs, VarData("c_0", shape=(n_t,), val=ones(n_t), units="m/s"))
    push!(inputs, VarData("rho_0", shape=(n_t,), val=ones(n_t), units="kg/m**3"))
    push!(inputs, VarData("mu_0", shape=(n_t,), val=ones(n_t), units="kg/m/s"))
    push!(inputs, VarData("T_0", shape=(n_t,), val=ones(n_t), units="K"))
    push!(inputs, VarData("theta", shape=(n_t,), val=ones(n_t), units="deg"))
    push!(inputs, VarData("phi", shape=(n_t,), val=ones(n_t), units="deg"))

    # Create outputs
    # Output
    outputs = Vector{VarData}()
    push!(outputs, VarData("msap_source", shape=((n_t, settings.N_f)), val=zeros((n_t, settings.N_f))))

    # Create partials
    partials = Vector{PartialsData}()
    ss_sizes = Dict(:i=>n_t, :j=>settings.N_f, :k=>1);
    rows_mv, cols_mv = get_rows_cols(ss_sizes=ss_sizes, of_ss=[:i, :j], wrt_ss=[:i])  # Partials of matrix (ntxNf) wrt vector (nt)
    if settings.jet_mixing && settings.jet_shock == false
        push!(partials, PartialsData("msap_source", "V_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "rho_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "A_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "Tt_j_star", rows=rows_mv, cols=cols_mv))
    elseif settings.jet_shock && settings.jet_mixing == false
        push!(partials, PartialsData("msap_source", "V_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "A_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "Tt_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "M_j", rows=rows_mv, cols=cols_mv))
    elseif settings.jet_shock && settings.jet_mixing
        push!(partials, PartialsData("msap_source", "V_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "rho_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "A_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "Tt_j_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "M_j", rows=rows_mv, cols=cols_mv))
    end

    if settings.core
        if settings.method_core_turb == "GE"
            push!(partials, PartialsData("msap_source", "mdoti_c_star", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "Tti_c_star", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "Ttj_c_star", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "Pti_c_star", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "DTt_des_c_star", rows=rows_mv, cols=cols_mv))
        elseif settings.method_core_turb == "PW"
            push!(partials, PartialsData("msap_source", "mdoti_c_star", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "Tti_c_star", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "Ttj_c_star", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "Pti_c_star", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "rho_te_c", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "c_te_c", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "rho_ti_c", rows=rows_mv, cols=cols_mv))
            push!(partials, PartialsData("msap_source", "c_ti_c", rows=rows_mv, cols=cols_mv))
        end
    end
    if settings.fan_inlet==true || settings.fan_discharge==true
        push!(partials, PartialsData("msap_source", "DTt_f_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "mdot_f_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "N_f_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "A_f_star", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "d_f_star", rows=rows_mv, cols=cols_mv))
    end
    if settings.airframe
        push!(partials, PartialsData("msap_source", "theta_flaps", rows=rows_mv, cols=cols_mv))
        push!(partials, PartialsData("msap_source", "I_landing_gear", rows=rows_mv, cols=cols_mv))
    end

    push!(partials, PartialsData("msap_source", "theta", rows=rows_mv, cols=cols_mv))
    push!(partials, PartialsData("msap_source", "phi", rows=rows_mv, cols=cols_mv))
    push!(partials, PartialsData("msap_source", "TS", rows=rows_mv, cols=cols_mv))
    push!(partials, PartialsData("msap_source", "M_0", rows=rows_mv, cols=cols_mv))
    push!(partials, PartialsData("msap_source", "c_0", rows=rows_mv, cols=cols_mv))
    push!(partials, PartialsData("msap_source", "mu_0", rows=rows_mv, cols=cols_mv))
    push!(partials, PartialsData("msap_source", "rho_0", rows=rows_mv, cols=cols_mv))
    push!(partials, PartialsData("msap_source", "T_0", rows=rows_mv, cols=cols_mv))

    return inputs, outputs, partials
end

function compute!(self::Source, inputs, outputs)
    # Loda options
    settings = self.settings
    data = self.data
    ac = self.ac
    shield = self.shield
    n_t = self.n_t
    idx_src = self.idx_src
    X = self.X
    Y = self.Y

    # Print start statement
    println("Computing noise sources")

    for i  in range(1, n_t, step=1)
        if settings.fan_inlet
            X[:DTt_f_star] = inputs["DTt_f_star"][i]
            X[:mdot_f_star] = inputs["mdot_f_star"][i]
            X[:N_f_star] = inputs["N_f_star"][i]
            X[:A_f_star] = inputs["A_f_star"][i]
            X[:d_f_star] = inputs["d_f_star"][i]
            X[:M_0] = inputs["M_0"][i]
            X[:c_0] = inputs["c_0"][i]
            X[:T_0] = inputs["T_0"][i]
            X[:rho_0] = inputs["rho_0"][i]
            X[:theta] = inputs["theta"][i]
        end
        if settings.fan_discharge
            X[:DTt_f_star] = inputs["DTt_f_star"][i]
            X[:mdot_f_star] = inputs["mdot_f_star"][i]
            X[:N_f_star] = inputs["N_f_star"][i]
            X[:A_f_star] = inputs["A_f_star"][i]
            X[:d_f_star] = inputs["d_f_star"][i]
            X[:M_0] = inputs["M_0"][i]
            X[:c_0] = inputs["c_0"][i]
            X[:T_0] = inputs["T_0"][i]
            X[:rho_0] = inputs["rho_0"][i]
            X[:theta] = inputs["theta"][i]
        end

        if settings.core
            if settings.method_core_turb == "GE"
                X[:mdoti_c_star] = inputs["mdoti_c_star"][i]
                X[:Tti_c_star] = inputs["Tti_c_star"][i]
                X[:Ttj_c_star] = inputs["Ttj_c_star"][i]
                X[:Pti_c_star] = inputs["Pti_c_star"][i]
                X[:DTt_des_c_star] = inputs["DTt_des_c_star"][i]
                X[:TS] = inputs["TS"][i]
                X[:theta] = inputs["theta"][i]
                X[:M_0] = inputs["M_0"][i]
            elseif settings.method_core_turb == "PW"
                X[:mdoti_c_star] = inputs["mdoti_c_star"][i]
                X[:Tti_c_star] = inputs["Tti_c_star"][i]
                X[:Ttj_c_star] = inputs["Ttj_c_star"][i]
                X[:Pti_c_star] = inputs["Pti_c_star"][i]
                X[:rho_te_c_star] = inputs["rho_te_c_star"][i]
                X[:c_te_c_star] = inputs["c_te_c_star"][i]
                X[:rho_ti_c_star] = inputs["rho_ti_c_star"][i]
                X[:c_ti_c_star] = inputs["c_ti_c_star"][i]
                X[:TS] = inputs["TS"][i]
                X[:theta] = inputs["theta"][i]
                X[:M_0] = inputs["M_0"][i]
            else
                throw(DomainError("Invalid method for turbine noise in core module. Specify: GE/PW."))
            end
        end

        if settings.jet_mixing
            X[:V_j_star] = inputs["V_j_star"][i]
            X[:rho_j_star] = inputs["rho_j_star"][i]
            X[:A_j_star] = inputs["A_j_star"][i]
            X[:Tt_j_star] = inputs["Tt_j_star"][i]
            X[:TS] = inputs["TS"][i]
            X[:M_0] = inputs["M_0"][i]
            X[:c_0] = inputs["c_0"][i]
            X[:theta] = inputs["theta"][i]
        end
        if settings.jet_shock
            X[:V_j_star] = inputs["V_j_star"][i]
            X[:M_j] = inputs["M_j"][i]
            X[:A_j_star] = inputs["A_j_star"][i]
            X[:Tt_j_star] = inputs["Tt_j_star"][i]
            X[:TS] = inputs["TS"][i]
            X[:M_0] = inputs["M_0"][i]
            X[:c_0] = inputs["c_0"][i]
            X[:theta] = inputs["theta"][i]
        end

        if settings.airframe
            X[:theta_flaps] = inputs["theta_flaps"][i]
            X[:I_landing_gear] = round(inputs["I_landing_gear"][i])
            X[:M_0] = inputs["M_0"][i]
            X[:c_0] = inputs["c_0"][i]
            X[:rho_0] = inputs["rho_0"][i]
            X[:mu_0] = inputs["mu_0"][i]
            X[:theta] = inputs["theta"][i]
            X[:phi] = inputs["phi"][i]
        end

        Y[:msap_source] = zeros(Float64, settings.N_f)

        self.source!(settings, data, ac, shield[i, 2:end], n_t, idx_src, Y, X)
        msap_source = Y[:msap_source]
        @. outputs["msap_source"][i,:] = msap_source

    end
end

function compute_partials!(self::Source, inputs, partials)
    # Load options
    n_t = self.n_t
    settings = self.settings
    data = self.data
    ac = self.ac
    shield = self.shield
    n_t = self.n_t
    idx_src = self.idx_src
    X = self.X
    Y = self.Y
    J = self.J

    # Print start statement
    println("Computing partials source")

    for i  in range(1, n_t, step=1)

        X[:M_0] = inputs["M_0"][i]
        X[:c_0] = inputs["c_0"][i]
        X[:T_0] = inputs["T_0"][i]
        X[:rho_0] = inputs["rho_0"][i]
        X[:mu_0] = inputs["mu_0"][i]
        X[:theta] = inputs["theta"][i]
        X[:phi] = inputs["phi"][i]
        X[:TS] = inputs["TS"][i]

        if settings.fan_inlet
            X[:DTt_f_star] = inputs["DTt_f_star"][i]
            X[:mdot_f_star] = inputs["mdot_f_star"][i]
            X[:N_f_star] = inputs["N_f_star"][i]
            X[:A_f_star] = inputs["A_f_star"][i]
            X[:d_f_star] = inputs["d_f_star"][i]
        end
        if settings.fan_discharge
            X[:DTt_f_star] = inputs["DTt_f_star"][i]
            X[:mdot_f_star] = inputs["mdot_f_star"][i]
            X[:N_f_star] = inputs["N_f_star"][i]
            X[:A_f_star] = inputs["A_f_star"][i]
            X[:d_f_star] = inputs["d_f_star"][i]
        end
        if settings.core
            if settings.method_core_turb == "GE"
                X[:mdoti_c_star] = inputs["mdoti_c_star"][i]
                X[:Tti_c_star] = inputs["Tti_c_star"][i]
                X[:Ttj_c_star] = inputs["Ttj_c_star"][i]
                X[:Pti_c_star] = inputs["Pti_c_star"][i]
                X[:DTt_des_c_star] = inputs["DTt_des_c_star"][i]
            elseif settings.method_core_turb == "PW"
                X[:mdoti_c_star] = inputs["mdoti_c_star"][i]
                X[:Tti_c_star] = inputs["Tti_c_star"][i]
                X[:Ttj_c_star] = inputs["Ttj_c_star"][i]
                X[:Pti_c_star] = inputs["Pti_c_star"][i]
                X[:rho_te_c_star] = inputs["rho_te_c_star"][i]
                X[:c_te_c_star] = inputs["c_te_c_star"][i]
                X[:rho_ti_c_star] = inputs["rho_ti_c_star"][i]
                X[:c_ti_c_star] = inputs["c_ti_c_star"][i]
            else
                throw(DomainError("Invalid method for turbine noise in core module. Specify: GE/PW."))
            end
        end

        if settings.jet_mixing
            X[:V_j_star] = inputs["V_j_star"][i]
            X[:rho_j_star] = inputs["rho_j_star"][i]
            X[:A_j_star] = inputs["A_j_star"][i]
            X[:Tt_j_star] = inputs["Tt_j_star"][i]
        end
        if settings.jet_shock
            X[:V_j_star] = inputs["V_j_star"][i]
            X[:M_j] = inputs["M_j"][i]
            X[:A_j_star] = inputs["A_j_star"][i]
            X[:Tt_j_star] = inputs["Tt_j_star"][i]
        end

        if settings.airframe
            X[:theta_flaps] = inputs["theta_flaps"][i]
            X[:I_landing_gear] = inputs["I_landing_gear"][i]
        end

        Y[:msap_source] = zeros(Float64, settings.N_f)
        source_fwd! = (y,x)->self.source!(settings, data, ac, shield[i, 2:end], n_t, idx_src, y, x)
        jacobian!(J, source_fwd!, Y, X)

        @. partials["msap_source", "TS"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :TS]
        @. partials["msap_source", "M_0"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :M_0]
        @. partials["msap_source", "c_0"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :c_0]
        @. partials["msap_source", "rho_0"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :rho_0]
        @. partials["msap_source", "mu_0"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :mu_0]
        @. partials["msap_source", "T_0"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :T_0]
        @. partials["msap_source", "theta"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :theta]
        @. partials["msap_source", "phi"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :phi]

        if settings.fan_inlet
            @. partials["msap_source", "DTt_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :DTt_f_star]
            @. partials["msap_source", "mdot_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :mdot_f_star]
            @. partials["msap_source", "N_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :N_f_star]
            @. partials["msap_source", "A_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :A_f_star]
            @. partials["msap_source", "d_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :d_f_star]
        end

        if settings.fan_discharge
            @. partials["msap_source", "DTt_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :DTt_f_star]
            @. partials["msap_source", "mdot_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :mdot_f_star]
            @. partials["msap_source", "N_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :N_f_star]
            @. partials["msap_source", "A_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :A_f_star]
            @. partials["msap_source", "d_f_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :d_f_star]
        end

        if settings.core
            if settings.method_core_turb == "GE"
                @. partials["msap_source", "mdoti_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :mdoti_c_star]
                @. partials["msap_source", "Tti_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :Tti_c_star]
                @. partials["msap_source", "Ttj_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :Ttj_c_star]
                @. partials["msap_source", "Pti_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :Pti_c_star]
                @. partials["msap_source", "DTt_des_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :DTt_des_c_star]

            elseif settings.method_core_turb == "GE"
                @. partials["msap_source", "mdoti_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :mdoti_c_star]
                @. partials["msap_source", "Tti_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :Tti_c_star]
                @. partials["msap_source", "Ttj_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :Ttj_c_star]
                @. partials["msap_source", "Pti_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :Pti_c_star]
                @. partials["msap_source", "rho_te_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :rho_te_c_star]
                @. partials["msap_source", "c_te_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :c_te_c_star]
                @. partials["msap_source", "rho_ti_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :rho_ti_c_star]
                @. partials["msap_source", "c_ti_c_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :c_ti_c_star]
            else
                throw(DomainError("Invalid method for turbine noise in core module. Specify: GE/PW."))
            end
        end

        if settings.jet_mixing == true && settings.jet_shock == false
            @. partials["msap_source", "V_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :V_j_star]
            @. partials["msap_source", "rho_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :rho_j_star]
            @. partials["msap_source", "A_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :A_j_star]
            @. partials["msap_source", "Tt_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :Tt_j_star]
        elseif settings.jet_shock == true && settings.jet_mixing == false
            @. partials["msap_source", "V_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :V_j_star]
            @. partials["msap_source", "M_j"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :M_j]
            @. partials["msap_source", "A_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :A_j_star]
            @. partials["msap_source", "Tt_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :Tt_j_star]
        elseif settings.jet_shock ==true && settings.jet_mixing == true
            @. partials["msap_source", "V_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :V_j_star]
            @. partials["msap_source", "rho_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :rho_j_star]
            @. partials["msap_source", "A_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :A_j_star]
            @. partials["msap_source", "Tt_j_star"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :Tt_j_star]
            @. partials["msap_source", "M_j"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :M_j]
        end

        if settings.airframe
            @. partials["msap_source", "theta_flaps"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :theta_flaps]
            @. partials["msap_source", "I_landing_gear"][(i-1)*settings.N_f+1:i*settings.N_f] = J[:msap_source, :I_landing_gear]
        end
    end
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


max_stars_repo_name: dillondaudert/JuliaMotifs

max_issues_repo_name: dillondaudert/JuliaMotifs

max_forks_repo_name: dillondaudert/JuliaMotifs


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# contributed by David Campbell
# based on the Go version
# modified by Jarrett Revels, Alex Arslan, Yichao Yu
#
# Bit-twiddle optimizations added by Kristoffer Carlsson

using Distributed
using Printf

addprocs(4)

@everywhere begin
const NucleotideLUT = zeros(UInt8, 256)
NucleotideLUT['A'%UInt8] = 0
NucleotideLUT['C'%UInt8] = 1
NucleotideLUT['G'%UInt8] = 2
NucleotideLUT['T'%UInt8] = 3
NucleotideLUT['a'%UInt8] = 0
NucleotideLUT['c'%UInt8] = 1
NucleotideLUT['g'%UInt8] = 2
NucleotideLUT['t'%UInt8] = 3

struct KNucleotides{L, T}
    i::T
end
Base.hash(kn::KNucleotides, h::UInt64) = hash(kn.i, h)
Base.isequal(kn1::KNucleotides, kn2::KNucleotides) = kn1.i == kn2.i
Base.show(io::IO, kn::KNucleotides) = print(io, '[', string(kn), ']')
function determine_inttype(l)
    l <= 4 && return UInt8
    l <= 8 && return UInt16
    l <= 16 && return UInt32
    l <= 32 && return UInt64
    error("invalid length")
end

function KNucleotides{L, T}(str::String) where {L, T}
    i = T(0)
    @inbounds for j in 1:L
        b = codeunit(str, j)
        i = (i << 2) | NucleotideLUT[b]
    end
    return KNucleotides{L, T}(i)
end

@inline function shift(kn::KNucleotides{L, T}, c::UInt8) where {L, T}
    i = kn.i
    i &= (~(3 << 2(L-1)) % T)
    KNucleotides{L, T}((i << T(2)) | @inbounds NucleotideLUT[c])
end

function Base.string(kn::KNucleotides{L}) where {L}
    sprint() do io
        for j in 1:L
            mask = 3 << (2(L-j))
            z = (kn.i & mask) >> 2(L-j)
            write(io,
                z == 0 ? 'A' :
                z == 1 ? 'C' :
                z == 2 ? 'G' :
                z == 3 ? 'T' :
                error())
        end
    end
end

function count_data(data::String, ::Type{KNucleotides{L, T}}) where {L, T}
    counts = Dict{KNucleotides{L, T}, Int}()
    kn = KNucleotides{L, T}(data)
    counts[kn] = 1
    @inbounds for offset = (L+1):sizeof(data)
        c = codeunit(data, offset)
        kn = shift(kn, c)
        token = Base.ht_keyindex2!(counts, kn)
        if token > 0
            counts.vals[token] += 1
        else
            Base._setindex!(counts, 1, kn, -token)
        end
    end
    return counts
end

function count_one(data::String, s::String)
    L = length(s)
    K = KNucleotides{L, determine_inttype(L)}
    d = count_data(data, K)
    return get(d, K(s), 0)
end

struct KNuc
    name::String
    count::Int
end

# sort down
function Base.isless(x::KNuc, y::KNuc)
    if x.count == y.count
        return x.name > y.name
    end
    x.count > y.count
end

function sorted_array(m)
    kn = Vector{KNuc}(undef, length(m))
    i = 1
    for (k, v) in m
        kn[i] = KNuc(string(k), v)
        i += 1
    end
    sort!(kn)
end

do_work(str::String, i::Int) = sorted_array(count_data(str, KNucleotides{i,determine_inttype(i)}))
do_work(str::String, i::String) = count_one(str, i)

end # @everywhere

function print_knucs(a::Array{KNuc, 1})
    sum = 0
    for kn in a
        sum += kn.count
    end
    for kn in a
        @printf("%s %.3f\n", kn.name, 100.0kn.count/sum)
    end
    println()
end

function perf_k_nucleotide(io = stdin)
    three = ">THREE "
    while true
        line = readline(io)
        if startswith(line, three)
            break
        end
    end
    data = read(io, String)
    str = filter(!isequal('\n'), data)

    vs = [1, 2, "GGT", "GGTA", "GGTATT", "GGTATTTTAATT", "GGTATTTTAATTTATAGT"]
    results = Vector{Any}(undef, length(vs))
    order = collect(enumerate(vs))

    @sync for w in workers()
        @async while !isempty(order)
            v = pop!(order)
            r = remotecall_fetch(do_work, w, str, v[2])
            results[v[1]] = r
        end
    end

    for (v, result) in zip(vs, results)
        if result isa Array
            print_knucs(result)
        end
        if result isa Int
            @printf("%d\t%s\n", result, v)
        end
    end
end

perf_k_nucleotide()