text:
using Jive # runtests
runtests(@__DIR__, skip=["Example", "errors", "jive/onlyonce/heavy.jl", "jive/__END__/included.jl", "jive/__REPL__"])


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Deep learning with Flux")

#

#
# ## Neural networks
#
# Now that we know what neurons are, we are ready for the final step: the neural network!. A neural network is literally made out of a network of neurons that are connected together.
#
# So far, we have just looked at single neurons, that only have a single output.
# What if we want multiple outputs?
#
#
# ### Multiple output models
#
# What if we wanted to distinguish between apples, bananas, *and* grapes? We could use *vectors* of `0` or `1` values to symbolize each output.
#
#
#
# The idea of using vectors is that different directions in the space of outputs encode information about different types of inputs.

#-

# Now we extend our previous model to give multiple outputs by repeating it with different weights. For the first element of the array we'd use:
#
# $$\sigma(x;w^{(1)},b^{(1)}) := \frac{1}{1 + \exp(-w^{(1)} \cdot x + b^{(1)})};$$
#
# then for the second we'd use
#
# $$\sigma(x;w^{(2)},b^{(2)}) := \frac{1}{1 + \exp(-w^{(2)} \cdot x + b^{(2)})};$$
#
# and if you wanted $n$ outputs, you'd have for each one
#
# $$\sigma(x;w^{(i)},b^{(i)}) := \frac{1}{1 + \exp(-w^{(i)} \cdot x + b^{(i)})}.$$

#-

# Notice that these equations are all the same, except for the parameters, so we can write this model more succinctly, as follows. Let's write $b$ in an array:
#
# $$b=\left[\begin{array}{c}
# b_{1}\\
# b_{2}\\
# \vdots\\
# b_{n}
# \end{array}\right]$$
#
# and put our array of weights as a matrix:
#
# $$ \mathsf{W}=\left[\begin{array}{c}
# \\
# \\
# \\
# \\
# \end{array}\begin{array}{cccc}
# w_{1}^{(1)} & w_{2}^{(1)} & \ldots & w_{n}^{(1)}\\
# w_{1}^{(2)} & w_{2}^{(2)} & \ldots & w_{n}^{(2)}\\
# \vdots & \vdots &  & \vdots\\
# w_{1}^{(n)} & w_{2}^{(n)} & \ldots & w_{n}^{(n)}
# \end{array}\right]
# $$
#
# We can write this all in one line as:
#
# $$\sigma(x;w,b)= \left[\begin{array}{c}
# \sigma^{(1)}\\
# \sigma^{(2)}\\
# \vdots\\
# \sigma^{(n)}
# \end{array}\right] = \frac{1}{1 + \exp(-\mathsf{W} x + b)}$$
#
# $\mathsf{W} x$ is the operation called "matrix multiplication"
#
# [Show small matrix multiplication]

W = [10 1;
     20 2;
     30 3]
x = [3;
     2]
W*x

# It takes each column of weights and does the dot product against $x$ (remember, that's how $\sigma^{(i)}$ was defined) and spits out a vector from doing that with each column. The result is a vector, which makes this version of the function give a vector of outputs which we can use to encode larger set of choices.
#
# Matrix multiplication is also interesting since **GPUs (Graphics Processing Units, i.e. graphics cards) are basically just matrix multiplication machines**, which means that by writing the equation this way, the result can be calculated really fast.

#-

# This "multiple input and multiple output" version of the sigmoid function is known as a *layer of neurons*.
#
# Previously we worked with a single neuron, which we visualized as
#
#
#
# where we have two pieces of data (green) coming into a single neuron (pink) that returned a single output. We could use this single output to do binary classification - to identify an image of a fruit as `1`, meaning banana or as `0`, meaning not a banana (or an apple).
#
# To do non-binary classification, we can use a layer of neurons, which we can visualize as
#
#
#
# We now have stacked a bunch of neurons on top of each other to hopefully work together and train to output results of more complicated features.
#
# We still have two input pieces of data, but now have several neurons, each of which produces an output for a given binary classification:
# * neuron 1: "is it an apple?"
# * neuron 2: "is it a banana?"
# * neuron 3: "is it a grape?"

#-

# # Multiple outputs with Flux.jl

#-

# First step: load the data.

using CSV, DataFrames, Flux, Plots
## Load apple data in CSV.read for each file
apples1 = DataFrame(CSV.File(datapath("data/Apple_Golden_1.dat"), delim='\t', allowmissing=:none, normalizenames=true))
apples2 = DataFrame(CSV.File(datapath("data/Apple_Golden_2.dat"), delim='\t', allowmissing=:none, normalizenames=true))
apples3 = DataFrame(CSV.File(datapath("data/Apple_Golden_3.dat"), delim='\t', allowmissing=:none, normalizenames=true))
## And then concatenate them all together
apples = vcat(apples1, apples2, apples3)
bananas = DataFrame(CSV.File(datapath("data/Banana.dat"), delim='\t', allowmissing=:none, normalizenames=true))
grapes1 = DataFrame(CSV.File(datapath("data/Grape_White.dat"), delim='\t', allowmissing=:none, normalizenames=true))
grapes2 = DataFrame(CSV.File(datapath("data/Grape_White_2.dat"), delim='\t', allowmissing=:none, normalizenames=true))
grapes = vcat(grapes1, grapes2)

#-

## Extract out the features and construct the corresponding labels
x_apples  = [ [apples[i, :red], apples[i, :blue]] for i in 1:size(apples, 1) ]
x_bananas  = [ [bananas[i, :red], bananas[i, :blue]] for i in 1:size(bananas, 1) ]
x_grapes = [ [grapes[i, :red], grapes[i, :blue]] for i in 1:size(grapes, 1) ]
xs = vcat(x_apples, x_bananas, x_grapes)
ys = vcat(fill([1,0,0], size(x_apples)),
          fill([0,1,0], size(x_bananas)),
          fill([0,0,1], size(x_grapes)))

# ### One-hot vectors

#-

# Recall:
#
#

#-

# `Flux.jl` provides an efficient representation for one-hot vectors, using advanced features of Julia so that it does not actually store these vectors, which would be a waste of memory; instead `Flux` just records in which position the non-zero element is. To us, however, it looks like all the information is being stored:

using Flux: onehot

onehot(2, 1:3)

#-

ys = vcat(fill(onehot(1, 1:3), size(x_apples)),
          fill(onehot(2, 1:3), size(x_bananas)),
          fill(onehot(3, 1:3), size(x_grapes)))

# ## The core algorithm from the previous lecture

## model = Dense(2, 1, σ)
## L(x,y) = Flux.mse(model(x), y)
## opt = SGD(params(model))
## Flux.train!(L, zip(xs, ys), opt)

# ### Visualization

using Plots
plot()

contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[1], levels=[0.5, 0.51], color = cgrad([:blue, :blue]))
contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[2], levels=[0.5,0.51], color = cgrad([:green, :green]))
contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[3], levels=[0.5,0.51], color = cgrad([:red, :red]))

scatter!(first.(x_apples), last.(x_apples), m=:cross, label="apples", color = :blue)
scatter!(first.(x_bananas), last.(x_bananas), m=:circle, label="bananas", color = :green)
scatter!(first.(x_grapes), last.(x_grapes), m=:square, label="grapes", color = :red)

ext: jl

max_stars_repo_path: Courses/Deep learning with Flux/1500.Intro-to-neural-networks.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
const θⁱ = π/2 # incidence angle, measured counter-clockwise from +x direction
const θᵗ = π/3 # desired transmitted field angle
const d = cos(θᵗ)-cos(θⁱ) # design parameter
const k₀ = 10.0 # wavenumber that attains the metasurface design goals
const λ₀ = 2π/k₀ # wavelength that attains the metasurface design goals
const L₀ = λ₀/abs(d) # period of metasurface

η(x) = sin(abs(θⁱ))*(1+exp.(im*k₀*d*x))
μ(x) = sin(abs(θⁱ))*(1-exp.(im*k₀*d*x))

const L=L₀ # length of domain (some integer multiple of period of metasurface)
const k=1.1k₀ # wavenumber of incident field

const α = k * cos(θⁱ)
const β = k * sin(θⁱ)
const i = 1 # == 1 for θⁱ = π/2
const exponents = 5:9 # defines the number of grid points, base 2


text:
@testset "incomplete beta" begin
#Compared with scipy.special.betainc
    ctr = 1
    ans1 = [0.063768560858519868, 0.090334470601733122, 0.11082468660445946, 0.12818843369794991, 0.14356629312870628, 0.15754242542404437, 0.17046342838319906, 0.18255489099382846, 0.19397336804135662, 0.2048327646991335, 0.21521902554364258, 0.22519890061842537, 0.23482546905318047, 0.24414177534616024, 0.253183311106635, 0.26197976086890934, 0.27055626226830082, 0.27893433611357099, 0.28713258625741256, 0.29516723530086664, 0.30305254089347, 0.31080112365240059, 0.3184242286319004, 0.32593193612447635, 0.33333333333333337, 0.34063665547711719, 0.34784940276303394, 0.35497843812588953, 0.36203006950165795, 0.36901011956554541, 0.375923985234415, 0.3827766887550389, 0.38957292183289544, 0.39631708397254195, 0.40301331597932183, 0.409665529398267, 0.41627743252839544, 0.42285255354095019, 0.42939426114224771, 0.43590578315102518, 0.44239022330290334, 0.44885057654807836, 0.45528974307049164, 0.46171054122585353, 0.46811571957074011, 0.47450796813452323, 0.48088992906936007, 0.48726420680024518, 0.49363337778673011, 0.50000000000000011, 0.50636662221327011, 0.51273579319975504, 0.5191100709306401, 0.52549203186547688, 0.53188428042926006, 0.53828945877414658, 0.54471025692950836, 0.55114942345192175, 0.55760977669709688, 0.56409421684897509, 0.57060573885775256, 0.57714744645905003, 0.58372256747160478, 0.59033447060173316, 0.59698668402067834, 0.60368291602745827, 0.61042707816710473, 0.61722331124496133, 0.62407601476558505, 0.63098988043445459, 0.63796993049834205, 0.64502156187411064, 0.65215059723696633, 0.65936334452288281, 0.66666666666666663, 0.67406806387552365, 0.68157577136809966, 0.68919887634759947, 0.69694745910653, 0.70483276469913336, 0.71286741374258733, 0.72106566388642879, 0.72944373773169924, 0.73802023913109061, 0.74681668889336494, 0.75585822465383967, 0.76517453094681931, 0.7748010993815746, 0.78478097445635719, 0.79516723530086619, 0.8060266319586431, 0.81744510900617107, 0.82953657161680039, 0.84245757457595483, 0.8564337068712925, 0.87181156630205003, 0.88917531339554046, 0.90966552939826684, 0.9362314391414801]
    ans2 = [0.087041654488676928, 0.12317839372423395, 0.15096439050993934, 0.1744376469592186, 0.19516124171507579, 0.21393635267239447, 0.23123866443564164, 0.24737796264883508, 0.26256933221933759, 0.27696931578459572, 0.29069605231718387, 0.30384130015177169, 0.31647801578481272, 0.32866534497198974, 0.34045203033834726, 0.35187880931248239, 0.36298014555988484, 0.37378550722502435, 0.38432032902829749, 0.39460674882922503, 0.40466418009101346, 0.41450976283569702, 0.42415872320430908, 0.43362466329673527, 0.44291979714289581, 0.4520551445675528, 0.46104069179265272, 0.46988552550762908, 0.47859794558683333, 0.48718556048052553, 0.49565536843945707, 0.50401382707500741, 0.51226691325208229, 0.52042017492130455, 0.52847877619206929, 0.53644753670811252, 0.54433096619706123, 0.55213329491367624, 0.55985850057457809, 0.56751033228373404, 0.57509233186792286, 0.58260785297598572, 0.59006007824198947, 0.59745203476815989, 0.60478660814681062, 0.61206655521005915, 0.61929451567078364, 0.62647302279711903, 0.63360451324511835, 0.64069133615944129, 0.64773576163960578, 0.65473998865910954, 0.66170615251627152, 0.66863633188876026, 0.67553255555825598, 0.68239680886742871, 0.68923103996827928, 0.69603716591885445, 0.70281707868434762, 0.70957265109866863, 0.71630574284372639, 0.72301820650601145, 0.72971189377369172, 0.73638866184254259, 0.74305038010582547, 0.74969893721203962, 0.75633624858569781, 0.76296426452045429, 0.76958497897174627, 0.77620043919849924, 0.78281275643163062, 0.78942411778262311, 0.79603679965050977, 0.80265318294305021, 0.80927577050162136, 0.81590720721472276, 0.82255030342945812, 0.82920806243432144, 0.83588371300491959, 0.84258074829830698, 0.84930297278272149, 0.85605455944435893, 0.86284012029219526, 0.86966479429526578, 0.876534358507564, 0.8834553705447894, 0.89043535424294118, 0.89748304605629159, 0.90460872897483335, 0.91182469611459804, 0.91914591281205937, 0.92659099459483496, 0.93418371187592941, 0.94195542513481079, 0.94994928861995065, 0.95822815489579449, 0.96689132900746566, 0.97611719389699281, 0.98631372146716156]
    ans3 = [0.012873483750773953, 0.024045778495719504, 0.034668868563965893, 0.044957919792848362, 0.05501104929272229, 0.064884481926849022, 0.074614610511183635, 0.084226884888052606, 0.093740088419712889, 0.10316864948642543, 0.11252399935922629, 0.12181542160322313, 0.13105061024399095, 0.14023605142518913, 0.14937729307408915, 0.15847914076041122, 0.16754580333263841, 0.1765810034326174, 0.18558806286120144, 0.19456996956034567, 0.20352943091060927, 0.21246891667645043, 0.22139069400685801, 0.23029685625993337, 0.23918934697038266, 0.24806997995710448, 0.2569404563342797, 0.26580237901716475, 0.27465726518533057, 0.28350655706915756, 0.29235163135146713, 0.30119380741925028, 0.31003435465621887, 0.31887449893226855, 0.32771542841860729, 0.33655829883559579, 0.34540423822300514, 0.3542543513084796, 0.36310972353877363, 0.37197142482927481, 0.38084051307999506, 0.38971803750029138, 0.39860504177981332, 0.40750256713936639, 0.41641165529238072, 0.42533335134536215, 0.43426870666399331, 0.44321878173037438, 0.45218464901620103, 0.46116739589644162, 0.47016812762827714, 0.47918797042070466, 0.48822807462128687, 0.49728961804808525, 0.50637380949686728, 0.51548189245630038, 0.5246151490670925, 0.53377490436500419, 0.5429625308524717, 0.55217945344936725, 0.56142715488038319, 0.57070718156489009, 0.58002115008515442, 0.58937075432090513, 0.59875777335283087, 0.60818408025528292, 0.61765165191996785, 0.62716258007870107, 0.63671908372557429, 0.6463235231787503, 0.65597841607162255, 0.66568645562500284, 0.67545053162996704, 0.68527375466991214, 0.6951594842368839, 0.70511136156043297, 0.71513334817976937, 0.72522977156950186, 0.73540537950094409, 0.74566540532108994, 0.75601564701303459, 0.76646256384385103, 0.77701339572859029, 0.78767631232750368, 0.79846060164096611, 0.80937691194811789, 0.82043756714268778, 0.83165698521328602, 0.84305224521672972, 0.85464387408595721, 0.86645696969790587, 0.87852285861624935, 0.89088164472968334, 0.90358633052105619, 0.91670992507941651, 0.93035879743560423, 0.94470095224765227, 0.96003789419656471, 0.97705809868034199]
    ans4 = [5.6598501882086638e-163, 1.2793886161240335e-138, 2.2576938967074611e-124, 2.8979946053659182e-114, 2.0091646864167303e-106, 5.1249028123107484e-100, 1.3384365465929538e-94, 6.5928463236275696e-90, 9.0834048016126672e-86, 4.5811711646109305e-82, 1.024712535445947e-78, 1.1712904358478081e-75, 7.61899830618238e-73, 3.0663999960291821e-70, 8.1593401517362853e-68, 1.5142302711017038e-65, 2.0475200955891378e-63, 2.0916713536958276e-61, 1.6640409038082517e-59, 1.057759838815703e-57, 5.4913266744369942e-56, 2.3725786887367975e-54, 8.6718335278405774e-53, 2.7197933177013249e-51, 7.4118251807120739e-50, 1.7744725579765609e-48, 3.7689144668459763e-47, 7.1639177579230802e-46, 1.2281691523716354e-44, 1.9124082898174377e-43, 2.7218251674439459e-42, 3.5610505888183843e-41, 4.3051319713293141e-40, 4.8320967998003304e-39, 5.0570559392737623e-38, 4.954351648811439e-37, 4.5601168072923055e-36, 3.9565121052116655e-35, 3.2458692695921472e-34, 2.5250201231776632e-33, 1.8674800424552049e-32, 1.31632743396051e-31, 8.862853181709709e-31, 5.71217922566864e-30, 3.5310511291576422e-29, 2.0973893941833731e-28, 1.1991528375926732e-27, 6.6098385509764045e-27, 3.5179140837714274e-26, 1.8104027085942885e-25, 9.0207305335112628e-25, 4.3574464645305935e-24, 2.0429719208160492e-23, 9.307222587504958e-23, 4.124448596486994e-22, 1.7796548980836483e-21, 7.4841565165421277e-21, 3.0702859704148574e-20, 1.2297464356858357e-19, 4.8128959875339663e-19, 1.8419917152708652e-18, 6.8988896144258444e-18, 2.5303817306725186e-17, 9.094925675449239e-17, 3.2055039961949189e-16, 1.1085187441597278e-15, 3.7635118688411301e-15, 1.255134520239119e-14, 4.1140267831674448e-14, 1.3260121388841605e-13, 4.2048092809091946e-13, 1.3124119186938038e-12, 4.0338344386888606e-12, 1.2214673719698131e-11, 3.6454138300247904e-11, 1.0727343115700782e-10, 3.1138059245591031e-10, 8.9189474914175862e-10, 2.5218693952601578e-09, 7.0417217701699225e-09, 1.9424040709629306e-08, 5.2949202359208791e-08, 1.4268930041189607e-07, 3.8026633846818863e-07, 1.0025399986934807e-06, 2.6156990305099061e-06, 6.756237677059056e-06, 1.7282977043885323e-05, 4.3803033456977776e-05, 0.00011004048858461744, 0.0002741427949192326, 0.00067767979451775295, 0.0016633950407366552, 0.0040576338171432955, 0.009848744173397107, 0.023828559323224537, 0.057638830150965534, 0.14020121761365209, 0.34813332055564178]
    ans5 = [0.013962049747018378, 0.042738447399963485, 0.080235508245346435, 0.12332127072222762, 0.16988388360150664, 0.21838766458819614, 0.26767957944384452, 0.31688247132531383, 0.36532771975839284, 0.41250897644193646, 0.45804848436293449, 0.50167153949877097, 0.54318655272778527, 0.58246914958502727, 0.61944929138313776, 0.65410072456419888, 0.68643226670884316, 0.71648056882768607, 0.74430408220004141, 0.76997801988794912, 0.79359014748064294, 0.8152372703488121, 0.83502230934244948, 0.85305187581929465, 0.86943427172671017, 0.88427785226197508, 0.89768969816152788, 0.90977455245768279, 0.92063398298239707, 0.93036573727946092, 0.93906326112402183, 0.94681535570471753, 0.95370595182410001, 0.95981398231488424, 0.96521333632950712, 0.96997288129947024, 0.97415654022800602, 0.97782341361450698, 0.9810279367445075, 0.98382006434176061, 0.98624547569157095, 0.98834579432593017, 0.99015881722713328, 0.99171874927107939, 0.993056439306095, 0.99419961485797559, 0.99517311297582334, 0.99599910519385948, 0.99669731498846048, 0.99728522646317797, 0.99777828330274299, 0.9981900773047403, 0.99853252602897213, 0.99881603930330098, 0.99904967449435056, 0.99924128059493644, 0.99939763130024395, 0.99952454734409613, 0.99962700844743568, 0.99970925529544441, 0.99977488200944087, 0.99982691961654058, 0.99986791104559958, 0.99989997819362175, 0.99992488161389415, 0.9999440733768028, 0.99995874364766302, 0.99996986151394762, 0.9999782105779107, 0.99998441981059161, 0.99998899014027587, 0.99999231722334625, 0.99999471081867697, 0.99999641115882809, 0.99999760268277504, 0.99999842546616124, 0.99999898465647652, 0.99999935819244778, 0.99999960305956725, 0.99999976030730964, 0.99999985902839927, 0.99999991947664768, 0.99999995547750509, 0.9999999762646713, 0.99999998785693356, 0.99999999407190587, 0.99999999725752275, 0.99999999880800261, 0.99999999951849838, 0.99999999982174348, 0.99999999994062116, 0.99999999998264033, 0.99999999999569866, 0.9999999999991418, 0.99999999999987266, 0.99999999999998768, 0.99999999999999944, 0.99999999999999989, 0.99999999999999989]
    ans6 = [3.3560143825993687e-143, 2.213518371644144e-113, 4.0145526832291446e-96, 5.2109131251766219e-84, 1.0203891990294805e-74, 3.2677379646480312e-67, 6.1015737258099584e-61, 1.4185802923209925e-55, 6.6985847646944376e-51, 8.9719726417923444e-47, 4.3251435416063278e-43, 8.9466478285748974e-40, 9.0742566193142145e-37, 5.0057011672920701e-34, 1.6303369052296386e-31, 3.3493886577780779e-29, 4.5813387723702873e-27, 4.3628327684456486e-25, 3.0028168459240925e-23, 1.5416424964057036e-21, 6.0647292131758379e-20, 1.8708373622428274e-18, 4.6165844488973474e-17, 9.2725535948439169e-16, 1.5390751973611367e-14, 2.1394034328449769e-13, 2.5200592719451178e-12, 2.5418366766245469e-11, 2.215842537778695e-10, 1.6834088907362684e-09, 1.122874450436064e-08, 6.6201731659826321e-08, 3.4707638415115294e-07, 1.6269341185131871e-06, 6.8526269592546404e-06, 2.6052419808499182e-05, 8.97713031375554e-05, 0.0002814332703651891, 0.00080554169029503023, 0.0021120197810150015, 0.0050879707571202049, 0.01129547597461162, 0.023174688420508269, 0.044064890651698461, 0.077870360082158169, 0.12827033082474237, 0.19756362321699061, 0.28548625344477141, 0.38850225052293164, 0.50000000000002465, 0.61149774947706836, 0.71451374655522859, 0.80243637678300939, 0.87172966917525962, 0.92212963991784314, 0.95593510934830217, 0.97682531157949204, 0.98870452402538833, 0.99491202924287991, 0.99788798021898495, 0.99919445830970499, 0.99971856672963477, 0.99991022869686241, 0.99997394758019154, 0.99999314737304079, 0.99999837306588146, 0.99999965292361581, 0.99999993379826835, 0.99999998877125551, 0.99999999831659114, 0.99999999977841569, 0.99999999997458167, 0.9999999999974799, 0.99999999999978606, 0.99999999999998457, 0.99999999999999911, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989, 0.99999999999999989]
    for x = 0.01:0.01:0.99
        @test beta_inc(0.5,0.5,x,1.0-x)[1] ≈ ans1[ctr]#scipy.special.betainc(0.5,0.5,x)
        ctr += 1
    end
    ctr = 1
    for x = 0.01:0.01:0.99
        @test beta_inc(0.5,0.8,x,1.0-x)[1] ≈ ans2[ctr]#scipy.special.betainc(0.5,0.8,x)
        ctr += 1
    end
    ctr = 1
    for x = 0.01:0.01:0.99
        @test beta_inc(0.9,0.8,x,1.0-x)[1] ≈ ans3[ctr]#scipy.special.betainc(0.9,0.8,x)
        ctr += 1
    end
    ctr = 1
    for x = 0.01:0.01:0.99
        @test beta_inc(80.9,0.8,x,1.0-x)[1] ≈ ans4[ctr]#scipy.special.betainc(80.9,0.8,x)
        ctr += 1
    end
    ctr = 1
    for x = 0.01:0.01:0.99
        @test beta_inc(1.7,10.5,x,1.0-x)[1] ≈ ans5[ctr]#scipy.special.betainc(1.7,10.5,x)
        ctr += 1
    end
    ctr = 1
    for x = 0.01:0.01:0.99
        @test beta_inc(100.5,100.5,x,1.0-x)[1] ≈ ans6[ctr]#scipy.special.betainc(100.5,100.5,x)
        ctr += 1
    end
    @test beta_inc(1.5,200.5,0.07,0.93)[1] ≈ 0.99999790408564
    @test SpecialFunctions.loggammadiv(13.89, 21.0001) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89)))
    @test SpecialFunctions.stirling_corr(11.99, 100.1) ≈ SpecialFunctions.stirling(11.99) + SpecialFunctions.stirling(100.1) - SpecialFunctions.stirling(11.99 + 100.1)
end
@testset "inverse of incomplete beta" begin
    f(a,b,p) = beta_inc_inv(a,b,p)[1]
    @test f(.5,.5,0.6376856085851985E-01) ≈ 0.01
    @test f(.5,.5,0.20483276469913355) ≈ 0.1
    @test f(.5,.5,1.0000) ≈ 1.0000
    @test f(1.0,.5,0.0) ≈ 0.00
    @test f(1.0,.5,0.5012562893380045E-02) ≈ 0.01
    @test f(1.0,.5,0.5131670194948620E-01) ≈ 0.1
    @test f(1.0,.5, 0.2928932188134525) ≈ 0.5
    @test f(1.0,1.0,.5) ≈ 0.5
    @test f(2.0,2.0,.028) ≈ 0.1
    @test f(2.0,2.0,0.104) ≈ 0.2
    @test f(2.0,2.0,.216) ≈ 0.3
    @test f(2.0,2.0,.352) ≈ 0.4
    @test f(2.0,2.0,.5) ≈ 0.5
    @test f(2.0,2.0,0.648) ≈ 0.6
    @test f(2.0,2.0,0.784) ≈ 0.7
    @test f(2.0,2.0,0.896) ≈ 0.8
    @test f(2.0,2.0,.972) ≈ 0.9
    @test f(5.5,5.0,0.4361908850559777) ≈ .5
    @test f(10.0,.5,0.1516409096347099) ≈ 0.9
    @test f(10.0,5.0,0.8978271484375000E-01) ≈ 0.5
    @test f(10.0,5.0,1.00) ≈ 1.00
    @test f(10.0,10.0,.5) ≈ .5
    @test f(20.0,5.0,0.4598773297575791) ≈ 0.8
    @test f(20.0,10.0,0.2146816102371739) ≈ 0.6
    @test f(20.0,10.0,0.9507364826957875) ≈ 0.8
    @test f(20.0,20.0,.5) ≈ .5
    @test f(20.0,20.0,0.8979413687105918) ≈ 0.6
    @test f(30.0,10.0,0.2241297491808366) ≈ 0.7
    @test f(30.0,10.0,0.7586405487192086) ≈ 0.8
    @test f(40.0,20.0,0.7001783247477069) ≈ 0.7
    @test f(1.0,0.5,0.5131670194948620E-01) ≈ 0.1
    @test f(1.0,0.5,0.1055728090000841) ≈ 0.2
    @test f(1.0,0.5,0.1633399734659245) ≈ 0.3
    @test f(1.0,0.5,0.2254033307585166) ≈ 0.4
    @test f(1.0,2.0,.36) ≈ 0.2
    @test f(1.0,3.0,.488) ≈ 0.2
    @test f(1.0,4.0,.5904) ≈ 0.2
    @test f(1.0,5.0,.67232) ≈ 0.2
    @test f(2.0,2.0,.216) ≈ 0.3
    @test f(3.0,2.0,0.837e-1) ≈ 0.3
    @test f(4.0,2.0,0.3078e-1) ≈ 0.3
    @test f(5.0,2.0,0.10935e-1) ≈ 0.3
    @test f(1.30625000,11.75620000,0.9188846846205182) ≈ 0.225609
    @test f(1.30625000,11.75620000,0.21053116418502513) ≈ 0.033557
    @test f(1.30625000,11.75620000,0.18241165418408148) ≈ 0.029522
end


text:
module ReusePatterns
using InteractiveUtils
using Combinatorics

# TODO
# - forward macro calls

export supertypes, forward, @forward,
    @copy_fields,
    @quasiabstract, concretetype, isquasiabstract, isquasiconcrete

"""
`supertypes(T::Type)`

Returns a vector with all supertypes of type `T` (excluding `Any`).

# Example
```julia-repl
julia> println(supertypes(Int))
Type[Signed, Integer, Real, Number]
```
"""
function supertypes(T::Type)::Vector{Type}
    out = Vector{Type}()
    st = supertype(T)
    if st != Any
        push!(out, st)
        push!(out, supertypes(st)...)
    end
    return out
end

"""
`forward(sender::Tuple{Type,Symbol}, receiver::Type, method::Method; withtypes=true, allargs=true)`

Return a `Vector{String}` containing the Julia code to properly forward `method` calls to from a `sender` type to a receiver type.

The `sender` tuple must contain a structure type, and a symbol with the name of one of its fields.

The `withtypes` keyword controls whether the forwarding method has type annotated arguments.  The `allargs` keyword controls wether all arguments should be used, or just the first ones up to the last containing the `receiver` type.

Both keywords are `true` by defult, but they can be set to `false` to decrease the number of forwarding methods.

# Example:
Implement a wrapper for an `Int` object, and forward the `+` method accepting `Int`:
```julia-repl
struct Wrapper{T}
    wrappee::T
    extra
    Wrapper{T}(args...; kw...) where T = new(T(args...; kw...), nothing)
end

# Prepare and evaluate forwarding methods:
m = forward((Wrapper, :wrappee), Int, which(+, (Int, Int)))
eval.(Meta.parse.(m))

# Instantiate two wrapped `Int`
i1 = Wrapper{Int}(1)
i2 = Wrapper{Int}(2)

# And add them seamlessly
println(i1 +  2)
println( 1 + i2)
println(i1 + i2)
```
"""
function forward(sender::Tuple{Type,Symbol}, receiver::Type, method::Method;
                 withtypes=true, allargs=true)
    function newmethod(sender_type, sender_symb, argid, method, withtypes, allargs)
        s = "p" .* string.(1:method.nargs-1)
        (withtypes)  &&  (s .*= "::" .* string.(fieldtype.(Ref(method.sig), 2:method.nargs)))
        s[argid] .= "p" .* string.(argid) .* "::$sender_type"
        if !allargs
            s = s[1:argid[end]]
            push!(s, "args...")
        end
        m = string(method.module.eval(:(parentmodule($(method.name))))) * "."
        l = "$m:(" * string(method.name) * ")(" * join(s,", ") * "; kw..."

        m = string(method.module) * "."
        l *= ") = $m:(" * string(method.name) * ")("
        s = "p" .* string.(1:method.nargs-1)
        if !allargs
            s = s[1:argid[end]]
            push!(s, "args...")
        end
        s[argid] .= "getfield(" .* s[argid] .* ", :$sender_symb)"
        l *= join(s, ", ") * "; kw...)"
        l = join(split(l, "#"))
        return l
    end

    sender_type = string(parentmodule(sender[1])) * "." * string(nameof(sender[1]))
    sender_symb = string(sender[2])
    code = Vector{String}()

    # Seacrh for receiver type in method arguments
    foundat = Vector{Int}()
    for i in 2:method.nargs
        argtype = fieldtype(method.sig, i)
        (sender[1] == argtype)  &&  (return code)
        if argtype != Any
            (typeintersect(receiver, argtype) != Union{})  &&  (push!(foundat, i-1))
        end
    end
    (length(foundat) == 0)  &&  (return code)
    if string(method.name)[1] == '@'
        @warn "Forwarding macros is not yet supported."
        display(method)
        println()
        return code
    end

    for ii in combinations(foundat)
        push!(code, newmethod(sender_type, sender_symb, ii, method, withtypes, allargs))
    end

    tmp = split(string(method.module), ".")[1]
    code = "@eval " .* tmp .* " " .* code .*
        " # " .* string(method.file) .* ":" .* string(method.line)
    if  (tmp != "Base")  &&
        (tmp != "Main")
        pushfirst!(code, "using $tmp")
    end
    code = unique(code)
    return code
end

"""
`forward(sender::Tuple{Type,Symbol}, receiver::Type, methods::Vector{Method}; kw...)`
"""
function forward(sender::Tuple{Type,Symbol}, receiver::Type, methods::Vector{Method}; kw...)
    code = Vector{String}()
    for m in methods
        append!(code, forward(sender, receiver, m; kw...))
    end
    code = unique(code)
    return code
end

"""
`forward(sender::Tuple{Type,Symbol}, receivers::Vector{T}, methods; kw...)`
"""
function forward(sender::Tuple{Type,Symbol}, receivers::Vector{T}, methods; kw...) where T <: Type
    code = Vector{String}()
    for t in receivers
        append!(code, forward(sender, t, methods; kw...))
    end
    code = unique(code)
    return code
end

"""
`forward(sender::Tuple{Type,Symbol}, receiver::Type; super=true, kw...)`

Wrapper for `forward(send, receiver, supertypes(receiver, super=super); kw...)`
"""
function forward(sender::Tuple{Type,Symbol}, receiver::Type; super=true, kw...)
    tt = [receiver]
    (super)  &&  (append!(tt, supertypes(receiver)))
    return forward(sender, tt, methodswith(receiver, supertypes=super); kw...)
end

"""
`@forward(sender, receiver, ekws...)`

Evaluate the Julia code to forward methods.  The syntax is exactly the same as the `forward` function.

# Example:
```julia-repl
julia> struct Book
    title::String
    author::String
end
julia> Base.show(io::IO, b::Book) = println(io, "\$(b.title) (by \$(b.author))")
julia> Base.print(b::Book) = println("In a hole in the ground there lived a hobbit...")
julia> author(b::Book) = b.author

julia> struct PaperBook
    b::Book
    number_of_pages::Int
end
julia> @forward((PaperBook, :b), Book)
julia> pages(book::PaperBook) = book.number_of_pages

julia> struct Edition
    b::PaperBook
    year::Int
end
julia> @forward((Edition, :b), PaperBook)
julia> year(book::Edition) = book.year

julia> book = Edition(PaperBook(Book("The Hobbit", "J.R.R. Tolkien"), 374), 2013)

julia> print(author(book), ", ", pages(book), " pages, Ed. ", year(book))
J.R.R. Tolkien, 374 pages, Ed. 2013

julia> print(book)
In a hole in the ground there lived a hobbit...
```
"""
macro forward(sender, receiver, ekws...)
    kws = Vector{Pair{Symbol,Any}}()
    methods = nothing
    for kw in ekws
        if isa(kw, Expr)  &&  (kw.head == :(=))
            push!(kws, Pair(kw.args[1], kw.args[2]))
        else
            @assert methods == nothing "Too many arguments"
            methods = kw
        end
    end
    out = quote
        counterr = 0
        if $methods == nothing
            list = forward($sender, $receiver; $kws...)
        else
            list = forward($sender, $receiver, $methods; $kws...)
        end
        for line in list
            try
                eval(Meta.parse("$line"))
            catch err
                global counterr += 1
                println()
                println("$line")
                @error err;
            end
        end
        println(length(list) - counterr, " method(s) forwarded")
        if counterr > 0
            println(counterr, " method(s) raised an error")
        end
    end
    return esc(out)
end

"""
`@copy_fields T`

Copy all field definitions from a structure into another.

# Example
```julia-repl

julia> struct First
    field1
    field2::Int
end
julia> struct Second
    @copy_fields(First)
    field3::String
end

julia> println(fieldnames(Second))
(:field1, :filed2, :field3)

julia> println(fieldtype.(Second, fieldnames(Second)))
(Any, Int64, String)
```
"""
macro copy_fields(T)
    out = Expr(:block)
    for name in fieldnames(__module__.eval(T))
        e = Expr(Symbol("::"))
        push!(e.args, name)
        push!(e.args, fieldtype(__module__.eval(T), name))
        push!(out.args, e)
    end
    return esc(out)
end

"""
`concretetype(T::Type)`

Return the concrete type associated with the *quasi abstract* type `T`. If `T` is not *quasi abstract* returns nothing.

See also: `@quasiabstract`
"""
concretetype(T::Type) = nothing

"""
`isquasiabstract(T::Type)`

Return `true` if `T` is *quasi abstract*, `false` otherwise.

See also: `@quasiabstract`
"""
isquasiabstract(T::Type) = false

"""
`isquasiconcrete(T::Type)`

Return `true` if `T` is a concrete type associated to a *quasi abstract* type, `false` otherwise.

See also: `@quasiabstract`
"""
isquasiconcrete(T::Type) = false

#  """
#  `concretesubtypes(T::Type)`
#
#  Return a vector whose elements are the concrete types associated with the *quasi abstract* subtypes of `T`.  If there are no such types returns nothing.
#
#  See also: `@quasiabstract`
#  """
#  concretesubtypes(T::Type) = something([concretetype(a) for a in subtypes(T)]...)
#
#
#  """
#  `quasiabstractsubtypes(T::Type)`
#
#  Return a vector whose elements are the *quasi abstract* subtypes of `T`.  If there are no such types returns an empty vector.
#
#  See also: `@quasiabstract`
#  """
#  quasiabstractsubtypes(T::Type) = filter(isquasiabstract, subtypes(T))

"""
`@quasiabstract expression`

Create a *quasi abstract* type.

This macro accepts an expression defining a (mutable or immutable) structure, and outputs the code for two new type definitions:
- an abstract type with the same name and (if given) supertype of the input structure;
- a concrete structure definition with name prefix `Concrete_`, subtyping the abstract type defined above.

The relation between the types ensure there will be a single concrete type associated to the abstract one, hence you can use the abstract type to annotate method arguments, and be sure to receive the associated concrete type, or one of its subtypes (which shares all field names and types of the ancestors).

The concrete type associated to an *quasi abstract* type can be retrieved with the `concretetype` function.  The `Concrete_` prefix can be customized passing a second symbol argument to the macro.

These newly defined types allows to easily implement concrete subtyping: simply use the *quasi abstract* type name for both object construction and to annotate method arguments.

# Example:
```julia-repl
julia> @quasiabstract struct Book
    title::String
    author::String
end
julia> Base.show(io::IO, b::Book) = println(io, "\$(b.title) (by \$(b.author))")
julia> Base.print(b::Book) = println("In a hole in the ground there lived a hobbit...")
julia> author(b::Book) = b.author

julia> @quasiabstract struct PaperBook <: Book
    number_of_pages::Int
end
julia> pages(book::PaperBook) = book.number_of_pages

julia> @quasiabstract struct Edition <: PaperBook
    year::Int
end
julia> year(book::Edition) = book.year

julia> println(fieldnames(concretetype(Edition)))
(:title, :author, :number_of_pages, :year)

julia> book = Edition("The Hobbit", "J.R.R. Tolkien", 374, 2013)

julia> print(author(book), ", ", pages(book), " pages, Ed. ", year(book))
J.R.R. Tolkien, 374 pages, Ed. 2013

julia> print(book)
In a hole in the ground there lived a hobbit...

julia> isquasiconcrete(typeof(book))
true

julia> isquasiabstract(supertype(typeof(book)))
true

julia> concretetype(typeof(book)) === supertype(typeof(book))
true
```
"""
macro quasiabstract(expr, prefix=:Concrete_)
    function change_symbol!(expr, from::Symbol, to::Symbol)
        if isa(expr, Expr)
            for i in 1:length(expr.args)
                if isa(expr.args[i], Symbol)  &&  (expr.args[i] == from)
                    expr.args[1] = to
                else
                    if isa(expr.args[i], Expr)
                        change_symbol!(expr.args[i], from, to)
                    end
                end
            end
        end
    end
    drop_parameter_types(expr::Symbol) = expr
    function drop_parameter_types(_aaa::Expr)
        aaa = deepcopy(_aaa)
        for ii in 1:length(aaa.args)
            if isa(aaa.args[ii], Expr)  &&  (aaa.args[ii].head == :<:)
                insert!(aaa.args, ii, aaa.args[ii].args[1])
                deleteat!(aaa.args, ii+1)
            end
        end
        return aaa
    end
    function prepare_where(expr::Expr, orig::Expr)
        @assert expr.head == :curly
        whereclause = Expr(:where)
        push!(whereclause.args, orig)
        for i in 2:length(expr.args)
            @assert (isa(expr.args[i], Symbol)  ||  (isa(expr.args[i], Expr)  &&  (expr.args[i].head == :<:)))
            push!(whereclause.args, expr.args[i])
        end
        return whereclause
    end

    @assert isa(expr, Expr)
    @assert expr.head == :struct "Expression must be a `struct`"

    # Ensure there is room for the super type in the expression
    if isa(expr.args[2], Symbol)
        name = deepcopy(expr.args[2])
        deleteat!(expr.args, 2)
        insert!(expr.args, 2, Expr(:<:, name, :Any))
    end
    if isa(expr.args[2], Expr)  &&  (expr.args[2].head != :<:)
        insert!(expr.args, 2, Expr(:<:, expr.args[2], :Any))
        deleteat!(expr.args, 3)
    end
    @assert isa(expr.args[2], Expr)
    @assert expr.args[2].head == :<:

    # Get name and super type
    tmp = expr.args[2].args[1]
    @assert (isa(tmp, Symbol)  ||  (isa(tmp, Expr)  &&  (tmp.head == :curly)))
    name = deepcopy(tmp)
    name_symb = (isa(tmp, Symbol)  ?  tmp  :  tmp.args[1])

    tmp = expr.args[2].args[2]
    @assert (isa(tmp, Symbol)  ||  (isa(tmp, Expr)  &&  (tmp.head == :curly)))
    super = deepcopy(tmp)
    super_symb = (isa(tmp, Symbol)  ?  tmp  :  tmp.args[1])

    # Output abstract type
    out = Expr(:block)
    push!(out.args, :(abstract type $name <: $super end))

    # Change name in input expression
    concrete_symb = Symbol(prefix, name_symb)
    if isa(name, Symbol)
        concrete = concrete_symb
    else
        concrete = deepcopy(name)
        change_symbol!(concrete, name_symb, concrete_symb)
    end
    change_symbol!(expr, name_symb, concrete_symb)

    # Drop all types from parameters in `name`, or they'll raise syntax errors
    name = drop_parameter_types(name)

    # Change super type in input expression to actual name
    deleteat!(expr.args[2].args, 2)
    insert!(  expr.args[2].args, 2, name)

    # If an ancestor type is quasi abstract retrieve the associated
    # concrete type and add its fields as members
    p = __module__.eval(super_symb)
    while (p != Any)
        if isquasiabstract(p)
            parent = concretetype(p)
            for i in fieldcount(parent):-1:1
                e = Expr(Symbol("::"))
                push!(e.args, fieldname(parent, i))
                push!(e.args, fieldtype(parent, i))
                pushfirst!(expr.args[3].args, e)
            end
            break
        end
        p = __module__.eval(supertype(p))
    end

    # Add modified input structure to output
    push!(out.args, expr)

    # Add a constructor whose name is the same as the abstract type
    if isa(name, Expr)
        whereclause = prepare_where(name, :($name(args...; kw...)))
        push!(out.args, :($whereclause = $concrete(args...; kw...)))
    else
        push!(out.args, :($name(args...; kw...) = $concrete(args...; kw...)))
    end

    # Add `concretetype`, `isquasiabstract` and `isquasiconcrete` methods
    push!(out.args, :(ReusePatterns.concretetype(::Type{$name_symb}) = $concrete_symb))
    push!(out.args, :(ReusePatterns.isquasiabstract(::Type{$name_symb}) = true))
    push!(out.args, :(ReusePatterns.isquasiconcrete(::Type{$concrete_symb}) = true))

    if isa(name, Expr)
        whereclause = prepare_where(name, :(ReusePatterns.concretetype(::Type{$name})))
        push!(out.args, :($whereclause = $concrete))

        whereclause = prepare_where(name, :(ReusePatterns.isquasiabstract(::Type{$name})))
        push!(out.args, :($whereclause = true))

        # Drop all types from parameters in `concrete`, or they'll raise syntax errors
        concrete = drop_parameter_types(concrete)

        whereclause = prepare_where(name, :(ReusePatterns.isquasiconcrete(::Type{$concrete})))
        push!(out.args, :($whereclause = true))
    end
    return esc(out)
end

end # module


text:
module ShaleDrillingLikelihood_ThreadUtils_Test

DOBTIME = false

using ShaleDrillingLikelihood

using Test
using StatsFuns
using Random
using BenchmarkTools
using Base.Threads

using Base.Threads: nthreads, @threads, Atomic, atomic_add!

using ShaleDrillingLikelihood: getrange,
    Mapper,
    Mapper2,
    AbstractThreadMapper,
    getrange,
    nextrange!,
    stop,
    step,
    add_iter!

function h2(x)
    n = length(x)
    mapper = Mapper2(n)
    ld = stop(mapper)
    batch_size = step(mapper)

    s = Atomic{Float64}(0.0)
    @threads for j in 1:nthreads()
        while true
            k = add_iter!(mapper)
            batch_start = 1 + (k-1) * batch_size
            batch_end = min(k * batch_size, ld)
            batch_start > ld && break
            y = 0.0
            @inbounds @simd for i in batch_start:batch_end
                y += exp(x[i])
            end
            atomic_add!(s, y)
        end
    end
    return s[]
end

function h3(x)
    mapper = Mapper(length(x))
    s = Atomic{Float64}(0.0)

    @threads for j in 1:nthreads()
        while true
            y = 0.0
            irng = nextrange!(mapper)
            isnothing(irng) && break
            @inbounds @simd for i in irng
                y += exp(x[i])
            end
            atomic_add!(s, y)
        end
    end

    return s[]
end

function f(x)
    s = Atomic{Float64}(0.0)
    n = length(x)
    @threads for k in 1:nthreads()
        y = 0.0
        @inbounds @simd for i in getrange(n)
            y += exp(x[i])
        end
        atomic_add!(s, y)
    end
    s[]
end

function g(x)
    n = length(x)
    s = 0.0
    @inbounds @simd for i in 1:n
        s += exp(x[i])
    end
    s
end

gg(x) = mapreduce(exp, +, x)

function ggg(x)
    n = length(x)
    s = 0.0
    for i in 1:n
        s += exp(x[i])
    end
    s
end

z = rand(10^3)
@testset "summation: getrange" begin
    println("--------")

    println("Mapper")
    _h3 = h3(z)
    DOBTIME && @btime h3($z)

    println("Mapper2")
    _h2 = h2(z)
    DOBTIME && @btime h2($z)

    println("getrange")
    _f = f(z)
    DOBTIME && @btime f($z)

    println("unthreaded")
    DOBTIME && @btime g($z)
    _g = g(z)

    println("mapreduce")
    DOBTIME && @btime gg($z)

    println("vanilla")
    DOBTIME && @btime ggg($z)

    @test _f ≈ _g
    @test _g ≈ _h2
    @test _g ≈ _h3
end

end


text:
import Base: +, -
#==
Power Log
==#
"""
      A powerlog can be used to perform a one-tailed test, as well as get the proper sample size for a testing population.\n
      --------------------\n
      ==PARAMETERS==\n
     p1 <- A Float64 percentage representing the probability of scenario one.\n
     p2 <- A Float64 percentage representing the probability of scenario two. These two probability values should follow these guidelines: p1 = p1 + x = p2\n
     alpha = 0.05 <- Sets an alpha value\n
     --------------------\n
     Returns power, sample_size
       """
function PowerLog(p1::Float64,p2::Float64; alpha::Float64 = 0.05, rsq::Real = 0)
    pd = p2 - p1
    l1 = p1/(1-p1)
    l2 = p2/(1-p2)
    θ = l2 / l1
    or = θ
    λ = log(θ)
    λ2 = λ^2
    za = quantile(normal(),1-alpha)
    println("One-tailed test: alpha = ",alpha,", p1 = ",p1,", p2 = ",p2,", rsq = ",rsq,", odds ratio = ",or)
    δ = (1 + (1 + λ2)*exp(5 * λ2/4))/(1 + exp(-1*λ2/4))
    pwr = zeros(Float64,8)
    nn = zeros(Int64,8)
    i = 1
    for power = 0.6:.05:.95
        zb = quantile(normal(),power)

        N = ((za + zb*exp(-1 * λ2/4))^2 * (1 + 2*p1*δ))/(p1*λ2)
        N /= (1 - rsq)
        pwr[i] = power
        nn[i] = ceil(Int64,N)
        i += 1
    end
    return(pwr, nn)
end
#==
Mean
    Baseline
==#
 # Model Type
 @doc """
       A mean baseline is great for getting a basic accuracy score in order
           to make a valid direction for your model.\n
         --------------------\n
         ==PARAMETERS==\n
        [y] <- Fill with your trainY values. Should be an array of shape (0,1) or (1,0)\n
        pipl = Pipeline([StandardScalar(),LinearRegression(trainX,trainy)])\n
        --------------------\n
        ==Functions==\n
        predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)
                     """
function MeanBaseline(y)
    m = mean(m.y)
    predict(xt) =
    xt = [v = m for v in xt]
    (test)->(m;predict)
end
#==
Majority
    Class
        Baseline
==#
@doc """
      Majority class baseline is used to find the most often interpreted
      classification in an array.\n
      --------------------\n
      ==PARAMETERS==\n
     [y] <- Fill with your trainY values. Should be an array of shape (0,1) or (1,0)\n
     --------------------\n
     ==Functions==\n
     predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)\n
     counts() <- Returns a dictionary with the counts of all inserted keys.\n
     highest() <- Will return a Dictionary key with the count as well as the value for the most interpreted classification.
       """
function ClassBaseline(y)
    u=unique(y)
    d=Dict([(i,count(x->x==i,y)) for i in u])
    d = sort(collect(d), by=x->x[2])
    maxkey = d[length(d)]
    predict(xt) = [p = maxkey[1] for p in xt]
    counts() = d
    highest() = maxkey
    (var)->(y;maxkey;d;predict;counts;highest)
end
      """
          ## Pipeline
          ### Description
            Rescales an array. Pipelines can contain a predictable Lathe model
            with preprocessing that occurs automatically. This is done by
            putting X array processing methods into the iterable steps, and
            then putting your Lathe model in.\n
            --------------------\n
          ### Input
            --------------------\n
            #### Positional Arguments
            LatheObject :: steps - An infinte argument of LatheObject types. These types
             are any Lathe model or preprocessor.\n
            --------------------\n
           ### Output
           Pipeline :: A Pipeline object.
           ---------------------\n
           ### Functions
           Pipeline.predict(xt) :: Applies the steps inside of the pipeline
           to xt.\n
           Pipeline.show() :: shows all of the steps inside of the pipeline
           in order, along with their respective count.\n
           ---------------------\n
           ### Data
           steps - An array of LatheObject types that are predicted with usng
           the predict() function call.\n
           ---------------------\n
           ### Methods
           Base.+ - The + operator can be used to add steps to a pipeline.\n
           Pipeline + LatheObject\n
           Base.- - The - operator can be used to remove steps from a pipeline.\n
           Pipeline - Int64(Position in steps)
             """
      mutable struct Pipeline{P, S} <: Tool
          steps::Array{LatheObject}
          predict::P
          show::S
          function Pipeline(steps::LatheObject ...)
              steps = [step for step in steps]
              show() = _show(steps)
              predict(xt::Array) = pipe_predict(xt, steps)
              predict(xt::DataFrame) = pipe_predict(xt, steps)
              new{typeof(predict), typeof(show)}(steps, predict, show)
          end
          function _show(steps::Array{LatheObject})
              count = 0
              for step in steps
                  count += 1
                  typetitle = string(typeof(step))
                  name = split(typetitle, '{')[1]
                  println(count, " = = > ", name)
              end
          end
          function pipe_predict(xt, steps)
              for step in steps
                  xt = step.predict(xt)
              end
              return(xt)
          end
end
-(p::Pipeline, n::Int64) = deleteat!(p.steps, n)
+(p::Pipeline, step::LatheObject) = push!(p.steps, step)
+(m1::LatheObject, m2::LatheObject) = Pipeline(m1, m2)
mutable struct Router{P} <: LatheObject
    fn::Function
    components::Array{LatheObject}
    predict::P
    function Router(components::LatheObject ...; fn::Function)
        components = Array([comp for comp in components])
        predict(xt) = router_load(xt, fn, components)
        new{typeof(predict)}(fn, components, predict)
    end
    function router_load(data, fn, components)
        returns = fn(data)
        preds = []
        count = 0
        for cp in components
            count += 1
            res = cp.predict(returns[count])
            push!(preds, res)
        end
        return([res for res in preds])
    end
end
function _compare_predCat(models, xbar)
    count = 0
    preddict = Dict()
    for model in models
        preddict[count] = model.predict(xbar)
        count += 1
    end
    count = 1
    n_features = length(preddict)
    encoder = OrdinalEncoder(preddict[1])
    y_hat = encoder.predict(preddict[1])
    for (key, value) in preddict
       encoded = encoder.predict(value)
        y_hat[count] = mean([y_hat[count], encoded[count]])
        count += 1
    end
    y_hat = Array{Int64}(y_hat)
    inv_lookup = Dict(value => key for (key, value) in encoder.lookup)
    for x in y_hat
        println(inv_lookup[x])
    end
end


text:
module test_julia_range

using Test

r = range(1, step=2, stop=10)

@test r == 1:2:9
@test first(r) == 1
@test step(r) == 2
@test last(r) == 9

end # module test_julia_range


text:
include("NMF.jl")

using NMF

# 通常版を使用する場合

# 分解対象の[2,4]の行列を用意
V = [1 2 3 4; 2 3 4 5]

# 因子分回数 K を５で設定する
k = 3

# NMF処理を実行
# @time W,H = nmf_euc(V,k)
# @time W,H = nmf_kl(V,k)
@time W,H = nmf_is(V,k)

# 結果を確認
println("Normal NMF Result : $(norm(V - W*H))")
println("-------------------------------------------")

# テンプレート版を使用する場合

# 分解対象の[2,4]の行列を用意
V = [1 2 3 4; 2 3 4 5]

# テンプレート行列を用意
W = [0. 1. 11.;2. 6. 10.]

# NMF処理を実行
# @time W,H = nmf_euc(V,W)
# @time W,H = nmf_kl(V,W)
@time H = nmf_is(V,W)

# 結果を確認
println("Template NMF Result : $(norm(V - W*H))")


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

name = "JuliaThriftBuilder"
version = v"0.2.0"

# Collection of sources required to build JuliaThriftBuilder
sources = [
    "https://github.com/tanmaykm/thrift.git" =>
    "ef6de66707eb6135402a73520deec3478d9e1ec7",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd thrift/
./bootstrap.sh
if [ $target != "x86_64-apple-darwin14" ] && [ $target != "x86_64-unknown-freebsd11.1" ]; then
    LDFLAGS="-static-libgcc -static-libstdc++"
    export LDFLAGS
fi
./configure --prefix=$prefix --host=$target --enable-tutorial=no --enable-tests=no --enable-libs=no --disable-werror
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "thrift", :thrift)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/tanmaykm/BisonBuilder/releases/download/v3.0.5/build_BisonBuilder.v1.0.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: tanmaykm/JuliaThriftBuilder

max_issues_repo_name: tanmaykm/JuliaThriftBuilder

max_forks_repo_name: tanmaykm/JuliaThriftBuilder


text:
"""
    abstract type AbstractEDM{T} <: AbstractMatrix{T} end
"""
abstract type AbstractEDM{T} <: AbstractMatrix{T} end

"""
    EDM{T} <: AbstractEDM{T} <: AbstractMatrix{T}

Euclidean Distance Matrix (EDM).

# Examples
Generate a random EDM:

    julia> D = EDM(2, 3)
    3×3 EDM{Float64}:
     0.0        0.0588086  0.361541
     0.0588086  0.0        0.635088
     0.361541   0.635088   0.0

---

Generate an EDM from a Matrix or TripletEmbeddings.Embedding:

    julia> X = Embedding(2,3)
    2×3 Embedding{Float64}:
     0.584104  0.372789  0.389338
     0.54629   0.547285  0.609064

    julia> D = EDM(X)
    3×3 EDM{Float64}:
     0.0        0.0446553   0.0418744
     0.0446553  0.0         0.00409058
     0.0418744  0.00409058  0.0
"""
mutable struct EDM{T} <: AbstractEDM{T}
    D::Matrix{Float64} # Euclidean distances matrix

    function EDM(d::Int, n::Int)
        d ≥ 1 || throw(DomainError("d < 1"))
        n ≥ 1 || throw(DomainError("n < 1"))
        @assert d ≤ n "We need d ≤ n (dimension ≤ # of points)"

        X = rand(d, n)

        # This is equivalent to D = diag(G) * ones(1,n) - 2G + ones(n) * diag(G)',
        # where G is the Gram matrix of (centered) X.
        D = pairwise(SqEuclidean(), X, dims=2)

        new{Float64}(D)
    end

    function EDM(D::Matrix{T}) where T <: Real
        isedm(D) || throw(DomainError(D, "D is not an EDM."))
        new{T}(D)
    end

    function EDM(X::TripletEmbeddings.Embedding{T}) where T
        D = pairwise(SqEuclidean(), X, dims=2)
        new{T}(D)
    end

end

mutable struct NoisyEDM{T} <: AbstractEDM{T}
    D::EDM{T} # Data
    N::Matrix{Float64}            # Noise

    function NoisyEDM(
        D::EDM{T},
        N::Matrix{T}) where T <: Float64

        size(N) == size(D) || throw(DimensionMismatch("D and N must be the same size."))

        new{T}(D, N)
    end

    function NoisyEDM(
        X::TripletEmbeddings.Embedding{T},
        N::Matrix{T}) where T <: Float64

        D = EDM(X)
        size(N) == size(D) || throw(DimensionMismatch("D and N must be the same size."))

        new{T}(D, N)
    end
end

mutable struct MaskedEDM{T} <: AbstractEDM{T}
    D::EDM{Float64} # Euclidean distances matrix
    mask::AbstractMatrix{<:Real}

    function MaskedEDM(
        X::TripletEmbeddings.Embedding{T},
        p::T;
        symmetric=false) where T <: Real

        @assert 0 < p ≤ 1 "p must be in the interval (0,1]"
        D = EDM(X)
        new{T}(D, randommask(D, Bernoulli(p); symmetric=symmetric))
    end

    function MaskedEDM(
        X::TripletEmbeddings.Embedding{T},
        mask::AbstractMatrix{S}) where {T <: Real, S <: Real}

        D = pairwise(SqEuclidean(), X, dims=2)
        size(D) == size(mask) || throw(DimensionMismatch("D and mask must be the same size"))
        new{T}(D, mask)
    end

    function MaskedEDM(
        D::EDM{T},
        p::T) where T <: Real

        @assert 0 < p ≤ 1 "p must be in the interval (0,1]"
        new{T}(D.D, rand(Bernoulli(p), size(D)))
    end

    function MaskedEDM(
        D::EDM{T},
        distribution::S) where {T <: Real, S <: Distribution{Univariate,Discrete}}

        new{T}(D.D, rand(distribution, size(D)))
    end

    function MaskedEDM(
        D::EDM{T},
        mask::AbstractMatrix{S}) where {T <: Real, S <: Real}

        size(D) == size(mask) || throw(DimensionMismatch("D and mask must be the same size"))
        new{T}(D.D, mask)
    end
end

mutable struct NoisyMaskedEDM{T} <: AbstractEDM{T}
    D::MaskedEDM{T}
    N::Matrix{<:Real}

    function NoisyMaskedEDM(
        X::TripletEmbeddings.Embedding{T},
        N::Matrix{P},
        mask::AbstractMatrix{S},) where {T <: Real, S <: Real, P <: Real}

        new{T}(MaskedEDM(X, mask), N)

    end

    function NoisyMaskedEDM(
        D::EDM{T},
        N::Matrix{P},
        mask::AbstractMatrix{S}) where {T <: Real, S <: Real, P <: Real}

        size(mask) == size(D) || throw(DimensionMismatch("D and mask must be the same size."))
        size(mask) == size(N) || throw(DimensionMismatch("N and mask must be the same size."))

        new{T}(MaskedEDM(D, mask), N)
    end

    function NoisyMaskedEDM(
        D::NoisyEDM{T},
        mask::AbstractMatrix{S}) where {T <: Real, S <: Real}

        size(mask) == size(D) || throw(DimensionMismatch("D and mask must be the same size."))
        size(mask) == size(D.N) || throw(DimensionMismatch("N and mask must be the same size."))

        new{T}(MaskedEDM(D, mask), D.N)
    end
end

masked(D::MaskedEDM) = D .* D.mask

Base.size(D::AbstractEDM) = size(D.D)
Base.getindex(D::AbstractEDM, inds...) = getindex(D.D, inds...)

nitems(D::AbstractEDM) = size(D.D, 1)

function Base.show(io::IO, ::MIME"text/plain", D::MaskedEDM)
    show(io, "text/plain", D.D .* D.mask)
end

function Base.:*(s::Number, D::AbstractEDM)
    D.D = s * D.D
    return D
end

LinearAlgebra.rank(D::AbstractEDM) = rank(D.D)

"""
    isedm(D[, atol=sqrt(eps())])

Check whether a matrix is an EDM. atol is passed to isapprox's atol.
"""
function isedm(D::Matrix{T}; atol=sqrt(eps())) where T <: Real
    n = try
        LinearAlgebra.checksquare(D)
    catch DimensionMismatch
        return false
    end

    n > 1 || throw(DomainError(D, "D must have dimensions n > 1"))

    λs = eigvals(-J(n) * D * J(n) / 2)

    # isposdef doesn't work in this case for eigvals that are slightly negative (i.e. -1e-16).
    return all((isapprox.(real(λs), zeros(n), atol=atol)) .| (real(λs) .≥ 0.0)) .& # Real part
           all(isapprox.(imag(λs), zeros(n), atol=atol))                           # Imag part
end

function J(n::Int)
    return Matrix{Float64}(I, n, n) - ones(n,n)/n
end


text:
using SpatialEcology
using Random
using Test

@testset "Grid" begin
    gr = SpatialEcology.creategrid([repeat(20.5 .+ (1:5), inner = 4) repeat(-30.5 .+ (2:2:8), outer = 5)])
    @test xmin(gr) == 21.5
    @test ymin(gr) == -28.5
    @test cellsize(gr) == (1.0, 2.0)
    @test cells(gr) == (5, 4)
    @test xrange(gr) == 21.5:1.0:25.5
    @test yrange(gr) == -28.5:2.0:-22.5
    @test boundingbox(gr).xmin == xmin(gr)
    @test sprint(show, boundingbox(gr)) == "xmin:\t21.5\nxmax:\t25.5\nymin:\t-28.5\nymax:\t-22.5\n"
end


text:
using StabTrussTopOpt
sto = StabTrussTopOpt

function TO_define(par::sto.TOProblemPar; E::Number=1.0, τ::Number=1.0, σ_c::Number = 10.0, σ_t::Number = 10.0)::sto.TOProblem

    ndim = par.ndim
    n = par.n
    m = par.m
    B = par.B
    edof_map = par.dof_map

    ndof = ndim * n
    # free dof
    mask_free = [v[1] for v in findall(i->i==0, reshape(par.Bound', ndof, 1))]
    nfree = length(mask_free)

    ff = reshape(par.F', ndof, 1)[mask_free]

    l = zeros(Float64, m)
    for i=1:m
        u_id = par.T[i,1]
        v_id = par.T[i,2]
        l[i] = norm(par.X[v_id,:] - par.X[u_id,:])
    end

    # linear objective coeff
    c = l

    g_I = Int[]
    g_J = Int[]
    g_V = Float64[]
    for i=1:nfree
        append!(g_I, i * ones(ndim))
        append!(g_J, edof_map[i, 1:ndim])
        append!(g_V, B[i, 1:ndim])

        append!(g_I, i * ones(ndim))
        append!(g_J, edof_map[i, ndim+1:2*ndim])
        append!(g_V, B[i, ndim+1:2*ndim])
    end
    Γ_full = sparse(g_I, g_J, g_V, m, ndof)
    Γ = Γ_full[:, mask_free] # ∈ m x nfree

    eK = Array{SparseMatrixCSC{Float64, Int}, 1}(undef, m)
    for i=1:m
        eK[i] = E / l[i] * (Γ[i, :] * Γ[i, :]')
    end

    # construct orthonormal basis for gamma_i, i ∈ [m]
    # aka. δ_i (and η_i in 3D)
    Bn = zeros(Float64, m, ndof)
    for i=1:m
        # TODO: general in 3D: Gram-Schmidt in the null space of gamma_i
        # here 2D specific construction to make orthonormal
        Bn[i,1] = B[i,2]
        Bn[i,2] = -B[i,1]
        Bn[i,3] = B[i,4]
        Bn[i,4] = -B[i,3]
    end

    # map Bn's local index 1..2xndim to global dof indices
    spBn_I = Int[]
    spBn_J = Int[]
    spBn_V = Float64[]
    for i=1:m
        append!(spBn_I, i * ones(ndim))
        append!(spBn_J, edof_map[i, 1:ndim])
        append!(spBn_V, Bn[i, 1:ndim])

        append!(spBn_I, i * ones(ndim))
        append!(spBn_J, edof_map[i, ndim+1:2*ndim])
        append!(spBn_V, Bn[i, ndim+1:2*ndim])
    end
    spBn = sparse(spBn_I, spBn_J, spBn_V, m, ndof)

    # construct Δ_i = (δ_i * δ_i^T + η_i * η_i^T), i ∈ [m]
    eDelta = Array{SparseMatrixCSC{Float64, Int}, 1}(undef, m)
    for i=1:m
       eDelta[i] = spBn[i, mask_free] * spBn[i, mask_free]'
    end

    a_init = 10 * ones(Float64, m)
    # TODO: solve the full comp-geometric stiffness matrix to get init_q
    q_init = 10 * ones(Float64, m)
    x_init  = vcat(a_init, q_init)

    return sto.TOProblem(m, c, τ, E, σ_c, σ_t, l, ff, Γ, eDelta, eK, x_init)
end


text:
using DelimitedFiles
import Base.size

# Bingo struct
mutable struct Bingo
  n ::Int
  nums :: Matrix{Int}
  mask :: Matrix{Bool}
end

Bingo(n::Int) = Bingo(n,zeros(n,n),zeros(n,n))
Bingo(n::Int, nums::Matrix{Int}) = Bingo(n,nums,zeros(n,n))
Bingo(nums::Matrix{Int}) = Bingo(size(nums)[1], nums)

size(b::Bingo) :: Int = b.n
row(b::Bingo, i::Int) :: Vector{Bool} = b.mask[i,:]
col(b::Bingo,i::Int) :: Vector{Bool} = b.mask[:,i]

# Returns true if the bingo is complete.
function bingo(b::Bingo) :: Bool
  res = false
  for i in 1:size(b)
    res = res || all(row(b,i)) || all(col(b,i))
  end
  return res
end

# Adds number appearance to bingo.
function newNum!(num::Int, b::Bingo)
  ind = findall(x -> x == num, b.nums)
  if length(ind) > 0
    b.mask[ind] .= true
  end
end

# Returns the score of a bingo board.
function score(b::Bingo)
  return sum(b.nums[map(x->!x,b.mask)])
end

# Returns the score of the first bingo to be complete.
function playBingo!(nums::Vector{Int}, bingos::Vector{Bingo}) :: Int
  for num in nums
    println("New num: ", num)
    map(x -> newNum!(num, x), bingos)
    complete = map(x -> bingo(x), bingos)
    if (any(complete))
      ind = findfirst(complete)
      println("BINGO!!!")
      println("Bingo sheet nº", ind, " has been completed: ")
      display(bingos[ind].mask)
      return score(bingos[ind]) * num
    end
  end
  return -1
end

# Returns the score of the last bingo to be complete.
function lastBingo!(nums::Vector{Int}, bingos::Vector{Bingo}) :: Int
  lastRound :: Vector{Bool} = zeros(length(bingos))
  for num in nums
    println("New num: ", num)
    map(x -> newNum!(num, x), bingos)
    complete = map(x -> bingo(x), bingos)
    if (all(complete))
      ind = findfirst(map(x -> !x, lastRound))
      println("BINGO!!!")
      println("Bingo sheet nº", ind, " has been completed: ")
      display(bingos[ind].mask)
      return score(bingos[ind]) * num
    end
    lastRound = deepcopy(complete)
  end
  return -1
end

function readInput()
  # Read numbers
  f = readdlm("./inputs/day4.txt", ',')
  nums :: Vector{Int} = f[1,:]

  # Read bingo sheets
  f = readdlm("./inputs/day4.txt", Int, skipstart=1, skipblanks=true)
  sBingos = size(f)[2]
  nBingos = size(f)[1] ÷ sBingos
  f = reshape(transpose(f),(sBingos,sBingos,nBingos))

  bingos :: Vector{Bingo} = []
  for i in 1:nBingos
    push!(bingos, Bingo(permutedims(f[:,:,i],(2,1))))
  end
  return nums, bingos
end

nums, bingos = readInput()

res =  playBingo!(nums,deepcopy(bingos))
println("Final score: ", res)

res = lastBingo!(nums,deepcopy(bingos))
println("Final score: ", res)


text:
using NumericalContinuation
using NumericalContinuation.FourierCollocation
using Test

function duffing!(res, u, p, t)
    res[1] = u[2]
    res[2] = p[2]*sin(p[1]*t) - 2*p[3]*u[2] - u[1] - u[1]^3
    return
end

@testset "FourierCollocation" begin
    @testset "Duffing" begin
        p0 = [1.0, 0.1, 0.05]
        # sol = solve(ODEProblem(duffing!, [1.0, 0.0], (0.0, 50*2π/p0[1]), p0), Tsit5())
        # sol = solve(ODEProblem(duffing!, sol[:, end], (0.0, 2π/p0[1]), p0), Tsit5())
        # u0 = sol(range(0, 2π/p0[1], length=21)[1:end-1])[:, :]
        u0 = [-0.233 -0.092 0.057 0.202 0.329 0.426 0.481 0.487 0.443 0.354 0.233 0.092 -0.057 -0.202 -0.329 -0.426 -0.481 -0.487 -0.443 -0.354 -0.233; 0.425 0.469 0.473 0.439 0.364 0.248 0.099 -0.063 -0.217 -0.34 -0.425 -0.468 -0.473 -0.439 -0.364 -0.248 -0.1 0.063 0.216 0.34 0.424]
        fourier = fourier_collocation("duffing", duffing!, [0, 2π/p0[1]], u0, p0, pnames=["ω", "Γ", "ξ"], phase=false, fix_t0=true, fix_t1=false)

    p1 = zero_problem("circle", (u, p) -> u[1]^2 + u[2]^2 - 1, [1, 0], [])
    p2 = zero_problem("plane", (u, p) -> u[1] + u[2] + u[3], [1, 0, -1], [])

    p1 = zero_problem("cubic", (u, p) -> u^3 - p, 1.5, 1.5^3)
    p2 = zero_problem("quad", (u, p) -> u^2 - p, 1.2, 1.2^2)
    p3 = push!(push!(Problem("group"), p1), p2)
    p4 = push!(Problem("main"), p3)
    mfunc = monitor_function("test", u -> u[1], )
    mfuncs = monitor_functions()
    flat = NumericalContinuation.flatten(p4)
end

using Pkg
Pkg.activate("nonlinear", shared=true)
using NumericalContinuation
using NumericalContinuation.FourierCollocation

function duffing!(res, u, p, t)
    res[1] = u[2]
    res[2] = p[2]*sin(p[1]*t) - 2*p[3]*u[2] - u[1] - u[1]^3
    return
end

p0 = [1.0, 0.1, 0.05]
u0 = [-0.233 -0.092 0.057 0.202 0.329 0.426 0.481 0.487 0.443 0.354 0.233 0.092 -0.057 -0.202 -0.329 -0.426 -0.481 -0.487 -0.443 -0.354 -0.233; 0.425 0.469 0.473 0.439 0.364 0.248 0.099 -0.063 -0.217 -0.34 -0.425 -0.468 -0.473 -0.439 -0.364 -0.248 -0.1 0.063 0.216 0.34 0.424]
fourier = fourier_collocation("duffing", duffing!, [0, 2π/p0[1]], u0, p0, pnames=["ω", "Γ", "ξ"], phase=false, fix_t0=true, fix_t1=false)
opt = Options()
closed = NumericalContinuation.ClosedProblem(fourier, opt)
NumericalContinuation.init!(closed)


text:
immutable JNINativeInterface #struct JNINativeInterface_ {

    reserved0::Ptr{Void} # void *reserved0;

    reserved1::Ptr{Void} #void *reserved1;
    reserved2::Ptr{Void} #void *reserved2;
    reserved3::Ptr{Void} #void *reserved3;

    GetVersion::Ptr{Void} #jint ( *GetVersion)(JNIEnv *env);
    DefineClass::Ptr{Void} #jclass ( *DefineClass) (JNIEnv *env, const char *name, jobject loader, const jbyte *buf, jsize len);
    FindClass::Ptr{Void} #jclass ( *FindClass) (JNIEnv *env, const char *name);
    FromReflectedMethod::Ptr{Void} #jmethodID ( *FromReflectedMethod) (JNIEnv *env, jobject method);
    FromReflectedField::Ptr{Void} #jfieldID ( *FromReflectedField)(JNIEnv *env, jobject field);
    ToReflectedMethod::Ptr{Void} #jobject ( *ToReflectedMethod) (JNIEnv *env, jclass cls, jmethodID methodID, jboolean isStatic);

    GetSuperClass::Ptr{Void}  #jclass ( *GetSuperclass) (JNIEnv *env, jclass sub);
    IsAssignableFrom::Ptr{Void} #jboolean ( *IsAssignableFrom) (JNIEnv *env, jclass sub, jclass sup);

    ToReflectedField::Ptr{Void} #jobject ( *ToReflectedField)(JNIEnv *env, jclass cls, jfieldID fieldID, jboolean isStatic);

    Throw::Ptr{Void} #jint ( *Throw) (JNIEnv *env, jthrowable obj);
    ThrowNew::Ptr{Void} #jint ( *ThrowNew)(JNIEnv *env, jclass clazz, const char *msg);
    ExceptionOccurred::Ptr{Void} #jthrowable ( *ExceptionOccurred) (JNIEnv *env);
    ExceptionDescribe::Ptr{Void} #void ( *ExceptionDescribe)(JNIEnv *env);
    ExceptionClear::Ptr{Void} #void ( *ExceptionClear) (JNIEnv *env);
    FatalError::Ptr{Void} #void ( *FatalError) (JNIEnv *env, const char *msg);

    PushLocalFrame::Ptr{Void} #jint ( *PushLocalFrame) (JNIEnv *env, jint capacity);
    PopLocalFrame::Ptr{Void} #jobject ( *PopLocalFrame) (JNIEnv *env, jobject result);

    NewGlobalRef::Ptr{Void} #jobject ( *NewGlobalRef) (JNIEnv *env, jobject lobj);
    DeleteGlobalRef::Ptr{Void} #void ( *DeleteGlobalRef) (JNIEnv *env, jobject gref);
    DeleteLocalRef::Ptr{Void} #void ( *DeleteLocalRef) (JNIEnv *env, jobject obj);
    IsSameObject::Ptr{Void} #jboolean ( *IsSameObject) (JNIEnv *env, jobject obj1, jobject obj2);
    NewLocalRef::Ptr{Void} #jobject ( *NewLocalRef) (JNIEnv *env, jobject ref);
    EnsureLocalCapacity::Ptr{Void} #jint ( *EnsureLocalCapacity) (JNIEnv *env, jint capacity);

    AllocObject::Ptr{Void} #jobject ( *AllocObject) (JNIEnv *env, jclass clazz);
    NewObject::Ptr{Void} #jobject ( *NewObject) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    NewObjectV::Ptr{Void} #jobject ( *NewObjectV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    NewObjectA::Ptr{Void} #jobject ( *NewObjectA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    GetObjectClass::Ptr{Void} #jclass ( *GetObjectClass) (JNIEnv *env, jobject obj);
    IsInstanceOf::Ptr{Void} #jboolean ( *IsInstanceOf) (JNIEnv *env, jobject obj, jclass clazz);

    GetMethodID::Ptr{Void} #jmethodID ( *GetMethodID) (JNIEnv *env, jclass clazz, const char *name, const char *sig);

    CallObjectMethod::Ptr{Void} #jobject ( *CallObjectMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallObjectMethodV::Ptr{Void} #jobject ( *CallObjectMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallObjectMethodA::Ptr{Void} #jobject ( *CallObjectMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue * args);

    CallBooleanMethod::Ptr{Void} #jboolean ( *CallBooleanMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallBooleanMethodV::Ptr{Void} #jboolean ( *CallBooleanMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallBooleanMethodA::Ptr{Void} #jboolean ( *CallBooleanMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue * args);

    CallByteMethod::Ptr{Void} #jbyte ( *CallByteMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallByteMethodV::Ptr{Void} #jbyte ( *CallByteMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallByteMethodA::Ptr{Void} #jbyte ( *CallByteMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue *args);

    CallCharMethod::Ptr{Void} #jchar ( *CallCharMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallCharMethodV::Ptr{Void} #jchar ( *CallCharMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallCharMethodA::Ptr{Void} #jchar ( *CallCharMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue *args);

    CallShortMethod::Ptr{Void} #jshort ( *CallShortMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallShortMethodV::Ptr{Void} #jshort ( *CallShortMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallShortMethodA::Ptr{Void} #jshort ( *CallShortMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue *args);

    CallIntMethod::Ptr{Void} #jint ( *CallIntMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallIntMethodV::Ptr{Void} #jint ( *CallIntMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallIntMethodA::Ptr{Void} #jint ( *CallIntMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue *args);

    CallLongMethod::Ptr{Void} #jlong ( *CallLongMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallLongMethodV::Ptr{Void} #jlong ( *CallLongMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallLongMethodA::Ptr{Void} #jlong ( *CallLongMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue *args);

    CallFloatMethod::Ptr{Void} #jfloat ( *CallFloatMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallFloatMEthodV::Ptr{Void} #jfloat ( *CallFloatMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallFloatMethodA::Ptr{Void} #jfloat ( *CallFloatMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue *args);

    CallDoubleMethod::Ptr{Void} #jdouble ( *CallDoubleMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallDoubleMethodV::Ptr{Void} #jdouble ( *CallDoubleMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallDoubleMethodA::Ptr{Void} #jdouble ( *CallDoubleMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue *args);

    CallVoidMethod::Ptr{Void} #void ( *CallVoidMethod) (JNIEnv *env, jobject obj, jmethodID methodID, ...);
    CallVoidMethodV::Ptr{Void} #void ( *CallVoidMethodV) (JNIEnv *env, jobject obj, jmethodID methodID, va_list args);
    CallVoidMethodA::Ptr{Void} #void ( *CallVoidMethodA) (JNIEnv *env, jobject obj, jmethodID methodID, const jvalue * args);

    CallNonvirtualObjectMethod::Ptr{Void} #jobject ( *CallNonvirtualObjectMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonvirtualObjectMethodV::Ptr{Void} #jobject ( *CallNonvirtualObjectMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualObjectMEthodA::Ptr{Void} #jobject ( *CallNonvirtualObjectMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, const jvalue * args);

    CallNonvirtualBooleanMethod::Ptr{Void} #jboolean ( *CallNonvirtualBooleanMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonvirtualBooleanMethodV::Ptr{Void} #jboolean ( *CallNonvirtualBooleanMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualBooleanMethodA::Ptr{Void} #jboolean ( *CallNonvirtualBooleanMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, const jvalue * args);

    CallNonvirtualByteMethod::Ptr{Void} #jbyte ( *CallNonvirtualByteMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonvirtualByteMethodV::Ptr{Void} #jbyte ( *CallNonvirtualByteMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualByteMethodA::Ptr{Void} #jbyte ( *CallNonvirtualByteMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, const jvalue *args);

    CallNonvirtualCharMethod::Ptr{Void} #jchar ( *CallNonvirtualCharMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonvirtualCharMethodV::Ptr{Void} #jchar ( *CallNonvirtualCharMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualCharMethodA::Ptr{Void} #jchar ( *CallNonvirtualCharMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, const jvalue *args);

    CallNonvirtualShortMethod::Ptr{Void} #jshort ( *CallNonvirtualShortMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonvirtualShortMethodV::Ptr{Void} #jshort ( *CallNonvirtualShortMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualShortMethodA::Ptr{Void} #jshort ( *CallNonvirtualShortMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, const jvalue *args);

    CallNonvirtualIntMethod::Ptr{Void} #jint ( *CallNonvirtualIntMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonvirtualIntMethodV::Ptr{Void} #jint ( *CallNonvirtualIntMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualIntMethodA::Ptr{Void} #jint ( *CallNonvirtualIntMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, const jvalue *args);

    CallNonvirtualLongMethod::Ptr{Void} #jlong ( *CallNonvirtualLongMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonVirtualLongMethodV::Ptr{Void} #jlong ( *CallNonvirtualLongMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualLongMethodA::Ptr{Void} #jlong ( *CallNonvirtualLongMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, const jvalue *args);

    CallNonVirtualFloatMethod::Ptr{Void} #jfloat ( *CallNonvirtualFloatMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonvirtualFloatMethodV::Ptr{Void} #jfloat ( *CallNonvirtualFloatMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualFloatMethodA::Ptr{Void} #jfloat ( *CallNonvirtualFloatMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, const jvalue *args);

    CallNonvirtualDoubleMethod::Ptr{Void} # jdouble ( *CallNonvirtualDoubleMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonvirtualDoubleMethodV::Ptr{Void} # jdouble ( *CallNonvirtualDoubleMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualDoubleMethodA::Ptr{Void} # jdouble ( *CallNonvirtualDoubleMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID,  const jvalue *args);

    CallNonvirtualVoidMethod::Ptr{Void} # void ( *CallNonvirtualVoidMethod) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, ...);
    CallNonvirtualVoidMethodV::Ptr{Void} # void ( *CallNonvirtualVoidMethodV) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, va_list args);
    CallNonvirtualVoidMethodA::Ptr{Void} # void ( *CallNonvirtualVoidMethodA) (JNIEnv *env, jobject obj, jclass clazz, jmethodID methodID, const jvalue * args);

    GetFieldID::Ptr{Void} # jfieldID ( *GetFieldID) (JNIEnv *env, jclass clazz, const char *name, const char *sig);

    GetObjectField::Ptr{Void} # jobject ( *GetObjectField) (JNIEnv *env, jobject obj, jfieldID fieldID);
    GetBooleanField::Ptr{Void} # jboolean ( *GetBooleanField) (JNIEnv *env, jobject obj, jfieldID fieldID);
    GetByteField::Ptr{Void} # jbyte ( *GetByteField) (JNIEnv *env, jobject obj, jfieldID fieldID);
    GetCharField::Ptr{Void} # jchar ( *GetCharField) (JNIEnv *env, jobject obj, jfieldID fieldID);
    GetShortField::Ptr{Void} # jshort ( *GetShortField) (JNIEnv *env, jobject obj, jfieldID fieldID);
    GetIntField::Ptr{Void} # jint ( *GetIntField)  (JNIEnv *env, jobject obj, jfieldID fieldID);
    GetLongField::Ptr{Void} # jlong ( *GetLongField) (JNIEnv *env, jobject obj, jfieldID fieldID);
    GetFloatField::Ptr{Void} # jfloat ( *GetFloatField) (JNIEnv *env, jobject obj, jfieldID fieldID);
    GetDoubleField::Ptr{Void} # jdouble ( *GetDoubleField) (JNIEnv *env, jobject obj, jfieldID fieldID);

    SetObjectField::Ptr{Void} # void ( *SetObjectField) (JNIEnv *env, jobject obj, jfieldID fieldID, jobject val);
    SetBooleanField::Ptr{Void} # void ( *SetBooleanField) (JNIEnv *env, jobject obj, jfieldID fieldID, jboolean val);
    SetByteField::Ptr{Void} # void ( *SetByteField)  (JNIEnv *env, jobject obj, jfieldID fieldID, jbyte val);
    SetCharField::Ptr{Void} # void ( *SetCharField) (JNIEnv *env, jobject obj, jfieldID fieldID, jchar val);
    SetShortField::Ptr{Void} # void ( *SetShortField) (JNIEnv *env, jobject obj, jfieldID fieldID, jshort val);
    SetIntField::Ptr{Void} # void ( *SetIntField) (JNIEnv *env, jobject obj, jfieldID fieldID, jint val);
    SetLongField::Ptr{Void} # void ( *SetLongField) (JNIEnv *env, jobject obj, jfieldID fieldID, jlong val);
    SetFloatField::Ptr{Void} # void ( *SetFloatField) (JNIEnv *env, jobject obj, jfieldID fieldID, jfloat val);
    SetDoubleField::Ptr{Void} # void ( *SetDoubleField) (JNIEnv *env, jobject obj, jfieldID fieldID, jdouble val);

    GetStaticMethodID::Ptr{Void} # jmethodID ( *GetStaticMethodID) (JNIEnv *env, jclass clazz, const char *name, const char *sig);

    CallStaticObjectMethod::Ptr{Void} # jobject ( *CallStaticObjectMethod) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    CallStaticObjectMethodV::Ptr{Void} #jobject ( *CallStaticObjectMethodV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    CallStaticObjectMethodA::Ptr{Void} # jobject ( *CallStaticObjectMethodA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    CallStaticBooleanMethod::Ptr{Void} # jboolean ( *CallStaticBooleanMethod) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    CallStaticBooleanMethodV::Ptr{Void} # jboolean ( *CallStaticBooleanMethodV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    CallStaticBooleanMethodA::Ptr{Void} # jboolean ( *CallStaticBooleanMethodA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    CallStaticByteMethod::Ptr{Void} #jbyte ( *CallStaticByteMethod) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    CallStaticByteMethodV::Ptr{Void} # jbyte ( *CallStaticByteMethodV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    CallStaticByteMethodA::Ptr{Void} #jbyte ( *CallStaticByteMethodA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    CallStaticCharMethod::Ptr{Void} # jchar ( *CallStaticCharMethod) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    CallStaticCharMethodV::Ptr{Void} # jchar ( *CallStaticCharMethodV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    CallStaticCharMethodA::Ptr{Void} # jchar ( *CallStaticCharMethodA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    CallStaticShortMethod::Ptr{Void} # jshort ( *CallStaticShortMethod) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    CallStaticShortMethodV::Ptr{Void} # jshort ( *CallStaticShortMethodV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    CallStaticShortMethodA::Ptr{Void} # jshort ( *CallStaticShortMethodA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    CallStaticIntMethod::Ptr{Void} #jint ( *CallStaticIntMethod) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    CallStaticIntMethodV::Ptr{Void} # jint ( *CallStaticIntMethodV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    CallStaticIntMethodA::Ptr{Void} # jint ( *CallStaticIntMethodA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    CallStaticLongMethod::Ptr{Void} # jlong ( *CallStaticLongMethod) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    CallStaticLongMethodV::Ptr{Void} # jlong ( *CallStaticLongMethodV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    CallStaticLongMethodA::Ptr{Void} # jlong ( *CallStaticLongMethodA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    CallStaticFloatMethod::Ptr{Void} # jfloat ( *CallStaticFloatMethod) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    CallStaticFloatMethodV::Ptr{Void} # jfloat ( *CallStaticFloatMethodV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    CallStaticFloatMethodA::Ptr{Void} # jfloat ( *CallStaticFloatMethodA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    CallStaticDoubleMethod::Ptr{Void} # jdouble ( *CallStaticDoubleMethod) (JNIEnv *env, jclass clazz, jmethodID methodID, ...);
    CallStaticDoubleMethodV::Ptr{Void} # jdouble ( *CallStaticDoubleMethodV) (JNIEnv *env, jclass clazz, jmethodID methodID, va_list args);
    CallStaticDoubleMethodA::Ptr{Void} # jdouble ( *CallStaticDoubleMethodA) (JNIEnv *env, jclass clazz, jmethodID methodID, const jvalue *args);

    CallStaticVoidMethod::Ptr{Void} # void ( *CallStaticVoidMethod) (JNIEnv *env, jclass cls, jmethodID methodID, ...);
    CallStaticVoidMethodV::Ptr{Void} # void ( *CallStaticVoidMethodV) (JNIEnv *env, jclass cls, jmethodID methodID, va_list args);
    CallStaticVoidMethodA::Ptr{Void} # void ( *CallStaticVoidMethodA) (JNIEnv *env, jclass cls, jmethodID methodID, const jvalue * args);

    GetStaticFieldID::Ptr{Void} # jfieldID ( *GetStaticFieldID) (JNIEnv *env, jclass clazz, const char *name, const char *sig);
    GetStaticObjectField::Ptr{Void} # jobject ( *GetStaticObjectField) (JNIEnv *env, jclass clazz, jfieldID fieldID);
    GetStaticBooleanField::Ptr{Void} # jboolean ( *GetStaticBooleanField) (JNIEnv *env, jclass clazz, jfieldID fieldID);
    GetStaticByteField::Ptr{Void} # jbyte ( *GetStaticByteField) (JNIEnv *env, jclass clazz, jfieldID fieldID);
    GetStaticCharField::Ptr{Void} # jchar ( *GetStaticCharField) (JNIEnv *env, jclass clazz, jfieldID fieldID);
    GetStaticShortField::Ptr{Void} # jshort ( *GetStaticShortField) (JNIEnv *env, jclass clazz, jfieldID fieldID);
    GetStaticIntField::Ptr{Void} # jint ( *GetStaticIntField) (JNIEnv *env, jclass clazz, jfieldID fieldID);
    GetStaticLongField::Ptr{Void} # jlong ( *GetStaticLongField) (JNIEnv *env, jclass clazz, jfieldID fieldID);
    GetStaticFloatField::Ptr{Void} # jfloat ( *GetStaticFloatField) (JNIEnv *env, jclass clazz, jfieldID fieldID);
    GetStaticDoubleField::Ptr{Void} # jdouble ( *GetStaticDoubleField) (JNIEnv *env, jclass clazz, jfieldID fieldID);

    SetStaticObjectField::Ptr{Void} # void ( *SetStaticObjectField) (JNIEnv *env, jclass clazz, jfieldID fieldID, jobject value);
    SetStaticBooleanField::Ptr{Void} # void ( *SetStaticBooleanField) (JNIEnv *env, jclass clazz, jfieldID fieldID, jboolean value);
    SetStaticByteField::Ptr{Void} # void ( *SetStaticByteField) (JNIEnv *env, jclass clazz, jfieldID fieldID, jbyte value);
    SetStaticCharField::Ptr{Void} # void ( *SetStaticCharField) (JNIEnv *env, jclass clazz, jfieldID fieldID, jchar value);
    SetStaticShortField::Ptr{Void} # void ( *SetStaticShortField) (JNIEnv *env, jclass clazz, jfieldID fieldID, jshort value);
    SetStaticIntField::Ptr{Void} # void ( *SetStaticIntField) (JNIEnv *env, jclass clazz, jfieldID fieldID, jint value);
    SetStaticLongField::Ptr{Void} # void ( *SetStaticLongField) (JNIEnv *env, jclass clazz, jfieldID fieldID, jlong value);
    SetStaticFloatField::Ptr{Void} # void ( *SetStaticFloatField) (JNIEnv *env, jclass clazz, jfieldID fieldID, jfloat value);
    SetStaticDoubleField::Ptr{Void} # void ( *SetStaticDoubleField) (JNIEnv *env, jclass clazz, jfieldID fieldID, jdouble value);

    NewString::Ptr{Void} #::Ptr{Void} # jstring ( *NewString) (JNIEnv *env, const jchar *unicode, jsize len);
    GetStringLength::Ptr{Void} # jsize ( *GetStringLength) (JNIEnv *env, jstring str);
    GetStringChars::Ptr{Void} # const jchar *( *GetStringChars) (JNIEnv *env, jstring str, jboolean *isCopy);
    ReleaseStringChars::Ptr{Void} # void ( *ReleaseStringChars) (JNIEnv *env, jstring str, const jchar *chars);

    NewStringUTF::Ptr{Void} # jstring ( *NewStringUTF) (JNIEnv *env, const char *utf);
    GetStringUTFLength::Ptr{Void} # jsize ( *GetStringUTFLength) (JNIEnv *env, jstring str);
    GetStringUTFChars::Ptr{Void} # const char* ( *GetStringUTFChars) (JNIEnv *env, jstring str, jboolean *isCopy);
    ReleaseStringUTFChars::Ptr{Void} # void ( *ReleaseStringUTFChars) (JNIEnv *env, jstring str, const char* chars);

    GetArrayLength::Ptr{Void} # jsize ( *GetArrayLength ) (JNIEnv *env, jarray array);

    NewObjectArray::Ptr{Void} # jobjectArray ( *NewObjectArray) (JNIEnv *env, jsize len, jclass clazz, jobject init);
    GetObjectArrayElement::Ptr{Void} # jobject ( *GetObjectArrayElement) (JNIEnv *env, jobjectArray array, jsize index);
    SetObjectArrayElement::Ptr{Void} # void ( *SetObjectArrayElement) (JNIEnv *env, jobjectArray array, jsize index, jobject val);

    NewBooleanArray::Ptr{Void} # jbooleanArray ( *NewBooleanArray) (JNIEnv *env, jsize len);
    NewByteArray::Ptr{Void} # jbyteArray ( *NewByteArray) (JNIEnv *env, jsize len);
    NewCharArray::Ptr{Void} # jcharArray ( *NewCharArray) (JNIEnv *env, jsize len);
    NewShortArray::Ptr{Void} # jshortArray ( *NewShortArray) (JNIEnv *env, jsize len);
    NewIntArray::Ptr{Void} # jintArray ( *NewIntArray) (JNIEnv *env, jsize len);
    NewLongArray::Ptr{Void} # jlongArray ( *NewLongArray) (JNIEnv *env, jsize len);
    NewFloatArray::Ptr{Void} # jfloatArray ( *NewFloatArray) (JNIEnv *env, jsize len);
    NewDoubleArray::Ptr{Void} # jdoubleArray ( *NewDoubleArray) (JNIEnv *env, jsize len);

    GetBooleanArrayElements::Ptr{Void} # jboolean * ( *GetBooleanArrayElements) (JNIEnv *env, jbooleanArray array, jboolean *isCopy);
    GetByteArrayElements::Ptr{Void} # jbyte * ( *GetByteArrayElements) (JNIEnv *env, jbyteArray array, jboolean *isCopy);
    GetCharArrayElements::Ptr{Void} # jchar * ( *GetCharArrayElements) (JNIEnv *env, jcharArray array, jboolean *isCopy);
    GetShortArrayElements::Ptr{Void} # jshort * ( *GetShortArrayElements) (JNIEnv *env, jshortArray array, jboolean *isCopy);
    GetIntArrayElements::Ptr{Void} # jint * ( *GetIntArrayElements) (JNIEnv *env, jintArray array, jboolean *isCopy);
    GetLongArrayElements::Ptr{Void} # jlong * ( *GetLongArrayElements ) (JNIEnv *env, jlongArray array, jboolean *isCopy);
    GetFloatArrayElements::Ptr{Void} # jfloat * ( *GetFloatArrayElements) (JNIEnv *env, jfloatArray array, jboolean *isCopy);
    GetDoubleArrayElements::Ptr{Void} # jdouble * ( *GetDoubleArrayElements) (JNIEnv *env, jdoubleArray array, jboolean *isCopy);

    ReleaseBooleanArrayElements::Ptr{Void} # void ( *ReleaseBooleanArrayElements ) (JNIEnv *env, jbooleanArray array, jboolean *elems, jint mode);
    ReleaseByteArrayElements::Ptr{Void} # void ( *ReleaseByteArrayElements ) (JNIEnv *env, jbyteArray array, jbyte *elems, jint mode);
    ReleaseCharArrayElements::Ptr{Void} # void ( *ReleaseCharArrayElements ) (JNIEnv *env, jcharArray array, jchar *elems, jint mode);
    ReleaseShortArrayElements::Ptr{Void} # void ( *ReleaseShortArrayElements) (JNIEnv *env, jshortArray array, jshort *elems, jint mode);
    ReleaseIntArrayElements::Ptr{Void} # void ( *ReleaseIntArrayElements) (JNIEnv *env, jintArray array, jint *elems, jint mode);
    ReleaseLongArrayElements::Ptr{Void} # void ( *ReleaseLongArrayElements ) (JNIEnv *env, jlongArray array, jlong *elems, jint mode);
    ReleaseFloatArrayElements::Ptr{Void} # void ( *ReleaseFloatArrayElements) (JNIEnv *env, jfloatArray array, jfloat *elems, jint mode);
    ReleaseDoubleArrayElements::Ptr{Void} # void ( *ReleaseDoubleArrayElements ) (JNIEnv *env, jdoubleArray array, jdouble *elems, jint mode);

    GetBooleanArrayRegion::Ptr{Void} # void ( *GetBooleanArrayRegion ) (JNIEnv *env, jbooleanArray array, jsize start, jsize l, jboolean *buf);
    GetByteArrayRegion::Ptr{Void} # void ( *GetByteArrayRegion) (JNIEnv *env, jbyteArray array, jsize start, jsize len, jbyte *buf);
    vGetCharArrayRegion::Ptr{Void} # void ( *GetCharArrayRegion) (JNIEnv *env, jcharArray array, jsize start, jsize len, jchar *buf);
    GetShortArrayRegion::Ptr{Void} # void ( *GetShortArrayRegion ) (JNIEnv *env, jshortArray array, jsize start, jsize len, jshort *buf);
    GetIntArrayRegion::Ptr{Void} # void ( *GetIntArrayRegion ) (JNIEnv *env, jintArray array, jsize start, jsize len, jint *buf);
    GetLongArrayRegion::Ptr{Void} # void ( *GetLongArrayRegion) (JNIEnv *env, jlongArray array, jsize start, jsize len, jlong *buf);
    GetFloatArrayRegion::Ptr{Void} # void ( *GetFloatArrayRegion) (JNIEnv *env, jfloatArray array, jsize start, jsize len, jfloat *buf);
    GetDoubleArrayRegion::Ptr{Void} # void ( *GetDoubleArrayRegion) (JNIEnv *env, jdoubleArray array, jsize start, jsize len, jdouble *buf);

    SetBooleanArrayRegion::Ptr{Void} # void ( *SetBooleanArrayRegion )(JNIEnv *env, jbooleanArray array, jsize start, jsize l, const jboolean *buf);
    SetByteArrayRegion::Ptr{Void} # void ( *SetByteArrayRegion) (JNIEnv *env, jbyteArray array, jsize start, jsize len, const jbyte *buf);
    SetCharArrayRegion::Ptr{Void} # void ( *SetCharArrayRegion) (JNIEnv *env, jcharArray array, jsize start, jsize len, const jchar *buf);
    SetShortArrayRegion::Ptr{Void} # void ( *SetShortArrayRegion ) (JNIEnv *env, jshortArray array, jsize start, jsize len, const jshort *buf);
    SetIntArrayRegion::Ptr{Void} # void ( *SetIntArrayRegion) (JNIEnv *env, jintArray array, jsize start, jsize len, const jint *buf);
    SetLongArrayRegion::Ptr{Void} # void ( *SetLongArrayRegion) (JNIEnv *env, jlongArray array, jsize start, jsize len, const jlong *buf);
    SetFloatArrayRegion::Ptr{Void} # void ( *SetFloatArrayRegion ) (JNIEnv *env, jfloatArray array, jsize start, jsize len, const jfloat *buf);
    SetDoubleArrayRegion::Ptr{Void} # void ( *SetDoubleArrayRegion) (JNIEnv *env, jdoubleArray array, jsize start, jsize len, const jdouble *buf);

    RegisterNatives::Ptr{Void} # jint ( *RegisterNatives) (JNIEnv *env, jclass clazz, const JNINativeMethod *methods, jint nMethods);
    UnregisterNatives::Ptr{Void} # jint ( *UnregisterNatives ) (JNIEnv *env, jclass clazz);

    MonitorEnter::Ptr{Void} # jint ( *MonitorEnter) (JNIEnv *env, jobject obj);
    MonitorExit::Ptr{Void} # jint ( *MonitorExit ) (JNIEnv *env, jobject obj);

    GetJavaVM::Ptr{Void} # jint ( *GetJavaVM) (JNIEnv *env, JavaVM **vm);

    GetStringRegion::Ptr{Void} # void ( *GetStringRegion ) (JNIEnv *env, jstring str, jsize start, jsize len, jchar *buf);
    GetStringUTFRegion::Ptr{Void} # void ( *GetStringUTFRegion) (JNIEnv *env, jstring str, jsize start, jsize len, char *buf);

    GetPrimitiveArrayCritical::Ptr{Void} # void * ( *GetPrimitiveArrayCritical) (JNIEnv *env, jarray array, jboolean *isCopy);
    ReleasePrimitiveArrayCritical::Ptr{Void} # void ( *ReleasePrimitiveArrayCritical ) (JNIEnv *env, jarray array, void *carray, jint mode);

    GetStringCritical::Ptr{Void} # const jchar * ( *GetStringCritical ) (JNIEnv *env, jstring string, jboolean *isCopy);
    ReleaseStringCritical::Ptr{Void} # void ( *ReleaseStringCritical ) (JNIEnv *env, jstring string, const jchar *cstring);

    NewWeakGlobalRef::Ptr{Void} # jweak ( *NewWeakGlobalRef) (JNIEnv *env, jobject obj);
    DeleteWeakGlobalRef::Ptr{Void} # void ( *DeleteWeakGlobalRef) (JNIEnv *env, jweak ref);

    ExceptionCheck::Ptr{Void} # jboolean ( *ExceptionCheck) (JNIEnv *env);

    NewDirectByteBuffer::Ptr{Void} # jobject ( *NewDirectByteBuffer ) (JNIEnv* env, void* address, jlong capacity);
    GetDirectBufferAddress::Ptr{Void} # void* ( *GetDirectBufferAddress ) (JNIEnv* env, jobject buf);
    GetDirectBufferCapacity::Ptr{Void} # jlong ( *GetDirectBufferCapacity) (JNIEnv* env, jobject buf);

    #/* New JNI 1.6 Features */

    GetObjectRefType::Ptr{Void} # jobjectRefType ( *GetObjectRefType) (JNIEnv* env, jobject obj);
end #};

immutable JNIEnv
    JNINativeInterface_::Ptr{JNINativeInterface}
end

immutable JNIInvokeInterface #struct JNIInvokeInterface_ {
    reserved0::Ptr{Void} #void *reserved0;
    reserved1::Ptr{Void} #vvoid *reserved1;
    reserved2::Ptr{Void} #vvoid *reserved2;

    DestroyJavaVM::Ptr{Void} #jint (JNICALL *DestroyJavaVM)(JavaVM *vm);

    AttachCurrentThread::Ptr{Void} #jint (JNICALL *AttachCurrentThread)(JavaVM *vm, void **penv, void *args);

    DetachCurrentThread::Ptr{Void} #jint (JNICALL *DetachCurrentThread)(JavaVM *vm);

    GetEnv::Ptr{Void} #jint (JNICALL *GetEnv)(JavaVM *vm, void **penv, jint version);

    AttachCurrentThreadAsDaemon::Ptr{Void} #jint (JNICALL *AttachCurrentThreadAsDaemon)(JavaVM *vm, void **penv, void *args);
end

immutable JavaVM
    JNIInvokeInterface_::Ptr{JNIInvokeInterface}
end


text:
using DynamicPolynomials
# include("D:\\Programs\\TSSOS\\src\\TSSOS.jl")
using TSSOS

# unconstrained optimization using the TSSOS hierarchy
@polyvar x[1:3]
f = 1+x[1]^4+x[2]^4+x[3]^4+x[1]*x[2]*x[3]+x[2]
opt,sol,data = tssos_first(f, x, newton=true, reducebasis=true, TS="block", solution=true, QUIET=true)
# optimum = 0.475274778453039

opt,sol,data = tssos_higher!(data, TS="MD", solution=true, QUIET=true)
# optimum = 0.4752747815682939

# constrained optimization using the TSSOS hierarchy
@polyvar x[1:3]
f = 1+x[1]^4+x[2]^4+x[3]^4+x[1]*x[2]*x[3]+x[2]
g = 1-x[1]^2-2*x[2]^2
h = x[2]^2+x[3]^2-1
pop = [f, g, h]
opt,sol,data = tssos_first(pop, x, 2, numeq=1, TS="block", solution=true, QUIET=true)
# optimum = 0.6861925732538557

opt,sol,data = tssos_higher!(data, TS="MD", solution=true, QUIET=true)
# optimum = 0.6861925611033903

# constrained optimization using the CS-TSSOS hierarchy
n = 6
@polyvar x[1:n]
f = 1+sum(x.^4)+x[1]*x[2]*x[3]+x[3]*x[4]*x[5]+x[3]*x[4]*x[6]+x[3]*x[5]*x[6]+x[4]*x[5]*x[6]
pop = [f, 1-sum(x[1:3].^2), 1-sum(x[3:6].^2)]
d = 2
opt,sol,data = cs_tssos_first(pop, x, d, numeq=1, TS="block", solution=true, QUIET=true)
# optimum = 0.7174253329322346

opt,sol,data = cs_tssos_higher!(data, TS="MD", solution=true, QUIET=true)
# optimum = 0.7174253409780793

println("Run successfully!")


text:
using DataFrames

mutable struct ISLogger
    logs::Dict{Any, DataFrame}
    nlogs::Int
end

function ISLogger(is::InventorySystem)
    logs = Dict{String, DataFrame}()
    for item in is.bom
        logs[item.name] = DataFrame()
    end
    for cons in is.constraints
        logs[cons.name] = DataFrame()
    end
    ISLogger(logs, 0)
end

function (logger::ISLogger)(is::InventorySystem; log_id = 0)
    logger.nlogs +=1
    for item in is.bom
        df = get_logs(item)
        df.simulation_id = fill(logger.nlogs,nrow(df))
        df.log_id = fill(log_id, nrow(df))
        append!(logger.logs[item.name], df)
    end
    for cons in is.constraints
        df = get_logs(cons)
        df.simulation_id = fill(logger.nlogs, nrow(df))
        df.log_id = fill(log_id, nrow(df))
        append!(logger.logs[cons.name], df)
    end
end

Base.getindex(logger::ISLogger, key) = logger.logs[key]

function RLBase.reset!(logger::ISLogger)
    logger.nlogs = 0
    for k in keys(logger.logs)
        logger.logs[k] = DataFrame()
    end
end

function get_logs(item::Item)
    df = hcat(get_logs(item.inventory), map(get_logs, item.sources)...)
    df.iteration = 1:nrow(df)
    return df
end

function get_logs(ep::EndProduct)
    df = hcat(get_logs(ep.market), get_logs(ep.inventory), map(get_logs, ep.sources)...)
    df.iteration = 1:nrow(df)
    return df
end

function get_logs(ass::Assembly)
    hcat(DataFrame(assembly_batchsize = ass.batchsize_log, assembly_cost = ass.cost_log), get_logs(ass.leadtime))
end

function get_logs(inv::Inventory)
    DataFrame(inventory_onhand = inv.on_hand_log, inventory_stockout = inv.stockout_log, inventory_fillrate = inv.fillrate_log, inventory_cost = inv.cost_log)
end

function get_logs(lt::LeadTime)
    DataFrame(leadtime_cost = lt.cost_log)
end

function get_logs(ma::Market)
    DataFrame(market_backorder = ma.backorder_log, market_fillrate = ma.fillrate_log, market_cost = ma.cost_log)
end

function get_logs(sup::Supplier)
    hcat(DataFrame(supplier_batchsize = sup.batchsize_log, supplier_cost = sup.cost_log), get_logs(sup.leadtime))
end

function get_logs(c::RessourceConstraint)
    df = DataFrame("utilization" => c.utilization_log, "setups" => c.setup_log)
    df.iteration = 1:nrow(df)
    return df
end