text:
# Anderson-Darling test

export OneSampleADTest, KSampleADTest

abstract type ADTest <: HypothesisTest end

## ONE SAMPLE AD-TEST
### http://www.itl.nist.gov/div898/handbook/eda/section3/eda35e.htm

function adstats(x::AbstractVector{T}, d::UnivariateDistribution) where T<:Real
    n = length(x)
    A² = convert(Float64, -n)
    for i = 1:n
        lcdfz = logcdf(d, x[i])
        lccdfz = logccdf(d, x[n - i + 1])
        A² -= (i+i - 1)/n * (lcdfz + lccdfz)
    end
    return A²
end

struct OneSampleADTest <: ADTest
    n::Int      # number of observations
    μ::Float64  # sample mean
    σ::Float64  # sample std
    A²::Float64 # Anderson-Darling test statistic
end

"""
    OneSampleADTest(x::AbstractVector{<:Real}, d::UnivariateDistribution)

Perform a one-sample Anderson–Darling test of the null hypothesis that the data in vector
`x` come from the distribution `d` against the alternative hypothesis that the sample
is not drawn from `d`.

Implements: [`pvalue`](@ref)
"""
function OneSampleADTest(x::AbstractVector{T}, d::UnivariateDistribution) where T<:Real
    n = length(x)
    μ, σ = mean_and_std(x)
    y = sort(x)
    if isa(d, Uniform)
        m = y[1]
        r = y[end]-m
        broadcast!(x->(x-m)/r, y, y)
        # to avoid -Inf when calculated logccdf(d, 1.0) or logcdf(d, 0.0) in the test.
        y[1] += eps()
        y[end] -= eps()
    elseif isa(d, Exponential)
        broadcast!(x->x/μ, y, y)
    else
        zscore!(y, μ, σ)
    end

    OneSampleADTest(n, μ, σ, adstats(y, d))
end

testname(::OneSampleADTest) = "One sample Anderson-Darling test"
default_tail(test::OneSampleADTest) = :right

function show_params(io::IO, x::OneSampleADTest, ident = "")
    println(io, ident, "number of observations:   $(x.n)")
    println(io, ident, "sample mean:              $(x.μ)")
    println(io, ident, "sample SD:                $(x.σ)")
    println(io, ident, "A² statistic:             $(x.A²)")
end

### G. and J. Marsaglia, "Evaluating the Anderson-Darling Distribution", Journal of Statistical Software, 2004
function pvalue(t::OneSampleADTest)
    g1(x) = sqrt(x)*(1.0-x)*(49.0x-102.0)
    g2(x) = -0.00022633 + (6.54034 - (14.6538 - (14.458 - (8.259 - 1.91864x)x)x)x)x
    g3(x) = -130.2137 + (745.2337 - (1705.091 - (1950.646 - (1116.360 - 255.7844x)x)x)x)x

    n = t.n
    z = t.A²
    y = if z < 2.0
        exp(-1.2337141/z)*(2.00012+(0.247105-(.0649821-(.0347962-(.0116720-.00168691*z)*z)*z)*z)*z)/sqrt(z)
    else
        exp(-exp(1.0776-(2.30695-(.43424-(.082433-(.008056-.0003146*z)*z)*z)*z)*z))
    end

    pv = y
    if y > 0.8
        pv += g3(y)/n
    else
        c = 0.01265 + 0.1757/n
        if y < c
            pv += (0.0037/n^3 + 0.00078/n^2 + 0.00006/n)*g1(y/c)
        else
            pv += (0.04213/n + 0.01365/n^2)*g2((y - c)/(0.8 - c))
        end
    end
    return 1.0 - pv
end

## K-SAMPLE ANDERSON DARLING TEST
struct KSampleADTest{T<:Real} <: ADTest
    k::Int             # number of samples
    n::Int             # number of observations
    σ::Float64         # variance A²k
    A²k::Float64       # Anderson-Darling test statistic
    modified::Bool     # whether the modified statistic is used
    nsim::Int          # number of simulations for P-value calculation (0 - for asymptotic calculation)
    samples::Vector{T} # pooled samples
    sizes::Vector{Int} # sizes of samples
end

"""
    KSampleADTest(xs::AbstractVector{<:Real}...; modified = true, nsim = 0)

Perform a ``k``-sample Anderson–Darling test of the null hypothesis that the data in the
``k`` vectors `xs` come from the same distribution against the alternative hypothesis that
the samples come from different distributions.

`modified` parameter enables a modified test calculation for samples whose observations
do not all coincide.

If `nsim` is equal to 0 (the default) the asymptotic calculation of p-value is used.
If it is greater than 0, an estimation of p-values is used by generating `nsim` random splits
of the pooled data on ``k`` samples, evaluating the AD statistics for each split, and
computing the proportion of simulated values which are greater or equal to observed.
This proportion is reported as p-value estimate.

Implements: [`pvalue`](@ref)

# References

  * F. W. Scholz and M. A. Stephens, K-Sample Anderson-Darling Tests, Journal of the
    American Statistical Association, Vol. 82, No. 399. (Sep., 1987), pp. 918-924.
"""
KSampleADTest(xs::AbstractVector{T}...; modified = true, nsim = 0) where T<:Real =
    a2_ksample(xs, modified, nsim)

testname(::KSampleADTest) = "k-sample Anderson-Darling test"
default_tail(test::KSampleADTest) = :right

function show_params(io::IO, x::KSampleADTest, ident = "")
    println(io, ident, "number of samples:        $(x.k)")
    println(io, ident, "number of observations:   $(x.n)")
    println(io, ident, "SD of A²k:                $(x.σ)")
    println(io, ident, "A²k statistic:            $(x.A²k)")
    println(io, ident, "standardized statistic:   $((x.A²k - x.k + 1) / x.σ)")
    println(io, ident, "modified test:            $(x.modified)")
    println(io, ident, "p-value calculation:      $(x.nsim == 0 ? "asymptotic" : "simulation" )")
    x.nsim != 0 && println(io, ident, "number of simulations:    $(x.nsim)")
end

"""Monte-Carlo simulation of the p-value for AD test"""
function pvaluesim(x::KSampleADTest)
    Z = sort(x.samples)
    Z⁺ = unique(Z)

    cn = cumsum(x.sizes)
    insert!(cn, 1, 0.0)
    Xr = [cn[i]+1:cn[i+1] for i in 1:x.k]
    idxs = collect(1:x.n)
    IV = [view(idxs, Xr[i]) for i in 1:x.k]
    Xr = [view(x.samples, IV[i]) for i in 1:x.k]

    pv = 0
    for j in 1:x.nsim
        shuffle!(idxs)
        A²k, A²km = adkvals(Z⁺, x.n, Xr)
        adv = x.modified ? A²km : A²k
        adv >= x.A²k && (pv += 1)
    end
    return pv/x.nsim
end

const M   = [1,2,3,4,6,8,10,20]
const PV  = [.00001,.00005,.0001,.0005,.001,.005,.01,.025,.05,.075,.1,.2,.3,.4,.5,.6,.7,.8,.9,.925,.95,.975,.99,.9925,.995,.9975,.999,.99925,.9995,.99975,.9999,.999925,.99995,.999975,.99999]
const TKM = [-1.1923472556232377 -1.1779755032374564 -1.1689824368072272 -1.1424999536072922 -1.1274927776533619 -1.079431992286445 -1.0504873030014337 -0.9989259798052196 -0.9425379414511752 -0.8984863629243318 -0.8603859837928499 -0.7262450719682598 -0.595805134143589 -0.4564972961242285 -0.2961017594430343 -0.10042062839743962 0.15779986189933617 0.5377525370208343 1.2270085983886296 1.5270251901768497 1.9597314709824634 2.7268894645107142 3.7805477044406355 4.121285194909089 4.59629487851025 5.4144038948723905 6.476499981200343 6.818209982829944 7.288086696921435 8.091812033522697 9.184790239278916 9.59709802719259 10.088416692202337 10.914846562007554 12.040789367454868;
             -1.5849548989018964 -1.5407036306002324 -1.5217650056531025 -1.468313040457585 -1.4383460114029611 -1.3510451980424227 -1.2994715185279744 -1.2107011085917234 -1.1182349910369536 -1.0480894352574122 -0.9888467602015968 -0.792282309179151 -0.6158884225411821 -0.4357973228266815 -0.24055211468641896 -0.013263022467496355 0.2679318847085987 0.6545019598726984 1.3035350042516591 1.5717238802812676 1.9472543460769594 2.5917105532062927 3.4417697301287804 3.708500478251771 4.08514071675675 4.729002778898586 5.553060449998679 5.8246988123153525 6.2044102945806 6.8092648417824355 7.564819358105057 7.8751443348359444 8.265612297736865 9.09572889834988 9.634634570153578;
             -1.8179842473942176 -1.7701061384541197 -1.7430115382228353 -1.6664153177996697 -1.6277529354269717 -1.509882048484231 -1.4423296156847323 -1.327021532216986 -1.2105162624023968 -1.1240541192874598 -1.0516241510661641 -0.8188455856721732 -0.6168586530804309 -0.4186963520248099 -0.20789666676568955 0.0295008713522725 0.31671631783726234 0.6992815580015851 1.3223749183602587 1.5732385052700641 1.9236845463814185 2.5085099219292766 3.2689645805547753 3.5044925007442225 3.841724525500276 4.399781356888528 5.114888256216994 5.346077131581993 5.659256532458322 6.257314173278837 6.986207519686372 7.219730636115486 7.592515097944003 8.120907844495369 8.723051996739617;
             -1.9989354254437397 -1.9372512947973601 -1.9085397775154465 -1.8141114782122008 -1.7646890249660294 -1.6175716157784974 -1.5372145799225008 -1.4011822939743703 -1.2672404489587346 -1.169926310440552 -1.0887002685066605 -0.8322081015222909 -0.6140517112328044 -0.404166231020204 -0.18335516739183028 0.06015392460374823 0.34855873939459553 0.7267832591857595 1.3282605791296271 1.5672982936172346 1.8987053911907181 2.4483775155779135 3.157299591873491 3.3747262669186058 3.6797595381737613 4.192286238183518 4.869623799822088 5.080242748754673 5.379051520938655 5.8747370994257535 6.5004003034881785 6.697073027583033 7.03600074408326 7.445733876827231 8.020320100696683;
             -2.260406253565197 -2.1753074859488235 -2.128991194356946 -2.0065848989178288 -1.945510672007325 -1.7624757482456306 -1.6635124607248777 -1.4990052836398609 -1.3392374211352471 -1.2252207360751253 -1.1319971726969635 -0.8445584556230463 -0.6076273028769662 -0.3834824380209529 -0.15483904255918668 0.09368829578808609 0.3831767280474241 0.7547390312514602 1.3312274974190232 1.5571379748618364 1.865766664385659 2.3730178158029442 3.0080178359348353 3.1978501989364085 3.472363926096708 3.931316296612841 4.536760955793221 4.724863565260161 4.962319348947046 5.408613188229584 6.007714186756578 6.198538388520974 6.457796867859084 6.849495066215343 7.4640951217005576;
             -2.440876951075476 -2.341332708197796 -2.2804325041192546 -2.1402371332779366 -2.06667689211804 -1.8548276337573943 -1.7393591986264094 -1.5563125073892916 -1.3811208567928928 -1.257037173358548 -1.1566367065963437 -0.8485596103567375 -0.5998647718737802 -0.36811584852000284 -0.13474953508678567 0.11602561466210987 0.40573641653024317 0.7727862506626891 1.3335441037600442 1.5496136443763067 1.8418545847979217 2.317705414653881 2.917031508317842 3.0961555736231836 3.350246199357208 3.7727674902587927 4.339241095752881 4.516143795536152 4.753763012043537 5.116810577647693 5.672738471706687 5.83139110084699 6.050071129419877 6.401072284239558 6.742573268315121;
             -2.551081192060893 -2.4547632386210463 -2.3949899123424965 -2.2354277972087395 -2.1584171777926153 -1.9275339742517925 -1.800538561713763 -1.5999251384582587 -1.4106960313169241 -1.2786101238180516 -1.171674441939073 -0.8505458546842916 -0.5938925384383626 -0.35689852586669696 -0.11989003705869453 0.13208631237916915 0.42063871492979277 0.7821012141619668 1.3302949325805677 1.5403968280367855 1.8241731782273896 2.2856162625735688 2.8570374013416835 3.0361153244189336 3.274335630236937 3.6852545534103474 4.195577479309726 4.35310621764469 4.569614918885558 4.92866257900564 5.385051459236165 5.55218650057524 5.765620478139155 6.140990692861778 6.68156168217425;
             -2.943481865323865 -2.793682011666347 -2.7120308613215105 -2.506772718605168 -2.39851440603294 -2.1014812457466228 -1.9452101045743528 -1.7041554289774754 -1.4829780613698689 -1.331579442365025 -1.2111003274077823 -0.8544988934261977 -0.5790017475022499 -0.3302088456189407 -0.08744838035713975 0.1673635508237027 0.45239685235540256 0.8039352789696216 1.3240150219349214 1.5194206287846455 1.782378056603685 2.1987886405752803 2.709552813734848 2.8617281816117868 3.0781111816058195 3.431554501112554 3.870563905935619 4.00537363494248 4.1845928150957805 4.503766651874711 4.929971363503398 5.076322471440581 5.27497591999295 5.52791651167518 5.848654461679204]
const TK  = [-1.1926243990804033 -1.1743913869790086 -1.165614653421677 -1.1392512344379109 -1.1246309930698315 -1.0778245813819973 -1.0492276925221506 -0.997028564406416 -0.9412882492496206 -0.89760798562343 -0.8592489056677113 -0.7249815774273046 -0.5958882137218224 -0.4570311457320469 -0.2966615591256979 -0.10064577139947004 0.1579123553929535 0.5374654181056042 1.226650205162352 1.525729224090128 1.9604816602713853 2.7300669502101567 3.7819466008446807 4.115630048318622 4.587277241771445 5.410120147245256 6.504734959614865 6.843435861344368 7.312351411407612 8.204052812351167 9.380908370226425 9.801397342146213 10.210787470021238 11.065700560614296 12.19099994622706;
             -1.5715146874157493 -1.5403073115879633 -1.5169776509733697 -1.4665221836828817 -1.4369102095850903 -1.3496290747813668 -1.298278486773026 -1.2090663275299973 -1.1171538249300548 -1.0478551970728252 -0.9887300447097083 -0.7930787562580826 -0.6164000565514812 -0.4373090523882904 -0.24235866195256756 -0.015768955906757698 0.266294430295106 0.651647823017675 1.3029824265648888 1.57130609192188 1.9479671937433227 2.5878519082111944 3.4484532844903657 3.7110094714746533 4.0901675267103235 4.742903325658023 5.610841153635776 5.869347075916751 6.233093442455383 6.861159495163239 7.605983240042428 7.922087425828266 8.226816741867589 8.961571218345487 9.728780112646144;
             -1.8183619044915538 -1.771831627984931 -1.7438028632206626 -1.6685917204535754 -1.6302858014857406 -1.5093422844644864 -1.4403419931324906 -1.3247658289009951 -1.2083281555973817 -1.1229821929800383 -1.0505714158465091 -0.8173179899230579 -0.6158754729097995 -0.4173931635603583 -0.20656854219311493 0.03167040057514411 0.3179328129178279 0.6988924030220678 1.322675757704026 1.571254330141332 1.9208235093185737 2.5057336598815785 3.2620767489026155 3.494827429906371 3.821580239468593 4.382580788046574 5.103728314536232 5.330640017501667 5.654153205869778 6.197044525573566 6.931220396914263 7.148576870375071 7.383472025041834 7.97419921438612 8.823816150239953;
             -1.993303035023817 -1.9303016465983012 -1.9011878493296313 -1.809597739589436 -1.7604555789343959 -1.6186992259552142 -1.5378205899525677 -1.4012654453161015 -1.2667859059777642 -1.1695165793416582 -1.088073813777657 -0.8302604711079156 -0.6128988275759992 -0.40341742777142625 -0.18391138560002934 0.05895056579716881 0.34774871252504863 0.725472516184709 1.3257863522605307 1.5639734688636224 1.8941870167732415 2.4378536760134804 3.1376765263251087 3.3597689516182268 3.6615766661639153 4.184503386991316 4.855271567764788 5.072227985025396 5.36726178896159 5.8624270357507 6.553832805641278 6.757858279731317 7.027061678867403 7.507509459299485 7.938773729580974;
             -2.2545517114287588 -2.1697075864885766 -2.12187201988376 -2.003627986840613 -1.9429252192242403 -1.7605849774379487 -1.660685147362639 -1.49762799934165 -1.3385380674472869 -1.2255478542919374 -1.132652337293431 -0.8441203984143166 -0.6065848532259017 -0.38245676699262365 -0.1541978126415042 0.09480190481288699 0.38367083538080354 0.7546803845107093 1.3294592722342449 1.5534810345825008 1.860970482439231 2.363961558156508 3.005369250953606 3.206591077178068 3.4759320563848055 3.936940884025735 4.516582521230445 4.709896224583401 4.965706936814709 5.380382889324348 5.902619104241206 6.033874108179298 6.3350817756964535 6.803415323229341 7.179770904792615;
             -2.4559420446701052 -2.34920084792834 -2.296976185073357 -2.1420394469796284 -2.065104153183162 -1.8532671437387496 -1.7404482818386349 -1.5569468645307465 -1.3799823219728318 -1.2560175975037071 -1.1548651680857345 -0.8475332427807808 -0.5986528521480281 -0.3671741298455163 -0.13421441960815286 0.11717930226165238 0.4055787330043358 0.7709422173846994 1.3331585127712215 1.5514648960874131 1.8454640848313635 2.3221910595050956 2.927199195920216 3.112307311387735 3.3678052821501323 3.8019172475687917 4.3616672339609766 4.535393895749832 4.762542611054845 5.174277473836703 5.730501341808227 5.919474059847239 6.235336726812749 6.723206431226167 7.231651147242072;
             -2.569928572144317 -2.4656196647138655 -2.4014758345752387 -2.236262068336355 -2.158516518380391 -1.923673308482446 -1.7979065748604046 -1.5997846679626186 -1.4103321257056918 -1.2784135432631765 -1.1717696031635132 -0.8499132658320339 -0.5943484194811337 -0.35694312707262416 -0.1214887880639173 0.13112100670750404 0.4191870332071281 0.7814077014565005 1.3316584559236748 1.5432210796311716 1.8266084594238228 2.2920622811804403 2.856304977260614 3.0300396084181944 3.2684205939333277 3.664624203263345 4.176983049634148 4.326393560377753 4.552857168246306 4.92111391208768 5.502216961106142 5.636662090993625 5.850215772428456 6.1557325597112955 6.609738032513054;
             -2.972701482938346 -2.812945948534419 -2.7269368807141112 -2.505073465280368 -2.3962139122991437 -2.099707039134466 -1.9451305240930674 -1.7042574729580975 -1.4825402524085212 -1.330796855221357 -1.2110926561356254 -0.8546019664189122 -0.5788952082877374 -0.33039963362646674 -0.08728908009829067 0.16763773343275448 0.4519318546878791 0.8027235952910927 1.3216636485933564 1.5182514190293643 1.7792399561819003 2.196289357383261 2.704436161637905 2.8602211164423332 3.071314876385227 3.4161591494006935 3.856086976310044 3.999063549670627 4.180960369206441 4.484057207941962 4.871734136391144 5.018239006694096 5.187665970098557 5.479057559707716 5.841992903950639]

"""Asymptotic evaluation of the p-value for Anderson-Darling test"""
function pvalueasym(x::KSampleADTest)
    # Computational Details:
    #
    # This function uses the upper Tₘ quantiles as obtained via simulation of
    # the Anderson-Darling test statistics (2*10^6) with sample from standard
    # Normal distribution with size n=500 for each sample.
    #
    # After standardization of AD statistics, p-quantiles were estimated
    # for p ∈ [.00001,.00005,.0001,.0005,.001,.005,.01,.025,.05,.075,
    # .1,.2,.3,.4,.5,.6,.7,.8,.9,.925,.95,.975,.99,.9925,.995,.9975,.999,
    # .99925,.9995,.99975,.9999,.999925,.99995,.999975,.99999]
    #
    # Particular p-quantiles are determined from the simulations for various
    # values of `m`, m ∈ [1,2,3,4,6,8,10,20]. By interpolation, for the given
    # value of `m` from the performed test is done by fitting a cubic polynomial
    # of 1/sqrt(ms) to the simulated Tₘ quantiles.
    #
    # Next, the quadratic polynomial is used to fit the log((1-p)/p) to above
    # interpolated quantiles and the value fitted to the estimated standardized
    # AD statistics Tₘ.
    #
    # The p-values from Table 1 of the original paper were reproduced with
    # relative error bounded bounded by 1% in 85% of cases (see the relative
    # error table below).
    #
    # m\α   .25      .10     .05     .025     .01
    #  1   0.0001   0.0018  0.0006  0.0112   0.0305
    #  2  -0.0089  -0.0064  0.0024  0.0121   0.0306
    #  3  -0.0068  -0.0026  0.0045  0.0082   0.0164
    #  4  -0.0042  -0.0009  0.0027  0.005    0.0079
    #  6  -0.0025  -0.0005  0.005   0.003   -0.0004
    #  8  -0.0024   0.0011  0.0027  0.0023  -0.0023
    #  10 -0.0038   0.0008  0.0023  0.0026  -0.0043
    #  ∞  -0.0038   0.0014  0.005   0.0037   0.0148
    #

    m = x.k - 1
    Tk = (x.A²k - m) / x.σ
    sqm = 1.0 ./ sqrt.(M)

    T = x.modified ? TKM : TK

    # Construct Tm curve
    n = length(PV)
    Tm = zeros(n)
    for i in 1:n
        A = [ sqm[i]^p for i = 1:length(sqm), p = 0:3 ] # fit in cubic
        C = A \ T[:,i]
        x = 1/sqrt(m)
        Tm[i] = C[1] + C[2]*x + C[3]*x^2 + C[4]*x^3
    end
    logP =  log.((1 .- PV) ./ PV )

    # locate curve area for extrapolation
    _, j = findmin(abs(tm - Tk) for tm in Tm)
    A = [ Tm[i]^p for i = j-1:j+1, p = 0:2 ] # fit in quadratic
    C = A \ logP[j-1:j+1]
    lp0 = C[1] + C[2]*Tk + C[3]*Tk^2
    return exp(lp0)/(1 + exp(lp0))
end

pvalue(x::KSampleADTest) = x.nsim == 0 ? pvalueasym(x) : pvaluesim(x)

function adkvals(Z⁺, N, samples)
    k = length(samples)
    n = map(length, samples)
    L = length(Z⁺)

    fij = zeros(Int, k, L)
    for i in 1:k
        for s in samples[i]
            fij[i, searchsortedfirst(Z⁺, s)] += 1
        end
    end
    ljs = sum(fij, 1)

    A²k = A²km = 0.
    for i in 1:k
        innerm = 0.
        inner = 0.
        Mij = 0.
        Bj = 0.
        for j = 1:L
            lj = ljs[j]
            Mij += fij[i, j]
            Bj += lj
            Maij = Mij - fij[i, j]/2.
            Baj = Bj - lj/2.
            innerm += lj/N * (N*Maij-n[i]*Baj)^2 / (Baj*(N-Baj) - N*lj/4.)
            if j < L
                inner += lj/N * (N*Mij-n[i]*Bj)^2 / (Bj*(N-Bj))
            end
        end
        A²km += innerm / n[i]
        A²k += inner / n[i]
    end
    return A²k, A²km * (N - 1.) / N
end

function a2_ksample(samples, modified, method)
    k = length(samples)
    k < 2 && error("Need at least two samples")

    n = map(length, samples)
    pooled = vcat(samples...)
    Z = sort(pooled)
    N = length(Z)
    Z⁺ = unique(Z)
    L = length(Z⁺)

    L < 2 && error("Need more then 1 observation")
    minimum(n) == 0 && error("One of the samples is empty")

    A²k, A²km = adkvals(Z⁺, N, samples)

    H = sum(map(i->1./i, n))
    h = sum(1./(1:N-1))
    g = 0.
    for i in 1:N-2
        for j in i+1:N-1
            g += 1. / ((N - i) * j)
        end
    end

    a = (4*g - 6)*(k - 1) + (10 - 6*g)*H
    b = (2*g - 4)*k^2 + 8*h*k + (2*g - 14*h - 4)*H - 8*h + 4*g - 6
    c = (6*h + 2*g - 2)*k^2 + (4*h - 4*g + 6)*k + (2*h - 6)*H + 4*h
    d = (2*h + 6)*k^2 - 4*h*k
    σ² = (a*N^3 + b*N^2 + c*N + d) / ((N - 1.) * (N - 2.) * (N - 3.))

    KSampleADTest(k, N, sqrt(σ²), (modified ? A²km : A²k), modified, method, pooled, [n...])
end


text:
export CrackMesh
@doc raw"""
    CrackMesh(m::Int64, n::Int64, h::Float64, k::Int64 = 1)

Creates a crack mesh.

```julia
mmesh = CrackMesh(20, 10, 0.1, 4)
visualize_mesh(mmesh)
```

![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/AdFem/crackmesh.PNG?raw=true)

To access the underlying [`Mesh`](@ref) object, use `mmesh.mesh`.
"""
mutable struct CrackMesh
    m::Int64
    n::Int64
    h::Float64
    k::Int64
    nslip::Int64
    lower::Array{Int64,1}
    upper::Array{Int64,1}
    bdedge::Array{Int64,2}
    mesh::Mesh
    function CrackMesh(m::Int64, n::Int64, h::Float64, k::Int64 = 1)
        @assert mod(n, 2)==0
        @assert 2*k<=m
        coords = zeros((m+1)*(n+1)+m+1-2k,2)
        elem = zeros(Int64, 2m*n,3)
        s = 1
        for j = 1:n+1
            for i = 1:m+1
                x = (i-1)*h
                y = (j-1)*h
                coords[s, :] = [x;y]
                s += 1
            end
        end
        lower = [n÷2*(m+1)+k]
        upper = [n÷2*(m+1)+k]
        for i = k+1:m+1-k
            x = (i-1)*h
            y = n/2*h
            coords[s, :] = [x;y]
            s += 1
            push!(lower, n÷2*(m+1)+i)
            push!(upper, i-k+(m+1)*(n+1))
        end
        push!(lower, lower[end]+1)
        push!(upper, lower[end]+1)

        s = 1
        for i = 1:m
            for j = 1:n
                if j==(n÷2)+1 && i==k
                    elem[s,:] = [j*(m+1)+i; (j-1)*(m+1)+i; 1 + (m+1)*(n+1)]
                    elem[s+1,:] = [j*(m+1)+i; 1 + (m+1)*(n+1); j*(m+1)+i+1]
                elseif j==(n÷2)+1 && i==m-k+1
                    elem[s,:] = [j*(m+1)+i; m+1-2k+(m+1)*(n+1); (j-1)*(m+1)+i+1]
                    elem[s+1,:] = [j*(m+1)+i; (j-1)*(m+1)+i+1; j*(m+1)+i+1]
                elseif j==(n÷2)+1 && (i>k) && (i<m-k+1)
                    elem[s,:] = [j*(m+1)+i; i-k + (m+1)*(n+1); i-k+1 + (m+1)*(n+1)]
                    elem[s+1,:] = [j*(m+1)+i; i-k+1 + (m+1)*(n+1); j*(m+1)+i+1]
                else
                    elem[s,:] = [j*(m+1)+i; (j-1)*(m+1)+i; (j-1)*(m+1)+i+1]
                    elem[s+1,:] = [j*(m+1)+i; (j-1)*(m+1)+i+1; j*(m+1)+i+1]
                end
                s += 2
            end
        end

        mmesh = Mesh(m, n, h)
        bdedge = bcedge(mmesh)
        mmesh = Mesh(coords, elem, -1, 1, -1)
        nslip = m+1-2k
        new(m, n, h, k, nslip, lower, upper, bdedge, mmesh)
    end
end

function visualize_mesh(mmesh::CrackMesh)
    m, n, h = mmesh.m, mmesh.n, mmesh.h
    mmesh2 = copy(mmesh.mesh)
    mmesh2.nodes[(m+1)*(n+1)+1:end,2] .+= 0.3*h
    visualize_mesh(mmesh2)
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 21346b42-29f2-11eb-37be-93e9e286032b
using Automata

# ╔═╡ 287d8030-29f2-11eb-11f1-57b2b9f1038f
using Plots

# ╔═╡ 81ce2f92-29f2-11eb-1255-f74ff4d43aa1
using Images

# ╔═╡ 447fe1d0-29f5-11eb-2a6e-a5c55163cd23
md"
# 1D Cellular Automata Function
"

# ╔═╡ 36ac5460-29f2-11eb-1fc8-1d02c14dfc7c
plot(Gray.(Automata.to_array(continuous(x -> (x / 2) % 1, 100))))

# ╔═╡ 2d272abe-29f5-11eb-14a2-33e070036822
a(x, θ) = continuous(x -> sum(fill(x / 3, length(θ)) .^ collect(0:(length(θ)-1)) .* θ) % 1, 3, Float32[x])[end][1][3]

# ╔═╡ 862fd4a0-29f5-11eb-262e-7d7b5bcb39ff
r = collect(0:0.01:1)

# ╔═╡ c27a8db0-29f5-11eb-2495-6f8b73b8acbd
anim = @animate for i = r
	plot(r, (x -> a(x, [i*3, i*3])).(r), ylim=(0, 1))
end

# ╔═╡ e38993b2-29f6-11eb-1c2b-2f961ace6f85
gif(anim, "tmp.gif")

# ╔═╡ Cell order:
# ╟─447fe1d0-29f5-11eb-2a6e-a5c55163cd23
# ╠═21346b42-29f2-11eb-37be-93e9e286032b
# ╠═287d8030-29f2-11eb-11f1-57b2b9f1038f
# ╠═81ce2f92-29f2-11eb-1255-f74ff4d43aa1
# ╠═36ac5460-29f2-11eb-1fc8-1d02c14dfc7c
# ╠═2d272abe-29f5-11eb-14a2-33e070036822
# ╠═862fd4a0-29f5-11eb-262e-7d7b5bcb39ff
# ╠═c27a8db0-29f5-11eb-2495-6f8b73b8acbd
# ╠═e38993b2-29f6-11eb-1c2b-2f961ace6f85


text:
using RDatasets, DataFrames, Gadfly

set_default_plot_size(6inch, 3inch)

df = dataset("plm", "Cigar")
df = df[df.State .== 1, :]
ymin = df.Sales .- 20*0.23603334566204692
ymax = df.Sales .+ 20*0.34651701419196046

plot(df, x="Year", y="Sales", ymin=ymin, ymax=ymax,
     Geom.bar, Geom.errorbar)


text:
abstract type AbstractParticle end

struct particle <: AbstractParticle
	x :: Int64
	y :: Int64
	pos :: Vector
	infected :: Bool
	function particle(x :: Int64,y :: Int64,infected :: Bool)
		typeof([x,y])
		new(x,y,[x,y],infected)
	end
	function particle(pos :: Vector, infected :: Bool)
		new(pos[1],pos[2],pos,infected)
	end
end

function createLattice(size :: Int64, parts :: Vector)
	lattice = initaliseLattice(size,parts)
	for part in parts
		if part.infected == true
			lattice[part.x,part.y] =  -1
	    else
			lattice[part.x,part.y] = 1
		end
	end
	return lattice
end

function createLattice(lattice :: Matrix, parts :: Vector)
	size = length(lattice[1:end,1]) -2
	lattice = initaliseLattice(size,parts)
	for part in parts
		if part.infected == true
			lattice[part.x,part.y] =  -1
	    else
			lattice[part.x,part.y] = 1
		end
	end
	return lattice
end

function collision(part:: particle, lattice:: Matrix)
	dirVectors = ([0,1],[1,1],[1,0],[1,-1],[0,-1],[-1,-1],[-1,0],[-1,1])
	collQuery = [false,false,false,false,false,false,false,false]
	infection = false
	NBH = [part.pos + dirVector for dirVector in dirVectors]
	allowedNBH = []
	for i in 1:8
		if lattice[NBH[i][1],NBH[i][2]] == 0
		elseif lattice[NBH[i][1],NBH[i][2]] == -1
			infection = true
			collQuery[i] = true
		else
			collQuery[i] = true
		end
	end

	for i in 1:8
		if collQuery[i] == false
			push!(allowedNBH,NBH[i])
		end
	end
	return (allowedNBH, infection)

end

function move!(part:: particle, lattice:: Matrix)
	coldata = collision(part,lattice)
	if coldata[2] == true
		part= particle(rand(coldata[1]),true)
	else
		part= particle(rand(coldata[1]),part.infected)
	end
return part
end

function dynamics!(lattice, parts)
	for i in 1:length(parts)
		parts[i] = move!(parts[i],lattice)
		lattice = createLattice(lattice, parts)
	end
return (lattice,parts)
end

function initaliseParticles(size,num,initPercent)
	parts = []
	infParts = []
	nonInfParts = []

	for i in  1:round(num * initPercent)
		push!(parts, particle(rand(2:size-1, (2)),true))
		push!(infParts, last(parts))
	end

	for i in  1:(num - round(num * initPercent))
		push!(parts, particle(rand(2:size-1, (2)),false))
		push!(nonInfParts, last(parts))
	end

	return (parts,nonInfParts,infParts)
end

function initaliseLattice(size,parts)
	lattice = zeros( Int64, size +2 , size +2)
	for i in 1:(size+2)
		lattice[1,i] = 1
		lattice[size+2,i] = 1
		lattice[i,1] = 1
		lattice[i,size+2] = 1
	end
return lattice
end

function sim!(lattice,parts,n)
	if  n == 0
		return (lattice,parts)
	else
		println(n)
		step = sim!(lattice,parts,n-1)
		return dynamics!(step[1],step[2])

	end
end

function AnimatedSim(lattice, parts)

end

#Initialising Packages
using Colors,Plots

#input variables
size = 100
num = 30
initPercent = .30

#initialisation
partLists = initaliseParticles(size,num,initPercent)
parts = partLists[1]
nonInfParts = partLists[2]
infParts = partLists[3]
lattice = createLattice(size,parts)

anim = @animate for i in 1:500
	global lattice
	global parts
	heatmap(lattice,color = :redsblues)
	dyn = dynamics!(lattice, parts)
	lattice = dyn[1]
	parts = dyn[2]
end
gif(anim, "lattice_100.gif", fps =30)

max_stars_repo_name: wintProg/LatticePy

max_issues_repo_name: wintProg/LatticePy

max_forks_repo_name: wintProg/LatticePy


text:
# rigid body with quaternion representation (sphere w/ radius r)
mutable struct RigidBody{T} <: ContactModel
    dim::Dimensions
    m::T # mass
	J::Vector{T} # inertia
    g::T # gravity
    μ_world::T # friction coefficient
	r::T # radius

	base::BaseMethods
	dyn::DynamicsMethods

	joint_friction::SVector
end

function kinematics(model::RigidBody, q)
	return q[1:3] - [0.0; 0.0; model.r]
end

function lagrangian(model::RigidBody, q, q̇)
	m = 1.0
	J = Diagonal([1.0, 1.0, 1.0])

	p = q[1:3]
	quat = q[4:7]
	v = q̇[1:3]
	ω = q̇[4:6]

	L = 0.0

	# linear
	L += 0.5 * m * transpose(v) * v
	L -= m * model.g * p[3]

	# angular
	L += 0.5 * transpose(ω) * J * ω

	return L
end

# mass matrix
function M_func(model::RigidBody, q)
    m = model.m
	J = model.J
    Diagonal(SVector{6}([m, m, m, J...]))
end

# gravity
function C_func(model::RigidBody, q, q̇)
    m = model.m
    g = model.g

	ω = q̇[4:6]

    SVector{6}([0.0, 0.0, m * g, cross(ω, Diagonal(model.J) * ω)...]) #TODO confirm cross product
end

# signed distance function
function ϕ_func(model::RigidBody, env::Environment, q)
	SVector{1}(q[3] - model.r - env.surf(q[1:2]))
end

# control Jacobian
function B_func(model::RigidBody, q)
    SMatrix{6, 6}([1.0 0.0 0.0 0.0 0.0 0.0;
                   0.0 1.0 0.0 0.0 0.0 0.0;
                   0.0 0.0 1.0 0.0 0.0 0.0;
				   0.0 0.0 0.0 1.0 0.0 0.0;
                   0.0 0.0 0.0 0.0 1.0 0.0;
                   0.0 0.0 0.0 0.0 0.0 1.0])
end

# disturbance Jacobian
function A_func(model::RigidBody, q)
	SMatrix{3, 6}([1.0 0.0 0.0 0.0 0.0 0.0
				   0.0 1.0 0.0 0.0 0.0 0.0
				   0.0 0.0 1.0 0.0 0.0 0.0])
end

# contact Jacobian
function J_func(model::RigidBody, env::Environment, q)
	r = [0.0; 0.0; -model.r]
    SMatrix{3, 6}([Diagonal(ones(3)) -1.0 * skew(r)])
end

function contact_forces(model::RigidBody, env::Environment{<:World, LinearizedCone}, γ1, b1, q2, k)
	m = friction_mapping(env)

	SVector{3}(transpose(rotation(env, k)) * [m * b1; γ1])
end

function contact_forces(model::RigidBody, env::Environment{<:World, NonlinearCone}, γ1, b1, q2, k)
	SVector{3}(transpose(rotation(env, k)) * [b1; γ1])
end

function velocity_stack(model::RigidBody, env::Environment{<:World, LinearizedCone}, q1, q2, k, h)
	p1 = q1[1:3]
	quat1 = q1[4:7]

	p2 = q2[1:3]
	quat2 = q2[4:7]

	v = J_func(model, env, q2) * [(p2 - p1) / h[1]; ω_finite_difference(quat1, quat2, h[1])]

	v1_surf = rotation(env, k) * v

	SVector{4}(transpose(friction_mapping(env)) * v1_surf[1:2])
end

function velocity_stack(model::RigidBody, env::Environment{<:World,NonlinearCone}, q1, q2, k, h)
	nc = model.dim.c
	ne = dim(env)

	p1 = q1[1:3]
	quat1 = q1[4:7]

	p2 = q2[1:3]
	quat2 = q2[4:7]

	v = J_func(model, env, q2) * [(p2 - p1) / h[1]; ω_finite_difference(quat1, quat2, h[1])]

	v_surf = [rotation(env, k[(i-1) * ne .+ (1:ne)]) * v[(i-1) * ne .+ (1:ne)] for i = 1:nc]
	vT_stack = vcat([v_surf[i][1:ne-1] for i = 1:nc]...)
end

function dynamics(model::RigidBody, h, q0, q1, u1, w1, Λ1, q2)

	p0 = q0[1:3]
	quat0 = q0[4:7]

	p1 = q1[1:3]
	quat1 = q1[4:7]

	p2 = q2[1:3]
	quat2 = q2[4:7]

	# evalutate at midpoint
	ω1 = ω_finite_difference(quat0, quat1, h)
	ω2 = ω_finite_difference(quat1, quat2, h)

	qm1 = [0.5 * (p0 + p1); zeros(4)]
    vm1 = [(p1 - p0) / h[1]; ω1]
    qm2 = [0.5 * (p1 + p2); zeros(4)]
    vm2 = [(p2 - p1) / h[1]; ω2]

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	d = [(0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2)[1:3];
		 -1.0 * (Diagonal(model.J) * ω2 * sqrt(4.0 / h[1]^2.0 - transpose(ω2) * ω2)
			+ cross(ω2, Diagonal(model.J) * ω2)
			- Diagonal(model.J) * ω1 * sqrt(4.0 / h[1]^2.0 - transpose(ω1) * ω1)
			+ cross(ω1, Diagonal(model.J) * ω1))]

	return (d
		+ transpose(B_fast(model, qm2)) * u1
		+ transpose(A_fast(model, qm2)) * w1
		+ Λ1)
end

function dynamics(model::RigidBody, env::Environment, h, q0, q1, u1, w1, λ1, q2)
	Λ1 = transpose(J_func(model, env, q2)) * λ1
	dynamics(model, h, q0, q1, u1, w1, Λ1, q2)
end

function G_func(::RigidBody, q)
	quat = q[4:7]
	[1.0 0.0 0.0 0.0 0.0 0.0;
     0.0 1.0 0.0 0.0 0.0 0.0;
	 0.0 0.0 1.0 0.0 0.0 0.0;
     zeros(4, 3) attitude_jacobian(quat)]
end

function Gz_func(model::RigidBody, env, z)
	nz = num_var(model, env)
	ny = nz - 1

	[G_func(model, z) zeros(7, ny - 6);
	 zeros(nz - 7, 6) I]
end

# Model
rigidbody = RigidBody(Dimensions(7, 6, 3, 1, 1),
	10.0, [0.01, 0.01, 0.01], 9.81, 1.0, 0.25,
	BaseMethods(), DynamicsMethods(),
	SVector{6}(zeros(6)))

#
# @variables x
# @variables y
# @variables z
# @variables λ
#
# @variables px
# @variables py
# @variables pz
#
# L = (px - x)^2.0 + (py - y)^2.0 + (pz - z)^2.0 + λ * (z - x^2.0 - y^2.0)
#
# dL = Symbolics.gradient(L, [x; y; z; λ])
# dL = simplify.(dL)
#
# ddL = Symbolics.hessian(L, [x; y; z; λ])


text:
# Helper types

struct PVTKArgs
  part::Int
  nparts::Int
  ismain::Bool
  ghost_level::Int
end

struct PVTKFile <: VTKFile
  pvtkargs::PVTKArgs
  xdoc::XMLDocument
  vtk::DatasetFile
  path::String
end

# Main functions

"""
    pvtk_grid(args...;
              part,
              nparts,
              ismain = (part == 1),
              ghost_level = 0,
              kwargs...)

Return a handler representing a parallel vtk file, which can be
eventually written to file with `vtk_save`.

Positional and keyword arguments in `args` and `kwargs`
are passed to `vtk_grid` verbatim (except file names that are augmented with the
corresponding part id).

The extra keyword arguments only apply for parallel vtk file formats.  Mandatory ones are:

- `part` current (1-based) part id
- `nparts` total number of parts

Optional ones are
- `ismain` True if the current part id `part` is the main (the only one that will write the .pvtk file).
- `ghost_level` Ghost level.
"""
function pvtk_grid(filename::AbstractString,
                   args...;
                   part,
                   nparts,
                   ismain = (part == 1),
                   ghost_level = 0,
                   kwargs...)

  bname = basename(filename)
  mkpath(filename)
  prefix = joinpath(filename,bname)
  fn = _serial_filename(part,nparts,prefix,"")
  pvtkargs = PVTKArgs(part,nparts,ismain,ghost_level)
  xdoc  = XMLDocument()
  vtk = vtk_grid(fn,args...;kwargs...)
  _,ext = splitext(vtk.path)
  path = filename*".p"*ext[2:end]
  pvtk = PVTKFile(pvtkargs,xdoc,vtk,path)
  _init_pvtk!(pvtk)
  pvtk
end

# Add point and cell data as usual

function Base.setindex!(pvtk::PVTKFile,
                        data,
                        name::AbstractString,
                        loc::AbstractFieldData)
  pvtk.vtk[name,loc]=data
end

function Base.setindex!(pvtk::PVTKFile, data, name::AbstractString)
  pvtk.vtk[name]=data
end

# Save as usual

function vtk_save(pvtk::PVTKFile)
  if isopen(pvtk)
    if pvtk.pvtkargs.ismain
      _update_pvtk!(pvtk)
      save_file(pvtk.xdoc, pvtk.path)
    end
    vtk_save(pvtk.vtk)
    close(pvtk)
  end
  return [pvtk.path]
end

# Helper functions

function _serial_filename(part,nparts,prefix,extension)
  p = lpad(part,ceil(Int,log10(nparts)),'0')
  fn = prefix*"_$p"*extension
end

function _init_pvtk!(pvtk::PVTKFile)

  # Recover some data
  vtk = pvtk.vtk
  pvtkargs = pvtk.pvtkargs
  pgrid_type = "P"*vtk.grid_type
  npieces = pvtkargs.nparts
  pref,_ = splitext(pvtk.path)
  _,ext = splitext(vtk.path)
  prefix = joinpath(pref,basename(pref))

  # VTKFile (root) node
  pvtk_root = create_root(pvtk.xdoc, "VTKFile")
  set_attribute(pvtk_root, "type", pgrid_type)
  set_attribute(pvtk_root, "version", "1.0")
  if IS_LITTLE_ENDIAN
    set_attribute(pvtk_root, "byte_order", "LittleEndian")
  else
    set_attribute(pvtk_root, "byte_order", "BigEndian")
  end

  # Grid node
  pvtk_grid = new_child(pvtk_root,pgrid_type)
  set_attribute(pvtk_grid, "GhostLevel", string(pvtkargs.ghost_level))

  # Pieces (i.e. Pointers to serial files)
  for piece in 1:npieces
    pvtk_piece = new_child(pvtk_grid,"Piece")
    fn = _serial_filename(piece,npieces,prefix,ext)
    set_attribute(pvtk_piece,"Source",fn)
  end

  # Recover point type and number of components
  vtk_root = root(vtk.xdoc)
  vtk_grid = find_element(vtk_root,vtk.grid_type)
  vtk_piece = find_element(vtk_grid,"Piece")
  vtk_points = find_element(vtk_piece,"Points")
  vtk_data_array = find_element(vtk_points,"DataArray")
  point_type = attribute(vtk_data_array,"type")
  Nc = attribute(vtk_data_array,"NumberOfComponents")

  ## PPoints node
  pvtk_ppoints = new_child(pvtk_grid,"PPoints")
  pvtk_pdata_array = new_child(pvtk_ppoints,"PDataArray")
  set_attribute(pvtk_pdata_array,"type",point_type)
  set_attribute(pvtk_pdata_array,"Name","Points")
  set_attribute(pvtk_pdata_array,"NumberOfComponents",Nc)

  pvtk
end

function _update_pvtk!(pvtk::PVTKFile)

  vtk = pvtk.vtk
  vtk_root = root(vtk.xdoc)
  vtk_grid = find_element(vtk_root,vtk.grid_type)
  vtk_piece = find_element(vtk_grid,"Piece")

  pgrid_type = "P"*vtk.grid_type
  pvtk_root = root(pvtk.xdoc)
  pvtk_grid = find_element(pvtk_root,pgrid_type)

  # Generate PPointData
  vtk_point_data = find_element(vtk_piece,"PointData")
  if vtk_point_data !== nothing
    pvtk_ppoint_data = new_child(pvtk_grid,"PPointData")
    for vtk_data_array in child_elements(vtk_point_data)
      t = attribute(vtk_data_array,"type")
      name = attribute(vtk_data_array,"Name")
      Nc = attribute(vtk_data_array,"NumberOfComponents")
      pvtk_pdata_array = new_child(pvtk_ppoint_data,"PDataArray")
      set_attribute(pvtk_pdata_array,"type",t)
      set_attribute(pvtk_pdata_array,"Name",name)
      set_attribute(pvtk_pdata_array,"NumberOfComponents",Nc)
    end
  end

  # Generate PCellData
  vtk_cell_data = find_element(vtk_piece,"CellData")
  if vtk_cell_data !== nothing
    pvtk_pcell_data = new_child(pvtk_grid,"PCellData")
    for vtk_data_array in child_elements(vtk_cell_data)
      t = attribute(vtk_data_array,"type")
      name = attribute(vtk_data_array,"Name")
      Nc = attribute(vtk_data_array,"NumberOfComponents")
      pvtk_pdata_array = new_child(pvtk_pcell_data,"PDataArray")
      set_attribute(pvtk_pdata_array,"type",t)
      set_attribute(pvtk_pdata_array,"Name",name)
      set_attribute(pvtk_pdata_array,"NumberOfComponents",Nc)
    end
  end

  pvtk
end


text:
"""
	checktargets(x::T where T<:AbstractArray, y::S where S<:AbstractArray)::S

Function that enforces the targets (e.g. **y** field) of a datacell to be aligned with the size of the data matrix.
"""
function checktargets(x::T where T<:AbstractArray, y::S where S<:AbstractArray)

	# If target list is empty, issue error
	@assert !isempty(y)

	# If it is not empty, enforce the fact that it should not contain fewer points than the data
	@assert nobs(x) <= nobs(y) "[checktargets] Data size larger than targets size. Resize accordingly."

	@inline _trunc_(y::T, n::Int) where T<:SubArray{S,1} where S = view(y,1:n)
	@inline _trunc_(y::T, n::Int) where T<:SubArray{S,2} where S = view(y,:,1:n)
	@inline _trunc_(y::T, n::Int) where T<:AbstractArray{S,1} where S = y[1:n]
	@inline _trunc_(y::T, n::Int) where T<:AbstractArray{S,2} where S = y[:,1:n]
	if (nobs(x) < nobs(y))
		warn("[checktargets] Data size size smaller than targets size; taking only the first $(nobs(x)) elements." )
		return _trunc_(y, nobs(x))                                                  	# Truncate label list to data size
	end
	return y
end

# Mask checking functions - takes as input arguments the column indices / mask and the number of columns
checkcolumnmask(x::Vector{T} where T<:Int, N::Integer) = begin
	if all( map( (x)-> x>=0 && x<=1, x ) )							# Check if all integer values are in [0 1] e.g. 0 or 1
		x = Vector{Bool}(x)								# If all values are 0 or 1, convert to boolean array and
		checkcolumnmask(x,N)								# apply booblean checks
	else											# else, apply integer checks
		@assert length(x) == length(unique(x))						# The mask must contain unique values only
		@assert all(x .> 0)								# No columns should be smaller than 0
		@assert all(x .<= N)	 							# No columns should be larger than the number of columns
		@assert ~isempty(x)								# Mask should not be empty
	end
	xnew = falses(N)
	xnew[x] = true
	return xnew
end

checkcolumnmask(x::Vector{T} where T<:Bool, N::Integer) = begin
	@assert length(x) == N									# Enforce size
	@assert ~isempty(x)									# Mask should not be empty
	@assert ~isequal(x, falses(size(x))) 							# At least one mask value should be 'true' e.g. to be transformed ...
	return BitVector(x)
end

checkcolumnmask(x::T where T<:BitVector, N::Integer) = begin
	@assert length(x) == N									# Enforce size
	@assert ~isempty(x)									# Mask should not be empty
	@assert ~isequal(x, T(falses(size(x))) ) 						# At least one mask value should be 'true' e.g. to be transformed ...i
	return x
end

"""
	checkpriors(priors::Vector{T} where T<:Real)

Function that checks the consistency of a priors vector (e.g. all elements between 0 and 1, they sum up to 1)
and returns true of false depending on wether the prior vector is usable or not.
"""
function checkpriors(priors::Vector{T} where T<:Real)
	if (abs(sum(priors) - 1) > 1e-3)							# Check that the priors summate to 1
		return false
	end

	if (any(priors .< 0) || any(priors .> 1))						# Check that all priors are between 0 and 1
		return false
	end
	return true
end


text:
using Stochy
import Stochy.support
using Base.Test

for f in ["inference", "variational", "erp", "cps", "store", "addressing", "mem", "rand"]
    include("$(f)tests.jl")
end

println("Passed!")


text:
@testset "merge PatchEnv instances" begin
    multiply(x::Number) = 2x
    multiply(x::Int) = 2x - 1
    add(x::Number) = x + 2
    add(x::Int) = x + 1

    patches = Patch[
        @patch multiply(x::Integer) = 3x
        @patch multiply(x::Int) = 4x
        @patch add(x::Int) = x + 4
    ]

    @testset "simple" begin
        pe1 = Mocking.PatchEnv(patches[1])
        pe2 = Mocking.PatchEnv(patches[2:3])
        pe = Mocking.PatchEnv(patches)

        @test pe == merge(pe1, pe2)
    end

    @testset "debug flag" begin
        pe1 = Mocking.PatchEnv(patches[1], true)
        pe2 = Mocking.PatchEnv(patches[2:3])
        pe = Mocking.PatchEnv(patches, true)

        @test pe == merge(pe1, pe2)
    end
end


text:
using StructuresKit

shape_name = "W14X90"

shape_info = StructuresKit.CrossSection.AISC(shape_name)

#A Ix Iy J Cw xc yc xs ys
section_properties = [(shape_info.A, shape_info.Ix, shape_info.Iy, shape_info.J, shape_info.Cw, 0.0, shape_info.d/2, 0.0, shape_info.d/2)]

#define column cross-section mesh
n_Wshape=(4, 2, 4, 4, 4)
mesh_size = 0.01
xcoords, ycoords = CrossSection.wshape_nodes(shape_info, n_Wshape)
section_geometry = Column.SectionGeometry(xcoords, ycoords, mesh_size)

#define material properties
Es = 29000.0
σy = 50.0                   #steel yield stress
σy1 = 50.0                    #steel yield stress
σu = 65.3                   #steel ultimate stress
σf= σu      #steel fracture stress#
ϵy = σy / Es              #steel yield strain
ϵy1 = σy / Es * 10    #steel strain at end of yield plateau
ϵu = 0.18             #steel ultimate strain
ϵf = 0.21             #steel fracture strain
n = [1000, 1000, 1000]

material_properties = [(29000,0.30, σy, σy1, σu, σf, ϵy, ϵy1, ϵu, ϵf, n)]

#kx ky kϕ hx hy
springs = [0.0, 0.0, 0.0, 0.0, 0.0]

#member information
#L dL SectionProperties MaterialProperties Springs
member_definitions = [(13.0*12, 1.0,1,1,1)]

num_nodes = Int(member_definitions[1][1]/member_definitions[1][2]) + 1

#P
loads = 10000.0 * ones(num_nodes)

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
end_boundary_conditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0, 13.0*12]

#imperfection
L = member_definitions[1][1]
dL = member_definitions[1][2]
z = 0: dL :L
Δo = L / 1000 * sin.(π*z/L)

uo = Δo  #in the weak axis
vo = zeros(Float64, num_nodes)
ϕo = zeros(Float64, num_nodes)

imperfections = [uo, vo, ϕo]

#Calculate weak axis flexural buckling load
E = material_properties[1][1]
Iy = shape_info.Iy
Pcry = π^2*E*Iy/L^2   #weak axis buckling load

#calculate yield load
Py = 50*shape_info.A

#define load range for a VR test
P = 0:Py/10:0.90*Py

P = [P; 0.901 * Py]

#inelasticity
#cross-section mesh level inelasticity is not working yet here, so keep it simple
#just reduce Iy at first yield for now
strain_limit = 0.0
solution_controls = Column.SolutionControls(" ", strain_limit)

def_coordinates = []
connectivity = []

u_total_max = zeros(Float64, length(P))

for i = 1:length(P)

    loads = P[i] * ones(num_nodes)

    #reduce Iy when column fails
    if i == length(P)

        #A Ix Iy J Cw xc yc xs ys
        section_properties = [(shape_info.A, shape_info.Ix, shape_info.Iy * 0.30, shape_info.J, shape_info.Cw, 0.0, shape_info.d/2, 0.0, shape_info.d/2)]

    end

    #define column information
    column_properties = Column.initialize(member_definitions, section_properties, section_geometry, material_properties, loads, springs, end_boundary_conditions, supports, imperfections, solution_controls)

#   if i == length(P)

#         #A Ix Iy J Cw xc yc xs ys
#         column_properties.Iy[76:80] .= shape_info.Iy * 0.30

#     end

    #solve for column deformation
    u, v, ϕ = Column.solve(column_properties, solution_controls)

    #add initial imperfection to solution
    u_total = uo .+ u
    v_total = vo .+ v
    ϕ_total = ϕo .+ ϕ

    #get max out-of-plane deflection
    u_total_max[i] = maximum(u_total)

    #create surface mesh of undeformed column

    #define discretization along member
    zcoords = 0:dL:L
    n=(4, 4, 4, 4, 4)

    xcoords, ycoords = StructuresKit.CrossSection.wshape_nodes(shape_info, n)

    #undeformed shape
    coordinates, connectivity = StructuresKit.Mesh.surface(xcoords, ycoords, zcoords)

    #deformed shape

    σ_total_normal = Column.normal_stresses(xcoords, ycoords, zcoords, u, v, column_properties)
    ϵ_total_normal = Column.normal_stains(column_properties, σ_total_normal)
    w = Column.warping_displacements(zcoords, ϵ_total_normal)

    scale_u = 1.0
    scale_v = 1.0
    scale_ϕ = 1.0
    scale_w = 1.0

    if i == length(P)
        scale_u = 1.0
        scale_w = 1.0
    end

    x_deformed, y_deformed, z_deformed = Column.deformed_shape(xcoords, ycoords, zcoords, u, v, ϕ, w, scale_u, scale_v, scale_ϕ, scale_w)

    filename = string(@__DIR__, "/", "W14x90", "_SS_", string(i-1), "_R1", ".ply")

    def_coordinates = [x_deformed y_deformed z_deformed]

    # Visualize.create_ply_file(def_coordinates, connectivity, filename)

end

using Plots
SSRC = plot(u_total_max * 25.4, P * 4.448, markershape = :o, legend = false, fmt = :pdf)
Plots.xlims!(0, 8)
Plots.ylims!(0, 6000)
Plots.xlabel!("column midheight deflection, mm")
Plots.ylabel!("column compressive load, kN")

display(SSRC)

savefig(SSRC, "SSRC.pdf")

using Makie
scene = poly(def_coordinates, connectivity, color=:lightgray, shading=:true, show_axis=:true, overdraw=:false, strokecolor = (:black, 0.6), strokewidth = 1)


text:
## SingFun stieltjes

export logabs
logabs(x)=log(abs2(x))/2

# sqrtx2 is analytic continuation of sqrt(z^2-1)
sqrtx2(z::Complex)=sqrt(z-1)*sqrt(z+1)
sqrtx2(x::Real)=sign(x)*sqrt(x^2-1)
@vectorize_1arg Number sqrtx2
function sqrtx2(f::Fun)
    B=Evaluation(first(domain(f)))
    A=Derivative()-f*differentiate(f)/(f^2-1)
    linsolve([B,A],[sqrtx2(first(f))];tolerance=ncoefficients(f)*10E-15)
end

# intervaloffcircle maps the slit plane to the interior(true)/exterior(false) disk
# intervaloncircle maps the interval to the upper(true)/lower(false) half circle

# it is more accurate near infinity to do 1/J_- as it avoids round off
intervaloffcircle(s::Bool,z)=s?1./intervaloffcircle(false,z):(z+sqrtx2(z))
intervaloncircle(s::Bool,x)=x+1.0im*(s?1:-1).*sqrt(1-x).*sqrt(x+1)

intervaloffcircle(s::Int,x)=intervaloffcircle(s==1,x)
intervaloncircle(s::Int,x)=intervaloncircle(s==1,x)

#TODO: These aren't quite typed correctly but the trouble comes from anticipating the unifying type without checking every element.

updownjoukowskyinverse{T<:Number}(s::Bool,x::T) = in(x,Interval(-one(T),one(T))) ? intervaloncircle(s,x) : intervaloffcircle(s,x)
updownjoukowskyinverse{T<:Number}(s::Bool,x::Vector{T}) = Complex{real(T)}[updownjoukowskyinverse(s,xk) for xk in x]
updownjoukowskyinverse{T<:Number}(s::Bool,x::Array{T,2}) = Complex{real(T)}[updownjoukowskyinverse(s,x[k,j]) for k=1:size(x,1),j=1:size(x,2)]

function hornersum{S<:Number,V<:Number}(cfs::AbstractVector{S},y::V)
    N,P = length(cfs),promote_type(S,V)
    ret = N > 0 ? convert(P,cfs[N]) : zero(P)
    for k=N-1:-1:1
        ret = muladd(y,ret,cfs[k])
    end
    y*ret
end

hornersum{S<:Number,V<:Number}(cfs::AbstractVector{S},y::AbstractVector{V}) = promote_type(S,V)[hornersum(cfs,y[k]) for k=1:length(y)]
hornersum{S<:Number,V<:Number}(cfs::AbstractVector{S},y::AbstractArray{V,2}) = promote_type(S,V)[hornersum(cfs,y[k,j]) for k=1:size(y,1),j=1:size(y,2)]

function divkhornersum{S<:Number,T<:Number,U<:Number,V<:Number}(cfs::AbstractVector{S},y::T,ys::U,s::V)
    N,P = length(cfs),promote_type(S,T,U,V)
    ret = N > 0 ? convert(P,cfs[N]/(N+s)) : zero(P)
    for k=N-1:-1:1
        ret=muladd(y,ret,cfs[k]/(k+s))
    end
    y*ys*ret
end

divkhornersum{S<:Number,T<:Number,U<:Number,V<:Number}(cfs::AbstractVector{S},y::AbstractVector{T},ys::AbstractVector{U},s::V) =
    promote_type(S,T,U,V)[divkhornersum(cfs,y[k],ys[k],s) for k=1:length(y)]
divkhornersum{S<:Number,T<:Number,U<:Number,V<:Number}(cfs::AbstractVector{S},y::AbstractArray{T,2},ys::AbstractArray{U,2},s::V) =
    promote_type(S,T,U,V)[divkhornersum(cfs,y[k,j],ys[k,j],s) for k=1:size(y,1),j=1:size(y,2)]

realdivkhornersum{S<:Real}(cfs::AbstractVector{S},y,ys,s) = real(divkhornersum(cfs,y,ys,s))
realdivkhornersum{S<:Complex}(cfs::AbstractVector{S},y,ys,s) = complex(real(divkhornersum(real(cfs),y,ys,s)),real(divkhornersum(imag(cfs),y,ys,s)))

function stieltjes{S<:PolynomialSpace,DD<:Interval}(sp::JacobiWeight{S,DD},u,zv::Array,s...)
    ret=similar(zv,Complex128)
    for k=1:length(zv)
        @inbounds ret[k]=stieltjes(sp,u,zv[k],s...)
    end
    ret
end

function stieltjes{S<:PolynomialSpace,DD<:Interval}(sp::JacobiWeight{S,DD},u,z)
    d=domain(sp)

    if sp.α == sp.β == .5
        cfs = coefficients(u,sp.space,Ultraspherical(1,d))
        π*hornersum(cfs,intervaloffcircle(true,mobius(sp,z)))
    elseif sp.α == sp.β == -.5
        cfs = coefficients(u,sp.space,ChebyshevDirichlet{1,1}(d))
        z=mobius(sp,z)

        sx2z=sqrtx2(z)
        sx2zi=1./sx2z
        Jm=1./(z+sx2z)  # intervaloffcircle(true,z)

        if length(cfs) ≥1
            ret = π*cfs[1]*sx2zi

            if length(cfs) ≥2
                ret += cfs[2]*π*(z.*sx2zi-1)
            end

            ret - 2π*hornersum(cfs[3:end],Jm)
        else
            zero(z)
        end
    elseif isapproxinteger(sp.α) && isapproxinteger(sp.β)
        stieltjes(coefficients(u,sp,sp.space),z)
    else
        if d==Interval()
            S2=JacobiWeight(sp.α,sp.β,Jacobi(sp.β,sp.α))  # convert and then use recurrence
            stieltjesintervalrecurrence(S2,coefficients(u,sp,S2),z)
        else
            # project to interval
            stieltjes(setdomain(sp,Interval()),u,mobius(sp,z))
        end
    end
end

function stieltjes{SS<:PolynomialSpace,DD<:Interval}(sp::JacobiWeight{SS,DD},u,x::Number,s::Bool)
    d=domain(sp)

    if sp.α == sp.β == .5
        cfs=coefficients(u,sp.space,Ultraspherical(1,d))
        π*hornersum(cfs,intervaloncircle(!s,mobius(sp,x)))
    elseif sp.α == sp.β == -.5
        cfs = coefficients(u,sp.space,ChebyshevDirichlet{1,1}(d))
        x=mobius(sp,x)

        if length(cfs) ≥1
            ret = -π*im*cfs[1]*(s?1:-1)/sqrt(1-x^2)

            if length(cfs) ≥2
                ret += cfs[2]*(-π*im*(s?1:-1)*x/sqrt(1-x^2)-π)
            end

            ret - 2π*hornersum(cfs[3:end],intervaloncircle(!s,x))
        else
            0.0+0.0im
        end
    else
        if d==Interval()
            S=JacobiWeight(sp.α,sp.β,Jacobi(sp.β,sp.α))
            cfs=coefficients(u,sp,S)
            cf=stieltjesforward(S,length(cfs),x,s)
            dotu(cf,cfs)
        else
            @assert isa(d,Interval)
            stieltjes(setdomain(sp,Interval()),u,mobius(sp,x),s)
        end
    end
end

##
#  hilbert on JacobiWeight space
#  hilbert is always equal to im*(C^+ + C^-)
#  hilbert(f,z)=im*(cauchy(true,f,z)+cauchy(false,f,z))
##

function hilbert{DD<:Interval}(sp::JacobiWeight{Chebyshev{DD},DD},u)
    d=domain(u)

    if sp.α == sp.β == .5
        # Corollary 5.7 of Olver&Trogdon
        cfs=coefficients(u,sp.space,Ultraspherical(1))
        Fun([0.;-cfs],d)
    elseif sp.α == sp.β == -.5
        # Corollary 5.11 of Olver&Trogdon
        cfs= coefficients(u,sp.space,ChebyshevDirichlet{1,1})
        if length(cfs)≥2
            Fun([cfs[2];2cfs[3:end]],d)
        else
            Fun(zeros(eltype(cfs),1),d)
        end
    else
        error("hilbert only implemented for Chebyshev weights")
    end
end

integratejin(c,cfs,y)=.5*(-cfs[1]*(log(y)+log(c))+divkhornersum(cfs,y,y,1)-divkhornersum(view(cfs,2:length(cfs)),y,zero(y)+1,0))
realintegratejin(c,cfs,y)=.5*(-cfs[1]*(logabs(y)+logabs(c))+realdivkhornersum(cfs,y,y,1)-realdivkhornersum(view(cfs,2:length(cfs)),y,zero(y)+1,0))

#########
# stieltjesintegral is an indefinite integral of stieltjes
# normalized so that there is no constant term
# logkernel is the real part of stieljes normalized by π.
#####

function logkernel{S<:PolynomialSpace,DD<:Interval}(sp::JacobiWeight{S,DD},u,z)
    d=domain(sp)
    a,b=d.a,d.b

    if sp.α == sp.β == .5
        cfs=coefficients(u,sp.space,Ultraspherical(1,d))
        z=mobius(sp,z)
        y = updownjoukowskyinverse(true,z)
        arclength(d)*realintegratejin(4/(b-a),cfs,y)/2
    elseif  sp.α == sp.β == -.5
        cfs = coefficients(u,sp.space,ChebyshevDirichlet{1,1}(d))
        z=mobius(sp,z)
        y = updownjoukowskyinverse(true,z)

        if length(cfs) ≥1
            ret = -cfs[1]*arclength(d)*(logabs(y)+logabs(4/(b-a)))/2

            if length(cfs) ≥2
                ret += -arclength(d)*cfs[2]*real(y)/2
            end

            if length(cfs) ≥3
                ret - arclength(d)*realintegratejin(4/(b-a),view(cfs,3:length(cfs)),y)
            else
                ret
            end
        else
            zero(z)
        end
    else
        error("logkernel not implemented for parameters "*string(sp.α)*","*string(sp.β))
    end
end

function stieltjesintegral{S<:PolynomialSpace,DD<:Interval}(sp::JacobiWeight{S,DD},u,z)
    d=domain(sp)
    a,b=d.a,d.b

    if sp.α == sp.β == .5
        cfs=coefficients(u,sp.space,Ultraspherical(1,d))
        y=intervaloffcircle(true,mobius(sp,z))
        π*complexlength(d)*integratejin(4/(b-a),cfs,y)/2
    elseif  sp.α == sp.β == -.5
        cfs = coefficients(u,sp.space,ChebyshevDirichlet{1,1}(d))
        z=mobius(sp,z)
        y=intervaloffcircle(true,z)

        if length(cfs) ≥1
            ret = -cfs[1]*π*complexlength(d)*(log(y)+log(4/(b-a)))/2

            if length(cfs) ≥2
                ret += -π*complexlength(d)*cfs[2]*intervaloffcircle(true,z)/2
            end

            if length(cfs) ≥3
                ret - π*complexlength(d)*integratejin(4/(b-a),view(cfs,3:length(cfs)),y)
            else
                ret
            end
        else
            zero(z)
        end
    else
        error("stieltjes integral not implemented for parameters "*string(sp.α)*","*string(sp.β))
    end
end

function stieltjesintegral{S<:PolynomialSpace,DD<:Interval}(sp::JacobiWeight{S,DD},u,z,s::Bool)
    d=domain(u)
    a,b=d.a,d.b     # TODO: type not inferred right now

    if sp.α == sp.β == .5
        cfs=coefficients(u,sp.space,Ultraspherical(1,d))
        y=intervaloncircle(!s,mobius(sp,z))
        π*complexlength(d)*integratejin(4/(b-a),cfs,y)/2
    elseif  sp.α == sp.β == -.5
        cfs = coefficients(u,sp.space,ChebyshevDirichlet{1,1}(d))
        z=mobius(sp,z)
        y=intervaloncircle(!s,z)

        if length(cfs) ≥1
            ret = -cfs[1]*π*complexlength(d)*(log(y)+log(4/(b-a)))/2

            if length(cfs) ≥2
                ret += -π*complexlength(d)*cfs[2]*intervaloncircle(!s,z)/2
            end

            if length(cfs) ≥3
                ret - π*complexlength(d)*integratejin(4/(b-a),view(cfs,3:length(cfs)),y)
            else
                ret
            end
        else
            zero(z)
        end
    else
        error("stieltjes integral not implemented for parameters "*string(sp.α)*","*string(sp.β))
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

"""
    interpolate(B, T, xi)

Given basis B, interpolate T at xi.

# Example
```jldoctest
B = Quad4()
X = Vec.([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)])
T = [1.0, 2.0, 3.0, 4.0]
interpolate(B, T, Vec(0.0, 0.0))

# output

2.5
```
"""
function interpolate(B::AbstractBasis{dim}, T::Vector, xi::Vec{dim}) where {dim}
    N = eval_basis(B, xi)
    return sum(b*t for (b, t) in zip(N, T))
end

"""
    jacobian(B, X, xi)

Given basis B, calculate jacobian at xi.

# Example
```jldoctest
B = Quad4()
X = Vec.([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)])
jacobian(B, X, Vec((0.0, 0.0)))

# output

2×2 Tensor{2,2,Float64,4}:
 0.5  0.0
 0.0  0.5

```
"""
jacobian(B::AbstractBasis{dim}, X::Vector{<:Vec{dim}}, xi::Vec{dim}) where {dim} = jacobian(B, X, xi, eval_dbasis(B, xi))

function jacobian(B::AbstractBasis{dim}, X::Vector{<:Vec{dim}}, xi::Vec{dim}, dB::Vector{<:Vec{dim}}) where {dim}
    @assert length(X) == length(B) == length(dB)
    J = zero(Tensor{2, dim})
    @inbounds for i in 1:length(X)
        J += otimes(dB[i], X[i]) # dB[i] ⊗ X[i]
    end
    return J
end

"""
    grad(B, X, xi)

Given basis B, calculate gradient dB/dX at xi.

# Example
```jldoctest
B = Quad4()
X = Vec.([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)])
grad(B, X, Vec(0.0, 0.0))

# output

4-element Array{Tensor{1,2,Float64,2},1}:
 [-0.5, -0.5]
 [0.5, -0.5]
 [0.5, 0.5]
 [-0.5, 0.5]

```
"""
grad(B::AbstractBasis{dim}, X::Vector{<:Vec{dim}}, xi::Vec{dim}) where {dim} =
    grad!(B, similar(X), X, xi, eval_dbasis(B, xi))

function grad!(B::AbstractBasis{dim}, dN::Vector{<:Vec{dim}}, X::Vector{<:Vec{dim}}, xi::Vec{dim}, dB::Vector{<:Vec{dim}}) where {dim}
    @assert length(dN) == length(dB)
    J = jacobian(B, X, xi, dB)
    @inbounds for i in 1:length(dN)
        dN[i] = inv(J) ⋅ dB[i]
    end
    return dN
end

"""
    grad(B, T, X, xi)

Calculate gradient of `T` with respect to `X` in point `xi` using basis `B`.

# Example
```jldoctest
B = Quad4()
X = Vec.([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)])
u = Vec.([(0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0)])
grad(B, u, X, Vec(0.0, 0.0))

# output

julia> grad(B, u, X, Vec(0.0, 0.0))
2×2 Tensor{2,2,Float64,4}:
 1.5  0.5
 1.0  2.0

```
"""
function grad(B::AbstractBasis{dim}, T::Vector{<:Vec{dim}}, X::Vector{<:Vec{dim}}, xi::Vec{dim}) where {dim}
    G = grad(B, X, xi) # <- allocates
    dTdX = sum(T[i] ⊗ G[i] for i=1:length(B))
    return dTdX
end
function grad(B::AbstractBasis{dim}, T::Vector{<:Number}, X::Vector{<:Vec{dim}}, xi::Vec{dim}) where {dim}
    G = grad(B, X, xi) # <- allocates
    dTdX = sum(T[i] * G[i] for i=1:length(B))
    return dTdX
end

"""
Data type for fast FEM.
"""
mutable struct BasisInfo{B<:AbstractBasis,dim, T, M}
    N::Vector{T}
    dN::Vector{Vec{dim, T}}
    grad::Vector{Vec{dim, T}}
    J::Tensor{2, dim, T, M}
    invJ::Tensor{2, dim, T, M}
    detJ::T
    basis::Type{B}
end

Base.length(B::BasisInfo{T}) where T<:AbstractBasis = length(T)
Base.size(B::BasisInfo{T})   where T<:AbstractBasis = size(T)

"""
Initialization of data type `BasisInfo`.

# Examples

```jldoctest

BasisInfo(Tri3)

# output

FEMBasis.BasisInfo{FEMBasis.Tri3,Float64}([0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0], 0.0)

```

"""
function BasisInfo(::Type{B}, T=Float64) where B <: AbstractBasis{dim} where dim
    nbasis = length(B)
    N = zeros(T, nbasis)
    dN = zeros(Vec{dim, T}, nbasis)
    grad = zeros(Vec{dim, T}, nbasis)
    J = zero(Tensor{2, dim, T})
    invJ = zero(Tensor{2, dim, T})
    detJ = zero(T)
    return BasisInfo(N, dN, grad, J, invJ, detJ, B)
end

"""
Evaluate basis, gradient and so on for some point `xi`.

# Examples

```jldoctest

b = BasisInfo(Quad4)
X = Vec.([(0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0)])
xi = Vec(0.0, 0.0)
eval_basis!(b, X, xi)

# output

BasisInfo{Quad4,2,Float64,4}([0.25, 0.25, 0.25, 0.25], Tensors.Tensor{1,2,Float64,2}[[-0.25, -0.25], [0.25, -0.25], [0.25, 0.25], [-0.25, 0.25]], Tensors.Tensor{1,2,Float64,2}[[-0.5, -0.5], [0.5, -0.5], [0.5, 0.5], [-0.5, 0.5]], [0.5 0.0; 0.0 0.5], [2.0 -0.0; -0.0 2.0], 0.25, Quad4)

```
"""
function eval_basis!(bi::BasisInfo{B},
                     X::Vector{<:Vec{dim}}, xi::Vec{dim}) where B <: AbstractBasis{dim} where dim
    # evaluate basis and derivatives
    eval_basis!(B, bi.N, xi)
    eval_dbasis!(B, bi.dN, xi)

    # calculate Jacobian
    bi.J = jacobian(B(), X, xi, bi.dN)

    # calculate determinant of Jacobian + gradient operator

    # TODO, fixup curve + manifold
 #   @assert dim[1] == dim[2]
    bi.invJ = inv(bi.J)
    @inbounds for i in 1:length(bi.dN)
        bi.grad[i] = bi.invJ ⋅ bi.dN[i]
    end
    bi.detJ = det(bi.J)
    #=
    elseif dim1 == 1 # curve
        bi.detJ = norm(bi.J)
    elseif dim1 == 2 # manifold
        bi.detJ = norm(cross(bi.J[1,:], bi.J[2,:]))
    end
    =#

    return bi
end

"""
    grad!(bi, gradu, u)

Evalute gradient ∂u/∂X and store result to matrix `gradu`. It is assumed
that `eval_basis!` has been already run to `bi` so it already contains
all necessary matrices evaluated with some `X` and `xi`.

# Example

First setup and evaluate basis using `eval_basis!`:
```jldoctest ex1
B = BasisInfo(Quad4)
X = Vec.([(0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0)])
xi = Vec(0.0, 0.0)
eval_basis!(B, X, xi)

# output

BasisInfo{Quad4,2,Float64}([0.25, 0.25, 0.25, 0.25], Tensors.Tensor{1,2,Float64,2}[[-0.25, -0.25], [0.25, -0.25], [0.25, 0.25], [-0.25, 0.25]], Tensors.Tensor{1,2,Float64,2}[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.5, 0.5]], [0.5 0.0; 0.0 0.5], [2.0 -0.0; -0.0 2.0], 0.25)

```

Next, calculate gradient of `u`:
```jldoctest ex1
u = Vec.([(0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0)])
grad(B, u)

# output

2×2 Tensors.Tensor{2,2,Float64,4}:
 1.5  0.5
 1.0  2.0

```
"""
function grad(bi::BasisInfo{B}, u::Vector{<:Vec{dim}}) where B <: AbstractBasis{dim} where dim
    gradu = zero(Tensor{2, dim})
    for k in 1:length(B)
        gradu += otimes(u[k], bi.grad[k])
    end
    return gradu
end


text:
module Diffusion

import MPI
using ImplicitGlobalGrid

@views interior(A) = A[2:end-1,2:end-1]

@views north(A) = A[1:end-2,2:end-1]
@views south(A) = A[3:end,  2:end-1]
@views east(A)  = A[2:end-1,1:end-2]
@views west(A)  = A[2:end-1,3:end]

const results_chan = Channel{Any}(1000)
const termination_event = Base.Event()

get_result() = take!(results_chan)
terminate() = notify(termination_event)

function diffusion(T, timesteps, comm)
    nx, ny = size(T)
    coords = MPI.Cart_coords(comm)
    me, dims = init_global_grid(nx, ny, 1; init_MPI=!MPI.Initialized(), comm=comm, reorder=0)  # Initialize the implicit global grid

    put!(results_chan, (coords, collect(interior(T))))

    for t in 1:timesteps
        interior(T) .-= ((north(T) .+ south(T) .+ east(T) .+ west(T)) ./ 4.0)
        update_halo!(T)

        if t % 2 == 0
            put!(results_chan, (coords, collect(interior(T))))
        end
    end

    @info "Finished diffusion" me
    wait(termination_event)
    @info "Worker terminated" me
end

end # module

max_stars_repo_name: JuliaLabs/DaggerDiffusion

max_issues_repo_name: JuliaLabs/DaggerDiffusion

max_forks_repo_name: JuliaLabs/DaggerDiffusion


text:
module DigilentWaveForms

import Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("DigilentWaveForms was not built properly. Please run Pkg.build(\"DigilentWaveForms\").")
end
include(depsjl_path)

include(joinpath(@__DIR__, "..", "gen", "ctypes.jl"))
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "libdwf_fixes.jl"))
include(joinpath(@__DIR__, "..", "gen", "libdwf_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "libdwf_consts.jl"))
include(joinpath(@__DIR__, "..", "gen", "libdwf_api.jl"))

include("utils.jl")
include("enumeration.jl")
include("system.jl")

include("powersupplies.jl")
include("analogin.jl")
include("analogout.jl")

include("digitalout.jl")

end # module


text:
"""
CloseStatusCode indicates a reason for closing the connection.
It is optionally sent as the first two bytes of a Close frames payload.
"""
struct CloseStatus
    code::UInt16
end

const CLOSE_STATUS_NORMAL                     = CloseStatus(1000)
const CLOSE_STATUS_GOING_AWAY                 = CloseStatus(1001)
const CLOSE_STATUS_PROTOCOL_ERROR             = CloseStatus(1002)
const CLOSE_STATUS_UNACCEPTABLE_DATA          = CloseStatus(1003)
const CLOSE_STATUS_RESERVED_1004              = CloseStatus(1004)
const CLOSE_STATUS_NO_STATUS                  = CloseStatus(1005)
const CLOSE_STATUS_ABNORMAL_CLOSE             = CloseStatus(1006)
const CLOSE_STATUS_INCONSISTENT_DATA          = CloseStatus(1007)
const CLOSE_STATUS_POLICY_VIOLATION           = CloseStatus(1008)
const CLOSE_STATUS_MESSAGE_TOO_BIG            = CloseStatus(1009)
const CLOSE_STATUS_EXPECTED_EXTENSION         = CloseStatus(1010)
const CLOSE_STATUS_FATAL_UNEXPECTED_CONDITION = CloseStatus(1011)
const CLOSE_STATUS_TLS_HANDSHAKE_FAILURE      = CloseStatus(1015)

struct CloseStatusAndReason
    status::CloseStatus
    reason::String
end

function readstatusandreason(frame::Frame)
    if length(frame.payload) >= 2
        buffer = IOBuffer(frame.payload)
        code = ntoh(read(buffer, UInt16))
        reason = ""
        if length(frame.payload) > 2
            reason = String(frame.payload[3:end])
        end
        CloseStatusAndReason(CloseStatus(code), reason)
    else
        CloseStatusAndReason(CLOSE_STATUS_NO_STATUS, "")
    end
end

# Requirement
# @5_5_1-4 No frames after Close frame
# @5_5_1-8 Endpoint closed
# @5_5_1-10 Client closes the connection
# @7_1_1-1 Closing the underlying connection
# @7_1_1-3 Close the connection non-cleanly when necessary
# @7_1_1-5 The client may close the TCP connection in abnormal circumstances
# @7_1_1-7 Client should wait for the server to close the TCP connection
# @7_1_2-2 Closing the connection
# @7_1_3-1 Closing handshake is started
# @7_1_4-1 Underlying connection is closed
# @7_1_4-2 Clean close
# @7_1_7-1 Clients fail the WebSocket connection
# @7_1_7-5 Connection established before fail
# @7_1_7-6 Omit the Close frame during fail
# @7_1_7-7 Processing data after WebSocket connection fail
# @7_2_1-1 Client initiated close on fail
# @7_4 Reason for closure
# @7_4_1-1 Pre-defined status codes
#
# By design, the connection behaviours do not send text or binary frames, even when requested by the
# user.

"""
Failing the WebSocket connection is an action taken at certain points in the protocol specification,
in response to error conditions.

The closing behaviour is to optionally send a Close frame, with an appropriate status code, and then
close the socket.
"""
struct FailTheConnectionBehaviour <: ClosingBehaviour
    framewriter::AbstractFrameWriter
    handler::WebSocketHandler
    status::CloseStatus
    issocketprobablyup::Bool
    reason::String

    FailTheConnectionBehaviour(w::AbstractFrameWriter,
                               handler::WebSocketHandler,
                               status::CloseStatus;
                               issocketprobablyup=true,
                               reason::String = "") = new(w, handler, status, issocketprobablyup, reason)
end

function closetheconnection(fail::FailTheConnectionBehaviour)
    if fail.issocketprobablyup
        sendcloseframe(fail.framewriter, fail.status; reason=fail.reason)
    end
    closesocket(fail.framewriter)
    state_closed(fail.handler)
end

clientprotocolinput(::FailTheConnectionBehaviour, ::ClientProtocolInput) = nothing

protocolstate(::FailTheConnectionBehaviour) = STATE_CLOSED

isclosedcleanly(::FailTheConnectionBehaviour) = false

closestatusandreason(::FailTheConnectionBehaviour) = CloseStatusAndReason(CLOSE_STATUS_ABNORMAL_CLOSE, "")

"""
Closing the WebSocket connection is a procedure for closing the connection during the normal course
the protocol lifetime.
"""
mutable struct ClientInitiatedCloseBehaviour <: ClosingBehaviour
    framewriter::AbstractFrameWriter
    handler::WebSocketHandler
    status::CloseStatus
    reason::String
    state::SocketState
    isclosereceived::Bool
    serverstatusandreason::Union{CloseStatusAndReason, Nothing}

    function ClientInitiatedCloseBehaviour(w::AbstractFrameWriter, handler::WebSocketHandler;
                                           status::CloseStatus = CLOSE_STATUS_NORMAL,
                                           reason::String = "")
        new(w, handler, status, reason, STATE_CLOSING, false, nothing)
    end
end

protocolstate(normal::ClientInitiatedCloseBehaviour) = normal.state

function closetheconnection(normal::ClientInitiatedCloseBehaviour)
    sendcloseframe(normal.framewriter, normal.status; reason = normal.reason)
    state_closing(normal.handler)
end

function clientprotocolinput(normal::ClientInitiatedCloseBehaviour, frame::FrameFromServer)
    if frame.frame.opcode == OPCODE_CLOSE
        if normal.state == STATE_CLOSING
            normal.isclosereceived = true
            if normal.serverstatusandreason == nothing
                normal.serverstatusandreason = readstatusandreason(frame.frame)
            end
        end
    end
end

function clientprotocolinput(normal::ClientInitiatedCloseBehaviour, ::SocketClosed)
    if normal.state == STATE_CLOSING
        normal.state = STATE_CLOSED
        state_closed(normal.handler)
    end
end

function clientprotocolinput(normal::ClientInitiatedCloseBehaviour, ::AbnormalNoCloseResponseReceived)
    closesocket(normal.framewriter)
end

clientprotocolinput(::ClientInitiatedCloseBehaviour, ::ClientProtocolInput) = nothing

isclosedcleanly(normal::ClientInitiatedCloseBehaviour) = normal.state == STATE_CLOSED && normal.isclosereceived

function closestatusandreason(normal::ClientInitiatedCloseBehaviour)
    if normal.serverstatusandreason == nothing
        CloseStatusAndReason(CLOSE_STATUS_ABNORMAL_CLOSE, "")
    else
        normal.serverstatusandreason
    end
end

"""
The server can initiate a Close, in which case this behaviour ensures a proper close.
"""
mutable struct ServerInitiatedCloseBehaviour <: ClosingBehaviour
    framewriter::AbstractFrameWriter
    handler::WebSocketHandler
    serverstatus::CloseStatus
    serverreason::String
    state::SocketState
    issocketclosedbyserver::Bool

    function ServerInitiatedCloseBehaviour(w::AbstractFrameWriter, h::WebSocketHandler, servercloseframe::Frame)
        statusandreason = readstatusandreason(servercloseframe)
        new(w, h, statusandreason.status, statusandreason.reason, STATE_CLOSING, false)
    end
end

protocolstate(b::ServerInitiatedCloseBehaviour) = b.state

function closetheconnection(behaviour::ServerInitiatedCloseBehaviour)
    state_closing(behaviour.handler)
    sendcloseframe(behaviour.framewriter, behaviour.serverstatus; reason = "")
end

function clientprotocolinput(behaviour::ServerInitiatedCloseBehaviour, ::SocketClosed)
    behaviour.state = STATE_CLOSED
    behaviour.issocketclosedbyserver = true
    state_closed(behaviour.handler)
end

clientprotocolinput(::ServerInitiatedCloseBehaviour, ::ClientProtocolInput) = nothing

function clientprotocolinput(b::ServerInitiatedCloseBehaviour, ::AbnormalSocketNotClosedByServer)
    b.state = STATE_CLOSED
    state_closed(b.handler)
    closesocket(b.framewriter)
end

isclosedcleanly(b::ServerInitiatedCloseBehaviour) = b.state == STATE_CLOSED && b.issocketclosedbyserver

closestatusandreason(b::ServerInitiatedCloseBehaviour) = CloseStatusAndReason(b.serverstatus, b.serverreason)


text:
using Test, PlotsPyPlot

plt = plot(sin)
@test typeof(plt) <: Plots.Plot{Plots.PyPlotBackend}


text:
#=

The "Logarithmic Method" in [wdht].

=#

export WdhtCoordinator, insert_node!, remove_node!, choose_nodes, num_nodes, node_weight, reseed!

# Make it mutable for reseeding
mutable struct WdhtNode
    id::String
    pos::Float64
    weight::Float64
end

"Weighted Distributed Hash Tables"
struct WdhtCoordinator <: RvCoordinator
    nodes::Dict{String, WdhtNode}

    WdhtCoordinator() = new(Dict{String, WdhtNode}())
end

function hash1f(x::String)::Float64
    h = hash(x)
    (h & MASK) / MAX_HASH
end

function insert_node!(coordinator::WdhtCoordinator, id::String, weight::Float64)
    node = WdhtNode(id, hash1f(id), weight)
    coordinator.nodes[id] = node
end

function remove_node!(coordinator::WdhtCoordinator, id::String)
    delete!(coordinator.nodes, id)
end

function choose_nodes(
    coordinator::WdhtCoordinator,
    data::String,
    n::Union{Int, Nothing}=nothing,
)::Vector{String}
    h = hash1f(data)
    nodes = [
        begin
            s = -log(mod(1.0 - (h - node.pos), 1.0)) / node.weight
            s => id
        end
        for (id, node) in coordinator.nodes
    ]

    n::Int = isnothing(n) ? length(nodes) : n
    last.(partialsort!(nodes, 1:n))
end

function num_nodes(coordinator::WdhtCoordinator)::Int
    length(coordinator.nodes)
end

function node_weight(coordinator::WdhtCoordinator, id::String)::Float64
    coordinator.nodes[id].weight
end

function reseed!(coordinator::WdhtCoordinator)
    for node in values(coordinator.nodes)
        node.pos = rand()
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using JSON3
using StructTypes
using Test

struct StringNumber
    value::String
end

@inline Base.getindex(x::StringNumber) = x.value::String

@inline StructTypes.StructType(::Type{StringNumber}) = JSON3.RawType()
@inline StructTypes.construct(::Type{StringNumber}, x::JSON3.RawValue) = StringNumber(unsafe_string(pointer(x.bytes, x.pos), x.len))
@inline JSON3.rawbytes(x::StringNumber) = codeunits(x.value)

@inline function test_roundtrip(json::String, object, ::Type{T}) where T
    # starting with object
    a = JSON3.write(object)
    @test typeof(a) === String
    b = JSON3.read(a, T)
    @test typeof(b) === typeof(object)
    c = JSON3.write(b)
    @test typeof(c) === String
    @test a === c

    # starting with json
    a = JSON3.read(json, T)
    @test typeof(a) == T
    @test typeof(a) == typeof(object)
    mutable_struct_equality(a, object)
    b = JSON3.write(a)
    @test typeof(b) === String
    @test b === json

    # @test false

    return nothing
end

@inline function mutable_struct_equality(a::T, b::T)::Bool where T
    if isstructtype(T) && ismutable(a)
        T_fieldnames = fieldnames(T)
        T_num_fieldnames = length(T_fieldnames)
        temp_vector = Vector{Bool}(undef, T_num_fieldnames)
        for i in 1:T_num_fieldnames
            name_i = T_fieldnames[i]
            temp_vector[i] = mutable_struct_equality(getfield(a, name_i), getfield(b, name_i))
        end
        return all(temp_vector)
    else
        return a == b
    end
end

mutable struct StringNumberTests_A{S}
    entries::Vector{S}
    StringNumberTests_A{S}()        where {S} = new{S}()
    StringNumberTests_A{S}(entries) where {S} = new{S}(entries)
end
mutable struct StringNumberTests_B{T}
    foo::T
    StringNumberTests_B{T}()    where {T} = new{T}()
    StringNumberTests_B{T}(foo) where {T} = new{T}(foo)
end
mutable struct StringNumberTests_C
    bar::StringNumber
    StringNumberTests_C() = new()
    StringNumberTests_C(bar) = new(bar)
end
StructTypes.StructType(::Type{<:StringNumberTests_A}) = StructTypes.Mutable()
StructTypes.StructType(::Type{<:StringNumberTests_B}) = StructTypes.Mutable()
StructTypes.StructType(::Type{<:StringNumberTests_C}) = StructTypes.Mutable()

@testset "JSON3.RawType: StringNumber" begin
    @testset "Basic tests" begin
        test_cases = ["1.2", "1.20", "1.200", "1.2000", "1.20000"]
        @testset "Reading" begin
            for a in test_cases
                b = JSON3.read(a, StringNumber)
                @test typeof(b) === StringNumber
                @test typeof(b[]) === String
                @test b[] === a
            end
        end
        @testset "Writing" begin
            for c in test_cases
                d = StringNumber(c)
                e = JSON3.write(d)
                @test typeof(e) === String
                @test e === c
            end
        end
        @testset "Roundtrip" begin
            for json in test_cases
                object = StringNumber(json)
                test_roundtrip(json, object, StringNumber)
            end
        end
    end
    @testset "More complex test" begin
        a = StringNumber("1.20")
        b = StringNumber("34.5600")
        c = StringNumberTests_C(a)
        d = StringNumberTests_C(b)
        e = StringNumberTests_B{StringNumberTests_C}(c)
        f = StringNumberTests_B{StringNumberTests_C}(d)
        g = [e, f]
        object = StringNumberTests_A{StringNumberTests_B{StringNumberTests_C}}(g)
        json = "{\"entries\":[{\"foo\":{\"bar\":1.20}},{\"foo\":{\"bar\":34.5600}}]}"
        T = StringNumberTests_A{StringNumberTests_B{StringNumberTests_C}}
        test_roundtrip(json, object, T)
    end
end


text:
any(pwd() .== LOAD_PATH) || push!(LOAD_PATH, pwd())
using Revise, SuperModule


text:
# Terra

function QuantumCircuit(x)
    qiskit.QuantumCircuit(x)
end

function QuantumCircuit(x, y)
    qiskit.QuantumCircuit(x, y)
end

function QuantumCircuit(x, y, n)
    qiskit.QuantumCircuit(x, y, n)
end


text:
"""

$(TYPEDEF)

Structures to contain the volumes obtained from calculations.

$(TYPEDFIELDS)

"""
@with_kw mutable struct Volume
  total::Float64
  bulk::Float64
  domain::Float64
  shell::Vector{Float64}
end

Volume(nbins::Int) = Volume( 0., 0., 0.,  zeros(Float64,nbins) )

function reset!(v::Volume)
  v.total = 0.
  v.bulk = 0.
  v.domain = 0.
  @. v.shell = 0.
  return nothing
end

#function Base.show(io::IO, v::Volume)
#  n = length(v.shell)
#  println(" Mean total box volume: $(v.total) ")
#  println(" Mean bulk volume: $(v.bulk) ")
#  println(" Mean solute domain volume: $(v.domain) ")
#  println(" Volumes of first, medium, and last solvation shells: $(v.shell[1]), $(v.shell[round(Int,n/2)]), $(v.shell[n])")
#end


text:
@testset "test_promotion" begin
G = tf(1)
sys = ss(1)
TG = typeof(G)
Tsys = typeof(sys)
@test promote_type(Tsys,TG) == Tsys

G = tf(1.)
sys = ss(1.)
TG = typeof(G)
Tsys = typeof(sys)
@test promote_type(Tsys,TG) == Tsys

G = tf(1.)
sys = ss(1)
TG = typeof(G)
Tsys = typeof(sys)
@test promote_type(Tsys,TG) == typeof(ss(1.))

G = zpk(1)
sys = ss(1)
TG = typeof(G)
Tsys = typeof(sys)
@test promote_type(Tsys,TG) == typeof(ss(1))

G1 = zpk(1)
G2 = zpk(1.)
TG1 = typeof(G1)
TG2 = typeof(G2)
@test promote_type(TG1,TG2) == TG2

@test promote_type(typeof(zpk(1)), typeof(zpk(1))) == typeof(zpk(1))
@test promote_type(typeof(tf(1)), typeof(tf(1))) == typeof(tf(1))
@test promote_type(typeof(ss(1)), typeof(ss(1))) == typeof(ss(1))

@test promote(ss(1), 1.0) == (ss(1.), ss(1.))
@test promote(ss(1), 1) == (ss(1.), ss(1.))
@test promote(ss(1), ss(1.)) == (ss(1.), ss(1.))
@test promote(ss(1), tf(1.)) == (ss(1.), ss(1.))
@test promote(ss(1), zpk(1.)) == (ss(1.), ss(1.))
@test promote(zpk(1), tf(1.)) == (zpk(1.), zpk(1.))
@test promote(tf([1 + im]), tf([1.])) == (tf([1. + im]), tf([1. + 0*im]))
@test promote(zpk(1), tf(1)) == (zpk(1.), zpk(1.))
@test promote(ss(1), ss([1 .*im])) == (ss([1+0*im]),ss([1 .*im]))
@test promote(ss(1), tf([1 .*im])) == (ss([1+0*im]),ss([1 .*im]))
@test promote(tf(1), 1) == (tf(1.), tf(1.))

end


text:
@testset "html partition" begin
    parts = """
        ABC1 <!-- DEF --> ABC2 <script> GHI </script> ABC3 {{ JKL }} ABC4
        """ |> FP.html_partition
    for (i, c) in enumerate(("ABC1", "DEF", "ABC2", "> GHI", "ABC3", "JKL", "ABC4"))
        @test isapproxstr(FP.content(parts[i]), c)
    end
end

@testset "script content" begin
    parts = """
        ABC <script> GHI </script> DEF
        """ |> FP.html_partition
    c = FP.content(parts[2])
    i = findfirst('>', c)
    @test strip(c[nextind(c, i):end]) == "GHI"
end

@testset "math content" begin
    parts = """
        \\(x\\)
        """ |> FP.html_partition
    parts[1] // "\\(x\\)"
    parts = """
        \\[x\\]
        """ |> FP.html_partition
    parts[1] // "\\[x\\]"
end


text:
@testset "Modularity" begin
    n = 10
    m = n*(n-1)/2
    c = ones(Int, n)
    gint = CompleteGraph(n)
    for g in testgraphs(gint)
      @test @inferred(modularity(g, c)) == 0
    end

    gint = Graph(n)
    for g in testgraphs(gint)
      @test @inferred(modularity(g, c)) == 0
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
module RandomExtensions

export Combine, Uniform, Normal, Exponential, CloseOpen, Rand

import Random: Sampler, rand, rand!

using Random
using Random: GLOBAL_RNG, SamplerTrivial, SamplerSimple, SamplerTag, Repetition

using SparseArrays: sprand, sprandn

include("distributions.jl")
include("sampling.jl")
include("containers.jl")
include("iteration.jl")

## updated rand docstring (TODO: replace Base's one)

"""
    rand([rng=GLOBAL_RNG], [S], [C...]) # RandomExtensions

Pick a random element or collection of random elements from the set of values specified by `S`;
`S` can be

* an indexable collection (for example `1:n` or `['x','y','z']`),
* an `AbstractDict` or `AbstractSet` object,
* a string (considered as a collection of characters), or
* a type: the set of values to pick from is then equivalent to `typemin(S):typemax(S)` for
  integers (this is not applicable to [`BigInt`](@ref)), and to ``[0, 1)`` for floating
  point numbers;
* a `Distribution` object, e.g. `Normal()` for a normal distribution (like `randn()`),
  or `CloseOpen(10.0, 20.0)` for uniform `Float64` numbers in the range ``[10.0, 20.0)``;
* a `Combine` object, which can be either `Combine(Pair, S1, S2)` or `Combine(Complex, S1, S2)`,
  where `S1` and `S2` are one of the specifications above; `Pair` or `Complex` can optionally be
  given as concrete types, e.g. `Combine(ComplexF64, 1:3, Int)` to generate `ComplexF64` instead
  of `Complex{Int}`.

`S` usually defaults to [`Float64`](@ref).

If `C...` is not specified, `rand` produces a scalar. Otherwise, `C...` can be:

* a set of integers, or a tuple of `Int`, which specify the dimensions of an `Array` to generate;
* `(p::AbstractFloat, m::Integer, [n::Integer])`, which produces a sparse array of dimensions `(m, n)`,
  in which the probability of any element being nonzero is independently given by `p`
* `(String, [n=8])`, which produces a random `String` of length `n`; the generated string consists of `Char`
  taken from a predefined set like `randstring`, and can be specified with the `S` parameter.
* `(Dict, n)`, which produces a `Dict` of length `n`; `S` must then specify the type of its elements,
  e.g. `Combine(Pair, Int, 2:3)`;
* `(Set, n)`, which produces a `Set` of length `n`;
* `(BitArray, dims...)`, which produces a `BitArray` with the specified dimensions.

# Examples
```julia-repl
julia> rand(Int, 2)
2-element Array{Int64,1}:
 1339893410598768192
 1575814717733606317

julia> rand(MersenneTwister(0), Dict(1=>2, 3=>4))
1=>2

julia> rand("abc", String, 12)
"bccaacaabaac"

julia> rand(1:10, Set, 3)
Set([3, 8, 6])
```

!!! note
    The complexity of `rand(rng, s::Union{AbstractDict,AbstractSet})`
    is linear in the length of `s`, unless an optimized method with
    constant complexity is available, which is the case for `Dict`,
    `Set` and `BitSet`. For more than a few calls, use `rand(rng,
    collect(s))` instead, or either `rand(rng, Dict(s))` or `rand(rng,
    Set(s))` as appropriate.
"""
rand

end # module


text:
"""
    getimports(imports::Dict)

Get a .wat string with any used functions that are not builtins or userdefined.

# Details
- a few basic are builtin to wasm. these can be translated.
- other basic functions are bultin to JavaScripts global Math object, these can be imported
- many more are builtin to julia... so they need to be implemented either in .jl, .wat or .js
- give warning if the function cant be imported from js Math.
"""
function getimports(imports::Dict)
    watimports = []
    jsimports = []
    for func in keys(imports)
        push!(watimports, imports[func][1])
        push!(jsimports, imports[func][2])
    end

    if length(jsimports) > 0
        println("INFO: Suggested import object when instantiating WebAssembly from JavaScript:")
        println(jsimportsstring(jsimports))
        println()
    end
    return join(watimports, "\n")
end

function jsimportsstring(jsimports)
    str = ["imports = {\n"]
    for func in jsimports
        push!(str, "  ")
        push!(str, func)
        push!(str, ",\n")
    end
    push!(str, "}")
    return join(str)
end

"""
    jsimportentry(func, argtypes)

Get a string of possible javascript Math module import, assuming it exists in the Math module.

# Details
from a func such as sin, return a string like
\"sin: (x) => Math.sin(x)\"
"""
function jsimportentry(func, argtypes)
    Nparams = length(argtypes)

    # the global Math object in js have the same function names as julia has, except for pow and random
    # ofcourse it may not exist in js, but this allows the ones that do to be imported with correct name
    if func == :(^)
        jsfunc = "pow"
    elseif func == :(rand)
        jsfunc = "random"
    else
        jsfunc = func
    end

    jskey = "\"$(func)\": "
    jsval = " => Math.$(jsfunc)"
    paramnames = []
    for i = 1:Nparams
        push!(paramnames, index2abc(i))
    end
    pn = join(["("; join(paramnames, ", "); ")"])

    return join([jskey,pn,jsval,pn])
end

index2abc(i) = '`' + i # unicode char, '`'+1 means a, '`'+2 means b

isgenericfunction(func) = length(methods(Base.eval(Evalscope, func)).ms) > 0
imports!(imports::Dict, ci::CodeInfo, funcs::Dict, builtinfuncs::Dict, item) = nothing
function imports!(imports::Dict, ci::CodeInfo, funcs::Dict, builtinfuncs::Dict, items::Array)
    jsglobalMath = ["^","rand","acos","acosh","asin","asinh","atan","atanh","atan2","cbrt","cos","cosh","exp","expm1","hypot","imul","log","log1p","log10","log2","sign","sin","sinh","tan","tanh","trunc"]

    if isimport(funcs, builtinfuncs, items[1])
        func = items[1]
        if isgenericfunction(func)
            argtypes = itemtype.((ci,), items[2:end])
            ct = code_typed(Base.eval(Evalscope, func), Tuple{argtypes...}; optimize=false, debuginfo=:none)[1]
            cinfo = ct[1]
            Rtype = ct[2]
            imports[func] = [importdeclaration(func.name, argtypes, Rtype), jsimportentry(func.name, argtypes)]

            isavailableinjavascript = string(func.name) in jsglobalMath
            if !isavailableinjavascript
                println("INFO: JavaScript does not have Math.$(func.name). Import your own function instead or write it in julia.")
            end
        end
    end
    imports!.((imports,), (ci,), (funcs,), (builtinfuncs,), items)
end

isimport(funcs, builtinfuncs, item) = false
isimport(funcs, builtinfuncs, item::GlobalRef) = !haskey(builtinfuncs, item.name) && !haskey(funcs, item.name) && !haskey(floatops, item.name) && !haskey(intops, item.name)

function importdeclaration(func, argtypes, Rtype)
    Nparams = length(argtypes)
    paramtypes = [paramtype <: AbstractFloat ? "f32" : "i32" for paramtype in argtypes]

    decl = ["func \$$func (import \"imports\" \"$func\")"]
    for i = 1:Nparams
        push!(decl, "(param \$$('`'+i) $(paramtypes[i]))")
    end
    if !(Rtype <: Nothing) && Rtype <: Number
        rt = Rtype <: AbstractFloat ? "f32" : "i32"
        push!(decl, "(result $rt)")
    end
    decl = join(decl," ")
    return "($decl)"
end


text:
using DataStructures
using Graphs
# using Statistics
# using Combinatorics

lines = readlines("inputs/11-15/12ex.txt")
# lines = readlines("inputs/11-15/12.txt")
# nums = map(str -> [parse(Int,str[l]) for l in findall(r"[[:digit:]]+", str)], lines)
# tallNums(lines) = [parse(Int64, x) for x in lines]
# A = zeros(Int, 23, 23)

global m = DefaultDict([])

global isBig = DefaultDict(false)

global caves = ["start","end"]

for line in lines
  (a, b) =  split(line,'-')
  for x in [a,b]
    if x ∉ caves
      push!(caves, x)
    end
    # determine big
    if x ∉ ["start","end"]
      if x == uppercase(x)
        isBig[x] = true
      else
        isBig[x] = false
      end
    end
  push!(m[a],b)
  push!(m[b],a)
  end
end

function cpaths()
  v = DefaultDict(false)
  pathCount = [0]
  cpathsuntil("start",v,pathCount)
  pathCount[1]
end

function cpathsuntil(a,v,pathCount)
  @show a, v, pathCount[1]
  if !isBig[a]
    v[a] = true
  end
  if a == "end"
    pathCount[1] += 1
  else
    for neighbor in m[a]
      if !v[neighbor]
        cpathsuntil(neighbor,v,pathCount)
      end
    end
  end
  v[a] = false
end

cpaths()

adj = DefaultDict([])
for l in lines
  a,b = split(l,'-')
  push!(adj[a], b)
  push!(adj[b], a)
end

small(x) = x != uppercase(x)

function cp(adj, start, goal, v=DefaultDict(false), twice = false)
  (start == goal) && return 1
  t = 0
  for ne in adj[start]

    if small(ne) && v[ne]
      if a in ["start",]
end