text:
function polrtest(nm::OrdinalMultinomialModel, Z::AbstractVecOrMat; test=:score)
    if test == :score
        polrtest(OrdinalMultinomialScoreTest(nm, reshape(Z, size(Z, 1), size(Z, 2))))
    elseif test == :LRT
        polrtest(OrdinalMultinomialLrtTest(nm, reshape(Z, size(Z, 1), size(Z, 2))))
    else
        throw(ArgumentError("unrecognized test $test"))
    end
end

polrtest(nm::StatsModels.TableRegressionModel{<:OrdinalMultinomialModel}, Z::AbstractVecOrMat; kwargs...) =
polrtest(nm.model, Z; kwargs...)

###########################
# Score test
###########################

struct OrdinalMultinomialScoreTest{TY<:Integer, T<:BlasReal, TL<:GLM.Link}
    "`q`: number of covariates to test significance"
    q::Int
    "`Z`: n-by-q covariate matrix to test significance"
    Z::Matrix{T}
    "`nm`: fitted null model"
    nm::OrdinalMultinomialModel{TY,T,TL}
    "`scoreγ`: score vector of γ"
    scoreγ::Vector{T}
    "`∂γ∂θβ`: information of γ vs (θ,β)"
    ∂γ∂θβ::Matrix{T}
    "`∂γ∂γ`: information of γ vs γ"
    ∂γ∂γ::Matrix{T}
    "`scratchv1`: working vector, same size as `scoreγ`"
    scratchv1::Vector{T}
    "`scratchm1`: working matrix, same size as `Z`"
    scratchm1::Matrix{T}
    "`scratchm2`: working matrix, same size as `∂γ∂θβ`"
    scratchm2::Matrix{T}
    "`scratchm3`: working matrix, same size as `∂γ∂γ`"
    scratchm3::Matrix{T}
end

function OrdinalMultinomialScoreTest(nm::OrdinalMultinomialModel, Z::Matrix)
    TY, T, TL = eltype(nm.Y), eltype(nm.X), typeof(nm.link)
    q = size(Z, 2)
    scoreγ = zeros(T, q)
    ∂γ∂θβ  = zeros(T, q, nm.npar)
    ∂γ∂γ   = zeros(T, q, q)
    scratchv1 = similar(scoreγ)
    scratchm1 = similar(Z)
    scratchm2 = similar(∂γ∂θβ)
    scratchm3 = similar(∂γ∂γ)
    OrdinalMultinomialScoreTest{TY, T, TL}(q, Z, nm, scoreγ, ∂γ∂θβ, ∂γ∂γ, scratchv1, scratchm1, scratchm2, scratchm3)
end

function polrtest(d::OrdinalMultinomialScoreTest)
    #
    # partition the informatrix matrix as
    # [P  W]
    # [W' Q],
    # test statistics is R^T (Q - W' P^{-1} W)^{-1} * R, where R is the score
    #
    # compute R = d.scoreγ, score of γ
    if isempty(d.nm.wts)
        mul!(d.scoreγ, transpose(d.Z), d.nm.wt∂β)
    else
        d.nm.wtwk .= d.nm.wts .* d.nm.wt∂β
        mul!(d.scoreγ, transpose(d.Z), d.nm.wtwk)
    end
    # compute W' = d.∂γ∂θβ
    if isempty(d.nm.wts)
        @views mul!(d.∂γ∂θβ[:, 1:d.nm.J-1], transpose(d.Z), d.nm.wt∂θ∂β) # wt∂θ∂β is n-by-(J-1)
    else
        copyto!(d.scratchm1, d.Z)
        lmul!(Diagonal(d.nm.wts), d.scratchm1)
        @views mul!(d.∂γ∂θβ[:, 1:d.nm.J-1], transpose(d.scratchm1), d.nm.wt∂θ∂β)
    end
    if isempty(d.nm.wts)
        copyto!(d.nm.wtwk, d.nm.wt∂β∂β)
    else
        d.nm.wtwk .= d.nm.wts .* d.nm.wt∂β∂β
    end
    copyto!(d.nm.scratchm1, d.nm.X)
    lmul!(Diagonal(d.nm.wtwk), d.nm.scratchm1)
    @views mul!(d.∂γ∂θβ[:, d.nm.J:end], transpose(d.Z), d.nm.scratchm1)
    # compute Q = ∂γ∂γ
    copyto!(d.scratchm1, d.Z)
    lmul!(Diagonal(d.nm.wtwk), d.scratchm1)
    mul!(d.∂γ∂γ, transpose(d.Z), d.scratchm1)
    # compute scratchm3 = Q - W' P^{-1} W
    mul!(d.scratchm2, d.∂γ∂θβ, d.nm.vcov)
    mul!(d.scratchm3, d.scratchm2, transpose(d.∂γ∂θβ))
    d.scratchm3 .= d.∂γ∂γ .- d.scratchm3
    # compute ts = R^T (Q - W' P^{-1} W)^{-1} R
    eigfact = eigen!(Symmetric(d.scratchm3))
    mul!(d.scratchv1, transpose(eigfact.vectors), d.scoreγ)
    T = eltype(d.scoreγ)
    atol = 1e-8 # tolerance for determining rank
    rk = 0 # rank
    ts = zero(T)
    for j in 1:d.q
        if eigfact.values[j] > atol
            ts += abs2(d.scratchv1[j]) / eigfact.values[j]
            rk += 1
        end
    end
    pval = ts ≤ 0 ? 1.0 : ccdf(Chisq(rk), ts)
end

###########################
# LRT test
###########################

struct OrdinalMultinomialLrtTest{TY<:Integer, T<:BlasReal, TL<:GLM.Link}
    "`q`: number of covariates to test significance"
    q::Int
    "`Z`: n-by-q covariate matrix to test significance"
    Z::Matrix{T}
    "`nm`: fitted null model"
    nm::OrdinalMultinomialModel{TY,T,TL}
    "`Xaug`: augmented X matrix for full model"
    Xaug::Matrix{T}
end

function OrdinalMultinomialLrtTest(nm::OrdinalMultinomialModel, Z::AbstractVecOrMat)
    TY, T, TL = eltype(nm.Y), eltype(nm.X), typeof(nm.link)
    q = size(Z, 2)
    Xaug = Matrix{T}(undef, nm.n, nm.p + q)
    @views copyto!(Xaug[:, 1:nm.p], nm.X)
    @views copyto!(Xaug[:, nm.p+1:end], Z)
    OrdinalMultinomialLrtTest{TY, T, TL}(q, Z, nm, Xaug)
end

function polrtest(
    d::OrdinalMultinomialLrtTest;
    solver::TSOLVER = NLopt.Optimizer()
    ) where {TSOLVER <: MOI.AbstractOptimizer}
    if typeof(solver) == NLopt.Optimizer
        algo = MOI.get(solver, MOI.RawOptimizerAttribute("algorithm"))
        if algo == :none
            MOI.set(solver, MOI.RawOptimizerAttribute("algorithm"), :LD_SLSQP)
            MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), 4000)
        end
    end
    am = polr(d.Xaug, d.nm.Y, d.nm.link, solver; wts = d.nm.wts)
    ccdf(Chisq(d.q), deviance(d.nm) - deviance(am))
end


text:
@testset "logseries dist" begin
  @test logseriescdf(0.01)[1:3] ≈ [0.0, 0.994992, 0.999967] atol=1.0e-5
  @test logseriescdf(0.000000000001)[1:5] ≈ [0., 1., 1., 1., 1.] atol=1.0e-2
  @test logseriesquantile(0.9, [0.25, 0.5, 0.75]) == [1, 2, 5]
  Random.seed!(43)
  v = logseriesquantile(0.4, rand(1000000))
  @test mean(v) ≈ 1.304 atol=1.0e-2
  @test std(v) ≈ 0.687 atol=1.0e-2
  @test skewness(v) ≈ 3.1 atol=1.0e-2
  @test kurtosis(v) ≈ 13.5 atol=1.0
end
@testset "stable levy dist" begin
  Random.seed!(43)
  @test levyel(2.) ≈ 0.18170339379413047
  Random.seed!(43)
  @test levygen(2., [0.2, 0.4, 0.6, 0.8]) ≈ [0.159748, 0.181703, 3.20539, 6.91497] atol=1.0e-4
  Random.seed!(43)
  ti = tiltedlevygen([0.2, 0.4, 0.6], 2.)
  @test ti ≈ [0.00409133, 0.0574922, 0.755548] atol=1.0e-5
end
@testset "nested copulas data generators" begin
  @test Ginv(0.5, 0.5) ≈ 1.2732395447351625
  @test InvlaJ(4, 0.5) ≈ 0.7265625
  @test sampleInvlaJ(0.5, 0.5) == 1
  @test sampleInvlaJ(0.5, 0.8) == 8
  Random.seed!(43)
  @test elInvlaF(4., 2.) == 7
  Random.seed!(43)
  @test elInvlaF(4., .5) == 16
  Random.seed!(43)
  @test nestedfrankgen(5., 3., [1, 1, 2]) == [3, 1, 49]
end

ext: jl


text:
using Twitter, DataFrames

#Authentication using consumer & oauth values
twitterauth("6nOtpXmf...",
            "sES5Zlj096S...",
            "98689850-Hj...",
            "UroqCVpWKIt...")

############################################## timelines.jl ##############################################

#Get mentions (@ replies to authenticated account), defaults to last 20
mentions_default20 = get_mentions_timeline();
DataFrame(mentions_default20) #Convert Array{TWEETS,1} to DataFrame

#Get tweets from any public user timeline, defaults to last 20
user_timeline_default20 = get_user_timeline(; options = Dict{Any, Any}("screen_name" => "randyzwitch"))

#Get tweets from authenticated user's timeline, defaults to 20
home_timeline_default20 = get_home_timeline();

#Get most recent re-tweets from authenticated user, defaults to 20
retweets_default20 = get_retweets_of_me()

############################################## tweets.jl ##########################################

#Get retweets by id number of tweet
get_retweets_by_id = get_retweets_id("434685122671939584")

#Get single tweet by id number
get_tweet_by_id = get_single_tweet_id("434685122671939584")

#Remove a tweet by a specific id
remove_tweet = post_destroy_single_tweet_id("439427946487742465")

#Post a status update to Twitter
status_update = post_status_update("I can't believe it snowed today, then is now 45 degrees and all the snow is gone. Strange times.")

#Retweet someone based on id
retweet = post_status_retweet_id("439428476106702848")

#post_status_update_media

#Get embedded tweet code: either "id" or "url" required by API, passed in through options Dict in Julia
get_embedded_tweet_code = get_oembed(options = Dict{Any, Any}("id" => "434685122671939584"))

#Get ids for people who retweeted a given status
retweeters_ids = get_retweeters_id("431782228826550272")

############################################## search.jl ##########################################

#Return tweets matching search term
duke_tweets = get_search_tweets("#Duke")

############################################## streaming.jl #######################################

#post_streaming_filtered

#get_streaming_sample

#get_firehose

#get_user_streaming

#get_site_streaming

############################################## dm.jl ##############################################

#Get direct messages received
dm_received = get_direct_messages()

#Get direct messages sent
dm_sent = get_direct_messages_sent()

#Get single direct message received
get_single_dm = get_direct_messages_show("435234003939905536")

#Delete a received DM
destroy = post_direct_messages_destroy("438517039385178113")

#Send DM to a user
send_dm = post_direct_messages_send("Testing from Julia, this might disappear later"; options = Dict{Any, Any}("screen_name" => "jowanza"))

############################################## friendsfollowers.jl ################################

#Get ids of friends that you don't want to see retweets from
friends_but_dont_see_rts = get_friendships_no_retweets()

#Get numeric ids of all friends
get_friends = get_friends_ids(options = Dict{Any, Any}("screen_name" => "randyzwitch"))

#Get numeric ids of followers
get_followers = get_followers_ids()

#Get authenticated user relationship to other users
friendship_status = get_friendships_lookup(options = Dict{Any, Any}("screen_name" => "usujason, hilamonster, kyrieirving"))

#Get list of pending requests to authenticated user
pending_friendships_incoming = get_friendships_incoming()

#Get list of pending requests by authenticated user outstanding
pending_friendships_outgoing = get_friendships_outgoing()

#Follow a user
add_friend = post_friendships_create(options = Dict{Any, Any}("screen_name" => "kyrieirving"))

#Unfollow a user
unfollow = post_friendships_destroy(options = Dict{Any, Any}("screen_name" => "kyrieirving"))

#Enable/disable retweets and device updates from specified user
friend_update = post_friendships_update(options = Dict{Any, Any}("screen_name" => "jowanza", "retweets" => "true", "device" => "false"))

#Get friendship details between any two arbitrary users
detailed_info_two_users = get_friendships_show(options = Dict{Any, Any}("source_screen_name" => "randyzwitch", "target_screen_name" => "usujason"))

#Get friendship list as names
friends_names = get_friends_list(options = Dict{Any, Any}("screen_name" => "randyzwitch"))

#Get followers names
followers_names = get_followers_list(options = Dict{Any, Any}("screen_name" => "randyzwitch"))

############################################## users.jl ###########################################

#Get settings for authenticated account
settings = get_account_settings();

#Verify credentials
verify_cred = get_verify_credentials()

#Change account settings
change_settings = post_account_settings(options = Dict{Any, Any}("lang" => "en", "sleep_time_enabled" => "false"))

#Allow for SMS or regular notifications from Twitter
update_delivery = post_update_delivery_device("sms")

#Update profile information
update_profile = post_update_profile(options = Dict{Any, Any}("name" => "Julia Twitter User"))

#Update picture background
update_background = post_update_profile_background(options = Dict{Any, Any}("tile" => "true"))

#Update profile colors
update_colors = post_update_profile_colors(options = Dict{Any, Any}("profile_background_color" => "3D3D3D", "profile_link_color" => "0000FF"))

#Update profile image
update_image = post_update_profile_image(options = Dict{Any, Any}("image" => "ABCDEFGH..."))

#Get list of users authenticated user is blocking
blocked_users = get_blocks_list()

#Get list of ids blocked by authenticated user
blocked_ids = get_blocks_ids()

#Block a user
block_a_user = post_blocks_create(options = Dict{Any, Any}("screen_name" => "BandonDunesGolf"))

#Unblock a user
unblock_a_user = post_blocks_destroy(options = Dict{Any, Any}("screen_name" => "BandonDunesGolf"))

#Get information about specific users
multiple_user_info = get_users_lookup(options = Dict{Any, Any}("screen_name" => "kyrieirving, jpiz1"))

#Get information about a single user
single_user_info = get_users_show(options = Dict{Any, Any}("screen_name" => "amilejefferson"))

#Search for users based on a topic
search_for_users = get_users_search("data science")

#Get users that authenticate user account can contribute to
contributees = get_users_contributees(options = Dict{Any, Any}("screen_name" => "randyzwitch"))

#Get users that can contribute to the authenticated account
contributors = get_users_contributors(options = Dict{Any, Any}("screen_name" => "randyzwitch"))

#Remove profile banner
remove_banner = post_account_remove_profile_banner()

#post_account_update_profile_banner
add_banner = post_account_update_profile_banner("base64_banner_data"; options = Dict{Any, Any}("width" => "1040", "height" => "520"))

#Get profile banner
profile_banner = get_profile_banner(options = Dict{Any, Any}("screen_name" => "randyzwitch"))

############################################## suggestedusers.jl ##################################

#Get user suggestions for a given category of suggested users
user_suggestions_twitter_slug = get_user_suggestions_slug("twitter")

#Get suggested user categories
suggested_categories = get_user_suggestions()

#Get users in suggested user group
members_suggested_category = get_user_suggestions_slug_members("sochi")

############################################## favorites.jl #######################################

#Get favorites list for authenticated user (when no arguments are specified)
favorites_auth_user = get_favorites_list()

#Unfavorite a tweet
unfavorite_tweet = post_favorites_destroy("439089844468588544")

#Favorite a tweet
favorite_tweet = post_favorites_create("439801028910010369")

############################################## lists.jl ###########################################

#Get lists - by default, no argument means lists subscribed to by authenticated user
lists_subscribed_to = get_lists()

#Get tweets for a given list id
tweets_given_list = get_lists_statuses(options = Dict{Any, Any}("list_id" => "4300970"))

#Remove a member from a list
remove_user = post_lists_members_destroy(options = Dict{Any, Any}("owner_screen_name" => "randyzwitch", "slug" => "julia-users", "screen_name" => "viralbshah"))

#Get lists which user has been added to
list_adds = get_lists_memberships(options = Dict{Any, Any}("screen_name" => "randyzwitch"))

#Get subscribers of a certain list id
list_subscribers = get_lists_subscribers(options = Dict{Any, Any}("list_id" => "4300970"))

#Subscribe authenticated user to list
subscribe_list = post_lists_subscribers_create(options = Dict{Any, Any}("list_id" => "4300970"))

#Check if user is subscribed to a list
check_user_subscriber = get_lists_subscribers_show(options = Dict{Any, Any}("list_id" => "4300970", "screen_name" => "randyzwitch"))

#Unsubscribe authenticated user from list
unsub_list = post_lists_subscribers_destroy(options = Dict{Any, Any}("owner_screen_name" => "nduley", "slug" => "omniture-webanalytics"))

#Add multiple users to list; no spaces allowed in list of screen_names
multiple_create = post_lists_members_createall(options = Dict{Any, Any}("owner_screen_name" => "randyzwitch", "slug" => "julia-users", "screen_name" => "stefankarpinski,viral_b_shah,jowanza"))

#Check if a member is a member of a specific list
mem_show = get_lists_members_show(options = Dict{Any, Any}("owner_screen_name" => "randyzwitch", "slug" => "julia-users", "screen_name" => "johnmyleswhite"))

#Get members of a list
list_members = get_lists_members(options = Dict{Any, Any}("list_id" => "4300970"))

#Add single user to list
add_one_user = post_lists_members_create(options = Dict{Any, Any}("owner_screen_name" => "randyzwitch", "slug" => "julia-users", "screen_name" => "johnmyleswhite"))

#Delete a list
destroy = post_lists_destroy(options = Dict{Any, Any}("owner_screen_name" => "randyzwitch", "slug" => "julia-folks"))

#Update characteristics of a list
update_list = post_lists_update(options = Dict{Any, Any}("list_id" => "106887912", "name" => "julia folks"))

#Create a new list owned by authenticated user
create_list = post_lists_create("julia users")

#Get specified list
get_list = get_lists_show(options = Dict{Any, Any}("list_id" => "4300970"))

#Get lists a user is subscribed to
get_list_subs = get_lists_subscriptions(options = Dict{Any, Any}("screen_name" => "randyzwitch"))

#Remove multiple members from a list
remove_multiple = post_lists_members_destroyall(options = Dict{Any, Any}("owner_screen_name" => "randyzwitch", "slug" => "julia-users", "screen_name" => "johnmyleswhite,jowanza"))

#Get lists user owns
owned_lists = get_lists_ownerships(options = Dict{Any, Any}("screen_name" => "nduley"))

############################################## savedsearches.jl ###################################

#Get saved searches for authenticated user
saved_searches = get_saved_searches_list()

#Get info about saved search by id
saved_search_info = get_saved_searches_show_id("23777812")

#Create a saved search
create_saved_search = post_saved_searches_create("#julialang")

#Remove a saved search by id
remove_saved_search = post_saved_searches_destroy_id("315945769")

############################################## geo.jl #############################################

#Get information about a specific place
geo_info = get_geo_id_place_id("df51dec6f4ee2b2c")

#Get places that can be used for a place_id from a set of lat/long coordinates
place_ids = get_geo_reverse_geocode("37.7821120598956", "-122.400612831116")

#Search for place_ids based on information that you know (IP, lat/long, etc.)
search_for_places = get_geo_search(options = Dict{Any, Any}("ip" => "74.125.19.104"))

#Get places near a given location
similar_places = get_geo_similar_places("37.7821120598956", "-122.400612831116", "Twitter HQ")

############################################## trends.jl ##########################################

#Get trends based on a certain WOEID
trends_place1 = get_trends_place("1")

#Get locations where Twitter has trend data
trend_locations = get_trends_available()

#Get trends closest to a specified location
close_trends = get_trends_closest("37.781157", "-122.400612831116")

############################################## spam.jl ############################################

#Report a user as a spammer (also blocks user for authenticated user)
report_spammer = post_users_report_spam(options = Dict{Any, Any}("screen_name" => "TEAMF0LL0W"))

############################################## oauth.jl ###########################################

#get_oauth_authenticate

#get_oauth_authorize

#post_oauth_access_token

#post_oauth_request_token

#post_oauth2_token

#post_oauth2_invalidate_token

############################################## help.jl ############################################

#Get configuration for authenticated user
user_configuration = get_help_configuration()

#Get languages/codes that Twitter supports
languages = get_help_languages()

#Get privacy policy
privacy_policy = get_help_privacy()

#Get terms of service
terms_of_service = get_help_tos()

#Get application rate limits
rate_limits = get_application_rate_limit_status()


text:
using MultipleViewGeometry, Test
using MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTypes
using MultipleViewGeometry.ModuleConstraints
using MultipleViewGeometry.ModuleConstruct
using MultipleViewGeometry.ModuleNoise
using LinearAlgebra, Random
using StaticArrays

# Fix random seed.
Random.seed!(1234)

𝒳 = [Point3D(x,y,z) for x=-1:1:10 for y=-1:1:10 for z=-1:1:10]
𝒟 = perturb(GaussianNoise(), 1.0, tuple(𝒳) )
𝒳ʹ = 𝒟[1]

N = length(𝒳ʹ)
for n = 1:N
    @test !isapprox(sum(abs.(𝒳[1][1:3]-𝒳ʹ[1][1:3])/4), 0.0; atol = 1e-12)
end


text:
module OpenGL

import OpenGLStd.glVertex2i, OpenGLStd.glVertex3i, OpenGLStd.glVertex2d, OpenGLStd.glVertex3d, OpenGLStd.glVertex4d,
       OpenGLStd.glTexCoord2i, OpenGLStd.glTexCoord3i, OpenGLStd.glTexCoord4i, OpenGLStd.glTexCoord2d, OpenGLStd.glTexCoord3d,
       OpenGLStd.glTexCoord4d, OpenGLStd.glNormal3d, OpenGLStd.glColor3f, OpenGLStd.glColor4f, OpenGLStd.glColor3b, OpenGLStd.glColor4b,
       OpenGLStd.glScaled, OpenGLStd.glTranslated, OpenGLStd.glRotated, OpenGLStd.glEnable, OpenGLStd.glDisable, OpenGLStd.glBegin,
       OpenGLStd.glEnd, OpenGLStd.glPushMatrix, OpenGLStd.glLoadIdentity

import Images.imread

#Handy stuff to make use of Julia features.
export glVertex, glColor, glColorb, glTexCoord, glNormal,
       glScale, glTranslate, glRotate, glRotateRad,
       glPrimitive, glPushed,
       unitFrame, unitFrameFrom, unitFrameTo,
       quad, centeredQuad, glimread #glimg

#TODO upgrade so glBegin/glPushMatrix can be used directly.
#glEnable (covered by autoFFI)

#Macro to conveniently also support tupled arguments.
#(TODO probably move elsewhere)
macro also_tuple(of, to)
    function getwhich(given::Expr)
        if given.head == symbol(":")
            assert( length(given.args) <=2 )
            return given.args[1]:given.args[2]
        end
        if given.head == :tuple
            return given.args
        end
        error("Invalid specification of the lengths of tuples: $to Should be a tuple or range.") #TODO this could be better.
    end
    getwhich(given::Integer) = {given}
    ret = {}
    #  nT(n) = map((i->(:T)), 1:n)
    for n = getwhich(to)
        push!(ret,
        Expr(:function,
            #TODO templated it earlier, but unequal argument types make that not work.
            #                   {Expr(:call, {Expr(:curly, {of,:T}, Any),
            #                                 Expr(symbol("::"),
            #                                      {:x,Expr(:tuple,nT(n),Any)},Any)},Any),
            {Expr(:call,
                {of, Expr(symbol("::"),
                    {:x,Expr(:tuple, map(i->(:Any), 1:n), Any)},
                    Any)},Any),
                Expr(:call,cat(1,{of},map((i)->:(x[$i]), 1:n)), Any)},
            Any))
    end
    return esc(Expr(:block,ret, Any))
end

#Overloading stuff
#Vertices
glVertex(i::Integer,j::Integer) = glVertex2i(i,j)
glVertex(i::Integer,j::Integer,k::Integer) = glVertex3i(i,j,k)
glVertex(i::Integer,j::Integer,k::Integer,l::Integer) = glVertex3i(i,j,k,l)

glVertex(x::Number,y::Number) = glVertex2d(x,y)
glVertex(x::Number,y::Number,z::Number) = glVertex3d(x,y,z)
glVertex(x::Number,y::Number,z::Number,w::Number) = glVertex4d(x,y,z,w)

@also_tuple glVertex 2:4

function glVertex{T}(v::Array{T,1})
    if length(v)==3
        return glVertex(v[1],v[2],v[3])
    end
    if length(v)==2
        return glVertex(v[1],v[2])
    end
    if length(v)==4
        return glVertex(v[1],v[2],v[3],v[4])
    end
end

#Texture coordinates
glTexCoord(i::Integer,j::Integer) = glTexCoord2i(i,j)
glTexCoord(i::Integer,j::Integer,k::Integer) = glTexCoord3i(i,j,k)
glTexCoord(i::Integer,j::Integer,k::Integer,l::Integer) = glTexCoord4i(i,j,k,l)

glTexCoord(x::Number,y::Number) = glTexCoord2d(x,y)
glTexCoord(x::Number,y::Number,z::Number) = glTexCoord3d(x,y,z)
glTexCoord(x::Number,y::Number,z::Number,w::Number) = glTexCoord4d(x,y,z,w)
@also_tuple glTexCoord 2:4

glNormal(x::Number,y::Number,z::Number) = glNormal3d(x,y,z)
#glNormal(i::Integer,j::Integer,k::Integer) = glNormal3b(i,j,k)
@also_tuple glNormal 2:4

glColor(r::Number,g::Number,b::Number) = glColor3f(r,g,b)
glColor(r::Number,g::Number,b::Number,a::Number) = glColor4f(r,g,b,a)
@also_tuple glColor 3:4

function glColor{T}(v::Array{T,1})
    if length(v)==3
        return glColor(v[1],v[2],v[3])
    end
    if length(v)==2
        return glColor(v[1],v[2])
    end
    if length(v)==4
        return glColor(v[1],v[2],v[3],v[4])
    end
end

function glColor{T}(v::Array{T,2})
    h, w = size(v)
    if h == 1
        return glColor(vec(v))
    else
        error("Array must have 1 row (i.e., it must be like a vector).")
    end
end

glColorb(r::Integer,g::Integer,b::Integer) = glColor3b(r,g,b)
glColorb(r::Integer,g::Integer,b::Integer,a::Integer) = glColor4b(r,g,b,a)
@also_tuple glColorb 3:4

function glColorb(v::Array{Integer,1})
    if length(v)==3
        return glColorb(v[1],v[2],v[3])
    end
    if length(v)==2
        return glColorb(v[1],v[2])
    end
    if length(v)==4
        return glColorb(v[1],v[2],v[3],v[4])
    end
end

function glColorb{T}(v::Array{T,2})
    h, w = size(v)
    if h == 1
        return glColorb(vec(v))
    else
        error("Array must have 1 row (i.e., it must be like a vector).")
    end
end

glScale(x::Number,y::Number,z::Number) = glScaled(x,y,z)
glScale(x::Number,y::Number) = glScaled(x,y,1)
glScale(s::Number) = glScaled(s,s,s)
@also_tuple glScale 1:3

glTranslate(x::Number,y::Number,z::Number) = glTranslated(x,y,z)
glTranslate(x::Number,y::Number) = glTranslated(x,y,0)
@also_tuple glTranslate 2:3

glRotate(angle::Number, nx::Number,ny::Number,nz::Number) = glRotated(angle, nx,ny,nz)
glRotate(angle::Number) = glRotated(angle,0,0,1)
glRotate(angle::Number, n::(Number,Number,Number)) = glRotate(angle, n[1],n[2],n[3])

#Damn degrees.
glRotateRad(angle::Number) = glRotate(angle*180/pi)
glRotateRad(angle::Number, nx::Number,ny::Number,nz::Number) = glRotate(angle*180/pi, nx,ny,nz)
glRotateRad(angle::Number, n::(Number,Number,Number)) = glRotate(angle*180/pi, n)

#Enabling lists of stuff.
type _GlEnable
    things::Vector
end

function glEnable(things::Vector)
    for thing in things
        glEnable(thing)
    end
    return _GlEnable(things)
end
glEnable(things...) = glEnable(things)

#Disabling lists of stuff
function glDisable(things::Vector)
    for thing in things
        glDisable(thing)
    end
end
glDisable(things...) = glDisable(things)

#The whole `begin` ... `end` structures are rather bad for the savings from the
# macros below..

#NOTE: if you `return` or something in the middle it won't end of course!
# (no `cl:unwind-protect`)

#glBegin for use _with_ @with
#TODO will want to upgrade that so glBegin can be used directly.
type _GlPrimitive
end

function glPrimitive(primitive)
    glBegin(primitive)
    return _GlPrimitive()
end

type _GlPushed
end

#glPushMatrix for with @with #TODO will want upgrade.
function glPushed()
    glPushMatrix()
    return _GlPushed()
end

#More functions
function unitFrame()
    glLoadIdentity()
    glTranslate(-1,-1)
    glScale(2)
end

#Map the given range to the unit range.
function unitFrameFrom(fx::Number,fy::Number,tx::Number,ty::Number)
    assert( fx!=tx && fy!=ty, "There might be a division by zero here.." )
    glScale(1/(tx-fx),1/(ty-fy))
    glTranslate(-fx,-fy)
end

typealias Vector2 Union((Number,Number),Vector) #(just for here)

unitFrameFrom(fr::Vector2, to::Vector2) = unitFrameFrom(fr[1],fr[2], to[1],to[2])

@also_tuple unitFrameFrom 2,4

#Map the unit range to the given range.
function unitFrameTo(fx::Number,fy::Number,tx::Number,ty::Number)
    glTranslate(fx,fy)
    glScale(tx-fx, ty-fy)
end
unitFrameTo(fr::Vector2, to::Vector2) = unitFrameTo(fr[1],fr[2], to[1],to[2])

@also_tuple unitFrameFrom 2,4

#Rectangle vertices (in QUADS, LINE_LOOP-able style)
function quad(fx::Number,fy::Number,tx::Number,ty::Number)
    glVertex(fx,fy)
    glVertex(fx,ty)
    glVertex(tx,ty)
    glVertex(tx,fy)
end
quad(fr::Vector2, to::Vector2) = quad(fr[1],fr[2], to[1],to[2])

@also_tuple quad 2,4

centeredQuad(x::Number,y::Number, r::Number) = quad(x-r,y-r, x+r, y+r)

centeredQuad(pos::Vector2, r::Number) = centeredQuad(pos[1],pos[2],r)

@also_tuple centeredQuad 2:3

# Going with a basic OpenGL RGB 1D-interleaved image-texture
function glimread(filename::String)
    img = imread(filename)
    w = size(img,2)
    h = size(img,3)
    img1D = Array(Uint8,w*h*3)
    for j=1:w
        for k=1:h
            for i=1:3
                img1D[3*w*(j-1)+3*(k-1)+i] = img[i,w-j+1,k]
            end
        end
    end
    return img1D, w, h
end

end


text:
module NahaJuliaLib

include("type_hierarchy.jl")
include("uri_utils.jl")
include("properties.jl")
include("trace.jl")
include("trace_analysis.jl")

end


text:
"""
    disrelation(A::AbstractMatrix)

Calculate a disrelation two-dimensional correlation spectrum.

`A` is a dataset with multiple spectra in rows.
They don't need to be at equally spaced measurements of the perturbation.
Normally, the data is given as the difference from a reference spectrum (see
below).

# Examples
```jldoctest
julia> A = rand(50, 1000);

julia> Λ = disrelation(A .- mean(A, dims=1));

julia> size(Λ)
(1000, 1000)
```
"""
function disrelation(A::AbstractMatrix)
    Φ = synchronous(A)
    Λ = similar(Φ)

    n = size(Λ, 1)
    for ν₁ in 1:n, ν₂ in 1:n
        Λ[ν₁, ν₂] = √(Φ[ν₁, ν₁] * Φ[ν₂, ν₂] - Φ[ν₁, ν₂]^2)
    end
    return Λ
end


text:
############################################################
# IJulia/ Jupyter output
############################################################
function inNotebook()
  return isdefined(Main, :IJulia) && Main.IJulia.inited
end

function math_html(io::IO, a::PolyElem)
  f = "$a"
  f = replace(f, "*" => "")
  f = replace(f, r"\^([0-9]*)" => s"^{\1}")
  print(io, f)
end

function Base.show(io::IO, ::MIME"text/html", a::PolyElem)
  print(io, "\$")
  math_html(io, a)
  print(io, "\$")
end

function math_html(io::IO, a::AnticNumberField)
  n = find_name(a)
  if n === nothing || !get(io, :compact, false)
    print(io, "\\text{Number field over Rational field with defining polynomial }")
    math_html(io, a.pol)
  else
    print(io, string(n))
  end
end

function Base.show(io::IO, ::MIME"text/html", a::AnticNumberField)
  print(io, "\$")
  math_html(io, a)
  print(io, "\$")
end

function Base.show(io::IO, ::MIME"text/html", A::Fac{T}) where {T}
  print(io, "\$")
  math_html(io, A)
  print(io, "\$")
end

function math_html(io::IO, a)
  try
    show(io, "text/html", a)
  catch e
    if isa(e, MethodError)
      show(io, a)
    else
      rethrow(e)
    end
  end
end

#=
function Base.show(io::IO, ::MIME"text/html", a)
  print(io, "\$")
  math_html(io, a)
  print(io, "\$")
end
=#

#=
function Base.show(io::IO, b::HTML{nf_elem})
  math_html(io, b.content)
end

function Base.show(io::IO, b::HTML{AnticNumberField})
  math_html(io, b.content)
end
=#

function Base.show(io::IO, ::MIME"text/html", T::Tuple)
  print(io, "(")
  for i =1:length(T)
    show(IOContext(io, :compact => true), "text/html", T[i])
    if i<length(T)
      print(io, ", ")
    end
  end
  print(io, ")")
end

function math_html(io::IO, a::nf_elem)
  s = parent(a).S
  if s in [:_a, Symbol("_\$")]
    parent(a).S = Symbol("\\alpha")
    fs = string(a)
    parent(a).S = s
  else
    fs = string(a)
  end
  f = replace(fs, "*" => "")
  f = replace(f, r"\^([0-9]*)" => s"^{\1}")
  f = replace(f, r"([0-9]*)//([0-9]*)" => s"\\frac{\1}{\2}")
  if true
    print(io, f)
  else
    print(io, "\\toggle{$f}{$fs}\\endtoggle")
  end
end

function Base.show(io::IO, ::MIME"text/html", a::nf_elem)
  print(io, "\$")
  math_html(io, a)
  print(io, "\$")
end

function math_html(io::IO, A::Fac{T}) where {T}
  empty = true
  if !isone(A.unit)
    math_html(io, A.unit)
    empty = false
  end
  for k = (collect(keys(A.fac)))
    if !empty
      print(io, "\\cdot")
    end
    empty = false
    AbstractAlgebra.needs_parentheses(k) && print(io, "(")
    math_html(io, k)
    AbstractAlgebra.needs_parentheses(k) && print(io, ")")
    if A.fac[k] != 1
      print(io, "^{")
      math_html(io, A.fac[k])
      print(io, "}")
    end
  end
end

Base.show(io::IO, ::MIME"text/html", a::Integer) = show(io, "text/html", fmpz(a))
math_html(io::IO, a::Integer) = math_html(io, fmpz(a))

function Base.show(io::IO, ::MIME"text/html", a::fmpz)
  print(io, "\$")
  math_html(io, a)
  print(io, "\$")
end

function math_html(io::IO, a::fmpz)
  nd = ndigits(a, 10)
  if nd < 20
    print(io, a)
  else
    d = string(abs(a) % 10^5)
    d = "00000"[1:5-length(d)] * d
    d = "("*string(div(a, fmpz(10)^(nd-5))) *  "..($nd \\text{ digits}).." * d *")"
#    print(io, "\$\\require{action}\$")
    print(io, "\\toggle{$d}{$a}\\endtoggle")
  end
end

function math_html(io::IO, M::MatElem)
  print(io, "\\begin{bmatrix}")
  for i=1:nrows(M)
    for j=1:nrows(M)
      math_html(io, M[i,j])
      if j<nrows(M)
        print(io, "&")
      else
        print(io, "\\\\")
      end
    end
  end
  print(io, "\\end{bmatrix}")
end

function Base.show(io::IO, ::MIME"text/html", M::MatElem)
  print(io, "\$")
  math_html(io, M)
  print(io, "\$")
end

function math_html(io::IO, a::fmpq)
  if denominator(a) == 1
    math_html(io, numerator(a))
  else
    print(io, "\\frac{")
    math_html(io, numerator(a))
    print(io, "}{")
    math_html(io, denominator(a))
    print(io, "}")
  end
end

function Base.show(io::IO, ::MIME"text/html", a::fmpq)
  print(io, "\$")
  math_html(io, a)
  print(io, "\$")
end

math_html(io::IO, a::Rational) = math_html(io, fmpq(a))
function Base.show(io::IO, ::MIME"text/html", a::Rational)
  print(io, "\$")
  math_html(io, a)
  print(io, "\$")
end

#= infinite recursion through generic math_html, so don't
function Base.show(io::IO, ::MIME"text/html", a)
  print(io, "\$")
  math_html(io, a)
  print(io, "\$")
end
=#

function math_html(io::IO, l::Array{T, 1}) where {T}
  print(io, "[")
  for i in 1:length(l)
    if i>1
      print(io, ", ")
    end
    math_html(IOContext(io, :compact => true), l[i])
  end
  print(io, "]")
end

function Base.show(io::IO, ::MIME"text/html", l::Array{T, 1}) where {T}
  print(io, "\$")
  math_html(io, l)
  print(io, "\$")
end

function math_html(io::IO, a::NfAbsOrdElem)
  math_html(io, elem_in_nf(a))
end

function math_html(io::IO, O::NfAbsOrd{AnticNumberField, nf_elem})
  c = get(io, :compact, false)
  if !c
    print(io, "\\text{Maximal order of }")
    math_html(io, nf(O))
    print(io, "\\text{ with basis }")
    math_html(io, basis(O))
    return
  end

  n = find_name(O)
  if !(n===nothing)
    print(io, string(n))
    return
  end
  n = find_name(nf(O))
  if n === nothing
    print(io, "\\text{Maximal order of }")
    math_html(io, nf(O))
  else
    print(io, "\\mathcal O_$(string(n))")
  end
end

function Base.show(io::IO, ::MIME"text/html", a::NfAbsOrd{AnticNumberField, nf_elem})
  print(io, "\$")
  math_html(io, a)
  print(io, "\$")
end

function math_html(io::IO, M::Map)
  n = find_name(M)
  cio = IOContext(io, :compact => true)
  if n === nothing
    print(io, "\\text{Map from }")
  else
    print(io, string(n))
    print(io, ": ")
  end
  n = find_name(domain(M))
  if n === nothing
    math_html(cio, domain(M))
  else
    print(io, string(n))
  end
  print(io, "\\to ")
  n = find_name(codomain(M))
  if n === nothing
    math_html(cio, codomain(M))
  else
    print(io, string(n))
  end
end

function Base.show(io::IO, ::MIME"text/html", M::Map)
  print(io, "\$")
  math_html(io, M)
  print(io, "\$")
end

function math_html(io::IO, I::NfAbsOrdIdlSet)
  print(io, "\\text{Set of ideals of }")
  n = find_name(order(I))
  if n === nothing || !get(io, :compact, false)
    math_html(IOContext(io, :compact => true), order(I))
  else
    print(io, string(n))
  end
end

function Base.show(io::IO, ::MIME"text/html", G::GrpAbFinGen)
  print(io, "\$")
  math_html(io, G)
  print(io, "\$")
end

function math_html(io::IO, G::GrpAbFinGen)
  n = find_name(G)
  if !(n === nothing) && get(io, :compact, false)
    print(io, string(n))
    return
  end
  if issnf(G)
    if get(io, :compact, false)
      print(io, "\\text{Abelian Group with Invariants: }")
    else
      print(io, "\\text{GrpAb: }")
    end
    show_snf_structure(io, G, "\\times")
  else
    print(io, G)
  end
end

function math_html(io::IO, R::PolyRing)
  n = find_name(R)
  if !(n === nothing) && get(io, :compact, false)
    print(io, string(n))
    return
  end
  print(io, "\\text{Polynomial ring over }")
  print(IOContext(io, :compact => true), base_ring(R))
end

function Base.show(io::IO, ::MIME"text/html", K::PolyRing)
  print(io, "\$")
  math_html(io, K)
  print(io, "\$")
end

function math_html(io::IO, K::NfRel)
  n = find_name(K)
  if !(n === nothing) && get(io, :compact, false)
    print(io, string(n))
    return
  end
  print(io, "\\text{Relative number field over }")
  math_html(IOContext(io, :compact => true), base_field(K))
  print(io, "\\text{ defined by }")
  math_html(io, K.pol)
end

function Base.show(io::IO, ::MIME"text/html", K::NfRel)
  print(io, "\$")
  math_html(io, K)
  print(io, "\$")
end

function Base.show(io::IO, ::MIME"text/html", b::Bool)
   print(io, b ? "true" : "false")
end

function find_name(A, M = Main)
  for a = names(Main)
    a === :ans && continue
    d = Meta.parse("$M.$a")
    try
      z = eval(d);
      if z === A
        return a
      end
    catch e
    end
  end
end


text:
# ---
# title: 623. Add One Row to Tree
# id: problem623
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/add-one-row-to-tree/description/>
# hidden: true
# ---
#
# Given the root of a binary tree, then value `v` and depth `d`, you need to add
# a row of nodes with value `v` at the given depth `d`. The root node is at
# depth 1.
#
# The adding rule is: given a positive integer depth `d`, for each NOT null tree
# nodes `N` in depth `d-1`, create two tree nodes with value `v` as `N's` left
# subtree root and right subtree root. And `N's` **original left subtree**
# should be the left subtree of the new left subtree root, its **original right
# subtree** should be the right subtree of the new right subtree root. If depth
# `d` is 1 that means there is no depth d-1 at all, then create a tree node with
# value **v** as the new root of the whole original tree, and the original tree
# is the new root's left subtree.
#
# **Example 1:**
#
#
#
#     Input:
#     A binary tree as following:
#            4
#          /   \
#         2     6
#        / \   /
#       3   1 5
#
#     **v = 1**
#
#     **d = 2**
#
#     Output:
#            4
#           / \
#          1   1
#         /     \
#        2       6
#       / \     /
#      3   1   5
#
#
#
# **Example 2:**
#
#
#
#     Input:
#     A binary tree as following:
#           4
#          /
#         2
#        / \
#       3   1
#
#     **v = 1**
#
#     **d = 3**
#
#     Output:
#           4
#          /
#         2
#        / \
#       1   1
#      /     \
#     3       1
#
#
# **Note:**
#
#   1. The given d is in range [1, maximum depth of the given tree + 1].
#   2. The given binary tree has at least one tree node.
#
#
## @lc code=start
using LeetCode

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


text:
using Documenter
using ParallelUtilities

DocMeta.setdocmeta!(ParallelUtilities, :DocTestSetup, :(using ParallelUtilities); recursive=true)

makedocs(;
    modules=[ParallelUtilities],
    authors="Jishnu Bhattacharya",
    repo="https://github.com/jishnub/ParallelUtilities.jl/blob/{commit}{path}#L{line}",
    sitename="ParallelUtilities.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jishnub.github.io/ParallelUtilities.jl",
        assets=String[],
    ),
    pages=[
        "ParallelUtilities" => "index.md",
        "Mapreduce" => "pmapreduce.md",
        "Examples" => [
            "pmapreduce" => "examples/pmapreduce.md",
            "SharedArrays" => "examples/sharedarrays.md",
            "Threads" => "examples/threads.md",
        ],
        "ClusterQueryUtils" => "clusterquery.md",
        "Reference" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/jishnub/ParallelUtilities.jl",
)


text:
# Important, for broadcast joins we cannot assume c and d have same number of columns:
# c could have more columns than d
rowcmp(::Tuple, i, ::Tuple{}, j) = 0

function rowcmp(tc::Tuple, i, td::Tuple, j)
    c, d = tc[1], td[1]
    let k = rowcmp(c, i, d, j)
        (k == 0) ? rowcmp(tail(tc), i, tail(td), j) : k
    end
end

function rowcmp(c::StructVector, i, d::StructVector, j)
    tc = Tuple(fieldarrays(c))
    td = Tuple(fieldarrays(d))
    return rowcmp(tc, i, td, j)
end

@inline function rowcmp(c::AbstractVector, i, d::AbstractVector, j)
    cmp(c[i], d[j])
end

struct GroupJoinPerm{LP<:GroupPerm, RP<:GroupPerm}
    left::LP
    right::RP
end

GroupJoinPerm(lkeys::AbstractVector, rkeys::AbstractVector, lperm=sortperm(lkeys), rperm=sortperm(rkeys)) =
    GroupJoinPerm(GroupPerm(lkeys, lperm), GroupPerm(rkeys, rperm))

function _pick(s, a, b)
    if a === nothing && b === nothing
        return nothing
    elseif a === nothing
        return (1:0, b[1]), (1, a, b)
    elseif b === nothing
        return (a[1], 1:0), (-1, a, b)
    else
        lp = sortperm(s.left)
        rp = sortperm(s.right)
        cmp = rowcmp(parent(s.left), lp[first(a[1])], parent(s.right), rp[first(b[1])])
        if cmp < 0
            return (a[1], 1:0), (-1, a, b)
        elseif cmp == 0
            return (a[1], b[1]), (0, a, b)
        else
            return (1:0, b[1]), (1, a, b)
        end
    end
end

function Base.iterate(s::GroupJoinPerm)
    l = iterate(s.left)
    r = iterate(s.right)
    _pick(s, l, r)
end

function Base.iterate(s::GroupJoinPerm, (select, l, r))
    (select <= 0) && (l = iterate(s.left, l[2]))
    (select >= 0) && (r = iterate(s.right, r[2]))
    _pick(s, l, r)
end

Base.IteratorSize(::Type{<:GroupJoinPerm}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:GroupJoinPerm}) = Base.HasEltype()
Base.eltype(::Type{<:GroupJoinPerm}) = Tuple{UnitRange{Int}, UnitRange{Int}}


text:
"""
spitchfork(u,lambda)

The nonlinearity f(u) = u^3 - lamba u. The dynamics for du/du = -f(u)
have a pitchfork bifurcation at lambda=0. The steady-state solution
u=0 is unique for lambda < 0 and there are three steady-state solutions
if lambda > 0. This is a simple-minded version of the buckling beam problem.

The function sptest(u) = spitchfork(u,.5) is the one I call in the testing.

"""
function spitchfork(u, lambda)
    fu = u^3 - lambda * u
    return fu
end

function sptest(u)
    lambda = 0.5
    spt = spitchfork(u, lambda)
    return spt
end

function spitchp(u, lambda)
    fp = 3 * u^2 - lambda
    return fp
end

function sptestp(u)
    lambda = 0.5
    sptp = spitchp(u, lambda)
    return sptp
end


text:
module Monads

export Monad, Maybe, Optional
export mreturn, mjoin, mbind, mmap
export maybe, ↣

using SimpleTraits
using ..Applicative: Callable
import Base: something, isnothing, convert

@traitdef Monad{T}

function mbind end
function mreturn end
function mjoin end

@traitfn @inline mreturn(::Type{M}, x) where {M; Monad{M}} = M(x)
@traitfn @inline mjoin(m::M) where {M; Monad{M}} = mbind(identity, m)
@traitfn @inline mbind(f::F, m::M) where {F,M; Callable{F},Monad{M}} =
    f(mjoin(m))
@traitfn @inline function mmap(f::F, m::M) where {F,M; Callable{F},Monad{M}}
    mbind(m) do x
        mreturn(M, f(x))
    end
end
@traitfn @inline ↣(m::M, f::F) where {M,F; Monad{M},Callable{F}} = mbind(f, m)

const Optional{T} = Union{Nothing,T}
const Maybe{T} = Optional{Some{T}}

@traitimpl Monad{Maybe}

Maybe{T}(x::T) where {T} = Some(x)
Maybe{T}(::Nothing) where {T} = nothing
Maybe(::Nothing) = nothing
Maybe(x::T) where {T} = Maybe{T}(x)
Maybe(m::Maybe) = m

@traitfn @inline mbind(f::F, m::Maybe) where {F; Callable{F}} =
    isnothing(m) ? nothing : f(something(m))

@inline maybe(b, m::Optional) = isnothing(m) ? b : m
@inline maybe(b) = m::Optional -> maybe(b, m)
@traitfn @inline maybe(f::F, b, m::Optional) where {F; Callable{F}} =
    maybe(b, Maybe(m) ↣ f)

end # module


text:
#Recursively traverses an expression and removes things like X^1, 1*X. Will not actually have any effect on the expression when used as a function, but will make it much easier to look at it for debugging, as well as if it is transformed to LaTeX code.
function recursive_clean!(expr)
    (expr isa Symbol) && (expr == :no___noise___scaling) && (return 1)
    (typeof(expr)!=Expr) && (return expr)
    for i = 1:length(expr.args)
        expr.args[i] = recursive_clean!(expr.args[i])
    end
    (expr.args[1] == :^) && (expr.args[3] == 1) && (return expr.args[2])
    if expr.args[1] == :*
        in(0,expr.args) && (return 0)
        i = 1
        while (i = i + 1) <= length(expr.args)
             if (typeof(expr.args[i]) == Expr) && (expr.args[i].head == :call) && (expr.args[i].args[1] == :*)
                 for arg in expr.args[i].args
                     (arg != :*) && push!(expr.args, arg)
                 end
             end
        end
        for i = length(expr.args):-1:2
            (typeof(expr.args[i]) == Expr) && (expr.args[i].head == :call) && (expr.args[i].args[1] == :*) && deleteat!(expr.args,i)
            (expr.args[i] == 1) && deleteat!(expr.args,i)
        end
        (length(expr.args) == 2) && (return expr.args[2])                   # We have a multiplication of only one thing, return only that thing.
        (length(expr.args) == 1) && (return 1)                              # We have only * and no real arguments.
        (length(expr.args) == 3) && (expr.args[2] == -1) && return :(-$(expr.args[3]))
        (length(expr.args) == 3) && (expr.args[3] == -1) && return :(-$(expr.args[2]))
    end
    if expr.head == :call
        (expr.args[1] == :/) && (expr.args[3] == 1) && (return expr.args[2])
        haskey(funcdict, expr.args[1]) && return funcdict[expr.args[1]](expr.args[2:end])
        in(expr.args[1],hill_name) && return hill(expr)
        in(expr.args[1],hillR_name) && return hillR(expr)
        in(expr.args[1],mm_name) && return mm(expr)
        in(expr.args[1],mmR_name) && return mmR(expr)
        (expr.args[1] == :binomial) && (expr.args[3] == 1) && return expr.args[2]
        #@isdefined($(expr.args[1])) || error("Function $(expr.args[1]) not defined.")
    end
    return expr
end

function chemical_arrows(rn::ModelingToolkit.ReactionSystem;
    expand = true, double_linebreak=false, mathjax=true, starred=false, kwargs...)
    str = starred ? "\\begin{align*}\n" : "\\begin{align}\n"
    eol = double_linebreak ? "\\\\\\\\\n" : "\\\\\n"

    mathjax && (str *= "\\require{mhchem}\n")
    backwards_reaction = false
    rxs = ModelingToolkit.equations(rn)
    @variables t

    # this should replace A(t) with A in equations however, currently substituter rewrites
    # things like x/y as inv(y)*x^1, which looks worse... for now we leave a stub that can
    # be updated when substitution preserves expressions better.
    # subber = ModelingToolkit.substituter([s(t) => s() for s in states(rn)])
    subber = x -> x

    for (i, r) in enumerate(rxs)
        if backwards_reaction
            backwards_reaction = false
            continue
        end
        str *= "\\ce{ "

        ### Expand functions to maths expressions
        rate = r.rate isa Operation ? Expr(subber(r.rate)) : r.rate
        expand && (rate = recursive_clean!(rate))
        expand && (rate = recursive_clean!(rate))

        ### Generate formatted string of substrates
        substrates = [latexraw("$(substrate[2]== 1 ? "" : "$(substrate[2]) * ") $(substrate[1].op.name)"; kwargs...) for substrate in zip(r.substrates,r.substoich)]
        isempty(substrates) && (substrates = ["\\varnothing"])

        str *= join(substrates, " + ")

        ### Generate reaction arrows
        prestr  = mathjax ? "[" : "[\$"
        poststr = mathjax ? "]" : "\$]"
        if i + 1 <= length(rxs) && issetequal(r.products,rxs[i+1].substrates) && issetequal(r.substrates,rxs[i+1].products)
            ### Bi-directional arrows
            rate_backwards = rxs[i+1].rate isa Operation ? Expr(subber(rxs[i+1].rate)) : rxs[i+1].rate
            expand && (rate_backwards = recursive_clean!(rate_backwards))
            expand && (rate_backwards = recursive_clean!(rate_backwards))
            str *= " &<=>"
            str *= prestr * latexraw(rate; kwargs...) * poststr
            str *= prestr * latexraw(rate_backwards; kwargs...) * poststr * " "
            backwards_reaction = true
        else
            ### Uni-directional arrows
            str *= " &->"
            str *= prestr * latexraw(rate; kwargs...) * poststr * " "
        end

        ### Generate formatted string of products
        products = [latexraw("$(product[2]== 1 ? "" : "$(product[2]) * ") $(product[1].op.name)"; kwargs...) for product in zip(r.products,r.prodstoich) ]
        isempty(products) && (products = ["\\varnothing"])
        str *= join(products, " + ")
        str *= "}$eol"
    end
    str = str[1:end-length(eol)] * "\n"

    str *= starred ? "\\end{align*}\n" : "\\end{align}\n"

    latexstr = Latexify.LaTeXString(str)
    Latexify.COPY_TO_CLIPBOARD && clipboard(latexstr)
    return latexstr
end

@latexrecipe function f(sys::ModelingToolkit.ReactionSystem)
    # Set default option values.
    env --> :chem
    cdot --> false

    return sys
end

function Latexify.infer_output(env, rs::ReactionSystem, args...)
    env in [:arrows, :chem, :chemical, :arrow] && return chemical_arrows

    error("The environment $env is not defined.")
    latex_function = Latexify.get_latex_function(rs, args...)

    return latex_function
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

"""
    initial_condition_rayleigh_taylor_instability(coordinates, t, equations::CompressibleEulerEquations2D)

Setup used for the Rayleigh-Taylor instability. Initial condition adapted from
- Shi, Jing, Yong-Tao Zhang, and Chi-Wang Shu (2003).
  Resolution of high order WENO schemes for complicated flow structures.
  [DOI](https://doi.org/10.1016/S0021-9991(03)00094-9).
- Remacle, Jean-François, Joseph E. Flaherty, and Mark S. Shephard (2003).
  An adaptive discontinuous Galerkin technique with an orthogonal basis applied to compressible
  flow problems.
  [DOI](https://doi.org/10.1137/S00361445023830)

The domain is [0, 0.25] x [0, 1]. Boundary conditions can be reflective wall boundary conditions on
all boundaries or
- periodic boundary conditions on the left/right boundaries
- Dirichlet boundary conditions on the top/bottom boundaries

This should be used together with `source_terms_rayleigh_taylor_instability`, which is
defined below.
"""
@inline function initial_condition_rayleigh_taylor_instability(x, t,
                                                               equations::CompressibleEulerEquations2D,
                                                               slope=1000)
  tol = 1e2*eps()

  if x[2] < 0.5
    p = 2*x[2] + 1
  else
    p = x[2] + 3/2
  end

  # smooth the discontinuity to avoid ambiguity at element interfaces
  smoothed_heaviside(x, left, right) = left + 0.5*(1 + tanh(slope * x)) * (right-left)
  rho = smoothed_heaviside(x[2] - 0.5, 2.0, 1.0)

  c = sqrt(equations.gamma * p / rho)
  # the velocity is multiplied by sin(pi*y)^6 as in Remacle et al. 2003 to ensure that the
  # initial condition satisfies reflective boundary conditions at the top/bottom boundaries.
  v = -0.025 * c * cos(8*pi*x[1]) * sin(pi*x[2])^6
  u = 0.0

  return prim2cons(SVector(rho, u, v, p), equations)
end

@inline function source_terms_rayleigh_taylor_instability(u, x, t,
                                                          equations::CompressibleEulerEquations2D)
  g = 1.0
  rho, rho_v1, rho_v2, rho_e = u

  return SVector(0.0, 0.0, g*rho, g*rho_v2)
end

# numerical parameters
volume_flux = flux_ranocha
solver = DGSEM(polydeg=3, surface_flux=flux_hll,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

# The domain is [0, 0.25] x [0, 1]
mapping(xi, eta) = SVector(0.25 * 0.5 * (1.0 + xi), 0.5 * (1.0 + eta))

num_elements_per_dimension = 32
cells_per_dimension = (num_elements_per_dimension, num_elements_per_dimension * 4)
mesh = StructuredMesh(cells_per_dimension, mapping)

initial_condition = initial_condition_rayleigh_taylor_instability
boundary_conditions = (
                       x_neg=boundary_condition_slip_wall,
                       x_pos=boundary_condition_slip_wall,
                       y_neg=boundary_condition_slip_wall,
                       y_pos=boundary_condition_slip_wall,
                      )

# # Alternative setup: left/right periodic BCs and Dirichlet BCs on the top/bottom.
# boundary_conditions = (
#                        x_neg=boundary_condition_periodic,
#                        x_pos=boundary_condition_periodic,
#                        y_neg=BoundaryConditionDirichlet(initial_condition),
#                        y_pos=BoundaryConditionDirichlet(initial_condition),
#                       )

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                                    source_terms=source_terms_rayleigh_taylor_instability,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.65)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback)

###############################################################################
# run the simulation

sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary


text:
function solve()
    n, x, t = [parse(Int, x) for x in split(readline())]
    t * ((n ÷ x) + (n % x > 0 ? 1 : 0))
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
#==============================================================================
    data wrangling tutorial I
        basic dataframe and wrangle excercises
        overview presentation video part 1

    project:    data-work
    author:     merlinr@
 =============================================================================#

using DesertIslandDisk, Markdown
using CSV, DataFrames, Pipe, Dates

# read / process / visualize #
# ========================== #
#=
recommended packages:
CSV, Parquet, JSON, XLSX, ExcelReaders, Taro, PDFIO
DataFrames, DataFramesMeta, Query, JuliaDB, JuliaDBMeta
Plots, StatsPlots, Makie, StatsMakie, Vegalite, Gadfly, Gaston, UnicodePlots
=#

# CSV.jl features #
# =============== #
md"""
## by default it:

* auto-delimiter detection
* skip leading/trail rows
* parse missing values
* parse date/time, boolean
* normalize col names
* select/drop cols on read
* define column types
* transpose
"""

# data cleaning #
# ============= #

# example data
df = DataFrame(CSV.File("data/youth_suicide.csv"))
dfcounty = DataFrame(CSV.File("data/youth_suicide.csv", drop=[2:18...]))

describe(df)

# SELECT
select(df, Not(), Between(), AsTable() )

# selecting and removing cols
select!(df, 1, 7:11)

select(df, [1, 3:7..., 19])

select(df,Not(["County", "Total (2012)", "Total (2011)", "Total (2009)"]))

# renaming
rename!(dfcounty, "County" => "my_nice_county", "Total (2008-2012)" => "all_muh_total")

# remove rows with missing data
dropmissing!(df)

dropmissing!(df, Between(3, 5))

# TIDY DATA #
# ========= #

md"""
# TIDY IS:
* every variable is a column
* every observation is a row
* every type of observed unit forms a table
"""

# tidy examples #
# ============= #

df = DataFrame(CSV.File("data/zhvi.csv"))
select!(df, :RegionName, 294:297)

# stack all grouping columns into rows
sdf = stack(df, Not(:RegionName); variable_name = :Date)

# refer to column name
sdf.ReportDate

describe(sdf, :eltype)

#=  need to retype that categorical value
│ Row │ variable   │ eltype                          │
│     │ Symbol     │ DataType                        │
├─────┼────────────┼─────────────────────────────────┤
│ 1   │ RegionName │ String                          │
│ 2   │ ReportDate │ CategoricalValue{String,UInt32} │
│ 3   │ value      │ Float64                         │
=#

# date handling #
# ============= #

sdf.ReportDate = [Date(get(x)) for x in sdf.ReportDate]
# 3660-element Array{Date,1}:
#  2020-01-31
#  2020-01-31
#  2020-01-31
# very nice

# turn back to wide format
df2 = unstack(sdf, :ReportDate, :value)

# retype from Union{Missing, T} if there
disallowmissing!(df2, 2:5)

# filter is part of Base - df argument is at the end
filter(:RegionName => ==("Abilene, TX"), df)

# filter from pattern regex match
filter(:RegionName => x -> occursin(r"Abil", x), df)

filter(:RegionName => x -> contains(x, "Abil"), df)

filter(2 => >(400_000), df)

sort(df, 4, rev=true)

# transforms #
# ========== #

# adding columns based on a transform
# ByRow applies a function to value in each row
transform(df, :RegionName => ByRow(length) => :RegionNameLength)

# more useful: extract the state code pattern from RegionName string
transform(df,
    :RegionName => ByRow(x -> begin
        m = match(r", (..)", x)
        if typeof(m) == Nothing ""
        else m[1]
    end end) => :State)

# alternately with function definiton, one line
tdf = transform(df, :RegionName => ByRow(
    function (x) m = match(r", (..)", x); if typeof(m) == Nothing "" else m[1] end end) => :State)

# grouping by #
# =========== #
# grouping allows aggregates on groups of a categorical or date column

# returns grouped dataframes
groupby(sdf, :ReportDate)

using Statistics

# agg over each group, union the groups back up
combine(groupby(sdf, :ReportDate), :value => median => :P50)

combine(groupby(sdf, :ReportDate), :value => x -> quantile(x, [0.25:0.25:0.75...]) => (:P25, :P50, :P75)
#=
4×2 DataFrame
│ Row │ ReportDate │ P50      │
│     │ Date       │ Float64  │
├─────┼────────────┼──────────┤
│ 1   │ 2020-01-31 │ 155163.0 │
│ 2   │ 2020-02-29 │ 155726.0 │
│ 3   │ 2020-03-31 │ 156410.0 │
│ 4   │ 2020-04-30 │ 156411.0 │
=#

# almost there, must explode final column
combine(groupby(sdf, :ReportDate), :value => x -> quantile(x, [0.25:0.25:0.75...]) => [:P25 :P50 :P75])
#=
4×2 DataFrame
│ Row │ ReportDate │ value_function                                     │
│     │ Date       │ Pair{Array{Float64,1},Array{Symbol,2}}             │
├─────┼────────────┼────────────────────────────────────────────────────┤
│ 1   │ 2020-01-31 │ [117126.0, 155163.0, 2.23714e5]=>[:P25 :P50 :P75]  │
│ 2   │ 2020-02-29 │ [117609.0, 155726.0, 224787.0]=>[:P25 :P50 :P75]   │
│ 3   │ 2020-03-31 │ [1.17992e5, 156410.0, 2.25734e5]=>[:P25 :P50 :P75] │
│ 4   │ 2020-04-30 │ [1.1849e5, 156411.0, 226593.0]=>[:P25 :P50 :P75]   │
=#

# pipe syntax #
# =========== #

@pipe sdf |>
    groupby(_, :ReportDate) |>
    combine(_, :value => mean => :avg)

# pipe syntax is nice because of the incremental stages of transform are easier to read

# syntax flows from inside-out of the function calls

# joins #
# ===== #

country = DataFrame(Name = ["United States"; "Zoldova"; "Aazonia"])
states = DataFrame(StateCode = ["CA"; "TX"; "GA"])

# joining is easy
innerjoin(df, country, on = [:RegionName => :Name])

innerjoin(tdf, states, on = [:State => :StateCode])

@pipe rightjoin(df, country, on = [:RegionName => :Name]) |>
    sort!(_, :RegionName)
#=
3×5 DataFrame
│ Row │ RegionName    │ 2020-01-31 │ 2020-02-29 │ 2020-03-31 │ 2020-04-30 │
│     │ String        │ Float64?   │ Float64?   │ Float64?   │ Float64?   │
├─────┼───────────────┼────────────┼────────────┼────────────┼────────────┤
│ 1   │ Aazonia       │ missing    │ missing    │ missing    │ missing    │
│ 2   │ United States │ 247060.0   │ 248046.0   │ 249140.0   │ 250271.0   │
│ 3   │ Zoldova       │ missing    │ missing    │ missing    │ missing    │
=#

@pipe first(tdf, 6) |>
    antijoin(_, states, on = [:State => :StateCode])

#=
│ Row │ RegionName       │ 2020-01-31 │ 2020-02-29 │ 2020-03-31 │ 2020-04-30 │ State    │
│     │ String           │ Float64    │ Float64    │ Float64    │ Float64    │ Abstrac… │
├─────┼──────────────────┼────────────┼────────────┼────────────┼────────────┼──────────┤
│ 1   │ United States    │ 247060.0   │ 248046.0   │ 249140.0   │ 250271.0   │          │
│ 2   │ New York, NY     │ 485111.0   │ 486070.0   │ 486979.0   │ 488002.0   │ NY       │
│ 3   │ Chicago, IL      │ 242376.0   │ 242743.0   │ 243262.0   │ 243627.0   │ IL       │
│ 4   │ Philadelphia, PA │ 252614.0   │ 253547.0   │ 254232.0   │ 255300.0   │ PA       │
=#

# TX and GA no longer members
# nice!


text:
# functions relating to parametric solutions of a single factor that is likely in need of consolidation

"""
    $SIGNATURES

Helper function to propagate a parametric estimate along a factor chain.
This function takes and returns variable values as coordinates.

Notes
- Not used during MM-iSAM inference.
- Expected uses are for user analysis of factors and estimates.
- real-time dead reckoning chain prediction.
- Parametric binary factor utility function, used by DRT

DevNotes
- TODO ensure type stability, likely returning types `Any` at this time.
- TODO MeanMaxPPE currently stored as coordinates, complicating fast calculation.

Related:

[`getMeasurementParametric`](@ref), [`approxConv`](@ref), [`MutablePose2Pose2Gaussian`](@ref)
"""
function solveFactorParameteric(dfg::AbstractDFG,
                                fct::DFGFactor,
                                # currval::P1,
                                srcsym_vals::AbstractVector{Pair{Symbol, P}},
                                trgsym::Symbol,
                                solveKey::Symbol=:default;
                                evaltmpkw...  ) where P
  #

  varLbls = getVariableOrder(fct)
  varTypes = tuple((getVariableType.(dfg, varLbls))...)
  sfidx = findfirst( varLbls .== trgsym )

  # get the measurement point
  fctTyp = getFactorType(fct)
  # this is definitely in coordinates, see JuliaRobotics/RoME.jl#465
  mea, _ = getMeasurementParametric(fctTyp)
  # must change measT to be a tangent vector
  M = getManifold(fctTyp)
  e0 = identity_element(M)
  mea_ = hat(M, e0, mea)
  measT = [mea_]

  # get variable points
  function _getParametric(vari::DFGVariable, key=:default)
    # hasp = haskey(getPPEDict(vari), key)
    # FIXME use PPE via Manifold points currently in coordinates
    # hasp ? getPPE(vari, key).suggested : calcMean(getBelief(vari, key))
    pt = calcMean(getBelief(vari, key))

    getCoordinates(getVariableType(vari),pt)
  end

  # overwrite specific src values from user
  coordVals = _getParametric.(getVariable.(dfg, varLbls), solveKey)
  for (srcsym, currval) in srcsym_vals
    coordVals[findfirst(varLbls .== srcsym)] = currval
  end
  crds = tuple(coordVals...)

  pts = tuple(map(t->getPoint(t...), zip(varTypes,crds))...)

  # do the calculation to find solvefor index using the factor, as manifold point
  pt = _evalFactorTemporary!( fctTyp, varTypes, sfidx, measT, pts; solveKey=solveKey, evaltmpkw... )[1]

  return getCoordinates(getVariableType(dfg, trgsym), pt)
end


text:
"""
    coef(MLMNet)

Extract all coefficients from Mlmnet object

# Arguments

- MLMNet = Mlmnet object

# Value

3d array of coefficients

"""
function coef(MLMNet::Mlmnet)

    return MLMNet.B
end

"""
    coef(MLMNet, lambda, alpha)

Extract coefficients from Mlmnet object at a given lambda

# Arguments

- MLMNet = Mlmnet object
- lambda = lambda penalty to use, a floating scalar
- alpha = alpha penalty to determine the mix of penalties between L1 and L2
  a floating scalar

# Value

2d array of coefficients

"""
function coef(MLMNet::Mlmnet, lambda::Float64, alpha::Float64)

    # Find the index corresponding to the lambda of interest
    idx = findall([isapprox(lam, lambda) for lam in MLMNet.lambdas])
    idy = findall([isapprox(alp, alpha) for alp in MLMNet.alphas])
    if (length(idx) == 0)
        error("This lambda was not used in the MLMNet object.")
    end

    if (length(idy) == 0)
        error("This alpha was not used in the MLMNet object.")
    end

    # Return the coefficient slice for the lambda of interest
    return MLMNet.B[:,:,idx[1],idy[1]]
end

"""
    coef_3d(MLMNet)

Extract coefficients from Mlmnet object as a flattened 2d array

# Arguments

- MLMNet = Mlmnet object

# Value

2d array of flattened coefficients, where each column corresponds to a
different lambda and alpha

"""
function coef_3d(MLMNet::Mlmnet)

    # Extract coefficients
    coeffs = coef(MLMNet) # this is a 4-dimensional array

    # Initialize 2d array to store flattened coefficients
    coeffs3d = Array{Float64}(undef, size(coeffs,1)*size(coeffs,2),
                              size(coeffs,3), size(coeffs,4))

    # Iterate through each coefficient slice
    for j in 1:size(coeffs, 4), i in size(coeffs, 3)
        coeffs3d[:,i,j] = vec(coeffs[:,:,i,j])
    end

    return coeffs3d
end

"""
    predict(MLMNet, lambda, alpha, newPredictors)

Calculate new predictions based on Mlmnet object and given a lambda

# Arguments

- MLMNet = Mlmnet object
- lambda = lambda penalty to use, a floating scalar
- alpha = alpha penalty to determine the mix of penalties between L1 and L2
  a floating scalar
- newPredictors = Predictors object. Defaults to the data.predictors field
  in the MLM object used to fit the model.

# Value

2d array of predicted values

"""
function predict(MLMNet::Mlmnet, lambda::Float64, alpha::Float64,
                 newPredictors::Predictors=MLMNet.data.predictors)

    # Include X and Z intercepts in new predictors if necessary
    if MLMNet.data.predictors.hasXIntercept==true &&
       newPredictors.hasXIntercept==false
        newPredictors.X = add_intercept(newPredictors.X)
        newPredictors.hasXIntercept = true
        println("Adding X intercept to newPredictors.")
  	end
    if MLMNet.data.predictors.hasZIntercept==true &&
       newPredictors.hasZIntercept==false
        newPredictors.Z = add_intercept(newPredictors.Z)
        newPredictors.hasZIntercept = true
        println("Adding Z intercept to newPredictors.")
  	end

    # Remove X and Z intercepts in new predictors if necessary
    if MLMNet.data.predictors.hasXIntercept==false &&
       newPredictors.hasXIntercept==true
        newPredictors.X = remove_intercept(newPredictors.X)
        newPredictors.hasXIntercept = false
        println("Removing X intercept from newPredictors.")
    end
    if MLMNet.data.predictors.hasZIntercept==false &&
       newPredictors.hasZIntercept==true
        newPredictors.Z = remove_intercept(newPredictors.Z)
        newPredictors.hasZIntercept = false
        println("Removing Z intercept from newPredictors.")
    end

    # Calculate new predictions
    return calc_preds(newPredictors.X, newPredictors.Z, coef(MLMNet, lambda, alpha))
end

"""
    predict(MLMNet, newPredictors)

Calculate new predictions based on Mlmnet object

# Arguments

- MLMNet = Mlmnet object
- newPredictors = Predictors object. Defaults to the data.predictors field
  in the MLM object used to fit the model.

# Value

4d array of predicted values

"""
function predict(MLMNet::Mlmnet,
                 newPredictors::Predictors=MLMNet.data.predictors)

    # Initialize 4d array for storing predictions
    all_preds = Array{Float64}(undef, size(newPredictors.X,1), size(newPredictors.Z,1),
                               length(MLMNet.lambdas), length(MLMNet.alphas))

	# Calculate new predictions for each lambda
    for j = 1:length(MLMNet.alphas), i = 1:length(MLMNet.lambdas)
        all_preds[:,:,i,j] = predict(MLMNet, MLMNet.lambdas[i], MLMNet.alphas[j],
                                        newPredictors)
    end

    return all_preds
end

"""
    fitted(MLMNet, lambda, alpha)

Calculate fitted values of an Mlmnet object, given a lambda

# Arguments

- MLMNet = Mlmnet object
- lambda = lambda penalty to use, a floating scalar
- alpha = alpha penalty to determine the mix of penalties between L1 and L2
  a floating scalar

# Value

2d array of fitted values

"""
function fitted(MLMNet::Mlmnet, lambda::Float64, alpha::Float64)

    return predict(MLMNet, lambda, alpha)
end

"""
    fitted(MLMNet)

Calculate fitted values of an Mlmnet object

# Arguments

- MLMNet = Mlmnet object

# Value

3d array of fitted values

"""
function fitted(MLMNet::Mlmnet)

    return predict(MLMNet)
end

"""
    resid(MLMNet, lambda, alpha, newData)

Calculate residuals of an MLMNet object, given a lambda

# Arguments

- MLMNet = Mlmnet object
- lambda = lambda penalty to use, a floating scalar
- alpha = alpha penalty to determine the mix of penalties between L1 and L2
  a floating scalar
- newData = RawData object. Defaults to the data field in the MLM object
  used to fit the model.

# Value

2d array of residuals

"""
function resid(MLMNet::Mlmnet, lambda::Float64, alpha::Float64, newData::RawData=MLMNet.data)

    # Include X and Z intercepts in new data if necessary
    if MLMNet.data.predictors.hasXIntercept==true &&
       newData.predictors.hasXIntercept==false
        newData.predictors.X = add_intercept(newData.predictors.X)
        newData.predictors.hasXIntercept = true
        newData.p = newData.p + 1
        println("Adding X intercept to newData.")
    end
    if MLMNet.data.predictors.hasZIntercept==true &&
       newData.predictors.hasZIntercept==false
        newData.predictors.Z = add_intercept(newData.predictors.Z)
        newData.predictors.hasZIntercept = true
        newData.q = newData.q + 1
        println("Adding Z intercept to newData.")
    end

    # Remove X and Z intercepts in new data if necessary
    if MLMNet.data.predictors.hasXIntercept==false &&
       newData.predictors.hasXIntercept==true
        newData.predictors.X = remove_intercept(newData.predictors.X)
        newData.predictors.hasXIntercept = false
        newData.p = newData.p - 1
        println("Removing X intercept from newData.")
    end
    if MLMNet.data.predictors.hasZIntercept==false &&
       newData.predictors.hasZIntercept==true
        newData.predictors.Z = remove_intercept(newData.predictors.Z)
        newData.predictors.hasZIntercept = false
        newData.q = newData.q - 1
        println("Removing Z intercept from newData.")
    end

    # Calculate residuals
    return calc_resid(get_X(newData), get_Y(newData), get_Z(newData),
                      coef(MLMNet, lambda, alpha))
end

"""
    resid(MLMNet, newData)

Calculate residuals of an MLMNet object

# Arguments

- MLMNet = Mlmnet object
- newData = RawData object. Defaults to the data field in the MLM object
  used to fit the model.

# Value

3d array of residuals

"""
function resid(MLMNet::Mlmnet, newData::RawData=MLMNet.data)

    # Initialize 4d array for storing residuals
    all_resid = Array{Float64}(undef, newData.n, newData.m,
                                length(MLMNet.lambdas), length(MLMNet.alphas))

    # Calculate residuals for each lambda
    for j = 1:length(MLMNet.alphas), i = 1:length(MLMNet.lambdas)
        all_resid[:,:,i,j] = resid(MLMNet, MLMNet.lambdas[i], MLMNet.alphas[j], newData)
    end

    return all_resid
end


text:
struct ExactFormulation <: AbstractRelaxation end

function reformulate_unary_quadratic_term!(
    model::MOI.ModelLike,
    reformulation::Reformulation{ExactFormulation},
    x::VI,
    y::VI,
    l::Float64,
    u::Float64,
)
    MOI.add_constraint(
        model,
        1.0 * SV(y) - 1.0 * SV(x) * SV(x),
        MOI.EqualTo(0.0),
    )
    return nothing
end


text:
function sameSpan(aVecs::Array{Float64,2},bVecs::Array{Float64,2})
    aRank=rank(aVecs)
    bRank=rank(bVecs)
    bothRank=rank(vcat(aVecs,bVecs),10e-12)
    aRank==bRank==bothRank
end


text:
Base.last(::OneToSRange{T,E}) where {T,E} = E

Base.last(r::OneToMRange) = getfield(r, :stop)

Base.last(::UnitSRange{T,F,L}) where {T,F,L} = L

Base.last(r::UnitMRange) = getfield(r, :stop)

Base.last(r::AbstractStepRangeLen) = unsafe_getindex(r, length(r))

Base.last(::LinSRange{T,B,E,L,D}) where {T,B,E,L,D} = E

Base.last(r::LinMRange) = getfield(r, :stop)

"""
    can_set_last(x) -> Bool

Returns `true` if the last element of `x` can be set. If `x` is a range then
changing the first element will also change the length of `x`.
"""
can_set_last(::T) where {T} = can_set_last(T)
can_set_last(::Type{T}) where {T} = can_setindex(T)
can_set_last(::Type{T}) where {T<:LinMRange} = true
can_set_last(::Type{T}) where {T<:StepMRange} = true
can_set_last(::Type{T}) where {T<:StepMRangeLen} = true
can_set_last(::Type{T}) where {T<:UnitMRange} = true
can_set_last(::Type{T}) where {T<:OneToMRange} = true

"""
    set_last!(x, val)

Set the last element of `x` to `val`.
"""
function set_last!(x::AbstractVector{T}, val::T) where {T}
    can_set_last(x) || throw(MethodError(set_last!, (x, val)))
    setindex!(x, val, lastindex(x))
    return x
end
set_last!(x::AbstractVector{T}, val) where {T} = set_last!(x, convert(T, val))
set_last!(r::LinMRange{T}, val::T) where {T} = (setfield!(r, :stop, val); r)
function set_last!(r::StepMRange{T,S}, val::T) where {T,S}
    setfield!(r, :stop, Base.steprange_last(first(r), step(r), val))
    return r
end
set_last!(r::UnitMRange{T}, val::T) where {T} = (setfield!(r, :stop, val); r)
function set_last!(r::OneToMRange{T}, val::T) where {T}
    setfield!(r, :stop, max(zero(T), T(val)))
    return r
end
function set_last!(r::StepMRangeLen{T}, val::T) where {T}
    len = unsafe_findvalue(val, r)
    len >= 0 || throw(ArgumentError("length cannot be negative, got $len"))
    1 <= r.offset <= max(1, len) || throw(ArgumentError("StepSRangeLen: offset must be in [1,$len], got $(r.offset)"))
    setfield!(r, :len, len)
    return r
end


text:
using JuMP
using Clp
using Cbc

m = Model(Cbc.Optimizer)
@variable(m, fa >= 0, Int)
@variable(m, fb >= 0, Int)
@variable(m, pa >= 0, Int)
@variable(m, pb >= 0, Int)
@variable(m, pc >= 0, Int)
@variable(m, pd >= 0, Int)

@objective(m, Min, 320 * (fa + fb) + 120 * (pa + pb + pc + pd))
@constraint(m, fa + pa >= 4)
@constraint(m, fa + pb >= 8)
@constraint(m, fb + pc >= 10)
@constraint(m, fb + pd >= 8)
@constraint(m, fa >= 2pa)
@constraint(m, fa >= 2pb)
@constraint(m, fb >= 2pc)
@constraint(m, fb >= 2pd)

print(m)
optimize!(m)
status = termination_status(m)
println("Solution status: ", status)

println("Objective value: ", objective_value(m))

println("Full-time consultants for shift 1: ", value(fa))
println("Full-time consultants for shift 2: ", value(fb))
println("Part-time consultants for shift 1: ", value(pa))
println("Part-time consultants for shift 2: ", value(pd))
println("Part-time consultants for shift 3: ", value(pc))
println("Part-time consultants for shift 4: ", value(pd))


text:
# test d2f zeros lookup functions

include("../src/overt_utils.jl")
using Test

@testset "Testing Second Derivative Zero Computation" begin
    # test get_regions_1arg
    # division
    # c/x, c > 0, c<0, interval all pos or all neg
    @test get_regions_1arg(:(5/x), :x, 1, 2) == ([], true)
    @test get_regions_1arg(:(5/x), :x, -2, -1) == ([], false)
    @test get_regions_1arg(:(-5/x), :x, 1, 2) == ([], false)
    @test get_regions_1arg(:(-5/x), :x, -2, -1) == ([], true)
    # exponent
    # c^x
    @test get_regions_1arg(:(5.1^x), :x, -1, 1) == ([], true)
    # x^c
    # fractional exponent
    @test get_regions_1arg(:(x^5.1), :x, 0, 2) == ([], true)
    # odd exponent
    @test get_regions_1arg(:(x^5), :x, 1, 2) == ([], true)
    @test get_regions_1arg(:(x^5), :x, -2, -1) == ([], false)
    @test get_regions_1arg(:(x^5), :x, -1, 2) == ([0], nothing)
    # even exponent
    @test get_regions_1arg(:(x^6), :x, -2, -1) == ([], true)

    # test get_regions_unary
    # sin, at boundary points
    @test all(get_regions_unary(:sin, -π, π)[1] .≈ [-π, 0, π])
    @test get_regions_unary(:sin, -π, π)[2] === nothing
    # at fractional points
    @test get_regions_unary(:sin, -1, 1) == ([0], nothing)
    @test all(get_regions_unary(:sin, 1, 4)[1] .≈ [π])
    @test get_regions_unary(:sin, 1, 4)[2] === nothing
    # cos
    @test get_regions_unary(:cos, -1, 1) == ([], false)
    @test all(get_regions_unary(:cos, -π, π)[1] .≈ [-π/2, π/2])
    @test get_regions_unary(:cos, -π, π)[2] === nothing
    # exp
    @test get_regions_unary(:exp, -1, 1) == ([], true)
    # log
    @test get_regions_unary(:log, .5, 1.5) == ([], false)
    # tanh
    @test get_regions_unary(:tanh, -2, 2) == ([0], nothing)
    @test get_regions_unary(:tanh, -2, -1) == ([], true)
    @test get_regions_unary(:tanh, 1, 2) == ([], false)
end


text:
export hamiltonian_cycle

"""
`hamiltonian_cycle(G::SimpleGraph)` returns a Hamiltonian cycle in the
graph (if one exists) or an empty array (otherwise). This works
reasonably well for small graphs.
"""
function hamiltonian_cycle(G::SimpleGraph)
    if cache_check(G, :hamiltonian_cycle)
        return cache_recall(G, :hamiltonian_cycle)
    end
    T = eltype(G)

    # rule out some simple non-Hamiltonian graphs
    if NV(G) < 3 || minimum(deg(G)) < 2 || !is_connected(G)
        return T[]
    end

    n = NV(G)
    path = Array{T,1}(undef, n)  # we'll put the answer here

    used = Dict{T,Bool}()
    for v in G.V
        used[v] = false
    end

    VV = vlist(G)
    v = VV[1]
    used[v] = true
    path[1] = v

    if ham_extend(G, VV, 1, used, path)
        cache_save(G, :hamiltonian_cycle, path)
        return path
    else
        path = T[]
        cache_save(G, :hamiltonian_cycle, path)
        return path
    end
end

function ham_extend(G::SimpleGraph, VV::Array, idx::Int, used::Dict, path::Array)
    # println(path[1:idx])  # debug
    n = NV(G)
    v = path[idx]

    if idx == n
        return has(G, path[1], path[n])
    end

    for w in G[v]
        if !used[w]
            path[idx+1] = w
            used[w] = true
            if ham_extend(G, VV, idx + 1, used, path)
                return true
            end
            used[w] = false
        end
    end
    return false
end

function ham_check(G, path::Array)
    n = NV(G)
    if length(path) != n
        return false
    end

    for k = 1:n-1
        if !has(G, path[k], path[k+1])
            return false
        end
    end
    return has(G, path[1], path[n])
end


text:
using ONNXNaiveNASflux, NaiveNASflux, Flux
import ONNXNaiveNASflux: ONNX
using Test

@testset "ONNXNaiveNASflux.jl" begin

    @info "Test BaseOnnx"
    @testset "BaseOnnx" begin
        include("baseonnx/readwrite.jl")
    end

    @info "Test Deserialization"
    @testset "Deserialize" begin
        include("deserialize/vertex.jl")
        include("deserialize/constraints.jl")
        include("deserialize/testdata.jl")
        include("deserialize/deserialize.jl")
    end

    @info "Test Serialization"
    @testset "Serialize" begin
        include("serialize/serialize.jl")
    end

    @info "Test validation"
    include("validate.jl")

    @info "Test README examples"
    include("examples.jl")
end


text:
# gRPCController only keeps a debug flag as of now
type gRPCController <: ProtoRpcController
    debug::Bool
end

debug_log(controller::gRPCController, msg) = controller.debug && println(msg)
error_log(controller::gRPCController, msg) = println(STDERR, msg)

# gRPCChannel implementation
type gRPCChannel <: ProtoRpcChannel
    session::HTTPConnection
    stream_id::UInt32

    gRPCChannel(session::HTTPConnection) = new(session, 0)
end
close(channel::gRPCChannel) = close(channel.session)

function to_delimited_message_bytes(msg)
    iob = IOBuffer()
    write(iob, UInt8(0))
    write(iob, hton(UInt32(0)))
    data_len = writeproto(iob, msg)
    seek(iob, 1)
    write(iob, hton(UInt32(data_len)))
    take!(iob)
end

function from_delimited_message_bytes{T}(data, t::T)
    iob = IOBuffer(data)
    compressed = read(iob, UInt8)
    datalen = ntoh(read(iob, UInt32))
    readproto(iob, t)
end

function check_grpc_status(headers)
    ("grpc-status" in keys(headers)) || (return (1, ""))

    grpc_status = parse(Int, headers["grpc-status"])
    (grpc_status == 0) && (return (0, ""))

    grpc_message = get(headers, "grpc-message", "unexpected grpc status: $grpc_status")
    (-1, grpc_message)
end

function check_http_status(headers)
    (":status" in keys(headers)) || (return (1, ""))
    status = parse(Int, headers[":status"])
    (status == 200) && (return (0, ""))
    (-1, "unexpected http status: $status")
end

assert_status(status) = (status[1] == -1) ? error(status[2]) : status
assert_http_status(headers) = assert_status(check_http_status(headers))
assert_grpc_status(headers) = assert_status(check_grpc_status(headers))


text:
importall ReinforcementLearning
using Knet, ArgParse

type ReinforceAgent <: AbsAgent
	weights
	prms
	states
	actions
	rewards
	γ
end

function ReinforceAgent(env::GymEnv; α=0.001, γ=0.99)
	w = Dict()
	D = length(getInitialState(env).data)
	O = length(env.actions)

	w["w"] = convert(Array{Float32}, randn(O, D) / sqrt(D))#xavier
	w["b"] = zeros(Float32, O, 1)

	prms = Dict()
	prms["w"] = Adam(w["w"]; lr=α)
	prms["b"] = Adam(w["b"]; lr=α)

	ReinforceAgent(w, prms; γ=γ)
end

ReinforceAgent(w, prms; γ=0.99) = ReinforceAgent(w, prms, Any[], Any[], Array{Float64, 1}(), γ)

function predict(w, b, x)
	linear = w * x .+ b
end

function sample_action(linear)
	probs = exp(linear) ./ sum(exp(linear), 1)
	c_probs = cumsum(probs)
	return indmax(c_probs .> rand())
end

function loss(w, x, actions)
	linear = predict(w["w"], w["b"], x)
	ynorm = logp(linear,1)
	-sum(actions .* ynorm) / size(actions, 2)
end

lossgradient = grad(loss)

function discount(rewards; γ=0.9)
	discounted = zeros(Float64, length(rewards), 1)
	discounted[end] = rewards[end]

	for i=(length(rewards)-1):-1:1
		discounted[i] = rewards[i] + γ * discounted[i+1]
	end
	return discounted
end

function play(agent::ReinforceAgent, state::GymState, env::GymEnv; learn=false)

	linear = predict(agent.weights["w"], agent.weights["b"], state.data)
	a = sample_action(linear)

	if learn
		push!(agent.states, convert(Array{Float32}, state.data))
		y = zeros(Float32, size(linear, 1), 1)
		y[a] = 1.0
		push!(agent.actions, y)
	end
	return env.actions[a]
end

function observe(agent::ReinforceAgent, state::AbsState, reward::Float64, env::AbsEnvironment; learn=true, terminal=false)
	if learn
		push!(agent.rewards, reward)
		if terminal
			if length(agent.rewards) < 195
				agent.rewards[end] = -100;
			end

			disc_rewards = discount(agent.rewards; γ=agent.γ)
			for i=1:length(agent.states)
				g = lossgradient(agent.weights, agent.states[i], agent.actions[i])
				for k in keys(g)
					update!(agent.weights[k], disc_rewards[i] * agent.γ * g[k], agent.prms[k])
				end
			end

			agent.states = Any[]
			agent.actions = Any[]
			agent.rewards = Array{Float64, 1}()
		end
	end
end

function main(args=ARGS)
	s = ArgParseSettings()
	s.description="Solution with REINFORCE algorithm"
	s.exc_handler=ArgParse.debug_handler
	@add_arg_table s begin
		("--lr"; arg_type=Float64; default=0.01; help="learning rate")
		("--gamma"; arg_type=Float64; default=0.9; help="discount rate")
		("--env"; default="CartPole-v0"; help="name of the environment")
		("--render"; action=:store_true; help="display the env")
		("--monitor"; default=""; help="path of the log for the experiment, empty if mointoring is disabled")
		("--epoch"; arg_type=Int; default=300; help="number of epochs")
		("--threshold"; arg_type=Int; default=1000; help="stop the episode even it is not terminal after number of steps exceeds the threshold")
	end
	srand(123)
	o = parse_args(args, s)
	env = GymEnv(o["env"])
	agent = ReinforceAgent(env; α=o["lr"], γ=o["gamma"])
	rewards = Array{Float64, 1}()

	o["monitor"] != "" && monitor_start(env, o["monitor"])
	for i=1:o["epoch"]
		totalRewards, numberOfStates = playEpisode(env, agent; learn=true, threshold = o["threshold"], render=o["render"])
		push!(rewards, totalRewards)
		msg = string("Episode ", i, " , total rewards: ", totalRewards)
		if i >= 100
			m = mean(rewards[(i-100+1):end])
			msg = string(msg, " , mean: ", m)
			if m >= 195
				break
			end
		end
		println(msg)
	end
	o["monitor"] != "" && monitor_close(env)
end

main()

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
import NMRDataSetup
import NMRSpectraSimulator

include("../src/NMRSpecifyRegions.jl")
import .NMRSpecifyRegions
#import NMRSpecifyRegions

using LinearAlgebra
using FFTW
import PyPlot
import BSON
import Statistics

PyPlot.close("all")
fig_num = 1

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

### user inputs.
projects_dir = "/home/roy/MEGAsync/outputs/NMR/calibrate/"

# project_name = "glucose-700"
# molecule_names = ["D-(+)-Glucose"; "DSS"]
# w = [20.0/4.6; 1.0] # BMRB-700 glucose: DSS is 0.0046 M = 4.6 mM.

project_name = "phenylalanine-700"
molecule_names = ["L-Phenylalanine"; "DSS"]
w = [20/0.5; 1.0] # BMRB-700 phenylalanine: DSS is 500 micro M.

# path to the GISSMO Julia storage folder.
base_path_JLD = "/home/roy/Documents/data/NMR/NMRData/src/input/molecules"

# proxy-related.
tol_coherence = 1e-2
α_relative_threshold = 0.05
λ0 = 3.4
Δcs_max = 0.2
κ_λ_lb = 0.5
κ_λ_ub = 2.5

### end inputs.

## load data.
load_folder_path = joinpath(projects_dir, project_name)
load_path = joinpath(load_folder_path, "$(project_name).bson")
dic = BSON.load(load_path)
s_t = dic[:s_t]
fs = dic[:fs]
SW = dic[:SW]
ν_0ppm = dic[:ν_0ppm]

# normalize.
s_t = s_t

hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)

offset_Hz = ν_0ppm - (ppm2hzfunc(0.3)-ppm2hzfunc(0.0))

N = length(s_t)
DFT_s = fft(s_t)
U_DFT, U_y, U_inds = NMRDataSetup.getwraparoundDFTfreqs(N, fs, offset_Hz)

Z = maximum(abs.(DFT_s))
y = DFT_s[U_inds] ./ Z

#
S_U = DFT_s[U_inds]
P_y = hz2ppmfunc.(U_y)

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P_y, real.(y), label = "data spectrum, y")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("data spectra, real part")

####### mixture proxy.

Δcs_max_mixture = collect( Δcs_max for i = 1:length(molecule_names))

dummy_SSFID = NMRSpectraSimulator.SpinSysFIDType1(0.0)
mixture_params = NMRSpectraSimulator.setupmixtureproxies(molecule_names,
    base_path_JLD, Δcs_max_mixture, hz2ppmfunc, ppm2hzfunc, fs, SW, λ0,
    ν_0ppm, dummy_SSFID;
    tol_coherence = tol_coherence,
    α_relative_threshold = α_relative_threshold)
As = mixture_params

ΩS_ppm = NMRSpectraSimulator.getPsnospininfo(As, hz2ppmfunc)
ΩS_ppm_sorted = sort(NMRSpectraSimulator.combinevectors(ΩS_ppm))

# u_min = ppm2hzfunc(-0.5)
# u_max = ppm2hzfunc(3.0)

u_offset = 0.5
u_min = ppm2hzfunc(ΩS_ppm_sorted[1] - u_offset)
u_max = ppm2hzfunc(ΩS_ppm_sorted[end] + u_offset)

NMRSpectraSimulator.fitproxies!(As;
    κ_λ_lb = κ_λ_lb,
    κ_λ_ub = κ_λ_ub,
    u_min = u_min,
    u_max = u_max,
    Δr = 1.0,
    Δκ_λ = 0.05)

### plot.

P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
U = ppm2hzfunc.(P)
#ΩS_ppm = collect( hz2ppmfunc.( NMRSpectraSimulator.combinevectors(A.Ωs) ./ (2*π) ) for A in mixture_params )

cs_config_path = "/home/roy/MEGAsync/inputs/NMR/configs/cs_config_reduced.txt"

cs_delta_group = NMRSpecifyRegions.extractinfofromconfig( cs_config_path, molecule_names)
Δsys_cs = NMRSpecifyRegions.condenseΔcsconfig(cs_delta_group)

ΩS0 = NMRSpectraSimulator.getΩS(As)
ΩS0_ppm = NMRSpectraSimulator.getPs(ΩS0, hz2ppmfunc)
exp_info = NMRSpecifyRegions.setupexperimentresults(molecule_names, ΩS0_ppm, Δsys_cs;
min_dist = 0.1)

U_cost0 = U_y
P_cost0 = hz2ppmfunc.(U_cost0)
y_cost0 = y

band_inds, band_inds_set = NMRSpecifyRegions.getcostinds(exp_info, P_cost0)

U_cost = U_cost0[band_inds]
P_cost = P_cost0[band_inds]

y_cost = y_cost0[band_inds]

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P_y, real.(y), label = "data spectrum")
PyPlot.plot(P_cost, real.(y_cost), "^", label = "positions")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("positions against data spectrum, real part")


text:
module GaramondTesting

using Test
using Logging
using Random
using SparseArrays
using JuliaDB
using EmbeddingsAnalysis
using Garamond
#using JSON

# TODO(Corneliu): Add more tests
test_logger = ConsoleLogger(stdout, Logging.Error)  # supresses warnings, infos
with_logger(test_logger) do
    # Necessary data/config generators
    include("data/datagenerator.jl")
    include("configs/configgenerator.jl")
    # Tests
    include("input_parsers.jl")
    include("db.jl")
    include("index.jl")
    include("embedder.jl")
    include("searcher.jl")
    include("env.jl")
    include("indexfilter.jl")
    include("config_parser.jl")
end

end


text:
using FastRationals, LinearAlgebra
using Test

include("base_fast.jl")
include("more_fast.jl")
include("compactify.jl")
include("round_etc.jl")


text:
using Symbolics, SparseArrays, LinearAlgebra, BenchmarkTools

nfact = 10
nitem = 5

nobs = nitem*nfact
ntot = nobs+nfact

npar = Int64(nobs + (nfact+1)*nfact/2 + nobs - nfact)

@Symbolics.variables x[1:npar]

x = [x[i] for i in 1:size(x, 1)]
#cov = [cov[i,j] for i in 1:size(cov, 1) for j in 1:size(cov, 2)]
S = zeros(Num, ntot, ntot)
S[diagind(S)[1:nobs]] .= x[1:nobs]
xind = nobs
for i in 1:nfact
    for j in i:nfact
        xind = xind+1
        ind = nobs+i
        jnd = nobs+j
        S[ind, jnd] = x[xind]
        S[jnd, ind] = x[xind]
    end
end

F = zeros(Num, nobs, ntot)

F[diagind(F)] .= 1

A = zeros(Num, ntot, ntot)

for i in 1:nfact
    for j in 2:nitem
        xind = xind+1
        jnd = i+nobs
        ind = j+(i-1)*nitem
        A[ind, jnd] = x[xind]
    end
end

for i in 1:nfact
    A[1+(i-1)*nitem, nobs+i] = 1
end

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

par = [x[i] for i = 1:npar]
invia = I + A + A^2 + A^3 + A^4
imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)
imp_cov_sym = Array(imp_cov_sym)
imp_cov_sym = Symbolics.simplify.(imp_cov_sym)

imp_fun = eval(Symbolics.build_function(imp_cov_sym, par)[2])

Σ = zeros(nobs, nobs)
randpar = rand(npar)

imp_fun(Σ, randpar)

Σ

using BenchmarkTools

@benchmark imp_fun(Σ, randpar)

########## Benchmark Gradients
using sem, Optim

obs_cov = rand(nobs, nobs); obs_cov = obs_cov*obs_cov'
semobserved = SemObsCommon(data = obs_cov)
start_val_ml = randpar
loss_ml = Loss([SemML(semobserved, [0.0], similar(start_val_ml))])

# imply
imply_ml = ImplySymbolic(A, S, F, x, startval)

grad_ml = sem.∇SemML_2()

diff_ana_ml =
    SemAnalyticDiff(
        LBFGS(),# P = one(rand(31, 31))),
        Optim.Options(
            ;f_tol = 1e-10,
            x_tol = 1.5e-8),
            (grad_ml,))

model_ml_ana = Sem(semobserved, imply_ml, loss_ml, diff_ana_ml)

solution_ana_ml = sem_fit(model_ml_ana)

startval = ones(size(randpar))

model_ml_ana(startval)

gradient = similar(randpar)
model_ml_ana(startval, gradient)

using FiniteDiff
gradient ≈ FiniteDiff.finite_difference_gradient(model_ml_ana, startval)

isposdef(Symmetric(model_ml_ana.imply.imp_cov))

@benchmark model_ml_ana(startval, gradient)

@benchmark model_ml_ana(startval)


text:
using TimeZones.WindowsTimeZoneIDs

xml_file = TimeZones.WindowsTimeZoneIDs.WINDOWS_XML_FILE
!isfile(xml_file) && error("Missing required XML file. Run Pkg.build(\"TimeZones\").")

trans = TimeZones.WindowsTimeZoneIDs.compile(xml_file)
@test trans["Central European Standard Time"] == "Europe/Warsaw"

mktempdir() do temp_dir
    xml_file = joinpath(temp_dir, "windowZones.xml")
    @test !isfile(xml_file)

    empty!(TimeZones.WindowsTimeZoneIDs.WINDOWS_TRANSLATION)
    @test isempty(TimeZones.WindowsTimeZoneIDs.WINDOWS_TRANSLATION)

    # Does not perform download
    TimeZones.WindowsTimeZoneIDs.build(xml_file)
    @test isfile(xml_file)
    @test !isempty(TimeZones.WindowsTimeZoneIDs.WINDOWS_TRANSLATION)
end


text:
## ABSTRACT TYPE

abstract type StoppingCriterion end

## FALL BACK METHODS

# initialization call is either:
update(::StoppingCriterion, loss) = nothing # state

# ... or:
update_training(::StoppingCriterion, loss) = nothing # state

# subsequent updating:
update(::StoppingCriterion, loss, state) = state
update_training(::StoppingCriterion, loss, state) = state

# returns whether it's time to stop:
done(::StoppingCriterion, state) = false

message(criterion::StoppingCriterion, state) = "Early stop triggered by "*
    "$criterion stopping criterion. "

# to be overloaded if the criterion requires one or more
# `update_training` calls per `update` call, ie, it requires both
# "training" and "out-of-sample" losses to work:
needs_in_and_out_of_sample(::Type) = false

needs_in_and_out_of_sample(c::StoppingCriterion) =
    needs_in_and_out_of_sample(typeof(c))


text:
# Autogenerated wrapper script for Qt6Base_jll for aarch64-linux-musl-cxx11
export libqt6concurrent, libqt6core, libqt6dbus, libqt6gui, libqt6network, libqt6opengl, libqt6printsupport, libqt6sql, libqt6test, libqt6widgets, libqt6xml

using Xorg_libXext_jll
using Xorg_libxcb_jll
using Xorg_xcb_util_wm_jll
using Xorg_xcb_util_image_jll
using Xorg_xcb_util_keysyms_jll
using Xorg_xcb_util_renderutil_jll
using Xorg_libXrender_jll
using xkbcommon_jll
using Libglvnd_jll
using Fontconfig_jll
using Glib_jll
using Zlib_jll
using CompilerSupportLibraries_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Qt6Base")
JLLWrappers.@declare_library_product(libqt6concurrent, "libQt6Concurrent.so.6")
JLLWrappers.@declare_library_product(libqt6core, "libQt6Core.so.6")
JLLWrappers.@declare_library_product(libqt6dbus, "libQt6DBus.so.6")
JLLWrappers.@declare_library_product(libqt6gui, "libQt6Gui.so.6")
JLLWrappers.@declare_library_product(libqt6network, "libQt6Network.so.6")
JLLWrappers.@declare_library_product(libqt6opengl, "libQt6OpenGL.so.6")
JLLWrappers.@declare_library_product(libqt6printsupport, "libQt6PrintSupport.so.6")
JLLWrappers.@declare_library_product(libqt6sql, "libQt6Sql.so.6")
JLLWrappers.@declare_library_product(libqt6test, "libQt6Test.so.6")
JLLWrappers.@declare_library_product(libqt6widgets, "libQt6Widgets.so.6")
JLLWrappers.@declare_library_product(libqt6xml, "libQt6Xml.so.6")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libXext_jll, Xorg_libxcb_jll, Xorg_xcb_util_wm_jll, Xorg_xcb_util_image_jll, Xorg_xcb_util_keysyms_jll, Xorg_xcb_util_renderutil_jll, Xorg_libXrender_jll, xkbcommon_jll, Libglvnd_jll, Fontconfig_jll, Glib_jll, Zlib_jll, CompilerSupportLibraries_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libqt6concurrent,
        "lib/libQt6Concurrent.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6core,
        "lib/libQt6Core.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6dbus,
        "lib/libQt6DBus.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6gui,
        "lib/libQt6Gui.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6network,
        "lib/libQt6Network.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6opengl,
        "lib/libQt6OpenGL.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6printsupport,
        "lib/libQt6PrintSupport.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6sql,
        "lib/libQt6Sql.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6test,
        "lib/libQt6Test.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6widgets,
        "lib/libQt6Widgets.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6xml,
        "lib/libQt6Xml.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## ScalarAffineFunction-in-Set

function _build_constraint(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
    s::MOI.GreaterThan{T},
) where {T <: Real}
    coeffs, vars, constant = _parse_to_coeffs_vars(model, f)
    return LinearInt(model.inner, coeffs, vars, ">=", s.constant - constant)
end

function _build_constraint(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
    s::MOI.LessThan{T},
) where {T <: Real}
    coeffs, vars, constant = _parse_to_coeffs_vars(model, f)
    return LinearInt(model.inner, coeffs, vars, "<=", s.constant - constant)
end

function _build_constraint(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
    s::MOI.EqualTo{T},
) where {T <: Real}
    coeffs, vars, constant = _parse_to_coeffs_vars(model, f)
    return LinearInt(model.inner, coeffs, vars, "==", s.constant - constant)
end

# No vector of constraints, there is no more efficient way to do it.
# No intervals.
# No constraint deletion.

## VectorOfVariables-in-SOS{I|II}
# Not available. Bridge it?

## VectorOfVariables-in-SecondOrderCone
# Not available.
# https://github.com/radsz/jacop/blob/9c206884e4501c5278ce588233bd171e107bd3c6/src/main/java/org/jacop/fz/constraints/ConstraintFncs.java#L359


text:
using BitTwiddlingConvenienceFunctions: intlog2, nextpow2, prevpow2
using Test

@testset "BitTwiddlingConvenienceFunctions.jl" begin

  @test all(i ->  intlog2(1 << i) == i, 0:(Int == Int64 ? 53 : 30))
  @test all(i -> nextpow2(i) == i, 0:2)
  for j in 1:10
    l, u = (1<<j)+1, 1<<(j+1)
    @test all(i -> nextpow2(i) == u, l:u)
  end

end


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
## Obtain disturbance or path disurb matrices from modif, rIndex and skeleton
## disturb matrix if user not comfortable setting it themselves
################################################################################

function disturb_mat(W::Array{Float64, 2}, rIndex::Array{Int64, 2},
                     modif::Array{Float64, 2})

    ## Get sizes
    snum = size(W)[1]
    dnum = size(modif)[2]
    rnum = size(modif)[1]
    result = zeros(Array{Float64}(snum, snum, dnum))

    ## Construct distubances
    for k = 1:dnum
        result[:, :, k] = W
        for i = 1:rnum
            result[rIndex[i, 1], rIndex[i, 2], k] = modif[i, k]
        end
    end

    return result

end

function path_mat(W::Array{Float64, 2}, rIndex::Array{Int64, 2},
                  modif::Array{Float64, 3})

    ## Get sizes
    snum = size(W)[1] # state
    pnum = size(modif)[2] # path
    rnum = size(modif)[1]
    tnum = size(modif)[3] # time
    result = zeros(Array{Float64}(snum, snum, pnum, tnum))

    ## Construct distubances
    for t = 1:tnum
        for k = 1:pnum
            result[:, :, k, t] = W
            for i = 1:rnum
                result[rIndex[i, 1], rIndex[i, 2], k, t] = modif[i, k, t]
            end
        end
    end

    return result

end

function subsim_mat(W::Array{Float64, 2}, rIndex::Array{Int64, 2},
                    modif::Array{Float64, 4})

    ## Get sizes
    snum = size(W)[1]
    pnum = size(modif)[2]
    rnum = size(modif)[1]
    tnum = size(modif)[4]
    subsimNum = size(modif)[3] # nested simulation
    result = zeros(Array{Float64}(snum, snum, pnum, subsimNum, tnum))

    ## Construct distubances
    for t = 1:tnum
        for k = 1:pnum
            for s = 1:subsimNum
                result[:, :, k, s, t] = W
                for i = 1:rnum
                    result[rIndex[i, 1], rIndex[i, 2], k, s, t] = modif[i, k, s, t]
                end
            end
        end
    end

    return result

end