text:
using OffsetArrays
#include("copyq1.jl")
include("step1.jl")

#
#     ================================================================================================================
function claw1(parms::CParam,meqn::Int, mwaves::Int, maux::Int, mbc::Int, mx::Int,
               q::OffsetArray{Float64,2}, aux::OffsetArray{Float64,2},
               xlower::Float64, dx::Float64, tstart::Float64, tend::Float64,
               dtv::Array{Float64,1}, cflv::Array{Float64,1},
               nv::Array{Int,1}, method::Array{Int,1}, mthlim::Array{Int,1}, mthbc::Array{Int,1},
               #work,mwork, # ???
               use_fwave::Bool,
               # higher-order functions
               bc1, rp1, src1, b4step1)

#     ================================================================================================================
#
#  Solves a hyperbolic system of conservation laws in one space dimension
#  of the general form
#
#     capa * q_t + A q_x = psi
#
#  The "capacity function" capa(x) and source term psi are optional
#  (see below).
#
#  For a more complete description see the documentation at
#      http://www.amath.washington.edu/~claw
#
#  Sample driver programs and user-supplied subroutines are available.
#  See the the directories claw/clawpack/1d/example* for some examples, and
#  codes in claw/applications for more extensive examples.
#
#  --------------------------------------------------------
#
#  The user must supply the following subroutines:
#
#    bc1, rp1        subroutines specifying the boundary conditions and
#                    Riemann solver.
#                    These are described in greater detail below.
#
#    b4step1            The routine b4step1 is called each time step and
#                       can be supplied by the user in order to perform
#                       other operations that are necessary every time
#                       step.  For example, if the variables stored in
#                       the aux arrays are time-dependent then these
#                       values can be set.
#
#  In addition, if the equation contains source terms psi, then the user
#  must provide:
#
#    src1               subroutine that solves capa * q_t = psi
#                       over a single time step.
#
#  These routines must be declared EXTERNAL in the main program.
#  For description of the calling sequences, see below.
#
#  Dummy routines b4step1.f and src1.f are available in
#       claw/clawpack/1d/lib
#
#
#
#  Description of parameters...
#  ----------------------------
#
#
#    meqn is the number of equations in the system of
#         conservation laws.
#
#    mwaves is the number of waves that result from the
#           solution of each Riemann problem.  Often mwaves = meqn but
#           for some problems these may be different.
#
#    mbc is the number of "ghost cells" that must be added on to each
#       side of the domain to handle boundary conditions.  The cells
#       actually in the physical domain are labelled from 1 to mx in x.
#       The arrays are dimensioned actually indexed from 1-mbc to mx+mbc.
#       For the methods currently implemented, mbc = 2 should be used.
#       If the user implements another method that has a larger stencil and
#       hence requires more ghost cells, a larger value of mbc could be used.
#       q is extended from the physical domain to the ghost cells by the
#       user-supplied routine bc1.
#
#    mx is the number of grid cells in the x-direction, in the
#       physical domain.  In addition there are mbc grid cells
#       along each edge of the grid that are used for boundary
#       conditions.
#
#    q(meqn, 1-mbc:mx+mbc)
#        On input:  initial data at time tstart.
#        On output: final solution at time tend.
#        q(m,i) = value of mth component in the i'th cell.
#        Values within the physical domain are in q(m,i)
#                for i = 1,2,...,mx
#        mbc extra cells on each end are needed for boundary conditions
#        as specified in the routine bc1.
#
#    aux(maux, 1-mbc:mx+mbc)
#        Array of auxiliary variables that are used in specifying the problem.
#        If method(7) = 0 then there are no auxiliary variables and aux
#                         can be a dummy variable.
#        If method(7) = maux > 0 then there are maux auxiliary variables
#                         and aux must be dimensioned as above.
#
#        Capacity functions are one particular form of auxiliary variable.
#        These arise in some applications, e.g. variable coefficients in
#        advection or acoustics problems.
#        See Clawpack Note # 5 for examples.
#
#        If method(6) = 0 then there is no capacity function.
#        If method(6) = mcapa > 0  then there is a capacity function and
#            capa(i), the "capacity" of the i'th cell, is assumed to be
#            stored in aux(mcapa,i).
#            In this case we require method(7).ge.mcapa.
#
#    dx = grid spacing in x.
#         (for a computation in ax <= x <= bx,  set dx = (bx-ax)/mx.)
#
#    tstart = initial time.
#
#    tend = Desired final time (on input).
#              If tend<tstart, then claw1 returns after a single successful
#                 time step has been taken (single-step mode).
#              Otherwise, as many steps are taken as needed to reach tend,
#                 up to a maximum of nv(1).
#         = Actual time reached (on output).
#
#    dtv(1:5) = array of values related to the time step:
#               (Note: method(1)=1 indicates variable size time steps)
#         dtv(1) = value of dt to be used in all steps if method(1) = 0
#                = value of dt to use in first step if method(1) = 1
#         dtv(2) = unused if method(1) = 0.
#                = maximum dt allowed if method(1) = 1.
#         dtv(3) = smallest dt used (on output)
#         dtv(4) = largest dt used (on output)
#         dtv(5) = dt used in last step (on output)
#
#    cflv(1:4) = array of values related to Courant number:
#         cflv(1) = maximum Courant number to be allowed.  With variable
#                   time steps the step is repeated if the Courant
#                   number is larger than this value.  With fixed time
#                   steps the routine aborts.  Usually cflv(1)=1.0
#                   should work.
#         cflv(2) = unused if method(1) = 0.
#                 = desired Courant number if method(1) = 1.
#                   Should be somewhat less than cflv(1), e.g. 0.9
#         cflv(3) = largest Courant number observed (on output).
#         cflv(4) = Courant number in last step (on output).
#
#    nv(1:2) = array of values related to the number of time steps:
#         nv(1) = unused if method(1) = 0
#               = maximum number of time steps allowed if method(1) = 1
#         nv(2) = number of time steps taken (on output).
#
#    method(1:7) = array of values specifying the numerical method to use
#         method(1) = 0 if fixed size time steps are to be taken.
#                       In this case, dt = dtv(1) in all steps.
#                   = 1 if variable time steps are to be used.
#                       In this case, dt = dtv(1) in the first step and
#                       thereafter the value cflv(2) is used to choose the
#                       next time step based on the maximum wave speed seen
#                       in the previous step.  Note that since this value
#                       comes from the previous step, the Courant number will
#                       not in general be exactly equal to the desired value
#                       If the actual Courant number in the next step is
#                       greater than 1, then this step is redone with a
#                       smaller dt.
#
#         method(2) = 1 if Godunov's method is to be used, with no 2nd order
#                       corrections.
#                   = 2 if second order correction terms are to be added, with
#                       a flux limiter as specified by mthlim.
#
#         method(3)  is not used in one-dimension.
#
#         method(4) = 0 to suppress printing
#                   = 1 to print dt and Courant number every time step
#
#         method(5) = 0 if there is no source term psi.  In this case
#                       the subroutine src1 is never called so a dummy
#                       parameter can be given.
#                   = 1 if there is a source term.  In this case
#                       the subroutine src1 must be provided and a
#                       fractional step method is used.
#                       In each time step the following sequence is followed:
#                            call bc to extend data to ghost cells
#                            call step1 to advance hyperbolic eqn by dt
#                            call src1 to advance source terms by dt
#                   = 2 if there is a source term and Strang splitting is to
#                       be used instead of the Godunov splitting above.
#                       In each time step the following sequence is followed:
#                            call bc to extend data to ghost cells
#                            call src1 to advance source terms by dt/2
#                            call step1 to advance hyperbolic equation by dt
#                            call src1 to advance source terms by dt/2
#                       For most problems 1 is recommended rather than 2
#                       since it is less expensive and works essentially as
#                       well on most problems.

#
#         method(6) = 0 if there is no capacity function capa.
#                   = mcapa > 0 if there is a capacity function.  In this case
#                       aux(i,mcapa) is the capacity of the i'th cell and you
#                       must also specify method(7) .ge. mcapa and set aux.
#
#         method(7) = 0 if there is no aux array used.
#                   = maux > 0  if there are maux auxiliary variables.
#
#
#         The recommended choice of methods for most problems is
#            method(1) = 1,  method(2) = 2.
#
#    mthlim(1:mwaves) = array of values specifying the flux limiter to be used
#                     in each wave family mw.  Often the same value will be used
#                     for each value of mw, but in some cases it may be
#                     desirable to use different limiters.  For example,
#                     for the Euler equations the superbee limiter might be
#                     used for the contact discontinuity (mw=2) while another
#                     limiter is used for the nonlinear waves.  Several limiters
#                     are built in and others can be added by modifying the
#                     subroutine philim.
#
#        mthlim(mw) = 0 for no limiter
#                   = 1 for minmod
#                   = 2 for superbee
#                   = 3 for van Leer
#                   = 4 for monotonized centered
#
#
#    work(mwork) = double precision work array of length at least mwork
#
#    mwork = length of work array.  Must be at least
#               (mx + 2*mbc) * (2 + 4*meqn + mwaves + meqn*mwaves)
#            If mwork is too small then the program returns with info = 4
#            and prints the necessary value of mwork to unit 6.
#
#
#    info = output value yielding error information:
#         = 0 if normal return.
#         = 1 if mbc.lt.2
#         = 2 if method(1)=0 and dt doesn't divide (tend - tstart).
#         = 3 if method(1)=1 and cflv(2) > cflv(1).
#         = 4 if mwork is too small.
#         = 11 if the code attempted to take too many time steps, n > nv(1).
#              This could only happen if method(1) = 1 (variable time steps).
#         = 12 if the method(1)=0 and the Courant number is greater than 1
#              in some time step.
#
#           Note: if info.ne.0, then tend is reset to the value of t actually
#           reached and q contains the value of the solution at this time.
#
#    User-supplied subroutines
#    -------------------------
#
#    bc1 = subroutine that specifies the boundary conditions.
#         This subroutine should extend the values of q from cells
#         1:mx to the mbc ghost cells along each edge of the domain.
#
#          The form of this subroutine is
#  -------------------------------------------------
#     subroutine bc1(meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt,mthbc)
#     implicit double precision (a-h,o-z)
#     dimension   q(meqn, 1-mbc:mx+mbc)
#     dimension aux(maux, 1-mbc:mx+mbc)
#     dimension mthbc(2)
#  -------------------------------------------------
#
#    The routine claw/clawpack/1d/lib/bc1.f can be used to specify
#    various standard boundary conditions.
#
#
#    rp1 = user-supplied subroutine that implements the Riemann solver
#
#          The form of this subroutine is
#  -------------------------------------------------
#     subroutine rp1(meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
#     implicit double precision (a-h,o-z)
#     dimension   ql(meqn, 1-mbc:mx+mbc)
#     dimension   qr(meqn, 1-mbc:mx+mbc)
#     dimension auxl(maux, 1-mbc:mx+mbc)
#     dimension auxr(maux, 1-mbc:mx+mbc)
#     dimension wave(meqn, mwaves, 1-mbc:mx+mbc)
#     dimension    s(mwaves, 1-mbc:mx+mbc)
#     dimension amdq(meqn, 1-mbc:mx+mbc)
#     dimension apdq(meqn, 1-mbc:mx+mbc)
#  -------------------------------------------------
#
#         On input, ql contains the state vector at the left edge of each cell
#                   qr contains the state vector at the right edge of each cell
#                 auxl contains auxiliary values at the left edge of each cell
#                 auxr contains auxiliary values at the right edge of each cell
#
#         Note that the i'th Riemann problem has left state qr(:,i-1)
#                                            and right state ql(:,i)
#         In the standard clawpack routines, this Riemann solver is
#         called with ql=qr=q along this slice.  More flexibility is allowed
#         in case the user wishes to implement another solution method
#         that requires left and rate states at each interface.

#         If method(7)=maux > 0 then the auxiliary variables along this slice
#         are passed in using auxl and auxr.  Again, in the standard routines
#         auxl=auxr=aux in the call to rp1.
#
#          On output,
#              wave(m,mw,i) is the m'th component of the jump across
#                              wave number mw in the ith Riemann problem.
#              s(mw,i) is the wave speed of wave number mw in the
#                              ith Riemann problem.
#              amdq(m,i) = m'th component of A^- Delta q,
#              apdq(m,i) = m'th component of A^+ Delta q,
#                     the decomposition of the flux difference
#                         f(qr(i-1)) - f(ql(i))
#                     into leftgoing and rightgoing parts respectively.
#
#           It is assumed that each wave consists of a jump discontinuity
#           propagating at a single speed, as results, for example, from a
#           Roe approximate Riemann solver.  An entropy fix can be included
#           into the specification of amdq and apdq.
#
#    src1 = subroutine for the source terms that solves the equation
#               capa * q_t = psi
#           over time dt.
#
#           If method(5)=0 then the equation does not contain a source
#           term and this routine is never called.  A dummy argument can
#           be used with many compilers, or provide a dummy subroutine that
#           does nothing (such a subroutine can be found in
#           claw/clawpack/1d/lib/src1.f)
#
#          The form of this subroutine is
#  -------------------------------------------------
#     subroutine src1(meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt)
#     implicit double precision (a-h,o-z)
#     dimension   q(meqn, 1-mbc:mx+mbc)
#     dimension aux(maux, 1-mbc:mx+mbc)
#  -------------------------------------------------
#      If method(7)=0  or the auxiliary variables are not needed in this solver,
#      then the latter dimension statement can be omitted, but aux should
#      still appear in the argument list.
#
#      On input, q(m,i) contains the data for solving the
#                source term equation.
#      On output, q(m,i) should have been replaced by the solution to
#                 the source term equation after a step of length dt.
#
#
#      b4step1 = subroutine that is called from claw1 before each call to
#                step1.  Use to set time-dependent aux arrays or perform
#                other tasks which must be done every time step.
#
#          The form of this subroutine is
#
#  -------------------------------------------------
#      subroutine b4step1(mbc,mx,meqn,q,xlower,dx,time,dt,maux,aux)
#      implicit double precision (a-h,o-z)
#      dimension   q(meqn, 1-mbc:mx+mbc)
#      dimension aux(maux, 1-mbc:mx+mbc)
#  -------------------------------------------------
#
#
#
#
# =========================================================================
#
#  Copyright 1994 -- 2002 R. J. LeVeque
#
#  This software is made available for research and instructional use only.
#  You may copy and use this software without charge for these non-commercial
#  purposes, provided that the copyright notice and associated text is
#  reproduced on all copies.  For all other uses (including distribution of
#  modified versions), please contact the author at the address given below.
#
#  *** This software is made available "as is" without any assurance that it
#  *** will work for your purposes.  The software may in fact have defects, so
#  *** use the software at your own risk.
#
#  --------------------------------------
#    CLAWPACK Version 4.1,  August, 2002
#    Webpage: http://www.amath.washington.edu/~claw
#  --------------------------------------
#    Author:  Randall J. LeVeque
#             Applied Mathematics
#             Box 352420
#             University of Washington,
#             Seattle, WA 98195-2420
#             rjl@amath.washington.edu
# =========================================================================
#
#
#
#
#    ======================================================================
#    Beginning of claw1 code
#    ======================================================================
#
#      implicit double precision (a-h,o-z)
#      external bc1,rp1,src1,b4step1
#      dimension q(meqn,1-mbc:mx+mbc)
#      dimension aux(maux,1-mbc:mx+mbc)
#      dimension work(mwork)
#      dimension mthlim(mwaves),method(7),dtv(5),cflv(4),nv(2)
#      dimension mthbc(2)
#      common /comxt/ dtcom,dxcom,tcom
#
#

    info::Int = 0
    q_save = OffsetArray{Float64}(undef, 1:meqn, 1-mbc:mx+mbc)

    info = 0
    t = tstart
    maxn = nv[1]
    dt = dtv[1]   # initial dt
    cflmax = 0.0
    dtmin = dt
    dtmax = dt
    nv[2] = 0

    # check for errors in data:

    if method[1] == 0
        # fixed size time steps.  Compute the number of steps:
        if tend < tstart
            # single step mode
            maxn = 1
        else
            maxn = convert(Uint, (tend - tstart)/dt)    # Round to nearest int
            if abs(maxn*dt - (tend-tstart)) >
                        1e-5*(tend-tstart)
                # dt doesn't divide time interval integer number of times
                info = 2
                # [local goto #101](https://github.com/JuliaLang/julia/issues/101)
                @goto l900
            end
        end
    end

    if method[1] == 1 && cflv[2] > cflv[1]
        info = 3
        @goto l900
    end

    #
    #     -----------
    #     # main loop
    #     -----------
    #

    if maxn == 0
        return
    end

    f    = OffsetArray{Float64}(undef, 1:meqn,              1-mbc:mx+mbc)
    s    = OffsetArray{Float64}(undef,            1:mwaves, 1-mbc:mx+mbc)
    wave = OffsetArray{Float64}(undef, 1:meqn,    1:mwaves, 1-mbc:mx+mbc)
    amdq = OffsetArray{Float64}(undef, 1:meqn,              1-mbc:mx+mbc)
    apdq = OffsetArray{Float64}(undef, 1:meqn,              1-mbc:mx+mbc)
    dtdx = OffsetArray{Float64}(undef,                      1-mbc:mx+mbc)

    cfl::Float64 = 0.0
    for n=1:maxn
        told = t   # time at beginning of time step.

        # adjust dt to hit tend exactly if we're near end of computation
        #  (unless tend < tstart, which is a flag to take only a single step)
        if told+dt > tend && tstart < tend
            dt = tend - told
        end

        if method[1] == 1
            # save old q in case we need to retake step with smaller dt:
            #copyq1(meqn,mbc,mx,q,work(i0qwork))
            q_save = q
        end

    #
    @label l40

        dt2 = dt / 2.0
        #println("dt2: $dt2")
        thalf = t + dt2  # midpoint in time for Strang splitting
        t = told + dt    # time at end of step

        # store dt and t in the common block comxt in case they are needed
        # in the Riemann solvers (for variable coefficients)
        tcom = told
        dtcom = dt
        dxcom = dx

        # ------------------------------------------------------------------
        # main steps in algorithm:
        # ------------------------------------------------------------------

        # extend data from grid to bordering boundary cells:
        bc1(meqn,mbc,mx,xlower,dx,q,maux,aux,told,dt,mthbc)

        # call user-supplied routine which might set aux arrays
        # for this time step, for example.

        b4step1(mbc,mx,meqn,q,xlower,dx,told,dt,maux,aux)

        if method[5] == 2
            # with Strang splitting for source term:
            src1(meqn,mbc,mx,xlower,dx,q,maux,aux,told,dt2)
        end

        # take a step on the homogeneous conservation law:

        cfl =
        step1(parms,meqn,mwaves,mbc,maux,mx,q,aux,dx,dt,method,mthlim,
              f,wave,s,amdq,apdq,dtdx,
              use_fwave,rp1)

        if method[5] == 2
            # source terms over a second half time step for Strang splitting:
            # Note it is not so clear what time t should be used here if
            # the source terms are time-dependent!
            src1(meqn,mbc,mx,xlower,dx,q,maux,aux,thalf,dt2)
        end

        if method[5] == 1
            # source terms over a full time step:
            src1(meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt)
        end

        #
        #        ------------------------------------------------------------------
        #
        if method[4] == 1
            @printf("CLAW1... Step %4d\n   Courant number =',%6.3f\n  dt =%12.4e,\n  t =',%12.4e\n",n,cfl,dt,t)
        end

        if method[1] == 1
            # choose new time step if variable time step
            if cfl > 0.0
                dt = min(dtv[2], dt * cflv[2]/cfl)
                dtmin = min(dt,dtmin)
                dtmax = max(dt,dtmax)
            else
                dt = dtv[2]
            end
        end

        # check to see if the Courant number was too large:

        if cfl <= cflv[1]
            # accept this step
            cflmax = max(cfl,cflmax)
            # as the step accepted exit the loop
            #println("ACCEPT: cfl : $cfl cflv[1] : $(cflv[1])")
        else
            #println("REJECT: cfl : $cfl cflv[1] : $(cflv[1])")
            # reject this step
            t = told
            #copyq1(meqn,mbc,mx,work(i0qwork),q)
            q = q_save

            if method[4] == 1
                println("CLAW1 rejecting step... ")
                println("Courant number too large")
            end

            if method[1] == 1
                # if variable dt, go back and take a smaller step
                @goto l40
            else # if fixed dt, give up and return
                cflmax = max(cfl,cflmax)
                break
            end

        end

        # see if we are done:
        nv[2] += 1

        println("t : $t tend : $tend")
        if t >= tend
            break
        end

    end

    #
@label l900

    # return information

    if method[1] == 1 && t < tend && nv[2] == maxn
        # too many timesteps
        info = 11
    end

    if method[1] == 0 && cflmax > cflv[1]
        # Courant number too large with fixed dt
        info = 12
    end

    tend = t
    cflv[3] = cflmax
    cflv[4] = cfl
    dtv[3] = dtmin
    dtv[4] = dtmax
    dtv[5] = dt

    return info

end # claw1


text:
@doc raw"""
Abstract type for indicator functions $\mathbb{R}^{2n} \rightarrow
\left\{0, 1\right\}$ where $n = 1, 2 \text{ or } 3$.
"""
abstract type AbstractIndicator end

@doc raw"""
    SeparableIndicator(χ₁, χ₂)

Type for separable indicator function, that is for such an indicator
function which can be written as $\chi(x,y) = \chi_1(x)\chi_2(y)$.

`χ1` and `χ2` must be functions of one argument which return a value
of `Bool` type.

**NB**: This indicator function is not symmetric (i.e. $\chi(x,y) \ne
\chi(y,x)$). This behaviour is intentional. For example you can write
such an indicator, so the corresponding correlation function is
sensitive to the spatial orientation of a system.

*"That one, too fat! This one, too tall! This one… too symmetrical!"*
"""
struct SeparableIndicator <: AbstractIndicator
    χ1 :: Function
    χ2 :: Function
end

"""
    SeparableIndicator(χ)

This is equvalent of writing `SeparableIndicator(χ, χ)`.
"""
SeparableIndicator(χ :: Function) = SeparableIndicator(χ, χ)

@doc raw"""
    InseparableIndicator(χ)

Type for inseparable indicator function, that is for such an indicator
function which cannot be written as $\chi(x,y) = \max(\Chi(x, y), \Chi(y,
x))$, where $\Chi(x, y) = \chi_1(x)\chi_2(y)$.

`χ` must be a function of two arguments which returns a value of `Bool`
type.
"""
struct InseparableIndicator <: AbstractIndicator
    χ :: Function
end

indicator_function(x :: InseparableIndicator) = x.χ
indicator_function(x :: SeparableIndicator)   = x.χ1, x.χ2

ext: jl


text:
"""
Clean up and extract MCMC_draws and MCMC_diagnostics from MetropolisHastings.
"""
function tidybayes(MCMC_draws, MCMC_diagnostics)
    param_df = DataFrame(MCMC_draws, :auto)
    K_terms = size(param_df, 1)
    param_df.term = "parameter_" .* string.(1:K_terms)
    param_df = DataFramesMeta.stack(param_df, Not(:term))
    transform!(param_df, :variable => ByRow(x -> parse(Int, split(x, "x")[2])) => :draw)
    select!(param_df, Not(:variable))

    draw_df = DataFrame(
        [
            :draw => 1:length(MCMC_draws),
            :accept => MCMC_diagnostics[:, 1]
        ])
        param_df
    all_draw_df = innerjoin(draw_df, param_df, on = :draw)
    return all_draw_df
end

max_stars_repo_name: EdJeeOnGitHub/Archibald

max_issues_repo_name: EdJeeOnGitHub/Archibald

max_forks_repo_name: EdJeeOnGitHub/Archibald


text:
# ---
# title: 171. Excel Sheet Column Number
# id: problem171
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/excel-sheet-column-number/description/>
# hidden: true
# ---
#
# Given a column title as appear in an Excel sheet, return its corresponding
# column number.
#
# For example:
#
#
#
#         A -> 1
#         B -> 2
#         C -> 3
#         ...
#         Z -> 26
#         AA -> 27
#         AB -> 28
#         ...
#
#
# **Example 1:**
#
#
#
#     Input: "A"
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: "AB"
#     Output: 28
#
#
# **Example 3:**
#
#
#
#     Input: "ZY"
#     Output: 701
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 7`
#   * `s` consists only of uppercase English letters.
#   * `s` is between "A" and "FXSHRXW".
#
#
## @lc code=start
using LeetCode

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


text:
"""
DuckDB vector
"""
mutable struct Vec
    handle::duckdb_vector

    function Vec(handle::duckdb_vector)
        result = new(handle)
        return result
    end
end

function get_array(vector::Vec, ::Type{T})::Vector{T} where {T}
    raw_ptr = duckdb_vector_get_data(vector.handle)
    ptr = Base.unsafe_convert(Ptr{T}, raw_ptr)
    return unsafe_wrap(Vector{T}, ptr, VECTOR_SIZE, own = false)
end

function get_validity(vector::Vec)::ValidityMask
    duckdb_vector_ensure_validity_writable(vector.handle)
    validity_ptr = duckdb_vector_get_validity(vector.handle)
    ptr = Base.unsafe_convert(Ptr{UInt64}, validity_ptr)
    validity_vector = unsafe_wrap(Vector{UInt64}, ptr, VECTOR_SIZE ÷ BITS_PER_VALUE, own = false)
    return ValidityMask(validity_vector)
end

function all_valid(vector::Vec)::Bool
    return duckdb_vector_get_validity(vector.handle) == C_NULL
end

function list_child(vector::Vec)::Vec
    return Vec(duckdb_list_vector_get_child(vector.handle))
end

function list_size(vector::Vec)::UInt64
    return duckdb_list_vector_get_size(vector.handle)
end

function struct_child(vector::Vec, index::UInt64)::Vec
    return Vec(duckdb_struct_vector_get_child(vector.handle, index))
end

function assign_string_element(vector::Vec, index, str::AbstractString)
    return duckdb_vector_assign_string_element_len(vector.handle, index, str, sizeof(str))
end

max_stars_repo_name: kimmolinna/duckdb

max_issues_repo_name: kimmolinna/duckdb

max_forks_repo_name: kimmolinna/duckdb


text:
"""
ChebyshevWeight{kind,T}()

is a quasi-vector representing the Chebyshev weight of the specified kind on -1..1.
That is, `ChebyshevWeight{1}()` represents `1/sqrt(1-x^2)`, and
`ChebyshevWeight{2}()` represents `sqrt(1-x^2)`.
"""
struct ChebyshevWeight{kind,T} <: AbstractJacobiWeight{T} end
ChebyshevWeight{kind}() where kind = ChebyshevWeight{kind,Float64}()
ChebyshevWeight() = ChebyshevWeight{1,Float64}()

getproperty(w::ChebyshevWeight{1,T}, ::Symbol) where T = -one(T)/2
getproperty(w::ChebyshevWeight{2,T}, ::Symbol) where T = one(T)/2

"""
Chebyshev{kind,T}()

is a quasi-matrix representing Chebyshev polynomials of the specified kind (1, 2, 3, or 4)
on -1..1.
"""
struct Chebyshev{kind,T} <: AbstractJacobi{T} end
Chebyshev{kind}() where kind = Chebyshev{kind,Float64}()

AbstractQuasiArray{T}(::Chebyshev{kind}) where {T,kind} = Chebyshev{kind,T}()
AbstractQuasiMatrix{T}(::Chebyshev{kind}) where {T,kind} = Chebyshev{kind,T}()

const WeightedChebyshev{kind,T} = WeightedBasis{T,<:ChebyshevWeight{kind},<:Chebyshev{kind}}

WeightedChebyshev{kind}() where kind = ChebyshevWeight{kind}() .* Chebyshev{kind}()
WeightedChebyshev{kind,T}() where {kind,T} = ChebyshevWeight{kind,T}(λ) .* Chebyshev{kind,T}(λ)

const ChebyshevTWeight = ChebyshevWeight{1}
const ChebyshevUWeight = ChebyshevWeight{2}
const ChebyshevT = Chebyshev{1}
const ChebyshevU = Chebyshev{2}
const WeightedChebyshevT = WeightedChebyshev{1}
const WeightedChebyshevU = WeightedChebyshev{2}

# conveniences...perhaps too convenient
Chebyshev() = Chebyshev{1}()
WeightedChebyshev() = WeightedChebyshevT()

chebyshevt() = ChebyshevT()
chebyshevt(d::AbstractInterval{T}) where T = ChebyshevT{float(T)}()[affine(d, ChebyshevInterval{T}()), :]
chebyshevt(d::Inclusion) = chebyshevt(d.domain)
chebyshevt(S::AbstractQuasiMatrix) = chebyshevt(axes(S,1))
chebyshevu() = ChebyshevU()
chebyshevu(d::AbstractInterval{T}) where T = ChebyshevU{float(T)}()[affine(d, ChebyshevInterval{T}()), :]
chebyshevu(d::Inclusion) = chebyshevu(d.domain)
chebyshevu(S::AbstractQuasiMatrix) = chebyshevu(axes(S,1))

"""
     chebyshevt(n, z)

computes the `n`-th Chebyshev polynomial of the first kind at `z`.
"""
chebyshevt(n::Integer, z::Number) = Base.unsafe_getindex(ChebyshevT{typeof(z)}(), z, n+1)
"""
     chebyshevt(n, z)

computes the `n`-th Chebyshev polynomial of the second kind at `z`.
"""
chebyshevu(n::Integer, z::Number) = Base.unsafe_getindex(ChebyshevU{typeof(z)}(), z, n+1)

chebysevtweight(d::AbstractInterval{T}) where T = ChebyshevTWeight{float(T)}[affine(d,ChebyshevInterval{T}())]
chebysevuweight(d::AbstractInterval{T}) where T = ChebyshevUWeight{float(T)}[affine(d,ChebyshevInterval{T}())]

==(a::Chebyshev{kind}, b::Chebyshev{kind}) where kind = true
==(a::Chebyshev, b::Chebyshev) = false
==(::Chebyshev, ::Jacobi) = false
==(::Jacobi, ::Chebyshev) = false
==(::Chebyshev, ::Legendre) = false
==(::Legendre, ::Chebyshev) = false

OrthogonalPolynomial(w::ChebyshevWeight{kind,T}) where {kind,T} = Chebyshev{kind,T}()
orthogonalityweight(P::Chebyshev{kind,T}) where {kind,T} = ChebyshevWeight{kind,T}()

function getindex(w::ChebyshevTWeight, x::Number)
    x ∈ axes(w,1) || throw(BoundsError())
    inv(sqrt(1-x^2))
end

function getindex(w::ChebyshevUWeight, x::Number)
    x ∈ axes(w,1) || throw(BoundsError())
    sqrt(1-x^2)
end

sum(::ChebyshevWeight{1,T}) where T = convert(T,π)
sum(::ChebyshevWeight{2,T}) where T = convert(T,π)/2

normalizationconstant(::ChebyshevT{T}) where T = Vcat(sqrt(inv(convert(T,π))), Fill(sqrt(2/convert(T,π)),∞))

Jacobi(C::ChebyshevT{T}) where T = Jacobi(-one(T)/2,-one(T)/2)
Jacobi(C::ChebyshevU{T}) where T = Jacobi(one(T)/2,one(T)/2)

#######
# transform
#######

factorize(L::SubQuasiArray{T,2,<:ChebyshevT,<:Tuple{<:Inclusion,<:OneTo}}) where T =
    TransformFactorization(grid(L), plan_chebyshevtransform(Array{T}(undef, size(L,2))))

# TODO: extend plan_chebyshevutransform
factorize(L::SubQuasiArray{T,2,<:ChebyshevU,<:Tuple{<:Inclusion,<:OneTo}}) where T<:FastTransforms.fftwNumber =
    TransformFactorization(grid(L), plan_chebyshevutransform(Array{T}(undef, size(L,2))))

########
# Jacobi Matrix
########

jacobimatrix(C::ChebyshevT{T}) where T =
    Tridiagonal(Vcat(one(T), Fill(one(T)/2,∞)), Zeros{T}(∞), Fill(one(T)/2,∞))

jacobimatrix(C::ChebyshevU{T}) where T =
    SymTridiagonal(Zeros{T}(∞), Fill(one(T)/2,∞))

# These return vectors A[k], B[k], C[k] are from DLMF.
recurrencecoefficients(C::ChebyshevT) = (Vcat(1, Fill(2,∞)), Zeros{Int}(∞), Ones{Int}(∞))
recurrencecoefficients(C::ChebyshevU) = (Fill(2,∞), Zeros{Int}(∞), Ones{Int}(∞))

# special clenshaw!
function copyto!(dest::AbstractVector{T}, v::SubArray{<:Any,1,<:Expansion{<:Any,<:ChebyshevT}, <:Tuple{AbstractVector{<:Number}}}) where T
    f = parent(v)
    (x,) = parentindices(v)
    P,c = arguments(f)
    clenshaw!(paddeddata(c), x, dest)
end

###
# Mass matrix
###

@simplify function *(Tc::QuasiAdjoint{<:Any,<:ChebyshevT}, wT::WeightedChebyshevT)
    V = promote_type(eltype(Tc), eltype(wT))
    Diagonal([convert(V,π); Fill(convert(V,π)/2,∞)])
end

@simplify function *(Tc::QuasiAdjoint{<:Any,<:ChebyshevU}, wT::WeightedChebyshevU)
    V = promote_type(eltype(Tc), eltype(wT))
    Diagonal(Fill(convert(V,π)/2,∞))
end

##########
# Derivatives
##########

# Ultraspherical(1)\(D*Chebyshev())
@simplify function *(D::Derivative{<:Any,<:ChebyshevInterval}, S::ChebyshevT)
    T = promote_type(eltype(D),eltype(S))
    A = _BandedMatrix((zero(T):∞)', ℵ₀, -1,1)
    ApplyQuasiMatrix(*, ChebyshevU{T}(), A)
end

#####
# Conversion
#####

\(::Chebyshev{kind,T}, ::Chebyshev{kind,V}) where {kind,T,V} = SquareEye{promote_type(T,V)}(ℵ₀)

function \(U::ChebyshevU, C::ChebyshevT)
    T = promote_type(eltype(U), eltype(C))
    _BandedMatrix(Vcat(-Ones{T}(1,∞)/2,
                        Zeros{T}(1,∞),
                        Hcat(Ones{T}(1,1),Ones{T}(1,∞)/2)), ℵ₀, 0,2)
end

function \(w_A::WeightedChebyshevT, w_B::WeightedChebyshevU)
    wA,A = w_A.args
    wB,B = w_B.args
    T = promote_type(eltype(w_A), eltype(w_B))
    _BandedMatrix(Vcat(Fill(one(T)/2, 1, ∞), Zeros{T}(1, ∞), Fill(-one(T)/2, 1, ∞)), ℵ₀, 2, 0)
end

\(w_A::WeightedChebyshevU, w_B::WeightedChebyshevT) = inv(w_B \ w_A)
\(T::ChebyshevT, U::ChebyshevU) = inv(U \ T)

####
# interrelationships
####

# (18.7.3)

function \(A::ChebyshevT, B::Jacobi)
    J = Jacobi(A)
    Diagonal(J[1,:]) * (J \ B)
end

function \(A::Jacobi, B::ChebyshevT)
    J = Jacobi(B)
    (A \ J) * Diagonal(inv.(J[1,:]))
end

function \(A::Chebyshev, B::Jacobi)
    J = Jacobi(A)
    Diagonal(A[1,:] .\ J[1,:]) * (J \ B)
end

function \(A::Jacobi, B::Chebyshev)
    J = Jacobi(B)
    (A \ J) * Diagonal(J[1,:] .\ B[1,:])
end

function \(A::Jacobi, B::ChebyshevU)
    T = promote_type(eltype(A), eltype(B))
    (A.a == A.b == one(T)/2) || throw(ArgumentError())
    Diagonal(B[1,:] ./ A[1,:])
end

# TODO: Toeplitz dot Hankel will be faster to generate
function \(A::ChebyshevT, B::Legendre)
    T = promote_type(eltype(A), eltype(B))
   UpperTriangular( BroadcastMatrix{T}((k,j) -> begin
            (iseven(k) == iseven(j) && j ≥ k) || return zero(T)
            k == 1 && return Λ(convert(T,j-1)/2)^2/π
            2/π * Λ(convert(T,j-k)/2) * Λ(convert(T,k+j-2)/2)
        end, 1:∞, (1:∞)'))
end

\(A::AbstractJacobi, B::Chebyshev) = ApplyArray(inv,B \ A)

function \(A::Jacobi, B::WeightedBasis{<:Any,<:JacobiWeight,<:Chebyshev})
    w, T = B.args
    J = Jacobi(T)
    wJ = w .* J
    (A \ wJ) * (J \ T)
end

function \(A::Chebyshev, B::WeightedBasis{<:Any,<:JacobiWeight,<:Jacobi})
    J = Jacobi(A)
    (A \ J) * (J \ B)
end

function \(A::Chebyshev, B::WeightedBasis{<:Any,<:JacobiWeight,<:Chebyshev})
    J = Jacobi(A)
    (A \ J) * (J \ B)
end

####
# sum
####

function _sum(A::WeightedBasis{T,<:ChebyshevUWeight,<:ChebyshevU}, dims) where T
    w, U = A.args
    @assert dims == 1
    Hcat(convert(T, π)/2, Zeros{T}(1,∞))
end

function _sum(A::WeightedBasis{T,<:ChebyshevWeight,<:Chebyshev}, dims) where T
    @assert dims == 1
    Hcat(convert(T, π), Zeros{T}(1,∞))
end

function cumsum(T::ChebyshevT{V}; dims::Integer) where V
    @assert dims == 1
    Σ = _BandedMatrix(Vcat(-one(V) ./ (-2:2:∞)', Zeros{V}(1,∞), Hcat(one(V), one(V) ./ (4:2:∞)')), ℵ₀, 0, 2)
    ApplyQuasiArray(*, T, Vcat((-1).^(0:∞)'* Σ, Σ))
end

cumsum(f::Expansion{<:Any,<:ChebyshevT}) = cumsum(f.args[1]; dims=1) * f.args[2]

####
# algebra
####

broadcastbasis(::typeof(+), ::ChebyshevT, U::ChebyshevU) = U
broadcastbasis(::typeof(+), U::ChebyshevU, ::ChebyshevT) = U


text:
module MultiPolynomials

import
Base: (*), /, ^, +, -, <, >, rem, divrem, diff, isless, lexless, string, show, write, convert, zero, one

export Field, Index, Term, Polynomial, PolyAsArray, @variables
export terms, polynomial, coeff, deg, index, divides, evaluate
export extract_univariate_polynomial, convert_base_field
export grlexless, grevlexless
export set_order_lex, set_order_grlex, set_order_grevlex
export grevlexelimless, set_order_grevlexelimless
export LT, LC, LM, multidegree
export LCM, S_polynomial, groebner
export groebner_fgb, call_fgb
export fgb_gbasis, fgb_gbasis_elim
export call_giac, giac_gbasis, giac_greduce

const Field = Number

Index{N} = NTuple{N,UInt16}

immutable Term{K<:Field,N}
    index::Index{N}
    coefficient::K
end

Polynomial{K<:Field,N} = Dict{Index{N},K}

PolyAsArray{K<:Field,N} = Array{Term{K,N},1}

function polynomial{K<:Field, N}(terms::PolyAsArray{K, N})
    Polynomial{K, N}([(t.index, t.coefficient) for t in terms])
end

function terms{K<:Field,N}(f::Polynomial{K, N})
    PolyAsArray{K,N}(Term{K,N}[Term{K,N}(key, value) for (key,value) in f])
end

convert{K<:Field,N}(::Type{Polynomial{K,N}}, t::Term{K,N}) =
    Polynomial{K,N}([(t.index, t.coefficient)])

zero{K,N}(::Type{Polynomial{K,N}}) = Polynomial{K,N}([])
one{K,N}(::Type{Polynomial{K,N}}) = Polynomial{K,N}([((zeros(K,N)...), 1)])

deg(i::Index) = Int(sum(i))
deg(t::Term) = deg(t.index)

function deg(f::Polynomial)
    d = 0
    for key in keys(f)
        d = max(d, deg(key))
    end
    d
end

function deg{K<:Field, N}(f::PolyAsArray{K, N})
    d = 0
    for t in f
        d = max(d, deg(t))
    end
    d
end

index(t::Term) = map(x->convert(Int,x),t.index)

function coeff{K<:Field,N}(t::Term{K,N}, f::Polynomial{K,N})
    get(f, t.index, zero(K))/t.coefficient
end

function divides{K<:Field,N}(t1::Term{K,N}, t2::Term{K,N})
    for i=1:N
        if t1.index[i]>t2.index[i]
            return false
        end
    end
    return true
end

function divides{K<:Field,N}(t::Term{K,N}, f::Polynomial{K,N})
    for (i,c) in f
        if !divides(t, Term{K,N}(i,c))
            return false
        end
    end
    return true
end

function divides{K<:Field,N}(t::Term{K,N}, tt::PolyAsArray{K,N})
    for s in tt
        if !divides(t, s)
            return false
        end
    end
    return true
end

+{K<:Field, N}(f::Polynomial{K, N}) = copy(f)
+{K<:Field, N}(t::Term{K, N}) = t
-{K<:Field, N}(t::Term{K, N}) = Term(t.index, -t.coefficient)

function -{K<:Field, N}(f::Polynomial{K, N})
    Polynomial{K,N}([(i,-c) for (i,c) in f])
end

function +{K<:Field, N}(f::Polynomial{K, N}, g::Polynomial{K, N})
    r = copy(f)
    for (key, val) in g
        if haskey(r, key)
            r[key] += val
        else
            r[key] = val
        end
        if r[key]==zero(K)
            delete!(r, key)
        end
    end
    r
end

function -{K<:Field, N}(f::Polynomial{K, N}, g::Polynomial{K, N})
    r = copy(f)
    for (key, val) in g
        if haskey(r, key)
            r[key] -= val
        else
            r[key] = -val
        end
        if r[key]==zero(K)
            delete!(r, key)
        end
    end
    r
end

function +{K<:Field, N}(f::Polynomial{K, N}, t::Term{K, N})
    r = copy(f)
    key, val = (t.index, t.coefficient)
    if haskey(r, key)
        r[key] += val
    else
        r[key] = val
    end
    if r[key]==zero(K)
        delete!(r, key)
    end
    r
end

function -{K<:Field, N}(f::Polynomial{K, N}, t::Term{K, N})
    r = copy(f)
    key, val = (t.index, t.coefficient)
    if haskey(r, key)
        r[key] -= val
    else
        r[key] = -val
    end
    if r[key]==zero(K)
        delete!(r, key)
    end
    r
end

function -{K<:Field, N}(t::Term{K, N}, f::Polynomial{K, N},)
    r = -f
    key, val = (t.index, t.coefficient)
    if haskey(r, key)
        r[key] += val
    else
        r[key] = val
    end
    if r[key]==zero(K)
        delete!(r, key)
    end
    r
end

+{K<:Field, N}(t::Term{K, N}, f::Polynomial{K, N}) = +(f,t)

+{K<:Field, N}(t::Term{K, N}, s::Term{K, N}) = +(Polynomial{K,N}(t.index=>t.coefficient), s)

+{K<:Field, N}(f::Polynomial{K, N}, a::Field) = +(f, Term{K,N}((zeros(K,N)...), convert(K, a)))

+{K<:Field, N}(a::Field, f::Polynomial{K, N}) = +(f, a)

+{K<:Field, N}(t::Term{K, N}, a::Field) = +(Polynomial{K,N}(t.index=>t.coefficient), convert(K, a))

+{K<:Field, N}(a::Field, t::Term{K, N}) = +(t, a)

-{K<:Field, N}(t::Term{K, N}, s::Term{K, N}) = -(Polynomial{K,N}(t.index=>t.coefficient), s)

-{K<:Field, N}(f::Polynomial{K, N}, a::Field) = -(f, Term{K,N}((zeros(K,N)...), convert(K, a)))

-{K<:Field, N}(a::Field, f::Polynomial{K, N}) = -(Term{K,N}((zeros(K,N)...), convert(K, a)), f)

-{K<:Field, N}(t::Term{K, N}, a::Field) = -(Polynomial{K,N}(t.index=>t.coefficient), convert(K,a))

-{K<:Field, N}(a::Field, t::Term{K, N}) = -(convert(K,a), Polynomial{K,N}(t.index=>t.coefficient))

function _mul{N}(i1::Index{N}, i2::Index{N})
    Index{N}(([i1[k]+i2[k] for k=1:N]...))
end

function *{K<:Field, N}(t::Term{K, N}, a::Field)
    Term{K, N}(t.index, t.coefficient*convert(K, a))
end

*{K<:Field, N}(a::Field, t::Term{K, N}) = *(t, a)

function *{K<:Field, N}(t::Term{K, N}, s::Term{K, N})
    Term{K, N}( _mul(t.index, s.index), t.coefficient*s.coefficient)
end

function *{K<:Field, N}(f::Polynomial{K, N}, a::Field)
    a1 = convert(K, a)
    if a1==zero(K)
        return Polynomial{K, N}([])
    end
    r = copy(f)
    for key in keys(r)
        r[key]*=a1
    end
    r
end

*{K<:Field, N}(a::Field, f::Polynomial{K, N}) = *(f, a)

function *{K<:Field, N}(f::Polynomial{K, N}, t::Term{K, N})
    Polynomial{K, N}([(_mul(i, t.index), c*t.coefficient) for (i, c) in f])
end

*{K<:Field, N}(t::Term{K, N}, f::Polynomial{K, N}) = *(f, t)

function *{K<:Field, N}(f::Polynomial{K, N}, g::Polynomial{K, N})
    tt = terms(g)
    if length(tt) == 0
        return copy(g)
    end
    r = f*tt[1]
    for k=2:length(tt)
        r += f*tt[k]
    end
    r
end

function _div{N}(i1::Index{N}, i2::Index{N})
    Index{N}(([ (i1[k]>=i2[k] ? i1[k]-i2[k] : -1) for k=1:N]...))
    #deliberately throws error if not divisible
end

function /{K<:Field, N}(t::Term{K, N}, s::Term{K, N})
    Term{K, N}( _div(t.index, s.index), t.coefficient/s.coefficient)
end

function /{K<:Field, N}(f::Polynomial{K, N}, t::Term{K, N})
    Polynomial{K, N}([(_div(i, t.index), c/t.coefficient) for (i, c) in f])
end

function /{K<:Field, N}(t::Term{K, N}, a::Field)
    Term{K, N}(t.index, t.coefficient/convert(K, a))
end

function /{K<:Field, N}(f::Polynomial{K, N}, a::Field)
    a1 = convert(K, a)
    r = copy(f)
    for key in keys(r)
        r[key]/=a1
    end
    r
end

function _pow{N}(i::Index{N}, e::Integer)
    Index{N}(([e*i[k] for k=1:N]...))
end

function ^{K<:Field, N}(t::Term{K, N}, e::Integer)
    Term{K, N}(_pow(t.index, e), t.coefficient^e)
end

function ^{K<:Field, N}(f::Polynomial{K, N}, e::Integer)
    if e<0
        error("Nonnegative exponent expected")
    elseif e==0
        return Term{K,N}((zeros(K,N)...), 1)
    elseif e==1
        return f
    else
        s = bits(e)
        s = reverse(s[first(search(s,"1")):end])
        r = 0
        q = 0
        isfirst = true
        for k=1:length(s)
            if k==1
                q = f
            else
                q = q*q
            end
            if s[k:k]=="1"
                if isfirst
                    r = q
                    isfirst = false
                else
                    r = r*q
                end
            end
        end
        return r
    end
end

_varnames = []

macro variables(K, x...)
    N = length(x)
    T = parse(string("Term{", K, ",", N, "}"))
    q = Expr(:block)
    if length(x) == 1 && isa(x[1],Expr)
        @assert x[1].head === :tuple "@variables expected a list of symbols"
        x = x[1].args
    end
    i = 1
    for s in x
        index = zeros(UInt16, N)
        index[i] = 1
        @assert isa(s,Symbol) "@variables expected a list of symbols"
        push!(q.args, Expr(:(=), s, Expr(:call, T, Expr(:tuple, index...), 1) ) )
        i = i + 1
    end
    push!(q.args, Expr(:tuple, x...))
    global _varnames = map(string, x)
    eval(Main, q)
end

function _str(i::Index; latex::Bool=false)
    s = ""
    f = false
    for k = 1:length(i)
        if i[k]>0
            n = ""
            if k<=length(_varnames)
                n = _varnames[k]
            else
                if latex
                    n = string("X_{", k, "}")
                else
                    n = string("X", k)
                end
            end
            if latex
                s = string(s, n, i[k]>1?string("^{", i[k], "}"):"")
            else
                s = string(s, f?"*":"", n, i[k]>1?string("^", i[k]):"")
            end
            f = true
        end
    end
    s
end

_str(x; latex::Bool=false) = string(x)

function _str(t::Term; latex::Bool=false)
    return string(sum(t.index)!=0 && abs(t.coefficient)==1?(t.coefficient<0?"-":""):string(_str(t.coefficient, latex=latex),(sum(t.index)==0||latex)?"":"*"), _str(t.index, latex=latex))
end

#function _str{K<:Field,N}(tt::PolyAsArray{K,N}; latex::Bool=false)
function _str{K<:Field,N}(tt::Array{Term{K,N},1}; latex::Bool=false)
#function _str(tt::PolyAsArray; latex::Bool=false)
    if length(tt)==0
        return "0"
    end
    s = ""
    f = false
    for t in tt
        if f && t.coefficient>0
            s = string(s, "+")
        end
        f = true
        s = string(s, _str(t, latex=latex))
    end
    s
end

function _str{K<:Field,N}(f::Polynomial{K,N}; latex::Bool=false)
#function _str(f::Polynomial; latex::Bool=false)
#    _str(terms(f), latex=latex)
    tt = terms(f)
    if length(tt)==0
        return "0"
    end
    s = ""
    f = false
    for t in tt
        if f && t.coefficient>0
            s = string(s, "+")
        end
        f = true
        s = string(s, _str(t, latex=latex))
    end
    s
end

function _str{T}(q::Rational{T}; latex::Bool=false)
    if latex
        if denominator(q)==1
            return string(numerator(q))
        end
        s = ""
        if numerator(q)<0
            s = "-"
        end
        return string(s, "\\frac{", abs(numerator(q)), "}{", denominator(q),"}")
    else
        return string(q)
    end
end

string(t::Term) = _str(t)
show(io::IO, t::Term) = print(io, _str(t))
#show(io::IO, ::MIME"application/x-latex", t::Term) = write(io, "\$", _str(t, latex=true), "\$")
#show(io::IO, ::MIME"text/latex",  t::Term) = write(io, "\$", _str(t, latex=true), "\$")

string{K<:Field,N}(tt::PolyAsArray{K,N}) = _str(tt)
show{K<:Field,N}(io::IO, tt::PolyAsArray{K,N}) = print(io, _str(tt))
#show{K<:Field,N}(io::IO, ::MIME"application/x-latex", tt::PolyAsArray{K,N}) = write(io, "\$", _str(tt, latex=true), "\$")
#show{K<:Field,N}(io::IO, ::MIME"text/latex",  tt::PolyAsArray{K,N}) = write(io, "\$", _str(tt, latex=true), "\$")

#string(f::Polynomial) = _str(f)
#show(io::IO, f::Polynomial) = print(io, _str(f))
string{K<:Field,N}(f::Polynomial{K,N}) = _str(f)
show{K<:Field,N}(io::IO, f::Polynomial{K,N}) = print(io, _str(f))
#show(io::IO, ::MIME"application/x-latex", f::Polynomial) = write(io, "\$", _str(f, latex=true), "\$")
#show(io::IO, ::MIME"text/latex",  f::Polynomial) = write(io, "\$", _str(f, latex=true), "\$")

#show(io::IO, ::MIME"application/x-latex", q::Rational) = write(io, "\$", _str(q, latex=true), "\$")
#show(io::IO, ::MIME"text/latex",  q::Rational) = write(io, "\$", _str(q, latex=true), "\$")

evaluate{K<:Field,N}(t::Term{K,N}, x) = t.coefficient * prod([x[k]^t.index[k] for k=1:N])
evaluate{K<:Field,N}(f::Polynomial{K,N}, x) = sum([evaluate(t, x) for t in terms(f)])

function diff{K<:Field,N}(t::Term{K,N}, k::Int)
    if t.index[k] > 0
        i = [t.index...]
        i[k] -= 1
        return Term{K,N}((i...), t.index[k]*t.coefficient)
    else
        return zero(K)
    end
end

function diff{K<:Field,N}(t::Term{K,N}, v::Term{K,N})
    @assert sum(v.index)==1
    diff(t, findfirst(v.index,1))
end

function diff{K<:Field,N}(f::Polynomial{K,N}, k::Int)
    r = Polynomial{K,N}()
    for (i,c) in f
        r = r + diff(Term{K,N}(i,c), k)
    end
    r
end

function diff{K<:Field,N}(f::Polynomial{K,N}, v::Term{K,N})
    @assert sum(v.index)==1
    diff(f, findfirst(v.index,1))
end

function convert_base_field{K<:Field, N}(K_new::DataType, p::Polynomial{K,N})
    Polynomial{K_new,N}([(i, convert(K_new, c))  for (i,c) in p])
end

function extract_univariate_polynomial{K<:Field, N}(p::Polynomial{K,N})
    t = terms(p)
    d = deg(t)
    up = zeros(K, d+1)
    v = 0
    for i = 1:length(t)
        if v==0 && deg(t[i])>0
            v = findfirst(t[i].index)
        end
        if v>0
            for j=1:N
                if j!=v && t[i].index[j]!=0
                    error("not a univariate polynomial")
                end
            end
            up[t[i].index[v]+1] = t[i].coefficient
        else
            up[1] = t[i].coefficient
        end
    end
    up
end

include("groebner_et_al.jl")
include("interfaces.jl")

end #MultiPolynomials


text:
include("filtering.jl");
gr(markerstrokecolor = :black,markersize = 8,)
##
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]
s = filter(r->r.Treatment in list &&
    r.Trial < 51 &&
    # r.Trial_Travel_to < 40 &&
    r.MouseID != "pc7",
    streaks)
b = filter(r->r.Treatment in list &&
    r.Trial < 51 &&
    # r.Trial_Travel_to < 40 &&
    r.MouseID != "pc7",
    bouts)
p = copy(pokes)
for df in [s,b, p]
    df[df.Treatment .== "PreVehicle",:Treatment] .= "Control"
    df[df.Treatment .== "Saline",:Treatment] = [o ? "Optogenetic" : "Control" for o in df[df.Treatment .== "Saline",:Stim]]
    df[df.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
    df[df.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in df[df.Treatment .== "SB242084_opt",:Stim]]
end
##
#look at survival function without censor data. That is looking only at last bout in a trial, not bout interrupted by a reward
f_bouts = filter(r-> r.Omissions_plus_one >= 2 &&
    r.Leave &&
    r.Treatment in ["Control", "Citalopram", "Optogenetic", "SB242084", "Altanserin"],
    b)
f_bouts.BoutOut
f_bouts.Bout_duration
res, plt = mediansurvival_analysis(f_bouts,:Bout_duration,:Treatment)
Drug_colors!(res)
plt = @df res scatter(string.(:Treatment), :Mean, yerror = :Sem, color = :color,
    xlabel = "Group", ylabel = "Median survival time", label = "")
savefig(joinpath(figs_loc,"ForagingMeeting/mediansurvival.pdf"))

f_res = function_analysis(f_bouts,:Bout_duration, survivalrate_algorythm; grouping = :Treatment, step =0.05)
Drug_colors!(f_res)
@df f_res plot(:Bout_duration,:Mean, ribbon = :Sem, group = :Treatment,
    color = :color, linecolor = :color, xlims = (0,15),
    xlabel = "Consecutive unrewarded poking time", ylabel = "Survival rate")
savefig(joinpath(figs_loc,"ForagingMeeting/SurvivalRate.pdf"))
##
b2 = filter(r -> r.Treatment in ["Control", "Citalopram", "Optogenetic", "SB242084", "Altanserin"], b)
KMdf = combine(groupby(b2,[:MouseID,:Treatment]), [:Bout_duration, :Leave] => ((d,l) -> fit(KaplanMeier,d,l)) => :KMfit)
transform!(KMdf, :KMfit => ByRow(SwitchingAnalysis.KM_median) => :KM_median)
dropmissing!(KMdf)
res2 = combine(groupby(KMdf, :Treatment), :KM_median .=> [mean, sem])
Drug_colors!(res2)
plt = @df res2 scatter(string.(:Treatment), :KM_median_mean, yerror = :KM_median_sem, color = :color,
    xlabel = "Group", ylabel = "Median Kaplan Meier survival time", label = "")
savefig(joinpath(figs_loc,"ForagingMeeting/KMmediansurvival.pdf"))

survdf = survival_analysis(b2,:Bout_duration, :Leave; grouping = :Treatment, step =0.05)
Drug_colors!(survdf)
@df survdf plot(:Bout_duration, :Mean, ribbon = :Sem, group = :Treatment,
    color = :color, linecolor = :color, xlims = (0,10),
    xlabel = "Consecutive unrewarded poking time", ylabel = "Survival rate")
savefig(joinpath(figs_loc,"ForagingMeeting/KMSurvivalRate.pdf"))
##
f_pokes = filter(r -> r.Treatment in ["Control", "Citalopram", "Optogenetic", "SB242084", "Altanserin"], p)
f_pokes.Treatment = categorical(f_pokes.Treatment)
levels!(f_pokes.Treatment,["Control", "Citalopram", "Optogenetic", "SB242084", "Altanserin"])
gd = groupby(f_pokes, [:MouseID, :Day, :Treatment, :Bout])
res3 = combine(gd, [:PokeIn, :PokeOut] => ((i,o) -> (In_bout = i .- first(i), Out_bout = o .- first(i))) => AsTable,
    :LastPoke => :Leave,
    :Protocol)
transform!(res3, :Out_bout => zscore)
verbFullBasic = @formula(Leave ~ 1 + Out_bout_zscore * Protocol +  (1+Out_bout_zscore+Protocol|MouseID));
FullBasic = fit(MixedModel,verbFullBasic, res3, Bernoulli(), fast = true)

verbSimpleBasic = @formula(Leave ~ 1 + Out_bout_zscore + Protocol +  (1+Out_bout_zscore+Protocol|MouseID));
SimpleBasic = fit(MixedModel,verbSimpleBasic, res3, Bernoulli(), fast = true)

SimpleFullBasicTest = MixedModels.likelihoodratiotest(LeaveSimpleBasic,LeaveBasic)

verbTreatment = @formula(Leave ~ 1 + Out_bout_zscore * Protocol * Treatment +  (1+Out_bout_zscore+Protocol|MouseID));
LeaveTreatment = fit(MixedModel,verbTreatment, res3, Bernoulli(), fast = true)

TreatmentTest = MixedModels.likelihoodratiotest(LeaveBasic,LeaveTreatment)


text:
using Test
using XLSX: readxlsx
using DataFrames
using Mimi
using MimiDICE2013
using CSVFiles

using MimiDICE2013: getparams

@testset "MimiDICE2013" begin

    # ------------------------------------------------------------------------------
    #   1. Run tests on the whole model
    # ------------------------------------------------------------------------------

    @testset "MimiDICE2013-model" begin

        m = MimiDICE2013.get_model();
        run(m)

        f = readxlsx(joinpath(dirname(@__FILE__), "..", "data", "DICE_2013_Excel.xlsm"))

        # Test Precision
        Precision = 1.0e-11

        # Time Periods
        T = 60

        # TATM Test (temperature increase)
        True_TATM = getparams(f, "B99:BI99", :all, "Base", T);
        @test maximum(abs, m[:climatedynamics, :TATM] .- True_TATM) ≈ 0. atol = Precision

        # MAT Test (carbon concentration atmosphere)
        True_MAT = getparams(f, "B87:BI87", :all, "Base", T);
        @test maximum(abs, m[:co2cycle, :MAT] .- True_MAT) ≈ 0. atol = Precision

        # DAMFRAC Test (damages fraction)
        True_DAMFRAC = getparams(f, "B105:BI105", :all, "Base", T);
        @test maximum(abs, m[:damages, :DAMFRAC] .- True_DAMFRAC) ≈ 0. atol = Precision

        # DAMAGES Test (damages $)
        True_DAMAGES = getparams(f, "B106:BI106", :all, "Base", T);
        @test maximum(abs, m[:damages, :DAMAGES] .- True_DAMAGES) ≈ 0. atol = Precision

        # E Test (emissions)
        True_E = getparams(f, "B112:BI112", :all, "Base", T);
        @test maximum(abs, m[:emissions, :E] .- True_E) ≈ 0. atol = Precision

        # YGROSS Test (gross output)
        True_YGROSS = getparams(f, "B104:BI104", :all, "Base", T);
        @test maximum(abs, m[:grosseconomy, :YGROSS] .- True_YGROSS) ≈ 0. atol = Precision

        # CPC Test (per capita consumption)
        True_CPC = getparams(f, "B126:BI126", :all, "Base", T);
        @test maximum(abs, m[:neteconomy, :CPC] .- True_CPC) ≈ 0. atol = Precision

        # FORC Test (radiative forcing)
        True_FORC = getparams(f, "B100:BI100", :all, "Base", T);
        @test maximum(abs, m[:radiativeforcing, :FORC] .- True_FORC) ≈ 0. atol = Precision

        True_UTILITY = getparams(f, "B129:B129", :single, "Base", T);
        @test maximum(abs, m[:welfare, :UTILITY] .- True_UTILITY) ≈ 0. atol = Precision

    end # MimiDICE2013-model testset

    # ------------------------------------------------------------------------------
    #   2. Run tests to make sure integration version (Mimi v0.5.0)
    #   values match Mimi 0.4.0 values
    # ------------------------------------------------------------------------------
    @testset "MimiDICE2013-integration" begin

        Precision = 1.0e-10
        nullvalue = -999.999

        m = MimiDICE2013.get_model();
        run(m)

        for c in map(nameof, Mimi.compdefs(m)), v in Mimi.variable_names(m, c)

            # load data for comparison
            filepath = joinpath(@__DIR__, "../data/validation_data_v040/$c-$v.csv")
            results = m[c, v]

            df = load(filepath) |> DataFrame
            if typeof(results) <: Number
                validation_results = df[1,1]

            else
                validation_results = Matrix(df)

                # remove NaNs and Missings
                results[ismissing.(results)] .= nullvalue
                results[isnan.(results)] .= nullvalue
                validation_results[ismissing.(validation_results)] .= nullvalue
                validation_results[isnan.(validation_results)] .= nullvalue

                # match dimensions

                if size(validation_results, 1) == 1
                    validation_results = validation_results'
                end

            end
            @test results ≈ validation_results atol = Precision

        end # for loop

    end # MimiDICE2013-integration testset

    # ------------------------------------------------------------------------------
    #   3. Standard API
    # ------------------------------------------------------------------------------
    @testset "Standard API" begin

        m = MimiDICE2013.get_model()
        run(m)

        # Test the errors
        @test_throws ErrorException MimiDICE2013.compute_scc()  # test that it errors if you don't specify a year
        @test_throws ErrorException MimiDICE2013.compute_scc(year=2021)  # test that it errors if the year isn't in the time index
        @test_throws ErrorException MimiDICE2013.compute_scc(last_year=2299)  # test that it errors if the last_year isn't in the time index
        @test_throws ErrorException MimiDICE2013.compute_scc(year=2105, last_year=2100)  # test that it errors if the year is after last_year

        # Test the SCC
        scc1 = MimiDICE2013.compute_scc(year=2020)
        @test scc1 isa Float64

        # Test that it's smaller with a shorter horizon
        scc2 = MimiDICE2013.compute_scc(year=2020, last_year=2200)
        @test scc2 < scc1

        # Test that it's smaller with a larger prtp
        scc3 = MimiDICE2013.compute_scc(year=2020, last_year=2200, prtp=0.02)
        @test scc3 < scc2

        # Test with a modified model
        m = MimiDICE2013.get_model()
        update_param!(m, :t2xco2, 5)
        scc4 = MimiDICE2013.compute_scc(m, year=2020)
        @test scc4 > scc1   # Test that a higher value of climate sensitivty makes the SCC bigger

        # Test compute_scc_mm
        result = MimiDICE2013.compute_scc_mm(year=2030)
        @test result.scc isa Float64
        @test result.mm isa Mimi.MarginalModel
        marginal_temp = result.mm[:climatedynamics, :TATM]
        @test all(marginal_temp[1:findfirst(isequal(2030), MimiDICE2013.model_years)] .== 0.)
        @test all(marginal_temp[findfirst(isequal(2035), MimiDICE2013.model_years):end] .!= 0.)

    end

    # ------------------------------------------------------------------------------
    #   4. Deterministic SCC values
    # ------------------------------------------------------------------------------
    @testset "SCC values" begin

        atol = 1e-4 # TODO what is a reasonable tolerance given we test on a few different machines etc.

        # Test several validation configurations against the pre-saved values from MimiDICE2013
        specs = Dict([
            :year => [2020, 2050],
            :eta => [0, 1.5],
            :prtp => [0.015, 0.03],
            :last_year => [2200, 2305],
        ])

        results = DataFrame(year = [], eta = [], prtp = [], last_year = [], SC = [])

        for year in specs[:year]
            for eta in specs[:eta]
                for prtp in specs[:prtp]
                    for last_year in specs[:last_year]
                        sc = MimiDICE2013.compute_scc(year=Int(year), eta=eta, prtp=prtp, last_year=Int(last_year))
                        push!(results, (year, eta, prtp, last_year, sc))
                    end
                end
            end
        end

        validation_results = load(joinpath(@__DIR__, "..", "data", "SC validation data", "deterministic_sc_values_v1-0-2.csv")) |> DataFrame
        # diffs = sort(results[!, :SC] - validation_results[!, :SC], rev = true)
        # println(diffs)
        @test all(isapprox.(results[!, :SC], validation_results[!, :SC], atol = atol))

    end # SCC values testset
end # MimiDICE2013 testset

nothing


text:
# helper functions and types to deal with the fact that we are discretizing the output variables.

struct MCEBHistogram{S<:StepRangeLen, H<:Histogram, T}
    grid::S
    hist::H
    infty_bound::T
    max_bias::T
    η_infl::T
end

function MCEBHistogram(grid, hist;
                       infty_bound = Inf,
                       max_bias = 0.0,
                       η_infl =  0.01)
    MCEBHistogram(grid, hist, infty_bound, max_bias, η_infl)
end

function MCEBHistogram(grid::StepRangeLen, Zs::AbstractArray=[])
    xs = [-Inf; collect(grid); +Inf]
    hist = fit(Histogram, Zs, xs; closed=:right)
    MCEBHistogram(grid, hist)
end

StatsBase.binindex(mh::MCEBHistogram, x) = StatsBase.binindex(mh.hist, x)

Base.first(mh::MCEBHistogram) = Base.first(mh.grid)
Base.last(mh::MCEBHistogram) = Base.last(mh.grid)
Base.step(mh::MCEBHistogram) = Base.step(mh.grid)
broadcastable(mh::MCEBHistogram) = Ref(mh)

getindex(mhist::MCEBHistogram, k) = getindex(mhist.hist.weights, k)
lastindex(mhist::MCEBHistogram) = lastindex(mhist.hist.weights)

pdf(mh::MCEBHistogram) = mh.hist.weights ./ sum(mh.hist.weights)
midpoints(mh::MCEBHistogram) = midpoints(mh.grid)

"""
	DiscretizedStandardNormalSamples(marginal_grid)

## Example usage:
```julia
DiscretizedStandardNormalSamples(-5:0.1:5)
```

## Description
This type represents discretization of Standard Normal Samples to the grid `marginal_grid`.
Let ``t_1 < \\dotsc, < t_K`` the elements of `marginal_grid`, then this
type implies that we map each sample ``Z \\sim \\mathcal{N}(0,1)`` to one of the intervals
``[t_j, t_{j+1}),j=0,\\dotsc,K``, with ``t_0 = -\\infty`` and ``t_{k+1} =  +\\infty``.
 Our inference will depend only on the discretized samples.

	DiscretizedStandardNormalSamples(Zs::AbstractVector{<:StandardNormalSample}, marginal_grid)

If we also provide the constructor with a vector of `StandardNormalSample`'s, then it also
stores the Histogram of the samples (according to the binning defined above) in the `mhist` slot.
This is the multinomial sufficient statistics of the discretized data.

The `DiscretizedStandardNormalSamples` type also contains three more slots:

* `f_min`,`f_max`: These can be either `nothing` or vectors of length equal to the number of
histogram bins. They contain bounds so that:
```math
f_{\\text{min}}[k] \\leq \\mathbb P_G[ Z \\in \\text{bin }k] \\leq f_{\\text{max}}[k]
```
* `var_proxy`: This is typically a vector of point estimates of ``\\mathbb P_G[ Z \\in \\text{bin }k]``.
"""
struct DiscretizedStandardNormalSamples{MH<:MCEBHistogram,
                                       IX<:Union{Nothing, Vector{Int}},
                                       Tvar,
                                       Tmin,
                                       Tmax}
    mhist::MH
    idx::IX
    var_proxy::Tvar
    f_min::Tmin
    f_max::Tmax
end

function DiscretizedStandardNormalSamples(Zs::AbstractVector{<:StandardNormalSample},
                                          grid::AbstractVector;
                                          kwargs...)
    mhist = MCEBHistogram(grid, response.(Zs))
    DiscretizedStandardNormalSamples(Zs, mhist; kwargs...)
end

# point is ... can think of above as Int{K} in MLJ-lang
function DiscretizedStandardNormalSamples(Zs::AbstractVector{<:StandardNormalSample},
                                          mhist::MCEBHistogram;
                                          keep_idx=true,
                                          var_proxy=1/sqrt(2*pi),
                                          f_min=nothing,
                                          f_max=nothing)
    Zs = response.(Zs) # rethink.
    if keep_idx
        idx = StatsBase.binindex.(mhist, Zs)
    else
        idx = nothing
    end
    DiscretizedStandardNormalSamples(mhist, idx, var_proxy, f_min, f_max)
end

# point is ... can think of above as Int{K} in MLJ-lang
function DiscretizedStandardNormalSamples(mhist::MCEBHistogram;
                                          var_proxy=1/sqrt(2*pi),
                                          f_min=nothing,
                                          f_max=nothing)
    DiscretizedStandardNormalSamples(mhist, nothing, var_proxy, f_min, f_max)
end

function DiscretizedStandardNormalSamples(marginal_grid::AbstractVector)
	DiscretizedStandardNormalSamples(MCEBHistogram(marginal_grid))
end

response(Z_discr::DiscretizedStandardNormalSamples) = Z_discr.mhist.hist # what I want to say.
nobs(Z_discr::DiscretizedStandardNormalSamples) = sum(Z_discr.mhist.hist.weights)
nbins(Z_discr::DiscretizedStandardNormalSamples) = length(Z_discr.mhist.grid) + 1
#eltype(Z_discr::DiscretizedStandardNormalSamples) = eltype(Z_discr.mhist.grid)
#length(Z_discr::DiscretizedStandardNormalSamples) = sum(Z_discr.mhist.hist.weights)

function support(Z_discr::DiscretizedStandardNormalSamples)
    mhist = Z_discr.mhist
    supp = zeros(eltype(mhist.grid), nbins(Z_discr))
    supp[2:(end-1)] = StatsBase.midpoints(mhist.grid)

    # replace by gap just adjacent
    h = step(Z_discr.mhist.grid)
    supp[1] = mhist.grid[1] - h/2
    supp[end] =  mhist.grid[end] + h./2
    supp
end

# TODO: a lot of overlap with code for set_neighborhood based on KDE
"""
	set_neighborhood(Zs_discr::DiscretizedStandardNormalSamples,
                     prior::Distribution;
				     C∞ = Inf)

This function is useful for demonstration and simulations.

It returnss a [`DiscretizedStandardNormalSamples`](@ref) instance identical to `Zs_discr` but with
`f_min`,`f_max`,`var_proxy` set as follows. `prior`, i.e. ``G`` is assumed to be the ground truth
generating the true effect sizes ``\\mu_i`` and so `var_proxy` corresponds to the marginal
density ``\\phi \\star G``. `f_min` and `f_max` correspond to ``\\phi \\star G \\pm C_{\\infty}``.

"""
function set_neighborhood(Zs_discr::DiscretizedStandardNormalSamples,
                          prior::Distribution;
						  C∞ = Inf,
						  varproxy_lower_bound = false)

    grid = Zs_discr.mhist.grid

	f_marginal = pdf(marginalize(prior, Zs_discr))

    f_max = zeros(length(f_marginal))
    f_min = zeros(length(f_marginal))
	var_proxy = zeros(length(f_marginal))
    # fill in intermediate points first
    for i in 2:length(grid)
        x_left = grid[i-1]
        x_right = grid[i]
        h = x_right - x_left  # step(grid)...
        f_max[i] = f_marginal[i] + C∞*h
        f_min[i] = max(f_marginal[i] - C∞*h, 0.0)
        var_proxy[i] = varproxy_lower_bound ? max(f_marginal[i], C∞*h) : f_marginal[i]
    end

    f_max[1] = Inf
    f_min[1] = 0.0
    var_proxy[1] = f_marginal[1]
    f_max[end] = Inf
    f_min[end] = 0.0
    var_proxy[end] = f_marginal[end]

    Zs_discr = @set Zs_discr.var_proxy = var_proxy
    Zs_discr = @set Zs_discr.f_min = f_min
    Zs_discr = @set Zs_discr.f_max = f_max
    Zs_discr
end

struct DiscretizedStandardNormalSample{DISC<:DiscretizedStandardNormalSamples}
    samples::DISC
    bin::Int
end

# (i) indexing -> give
function (Zs_disc::DiscretizedStandardNormalSamples)(i::Int)
    DiscretizedStandardNormalSample(Zs_disc, i)
end

# FIXME: this will have to go once implement iterator protocol for
# DiscretizedStandardNormalSample
broadcastable(Zs::DiscretizedStandardNormalSamples) = Ref(Zs)

# convert to discretenonparametric dbn
function DiscreteNonParametric(Zs_discr::DiscretizedStandardNormalSamples)
    supp = support(Zs_discr)
    probs = pdf(Zs_discr.mhist)
    DiscreteNonParametric(supp, probs)
end

function DiscreteNonParametric(Zs_discr::MCEBHistogram)
   DiscreteNonParametric(DiscretizedStandardNormalSamples(Zs_discr))
end

# AffineEstimator
struct DiscretizedAffineEstimator{MH<:MCEBHistogram}
    mhist::MH
    Q::Vector{Float64}
    Qo::Float64 #offset
end

DiscretizedAffineEstimator(Z::DiscretizedStandardNormalSamples, args...) = DiscretizedAffineEstimator(Z.mhist, args...)

DiscretizedAffineEstimator(mhist, Q) = DiscretizedAffineEstimator(mhist, Q, Base.zero(Float64))

function DiscretizedAffineEstimator(mhist::MCEBHistogram, F::Function)
    Q = zeros(Float64, length(mhist.grid) + 1)
    Q[2:(end-1)] = F.(midpoints(mhist))
    DiscretizedAffineEstimator(mhist, Q)
end

function (calib::DiscretizedAffineEstimator)(x)
    idxs = StatsBase.binindex(calib.mhist, x)
    calib.Q[idxs] + calib.Qo
end

function (Q::DiscretizedAffineEstimator)(Zs::DiscretizedStandardNormalSamples)
	# TODO: what if idx is nothing?
	Q.Q[Zs.idx] .+ Q.Qo
end

function (Q::DiscretizedAffineEstimator)(x::StandardNormalSample)
    Q(response(x))
end

function std_proxy(estimator::DiscretizedAffineEstimator, Zs::DiscretizedStandardNormalSamples, n)
	first_moment = sum(estimator.Q .* Zs.var_proxy)
	second_moment = sum(estimator.Q.^2 .* Zs.var_proxy)
	sqrt(second_moment - first_moment^2)/sqrt(n)
end

#  function (calib::DiscretizedAffineEstimator)(x)
#    idxs = StatsBase.binindex(calib.mhist, x)
#   calib.Q[idxs] + calib.Qo
#end

function _get_plot_x(mhist::MCEBHistogram)
    x = midpoints(mhist.grid)
end

function _get_plot_y(mhist::MCEBHistogram)
    y = pdf(mhist)[2:(end-1)]
end

function _get_plot_x(affine::DiscretizedAffineEstimator)
    x = midpoints(affine.mhist.grid)
end

function _get_plot_y(affine::DiscretizedAffineEstimator)
    y = affine.Q[2:(end-1)] .+ affine.Qo
end

@recipe function f(affine::DiscretizedAffineEstimator)
    seriestype  -->  :path
    _get_plot_x(affine), _get_plot_y(affine)
end

@recipe function f(affines::Vector{<:DiscretizedAffineEstimator})
    seriestype  -->  :path
    # TODO: Check this is the same for all of them.
    x = _get_plot_x(affines[1])
    y = hcat(_get_plot_y.(affines)...)
    x,y
end

@recipe function f(mhist::MCEBHistogram;
                        show_bands=true,
                        as_density=false)
    x = MinimaxCalibratedEBayes._get_plot_x(mhist)
    y = MinimaxCalibratedEBayes._get_plot_y(mhist)
    infty_bound = mhist.infty_bound
    h = step(mhist)

    if as_density
        x = x
        y = y ./ h #well... somewhat hacky.
        infty_bound = infty_bound ./ h
        ylab --> "Density"
    else
        ylab --> "Probability"
    end

    seriestype  -->  :path

    if show_bands
        lower_lims = y .- max.(y .- infty_bound, 0)
        upper_lims = infty_bound
        ribbons --> (lower_lims, upper_lims)
        fillalpha --> 0.2
    end

    xlab --> "x"
    x,y
end

#TODO IDEA

# AbstractSummarizedEBayesSamoles

# dispatch on the below....
# HistogramSamples {BaseSamples + MCEBHistogram + optional idx mapping everyone to their bin}
# through kwarg, retain idx
# NeighborhoodHistogramSamples {BaseSamples + MCEBHistogram}


text:
using Pkg
using Plots
using StatsBase
using Distributions
using StatsPlots
using LinearAlgebra
using Random
using Plots.PlotMeasures

Pkg.activate(".")
using CanEvolve

function plotter(sim_object; bins::Int64 = 100)
    VAF = sim_object[2][1]
    scaler = fit(Histogram, VAF, nbins = bins)
    scaler = normalize(scaler, mode = :pdf)
    plot(scaler, fillcolor = "#008EA0", linecolor = "#008EA0", linealpha = 0.5, legend = false, xlims = (0,1), ytickfont = font(8, "Helvetica"), xtickfont = font(8, "Helvetica"), xtickfontsize=8,ytickfontsize=8)
    if sum(sim_object[3][1]) != 0
        vline!(sim_object[3][1], color = "#FF6F00", alpha = 1, linewidth = 3)
    end
    xlabel!("Variant allele frequency (VAF)")
    ylabel!("Density")
    #density!(VAF, linecolor = RGBA(249/255,101/255,103/255,0), linewidth = 5, linealpha = 0.95)
end

anim = Animation()
function plotter_two(x, y)
    gr(display_type=:inline)
    pl = plot(plotVAF2(x, bins = 100), plotVAF2(y, bins = 100), size=(1000,400), bottom_margin = 10mm)
end
for j in 1:20
    p, n = autoSimulation_fixed(log(2), 0.0, u = 100, depth = 120, rho = 0.001, lower_cutoff = 0.2, upper_cutoff = 0.35)
    push!(p, j, plotter_two(p, n))
    frame(anim)
end

gif(anim, "img/autosimulation.gif", fps=1)


text:
# ---
# title: 447. Number of Boomerangs
# id: problem447
# author: Qling
# date: 2021-05-30
# difficulty: Medium
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/number-of-boomerangs/description/>
# hidden: true
# ---
#
# You are given `n` `points` in the plane that are all **distinct** , where
# `points[i] = [xi, yi]`. A **boomerang**  is a tuple of points `(i, j, k)` such
# that the distance between `i` and `j` equals the distance between `i` and `k`
# **(the order of the tuple matters)**.
#
# Return  _the number of boomerangs_.
#
#
#
# **Example 1:**
#
#
#
#     Input: points = [[0,0],[1,0],[2,0]]
#     Output: 2
#     Explanation: The two boomerangs are [[1,0],[0,0],[2,0]] and [[1,0],[2,0],[0,0]].
#
#
# **Example 2:**
#
#
#
#     Input: points = [[1,1],[2,2],[3,3]]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: points = [[1,1]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `n == points.length`
#   * `1 <= n <= 500`
#   * `points[i].length == 2`
#   * `-104 <= xi, yi <= 104`
#   * All the points are **unique**.
#
#
## @lc code=start
using LeetCode
using DataStructures

function number_of_boomerangs(points::Vector{Vector{Int}})::Int
    res = 0
    for p in points
        distance_frequency = DefaultDict(0)
        for q in points
            p == q && continue
            dis = (p[1] - q[1])^2 + (p[2] - q[2])^2
            distance_frequency[dis] += 1
        end

        for k in values(distance_frequency)
            k > 1 && (res += k * (k - 1))
        end
    end

    return res
end
## @lc code=end


text:
# Convert a prolongation matrix P to a tentative operator S.
# P is piecewise constant over some number of aggregates.
# S has a n x n for each aggregate where n is the size of the aggregate.

using PETScBinaryIO

P = readPETSc(ARGS[1])

rows = rowvals(P)
m, n = size(P)
is = Vector{Int}()
js = Vector{Int}()
ks = Vector{Float64}()
j = 0
for i = 1:n
    nzrng = nzrange(P,i)
    agg_size = length(nzrng)
    if agg_size > 10000
        println("Aggregate is too large ($agg_size)")
        exit(-1)
    end
    block = fill(-1/agg_size, (agg_size, agg_size)) # B(BᵀB)⁻¹Bᵀ
    block[diagind(block)] += 1 # I - B(BᵀB)⁻¹Bᵀ

    U, _, _ = svd(block)
    for x = 1:agg_size
        for y = 1:agg_size-1
            push!(is, rows[nzrng[x]])
            push!(js, j + y)
            push!(ks, U[x,y])
        end
    end
    j += agg_size-1
end
S = sparse(is,js,ks,m,j)

writePETSc(ARGS[2], S)

max_stars_repo_name: ligmg/ligmg

max_issues_repo_name: ligmg/ligmg

max_forks_repo_name: ligmg/ligmg


text:
include("header.jl")
using Statistics, LinearAlgebra

sizes = [((2,4,3),(4,1,3)),((2,4,5),(4,8,5)),((2,8,4,3),(8,2,4,3))]
@testset "bmm" begin
    for t in (Float32, Float64)
        for s in sizes
            a = t(0.1)*randn(t, s[1]...)
            b = t(0.1)*randn(t, s[2]...)
            bmmul(w)=bmm(w[1],w[2])
            @test gradcheck(bmmul, (a,b))
            if gpu() >= 0
                c = bmm(a, b)
                ka = KnetArray(a)
                kb = KnetArray(b)
                kc = bmm(ka, kb)
                @test isapprox(c, Array(kc))
                @test gradcheck(bmmul, (ka,kb))
            end
        end
    end
end #end of testset

# suppress the return
nothing


text:
"""
    Either <: Augmentor.ArrayOperation

Description
--------------

Chooses between the given `operations` at random when applied.
This is particularly useful if one for example wants to first
either rotate the image 90 degree clockwise or anticlockwise (but
never both), and then apply some other operation(s) afterwards.

When compiling a pipeline, `Either` will analyze the provided
`operations` in order to identify the preferred formalism to use
when applied. The chosen formalism is chosen such that it is
supported by all given `operations`. This way the output of
applying `Either` will be inferable and the whole pipeline will
remain type-stable (even though randomness is involved).

By default each specified image operation has the same
probability of occurrence. This default behaviour can be
overwritten by specifying the `chance` manually.

Usage
--------------

    Either(operations, [chances])

    Either(operations...; [chances])

    Either(pairs...)

    *(operations...)

    *(pairs...)

Arguments
--------------

- **`operations`** : `NTuple` or `Vararg` of `Augmentor.ArrayOperation`
    that denote the possible choices to sample from when applied.

- **`chances`** : Optional. Denotes the relative chances for an
    operation to be sampled. Has to contain the same number of
    elements as `operations`. Either an `NTuple` of numbers if
    specified as positional argument, or alternatively a
    `AbstractVector` of numbers if specified as a keyword
    argument. If omitted every operation will have equal
    probability of occurring.

- **`pairs`** : `Vararg` of `Pair{<:Real,<:Augmentor.ArrayOperation}`.
    A compact way to specify an operation and its chance of
    occurring together.

See also
--------------

[`NoOp`](@ref), [`augment`](@ref)

Examples
--------------

```julia
using Augmentor
img = testpattern()

# all three operations have equal chance of occuring
augment(img, Either(FlipX(), FlipY(), NoOp()))
augment(img, FlipX() * FlipY() * NoOp())

# NoOp is twice as likely as either FlipX or FlipY
augment(img, Either(1=>FlipX(), 1=>FlipY(), 2=>NoOp()))
augment(img, Either(FlipX(), FlipY(), NoOp(), chances=[1,1,2]))
augment(img, Either((FlipX(), FlipY(), NoOp()), (1,1,2)))
augment(img, (1=>FlipX()) * (1=>FlipY()) * (2=>NoOp()))
```
"""
struct Either{N,T<:Tuple} <: ArrayOperation
    operations::T
    chances::SVector{N,Float64}
    cum_chances::SVector{N,Float64}

    function Either(operations::NTuple{N,ArrayOperation}, chances::SVector{N}) where N
        all(c->c>=0, chances) || throw(ArgumentError("All provided \"chances\" must be positive"))
        length(operations) > 0 || throw(ArgumentError("Must provide at least one operation in the constructor of \"Either\""))
        sum_chances = sum(chances)
        sum_chances > 0 || throw(ArgumentError("The sum of all provided \"chances\" must be strictly positive"))
        norm_chances = map(x -> Float64(x/sum_chances), chances)
        cum_chances = SVector(cumsum(norm_chances))
        new{N,typeof(operations)}(operations, norm_chances, cum_chances)
    end
end

Either() = throw(ArgumentError("Must provide at least one operation in the constructor of \"Either\""))

function Either(operations::NTuple{N,ArrayOperation}, chances::NTuple{N,Real} = map(op -> 1/length(operations), operations)) where N
    Either(operations, SVector{N}(chances))
end

function Either(operations::Vararg{ArrayOperation,N}; chances = map(op -> 1/length(operations), operations)) where N
    Either(operations, SVector{N}(map(Float64, chances)))
end

function Either(operations::Pair...)
    Either(map(last, operations), map(first, operations))
end

function Either(op::ArrayOperation, p::Real = .5)
    0 <= p <= 1. || throw(ArgumentError("The propability \"p\" has to be in the interval [0, 1]"))
    p1 = Float64(p)
    p2 = 1 - p1
    Either((op, NoOp()), (p1, p2))
end

Base.:*(op1::Pair{<:Number,<:Operation}, ops::Pair...) =
    Either(op1, ops...)
Base.:*(op1::Operation, ops::Operation...) = Either((op1, ops...))

for FUN in (:supports_view,
            :supports_stepview,
            :supports_permute,
            :supports_affine,
            :supports_affineview)
    # A predicate must be true for all the contained operations,
    # in order for it to be true for the "Either" containing them.
    @eval @inline ($FUN)(::Type{Either{N,T}}) where {N,T} =
        all($FUN, T.types)
end

function randparam(op::Either, img)
    p = safe_rand()
    for (i, p_i) in enumerate(op.cum_chances)
        if p <= p_i
            return i
        end
    end
    error("unreachable code reached")
end

# Choose lazy strategy based on shared support of operations.
# Note: We prefer "affine" only if "img" already is some
#   "InvWarpedView", otherwise the preference is
#   view > stepview > permute > affine > affineview
@generated function applylazy(op::Either, img::AbstractArray, idx)
    if isinvwarpedview(img) && supports_affine(op)
        :(applyaffine(op, img, idx))
    elseif isinvwarpedview(img) && supports_affineview(op)
        :(applyaffineview(op, img, idx))
    elseif supports_view(op)
        :(applyview(op, img, idx))
    elseif supports_stepview(op)
        :(applystepview(op, img, idx))
    elseif supports_permute(op)
        :(applypermute(op, img, idx))
    elseif supports_affine(op)
        :(applyaffine(op, prepareaffine(img), idx))
    elseif supports_affineview(op)
        :(applyaffineview(op, prepareaffine(img), idx))
    else
        :(throw(MethodError(applylazy, (op, img, idx))))
    end
end

@inline isinvwarpedview(::Type{SubArray{T,N,P,I,L}}) where {T,N,P<:InvWarpedView,I,L} = true
@inline isinvwarpedview(::Type{<:InvWarpedView}) = true
@inline isinvwarpedview(::Type) = false

# Only provided for the sake of interface completeness (not used)
function toaffinemap(op::Either, img, idx)
    supports_affine(typeof(op)) || throw(MethodError(toaffinemap, (op, img, idx)))
    op_inner = op.operations[idx]
    toaffinemap_common(op_inner, img, randparam(op_inner, img))
end

# To preserve type stability of "applyeager", always promote
# arrays to offset arrays
@inline _offsetarray(A::Array) = OffsetArray(A, 0, 0)
@inline _offsetarray(A::OffsetArray) = A
function applyeager(op::Either, img::AbstractArray, idx)
    _offsetarray(applyeager(op.operations[idx], img))
end

# Sample a random operation and pass the function call along.
# Note: "applyaffine" needs to map to "applyaffine_common" for
#   type stability, because otherwise the concrete type of the
#   "AffineMap" may differ from one operation to the next
#   (e.g. "Rotate" uses a "RotMatrix" by default, while "Scale"
#   for obvious reasons does not)
for KIND in (:permute, :view, :stepview, :affine, :affineview)
    FUN = Symbol(:apply, KIND)
    SUP = Symbol(:supports_, KIND)
    APP = startswith(String(KIND),"affine") ? Symbol(FUN, :_common) : FUN
    @eval function ($FUN)(op::Either, img::AbstractArray, idx)
        ($SUP)(typeof(op)) || throw(MethodError($FUN, (op, img, idx)))
        ($APP)(op.operations[idx], img)
    end
end

function showconstruction(io::IO, op::Either)
    chances_float = map(c->round(c,  digits=3), op.chances)
    if all(x->x≈chances_float[1], chances_float)
        for (i, op_i) in enumerate(op.operations)
            showconstruction(io, op_i)
            i < length(op.operations) && print(io, " * ")
        end
    else
        for (i, (op_i, p_i)) in enumerate(zip(op.operations, chances_float))
            print(io, '(', p_i, "=>")
            showconstruction(io, op_i)
            print(io, ')')
            i < length(op.operations) && print(io, " * ")
        end
    end
end

function Base.show(io::IO, op::Either)
    if get(io, :compact, false)
        print(io, "Either:")
        for (op_i, p_i) in zip(op.operations, op.chances)
            print(io, " (", round(Int, p_i*100), "%) ")
            show(IOContext(io, :compact => true), op_i)
            print(io, '.')
        end
    else
        print(io, "Augmentor.", typeof(op).name, " (1 out of ", length(op.operations), " operation(s)):")
        percent_int   = map(c->round(Int, c*100), op.chances)
        percent_float = map(c->round(c*100, digits=1), op.chances)
        percent = if any(i != f for (i,f) in zip(percent_int,percent_float))
            percent_float
        else
            percent_int
        end
        k = maximum(length(string(c)) for c in percent)
        for (op_i, p_i) in zip(op.operations, percent)
            println(io)
            print(io, "  - ", lpad(string(p_i), k, " "), "% chance to: ")
            show(IOContext(io, :compact => true), op_i)
        end
    end
end


text:
# This file was generated, do not modify it. # hide
opts = (resampling=Holdout(fraction_train=0.9), measure=cross_entropy)
res = evaluate!(knn; opts...)
round(res.measurement[1], sigdigits=3)


text:
function build_MultiPPOManager(env::NatureEnv, num_starting_players::Int, update_freq::Int, clip::Float32)
    ns, na = size(state(env, 1)), length(action_space(env,1))

    cnn_output_shape = Int.(floor.([ns[1], ns[2], 32]))
    create_actor() = Chain(
        Conv((3,3), ns[3]=>64, pad=(1,1), relu),
        Conv((3, 3), 64=>32, pad=(1,1), relu),
        Conv((3, 3), 32=>32, pad=(1,1), relu),
        flatten,
        Dense(prod(cnn_output_shape), 64),
        Dense(64, na)
    )

    create_critic() = Chain(
        Conv((3,3), ns[3]=>64, pad=(1,1), relu),
        Conv((3, 3), 64=>32, pad=(1,1), relu),
        Conv((3, 3), 32=>32, pad=(1,1), relu),
        flatten,
        Dense(prod(cnn_output_shape), 64),
        Dense(64, 1)
    )

    create_trajectory() = PPOTrajectory(;
        capacity = update_freq,
        state = Array{Float32, 4} => (ns..., num_starting_players),
        action = Vector{Int} => (num_starting_players,),
        action_log_prob = Vector{Float32} => (num_starting_players,),
        reward = Vector{Float32} => (num_starting_players,),
        terminal = Vector{Bool} => (num_starting_players,),
    )

    # TODO CHANGE THIS TO BE DYNAMIC UPDATE FREQ AS NUM AGENTS CHANGES
    create_policy() = PPOPolicy(
        approximator = ActorCritic(
            actor = create_actor() |> gpu,
            critic = create_critic() |> gpu,
            optimizer = ADAM(1e-3),
        ),
        γ = 0.99f0,
        λ = 0.95f0,
        clip_range = clip,
        max_grad_norm = 0.5f0,
        n_epochs = 4,
        n_microbatches = 4,
        actor_loss_weight = 1.0f0,
        critic_loss_weight = 0.5f0,
        entropy_loss_weight = 0.001f0,
        update_freq = update_freq,
    )

    create_agent() = Agent(
        policy = create_policy(),
        trajectory = create_trajectory()
    )
    agent1 = create_agent()
    agent_map = Dict(player => agent1 for player in 1:num_starting_players)
    agent_inv_map = Dict(agent1 => [player for player in 1:num_starting_players])

    agents = MultiPPOManager(
        agent_map,
        agent_inv_map,
        PPOTrajectory, # trace's type
    )
end