text:
# Model
include_model("double_pendulum")
include_constraints("free_time")

model = free_time_model(model)

function fd(model::DoublePendulum, x⁺, x, u, w, h, t)
    midpoint_implicit(model, x⁺, x, u, w, u[end])
end

# Horizon
T = 50
Tm = convert(Int, floor(T / 2))

# Time step
tf = 5.0
h = tf / (T - 1)

# Bounds

# ul <= u <= uu
ul, uu = control_bounds(model, T,
        [-50.0 * ones(model.m - 1); 0.01],
        [50.0 * ones(model.m - 1); 2.0 * h])

# Initial and final states
θ1 = pi / 3.0
θM = pi / 6.0
x1 = [pi / 2.0 + θ1; -2.0 * θ1; 0.0; 0.0]
xM = [pi / 2.0 + θM; -2.0 * θM; 0.0; 0.0]
xl, xu = state_bounds(model, T,
        x1 = [x1[1:2]; Inf * ones(2)])
xl[Tm][1:2] = copy(xM[1:2])
xu[Tm][1:2] = copy(xM[1:2])

# Trajectory initialization
x0 = [linear_interpolation(x1, xM, Tm)..., linear_interpolation(xM, x1, Tm)...] # linear interpolation on state
u0 = [[0.001 * rand(model.m - 1); h] for t = 1:T-1] # random controls

# Objective
include_objective("nonlinear_stage")
function lz_stage(x, u, t)
    if t > 1
        return 100.0 * (kinematics_ee(model, x)[2])^2.0
    else
        return 0.0
    end
end
lz_terminal(x) = 0.0
obj_z = nonlinear_stage_objective(lz_stage, lz_terminal)

obj_track = quadratic_time_tracking_objective(
        [Diagonal([1.0 * ones(2); zeros(2)]) for t = 1:T],
        [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1],
        [x0[t] for t = 1:T],
        [zeros(model.m) for t = 1:T],
        1.0)
obj = MultiObjective([obj_z, obj_track])

# Constraints
include(joinpath(pwd(), "src/constraints/loop.jl"))
con_loop = loop_constraints(model, 1, T)
con_free_time = free_time_constraints(T)
con = multiple_constraints([con_free_time, con_loop])

# Problem
prob = trajectory_optimization_problem(model,
           obj,
           T,
           # h = h,
           xl = xl,
           xu = xu,
           ul = ul,
           uu = uu,
           con = con)

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
@time z , info = solve(prob, copy(z0))

# Visualize
using Plots
x, u = unpack(z, prob)
_u = [u[t][1:2] for t = 1:T-1]
tf, t, _h = get_time(u)

plot(hcat(x...)', width = 2.0)
plot(hcat(u...)', width = 2.0, linetype = :steppost)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, x, Δt = h)

# LQR
Q = [Diagonal(ones(model.n)) for t = 1:T]
R = [Diagonal(0.1 * ones(model.m)) for t = 1:T-1]
K, P = tvlqr(model, x, u, Q, R, 0.0)

plot(hcat([vec(K[t]) for t = 1:T-1]...)',
    linetype = :steppost,
    width = 2.0,
    labels = "",
    title = "LQR gains",
    xlabel = "time step t")

plot(hcat([vec(P[t]) for t = 1:T-1]...)',
    linetype = :steppost,
    width = 2.0,
    labels = "",
    title = "LQR cost-to-go",
    xlabel = "time step t")

N = 5
TN = N * T - (N - 1)
xN = deepcopy(x)
uN = deepcopy(u)
for i = 1:(N - 1)
    xN = [xN..., x[2:end]...]
    uN = [uN..., u...]
end
QN = [Diagonal(ones(model.n)) for t = 1:TN]
RN = [Diagonal(0.1 * ones(model.m)) for t = 1:TN-1]
KN, PN = tvlqr(model, xN, uN, QN, RN, 0.0)

plot(hcat([vec(KN[t]) for t = 1:(TN - 1)]...)',
    linetype = :steppost,
    width = 2.0,
    title = "LQR gains ($N cycles)",
    xlabel = "time step t",
    labels = "")

plot(hcat([vec(PN[t]) for t = 1:(TN - 1)]...)',
    linetype = :steppost,
    width = 2.0,
    title = "LQR cost-to-go ($N cycles)",
    xlabel = "time step t",
    labels = "")

ctg = PN[1:T]
@save joinpath(@__DIR__, "double_pendulum_limit_cycle.jld2") x _u _h ctg


text:
#ExprC Stuff
abstract type exprC end
abstract type Value end

#Environment Stuff
struct Binding
    name::String
    val::Value
end

struct Env
    bindings::Array
end

struct numC <: exprC
    val::Number
end

struct strC <: exprC
    str::String
end

struct idC <: exprC
    id::String
end

struct lamC <: exprC
    parms::Array
    body::exprC
end

struct ifC <: exprC
    val::exprC
    then::exprC
    elsee::exprC
end

struct appC <: exprC
    fun::exprC
    args::Array
end

function valid_id(id)
    if id == "if"
        false
    elseif id == "let"
        false
    elseif id == "in"
        false
    elseif id == "fn"
        false
    else
        true
    end
end

#checks for no duplicates in function variables
function check_duplicates(parms)
    for i in range(1, length(parms))
        for j in range(i + 1, length(parms))
            if parms[i] == parms[j]
                return true
            end
        end
    end
    return false
end

#checks to sure all parameters are strings
function check_all_strings(parms)
    if !isa(parms, Array)
        return false
    end
    for s in parms
        if isa(s, String)
            return false
        end
    end
    return true
end

function parse_array(list)
    tr = []
    for l in list
        append!(tr, parse(l))
    end
    return tr
end

#input is in format ['var', ['z','=',14], ['+', 'z', 'z']]
#parses concrete syntax into an exprC
function parse(input)
    if isa(input, Number)
        return numC(input)
    elseif length(input) == 1 && isa(input, Array) && isa(input[1], Number)
        return numC(input[1])
    elseif typeof(input) == String
        if valid_id(input)
            return idC(input)
        else
            error("Invalid TULI5 id ", input, " in parse")
        end
    elseif length(input) == 4 && input[1] == "if"
        return ifC(parse(input[2]), parse(input[3]), parse(input[4]))
    elseif length(input) == 3 && input[1] == "fn" && check_all_strings(input[2])
        if check_duplicates(input[2])
            error("TULI5 duplicate parms parse")
        else
            print("here")
            return lamC(parse_array(input[2]), parse[input[3]])
        end
    elseif length(input) >= 2
        for i in range(2, length(input))
            append!(args, [parse(input[i])])
        end
        return appC(parse(input[1]), args)
    else
        error("Invlid TULI5 expression reading ", input, " in parse")
    end
end

parse(["fn", ["x"], ["+", "x", 3]])

ext: jl

max_stars_repo_name: nhausmna/julia-hello-world

max_issues_repo_name: nhausmna/julia-hello-world

max_forks_repo_name: nhausmna/julia-hello-world


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

module BB

frac(x) = x - floor(x)

include("structs.jl")

include("nodepool.jl")
include("optimize.jl")
include("log.jl")
include("lp.jl")
include("varbranch/hybrid.jl")
include("varbranch/infeasibility.jl")
include("varbranch/pseudocost.jl")
include("varbranch/random.jl")
include("varbranch/reliability.jl")
include("varbranch/strong.jl")

end # module


text:
"""
Holds collection of functions - mainly macros that when invoked from inside of the kernel will lead to iteration
Generally by convention we will use UInt32 as this is a format of threadid.x and other basic CUDA variables
arr - 3 dimensional data we analyze
maxX maxY maxZ- we will need also to supply data for variables called maxX maxY maxZ (if we have 3 dimensional loop) - those will create boundary checks  - for 2 dimensional case we obviously need only max X and max Y
minX, minY and minZ - this are assumed to be 1 if not we can supply them in the macro
loopX, loopY,loopZ - how many iterations should be completed while iterating over given dimension

we need also to supply functions of iterating the 3 dimensional data but with chosen dimension fixed -
  so for example if we want to analyze top padding from shared memory we will set the dimension 3 at 1
"""
module IterationUtils
using CUDA,Logging
export generalizedItermultiDim
export generalizedItermultiDimZdeepest,@iterateLinearlyCheckAll,@iterateLinearlyMultipleBlocks,@iterateLinearlyWithStart, @iterateLinearly,@exOnWarp,@exOnWarpIfBool, @iter3d, @iter3dAdditionalxyzActsAndZcheck, @iter3dAdditionalxyzActs, @iter3dAdditionalzActs,@iter3dWithVal

"""
arrDims- dimensions of main arrya
loopIterNameX ,loopIterNameY, loopIterNameZ  - variable name that will be used in a left hand side of loop
xname,yname,zname - symbols representing the x,y and z that are calculated for currnt position
xDimName,yDimName,zDimName - symbols representing the left hand sides of the loop
loopDims- information how many times we need to loop over given dimension
zOffset,zAdd - calculate offset and thread/block dependent add
Offset,xOffset,xAdd, yAdd - offsets for x and y  and what to add to them
xCheck,yCheck,zCheck - checks performed just after for - and determining wheather to continue
additionalActionAfterZ,additionalActionAfterY,additionalActionAfterX  - gives possibility of invoking more actions
    - invoked after the checks (checks are avoided to prevent warp stall if warp sync will be invoked)
is3d - if true we use 3 dimensional loop if not we iterate only over x and y
isFullBoundaryCheckX, isFullBoundaryCheckY, isFullBoundaryCheckZ - indicates wheather we want to check boundaries on all iterations if false it will be done only on last iteration if not stated explicitely to avoid all boundary checks
nobundaryCheckX, nobundaryCheckY, nobundaryCheckZ - true if we want to avoid completely boundary checks

ex - main expression around which we build loop
"""
macro iter3d(arrDims,loopXdim,loopYdim,loopZdim   ,ex   )
  mainExp = generalizedItermultiDim(; arrDims,loopXdim ,loopYdim,loopZdim, ex = ex)
  return esc(:( $mainExp))
  end#iter3d

"""
modification where value will be loaded from the suppplied array in x,y,z position
"""
  macro iter3dWithVal(arrMain,arrDims,loopXdim,loopYdim,loopZdim   ,ex   )
    mainExp = generalizedItermultiDim(; arrMain= arrMain,isVal=true, arrDims=arrDims,loopXdim ,loopYdim,loopZdim, ex = ex)
    return esc(:( $mainExp))
    end#iter3d

"""
modification of iter3d loop  where wa allow additional action after z check
"""
macro iter3dAdditionalxyzActsAndZcheck(arrDims,loopXdim,loopYdim,loopZdim,zCheck
  ,ex,additionalActionAfterX,additionalActionAfterY,additionalActionAfterZ)
  mainExp = generalizedItermultiDim(; arrDims=arrDims,loopXdim=loopXdim ,loopYdim=loopYdim,loopZdim=loopZdim,zCheck=zCheck, ex=ex
  ,additionalActionAfterX=additionalActionAfterX ,additionalActionAfterY=additionalActionAfterY
  ,additionalActionAfterZ=additionalActionAfterZ )

  return esc(:( $mainExp))
end#iter3dAdditionalxyzActs

"""
modification of iter3d loop  where wa allow additional actions to be performed after each loop check
"""
macro iter3dAdditionalzActs(arrDims, loopXdim,loopYdim,loopZdim,ex,additionalActionAfterZ)

  mainExp = generalizedItermultiDim(; arrDims=arrDims,loopXdim=loopXdim ,loopYdim=loopYdim,loopZdim=loopZdim
  ,additionalActionAfterZ=additionalActionAfterZ, ex = ex )

  return esc(:( $mainExp))
end#iter3dAdditionalxyzActs

"""
generalized version of iter3d we will specilize it in the macro on the basis of multiple dispatch
  arrDims- dimensions of main arrya
  loopIterNameX ,loopIterNameY, loopIterNameZ  - variable name that will be used in a left hand side of loop
  xname,yname,zname - symbols representing the x,y and z that are calculated for currnt position
  xDimName,yDimName,zDimName - symbols representing the left hand sides of the loop
  loopXdim, loopYdim,loopZdim - information how many times we need to loop over given dimension
  zOffset,zAdd - calculate offset and thread/block dependent add
  Offset,xOffset,xAdd, yAdd - offsets for x and y  and what to add to them
  xCheck,yCheck,zCheck - checks performed just after for - and determining wheather to continue
  additionalActionAfterZ,additionalActionAfterY,additionalActionAfterX  - gives possibility of invoking more actions
      - invoked after the checks (checks are avoided to prevent warp stall if warp sync will be invoked)
  is3d - if true we use 3 dimensional loop if not we iterate only over x and y
  isFullBoundaryCheckX, isFullBoundaryCheckY, isFullBoundaryCheckZ - indicates wheather we want to check boundaries on all iterations if false it will be done only on last iteration if not stated explicitely to avoid all boundary checks
  nobundaryCheckX, nobundaryCheckY, nobundaryCheckZ - true if we want to avoid completely boundary checks
  arrMain - array from which we will take value in case isVal is true
  isVal  - if true will indicate that we want  to have the value of the supplied array arrMain in x,y,z point

  ex - main expression around which we build loop
"""
function generalizedItermultiDim(; #we keep all as keyword arguments
  arrDims = (UInt32(1),UInt32(1),UInt32(1) )
   ,xname = :x
   ,yname = :y
   ,zname = :z
   ,loopIterNameX = :xdim
   ,loopIterNameY = :ydim
   ,loopIterNameZ = :zdim
   ,loopXdim = 1
   ,loopYdim= 1
   ,loopZdim= 1
  ,zOffset= :(zdim*gridDim().x)
   ,zAdd =:(blockIdxX())
  ,yOffset = :(ydim* blockDimY())
  ,yAdd= :(threadIdxY())
  ,xOffset= :(xdim * blockDimX())
   ,xAdd= :(threadIdxX())
  ,xCheck = :($xname <= $arrDims[1])
  ,yCheck = :($yname <= $arrDims[2])
  ,zCheck = :($zname <= $arrDims[3])
  ,additionalActionAfterZ= :()
   ,additionalActionAfterY= :()
   ,additionalActionAfterX = :()
   ,additionalActionBeforeZ= :()
   ,additionalActionBeforeY= :()
   ,additionalActionBeforeX = :()

   , is3d = true
   ,ex
   ,isFullBoundaryCheckX =false
   , isFullBoundaryCheckY=false
   , isFullBoundaryCheckZ=true
   ,nobundaryCheckX=false
   , nobundaryCheckY=false
   , nobundaryCheckZ =false
   ,isVal=false
   ,arrMain=:()
   ,loopStartX = 0
   ,loopStartY= 0
   ,loopStartZ= 0

   )
#we will define expressions from deepest to most superficial

# ,zOffset= :($loopIterNameZ*gridDim().x)
# ,zAdd =:(blockIdxX())
# ,yOffset = :($loopIterNameY* blockDimY())
# ,yAdd= :(threadIdxY())
# ,xOffset= :($loopIterNameX * blockDimX())

  # xState = :(x= $xOffset +$xAdd)
  # yState = :(y= $yOffset +$yAdd)
  # zState = :(z= $zOffset +$zAdd)
  # xState = :($xname= $xOffset +$xAdd)
  # yState = :($yname= $yOffset +$yAdd)
  # zState = :($zname= $zOffset +$zAdd)

  valExp= :()
  if(isVal)
    valExp=:(value = $arrMain[$xname,$yname,$zname ])
  end

  exp1= :()
  if(isFullBoundaryCheckX)
    exp1 = quote
      @unroll for $loopIterNameX in $loopStartX:$loopXdim
          $xname= $xOffset +$xAdd
          $additionalActionBeforeX
          if( $xCheck)
            $valExp
            $ex
          end#if x
        $additionalActionAfterX
        end#for x dim
      end#quote
  elseif(nobundaryCheckX)
    exp1 = quote
      @unroll for $loopIterNameX in $loopStartX:$loopXdim
          $xname= $xOffset +$xAdd
          $additionalActionBeforeX
          $valExp
          $ex
        $additionalActionAfterX
        end#for x dim
      end#quote

  else
    exp1 = quote
      @unroll for $loopIterNameX in $loopStartX:$loopXdim-1
          $xname= $xOffset +$xAdd
          $additionalActionBeforeX
          $valExp
          $ex
          $additionalActionAfterX
      end#for x dim

      $loopIterNameX=$loopXdim
      $xname= $xOffset +$xAdd
      $additionalActionBeforeX
        if( $xCheck)
          $valExp
          $ex
        end#if x
      $additionalActionAfterX
      end#quote
  end

  exp2= :()
if(isFullBoundaryCheckY)
    exp2= quote
      @unroll for $loopIterNameY in  $loopStartY:$loopYdim
        $yname = $yOffset +$yAdd
        $additionalActionBeforeY
          if($yCheck)
            $exp1
          end#if y
          $additionalActionAfterY
          end#for  yLoops
        end#quote
elseif(nobundaryCheckY)
  exp2= quote
    @unroll for $loopIterNameY in  $loopStartY:$loopYdim
      $yname = $yOffset +$yAdd
      $additionalActionBeforeY
          $exp1
        $additionalActionAfterY
        end#for  yLoops
      end#quote
else
  exp2= quote
    @unroll for $loopIterNameY in  $loopStartY:$loopYdim-1
      $yname = $yOffset +$yAdd
      $additionalActionBeforeY
        $exp1
        $additionalActionAfterY
    end#for  yLoops
        $loopIterNameY=$loopYdim
        $yname = $yOffset +$yAdd
        $additionalActionBeforeY

        if($yCheck)
          $exp1
        end#if y
        $additionalActionAfterY
      end#quote
end

############################### z
if(isFullBoundaryCheckZ)
        exp3= quote
          @unroll for $loopIterNameZ in $loopStartZ:$loopZdim
            $zname = $zOffset + $zAdd#multiply by blocks to avoid situation that a single block of threads would have no work to do
            $additionalActionBeforeZ
            if($zCheck)
              $exp2
            end#if z
            $additionalActionAfterZ
        end#for z dim
        end

elseif(nobundaryCheckZ)
        exp3= quote
          @unroll for $loopIterNameZ in $loopStartZ:$loopZdim
            $zname = $zOffset + $zAdd#multiply by blocks to avoid situation that a single block of threads would have no work to do
            $additionalActionBeforeZ
            $exp2
            $additionalActionAfterZ
        end#for z dim
        end
 else
        exp3= quote
          @unroll for $loopIterNameZ in $loopStartZ:$loopZdim-1
            $zname = $zOffset + $zAdd#multiply by blocks to avoid situation that a single block of threads would have no work to do
            $additionalActionBeforeZ
            $exp2
            $additionalActionAfterZ
        end#for z dim
        $loopIterNameZ=$loopZdim
        $additionalActionBeforeZ
        $zname = $zOffset + $zAdd
        if($zCheck)
          $exp2
        end#if z

        end
end

if(is3d)
  return exp3
end
#if 2d
return exp2

end#generalizedIter3d

"""
alternative with z dimension as the most inner iteration
"""

function generalizedItermultiDimZdeepest(; #we keep all as keyword arguments
  arrDims = (UInt32(1),UInt32(1),UInt32(1) )
   ,xname = :x
   ,yname = :y
   ,zname = :z
   ,loopIterNameX = :xdim
   ,loopIterNameY = :ydim
   ,loopIterNameZ = :zdim
   ,loopXdim = 1
   ,loopYdim= 1
   ,loopZdim= 1
  ,zOffset= :(zdim*gridDim().x)
   ,zAdd =:(blockIdxX())
  ,yOffset = :(ydim* blockDimY())
  ,yAdd= :(threadIdxY())
  ,xOffset= :(xdim * blockDimX())
   ,xAdd= :(threadIdxX())
  ,xCheck = :($xname <= $arrDims[1])
  ,yCheck = :($yname <= $arrDims[2])
  ,zCheck = :($zname <= $arrDims[3])
  ,additionalActionAfterZ= :()
   ,additionalActionAfterY= :()
   ,additionalActionAfterX = :()
   ,additionalActionBeforeZ= :()
   ,additionalActionBeforeY= :()
   ,additionalActionBeforeX = :()

   , is3d = true
   ,ex
   ,isFullBoundaryCheckX =false
   , isFullBoundaryCheckY=false
   , isFullBoundaryCheckZ=true
   ,nobundaryCheckX=false
   , nobundaryCheckY=false
   , nobundaryCheckZ =false
   ,isVal=false
   ,arrMain=:()
   ,loopStartX = 0
   ,loopStartY= 0
   ,loopStartZ= 0

   )
#we will define expressions from deepest to most superficial

# ,zOffset= :($loopIterNameZ*gridDim().x)
# ,zAdd =:(blockIdxX())
# ,yOffset = :($loopIterNameY* blockDimY())
# ,yAdd= :(threadIdxY())
# ,xOffset= :($loopIterNameX * blockDimX())

  # xState = :(x= $xOffset +$xAdd)
  # yState = :(y= $yOffset +$yAdd)
  # zState = :(z= $zOffset +$zAdd)
  # xState = :($xname= $xOffset +$xAdd)
  # yState = :($yname= $yOffset +$yAdd)
  # zState = :($zname= $zOffset +$zAdd)

  valExp= :()
  if(isVal)
    valExp=:(value = $arrMain[$xname,$yname,$zname ])
  end

############################### z
exp1= :()

if(isFullBoundaryCheckZ)
  exp1= quote
    @unroll for $loopIterNameZ in $loopStartZ:$loopZdim
      $zname = $zOffset + $zAdd#multiply by blocks to avoid situation that a single block of threads would have no work to do
      $additionalActionBeforeZ
      if($zCheck)
        $valExp
        $ex
      end#if z
      $additionalActionAfterZ
  end#for z dim
  end

elseif(nobundaryCheckZ)
  exp1= quote
    @unroll for $loopIterNameZ in $loopStartZ:$loopZdim
      $zname = $zOffset + $zAdd#multiply by blocks to avoid situation that a single block of threads would have no work to do
      $additionalActionBeforeZ
      $valExp
      $ex
      $additionalActionAfterZ
  end#for z dim
  end
else
  exp1= quote
    @unroll for $loopIterNameZ in $loopStartZ:$loopZdim-1
      $zname = $zOffset + $zAdd#multiply by blocks to avoid situation that a single block of threads would have no work to do
      $additionalActionBeforeZ
      $valExp
      $ex
      $additionalActionAfterZ
  end#for z dim
  $loopIterNameZ=$loopZdim
  $additionalActionBeforeZ
  $zname = $zOffset + $zAdd
  if($zCheck)
    $valExp
    $ex
  end#if z

  end
end

  ##############  x
  exp2= :()

  if(isFullBoundaryCheckX)
    exp2 = quote
      @unroll for $loopIterNameX in $loopStartX:$loopXdim
          $xname= $xOffset +$xAdd
          $additionalActionBeforeX
          if( $xCheck)
            $exp1
          end#if x
        $additionalActionAfterX
        end#for x dim
      end#quote
  elseif(nobundaryCheckX)
    exp2 = quote
      @unroll for $loopIterNameX in $loopStartX:$loopXdim
          $xname= $xOffset +$xAdd
          $additionalActionBeforeX
          $exp1
        $additionalActionAfterX
        end#for x dim
      end#quote

  else
    exp2 = quote
      @unroll for $loopIterNameX in $loopStartX:$loopXdim-1
          $xname= $xOffset +$xAdd
          $additionalActionBeforeX
          $exp1
          $additionalActionAfterX
      end#for x dim

      $loopIterNameX=$loopXdim
      $xname= $xOffset +$xAdd
      $additionalActionBeforeX
        if( $xCheck)
          $exp1
        end#if x
      $additionalActionAfterX
      end#quote
  end

#######y
  exp3= :()
if(isFullBoundaryCheckY)
    exp3= quote
      @unroll for $loopIterNameY in  $loopStartY:$loopYdim
        $yname = $yOffset +$yAdd
        $additionalActionBeforeY
          if($yCheck)
            $exp2
          end#if y
          $additionalActionAfterY
          end#for  yLoops
        end#quote
elseif(nobundaryCheckY)
  exp3= quote
    @unroll for $loopIterNameY in  $loopStartY:$loopYdim
      $yname = $yOffset +$yAdd
      $additionalActionBeforeY
          $exp2
        $additionalActionAfterY
        end#for  yLoops
      end#quote
else
  exp3= quote
    @unroll for $loopIterNameY in  $loopStartY:$loopYdim-1
      $yname = $yOffset +$yAdd
      $additionalActionBeforeY
        $exp2
        $additionalActionAfterY
    end#for  yLoops
        $loopIterNameY=$loopYdim
        $yname = $yOffset +$yAdd
        $additionalActionBeforeY

        if($yCheck)
          $exp2
        end#if y
        $additionalActionAfterY
      end#quote
end

if(is3d)
  return exp3
end
#if 2d
return exp2

end#generalizedIter3d

"""
macro will know about the number of available warps as this will equall the  y dimension of thread block
    now we will supply on what warp we want to execute the function  if the number will be smaller than number of warps
    it will be executed on chosen warp otherwise macro will perform modulus operation to establish index that is indicating some
    warp that exists
"""
macro exOnWarp(numb,ex)
    return  esc(quote
    if( ( (threadIdxY()==$numb) ||  (rem($numb,blockDimY()))==threadIdxY()) ||  (rem($numb,blockDimY())+blockDimY() )==threadIdxY())
      $ex
    end
        # if($numb<=blockDimY())
        #   if (threadIdxY()== $numb)
        #      $ex
        #   end
        # else
        #   if (threadIdxY()==(mod($numb,blockDimY())+1) )
        #     $ex
        #   end
        # end
    end)
 end

"""
macro will know about the number of available warps as this will equall the  y dimension of thread block
    now we will supply on what warp we want to execute the function  if the number will be smaller than number of warps
    it will be executed on chosen warp otherwise macro will perform modulus operation to establish index that is indicating some
    warp that exists
    we execute only if isMaskFull is false what indicates that there is a metadata block that is associated with this idX
        tobeEx - indicates is it to be executed
"""
macro exOnWarpIfBool(tobeEx,numb, ex)
    return  esc(quote
        if($tobeEx)
            if($numb<=blockDimY())
                @ifY $numb $ex
            else
                @ifY (mod($numb,blockDimY())+1) $ex
            end
        end
    end)
end

"""
iterate over array that is treated as one dimensional with given length lengthh as argument
amount of iterations needed is also passed as an argument - iterLoop
"""
macro iterateLinearly(iterLoop,lengthh, ex)
  return  esc(quote
  i = UInt32(0)
  @unroll for j in 0:($iterLoop-1)
    i= threadIdxX()+(threadIdxY()-1)*blockDimX()+ j* blockDimX()*blockDimY()
    $ex
  end
    i= threadIdxX()+(threadIdxY()-1)*blockDimX()+ $iterLoop* blockDimX()*blockDimY()
      if(i<=$lengthh)
      $ex
    end
end)

end

"""
iterate over array that is treated as one dimensional with given length lengthh as argument
amount of iterations needed is also passed as an argument - iterLoop
"""
macro iterateLinearlyCheckAll(iterLoop,lengthh, ex)
  return  esc(quote
  i = UInt32(0)
  @unroll for j in 0:($iterLoop)
    i= threadIdxX()+(threadIdxY()-1)*blockDimX()+ j* blockDimX()*blockDimY()
    if(i<=$lengthh)
      $ex
    end
    end
end)

end

"""
as x,y,z positions are not important we can just treat the array as it would be linear
loop will be divided into two parts in order to reduce number of if statements
each thread block iterates over just one slice in order to enable also slice wise results
iterLoop - how many times block needs to loop through slice to cover all
pixPerSlice - pixels per slice
"""
macro iterateLinearlyMultipleBlocks(iterLoop,pixPerSlice,totalNumbOfVoxels, ex)
    return  esc(quote
    i = UInt32(0)
    @unroll for j in 0:($iterLoop-1)
        i= threadIdxX()+(threadIdxY()-1)*blockDimX()+ j* blockDimX()*blockDimY()+ ((blockIdxX()-1) *$pixPerSlice)
            if(i<=$totalNumbOfVoxels)
                $ex
            end
    end
      i= threadIdxX()+(threadIdxY()-1)*blockDimX()+ $iterLoop* blockDimX()*blockDimY()
      if(i<=$pixPerSlice )
        i= i+ ((blockIdxX()-1) *$pixPerSlice)
        if(i<=$totalNumbOfVoxels)
            $ex
        end
      end
  end)

  end

"""
iterate over array that is treated as one dimensional with given length lengthh as argument
amount of iterations needed is also passed as an argument - iterLoop

modification that wil be used when very short arrays are possibility
"""
macro iterateLinearlyWithStart(start,iterLoop,lengthh, ex)
  return  esc(quote
  i = UInt32(0)
  @unroll for j in $start:($iterLoop-1)
    i= threadIdxX()+(threadIdxY()-1)*blockDimX()+ j* blockDimX()*blockDimY()
    $ex
  end
    i= threadIdxX()+(threadIdxY()-1)*blockDimX()+ $iterLoop* blockDimX()*blockDimY()
      if(i<=$lengthh)
      $ex
    end
end)

end

"""
iteration through 3 dimensional data but  with one dimension fixed - ie we will analyze plane or part of the plane where value of given dimension is as we supply it
this can be used to analyze for example padding in the stencil - of course it may lead to non coalesced memory access
chosenDim (1,2 or 3) - the dimension of choice through which we are NOT iterating
dimValue - value of this dimension that we will keep fixed
        DimA, DimB - those dimensions will be chosen to iterate over (macro will analyze those on its own - no need to supply them explicitely)
              - in case chosenDim = 1  DimA, DimB will be 2,3
              - in case chosenDim = 2  DimA, DimB will be 1,3
              - in case chosenDim = 3  DimA, DimB will be 1,2
maxDimA, maxDimB - maximum values available in source array
loopDimA, loopDimB - how many times we should iterate in those dimensions

"""
macro iterDimFixed(chosenDim, dimValue,maxDimA, maxDimB,loopDimA, loopDimB  )

end#iterDimFixed

end#module

# function generalizedIter3d( arrDims, loopXdim, loopYdim,loopZdim,zOffset,zAdd ,ex  )
#   quote
#     @unroll for zdim::UInt32 in 0:$loopZdim
#       z::UInt32 = (zdim*gridDim().x) + blockIdxX()#multiply by blocks to avoid situation that a single block of threads would have no work to do
#       if( (z<= $arrDims[3]) )
#           @unroll for ydim::UInt32 in  0:$loopYdim
#               y::UInt32 = ((ydim* blockDimY()) +threadIdxY())
#                 if( (y<=$arrDims[2])  )
#                   @unroll for xdim::UInt32 in 0:$loopXdim
#                       x::UInt32=(xdim* blockDimX()) +threadIdxX()
#                       if((x <=$arrDims[1]))
#                         $ex
#                       end#if x
#                  end#for x dim
#             end#if y
#           end#for  yLoops
#       end#if z
#   end#for z dim
#   end

# function generalizedIter3dFullChecks( arrDims, loopXdim, loopYdim,loopZdim,zOffset,zAdd,yOffset,xOffset,ex  )
#   quote
#     @unroll for zdim::UInt32 in 0:$loopZdim
#       z::UInt32 = $zOffset + $zAdd#multiply by blocks to avoid situation that a single block of threads would have no work to do
#       if( (z<= $arrDims[3]) )
#           @unroll for ydim::UInt32 in  0:$loopYdim
#               y::UInt32 = $yOffset +threadIdxY()
#                 if( (y<=$arrDims[2])  )
#                   @unroll for xdim::UInt32 in 0:$loopXdim
#                       x::UInt32= $xOffset +threadIdxX()
#                       if((x <=$arrDims[1]))
#                         $ex
#                       end#if x
#                  end#for x dim
#             end#if y
#           end#for  yLoops
#       end#if z
#   end#for z dim
#   end
# end#generalizedIter3d

# quote
#   @unroll for xdim::UInt32 in 0:($loopXdim-1) # we subtract one as we are intrsted only for those iterations that will not need to be bound checked
#     $xname::UInt32= $xOffset +threadIdxX()
#      $ex
#       $additionalActionAfterX
#     end#for x dim
#     if( $xCheck)
#         $ex
#     end#if x
#      $additionalActionAfterX

#     end#quote

# exp2= quote
#     @unroll for ydim::UInt32 in  0:($loopYdim-1)  # we subtract one as we are intrsted only for those iterations that will not need to be bound checked
#      $yname::UInt32 = $yOffset +threadIdxY()

#         $exp1
#         $additionalActionAfterY
#     end#for  yLoops

#     if($yCheck)
#           $exp1
#     end#if y
#     $additionalActionAfterY
#   end

# exp3= quote
# @unroll for zdim::UInt32 in 0:$loopZdim # we subtract one as we are intrsted only for those iterations that will not need to be bound checked
#   $zname::UInt32 = $zOffset + $zAdd#multiply by blocks to avoid situation that a single block of threads would have no work to do
#     $exp2
#     $additionalActionAfterZ
#  end#for z dim
# if($zCheck)
#     $exp2
#   end#if z
#   $additionalActionAfterZ

# function generalizedItermultiDim(; #we keep all as keyword arguments
#   arrDims = (UInt32(1),UInt32(1),UInt32(1) )
#    ,xname = :x
#    ,yname = :y
#    ,zname = :z
#    ,loopIterNameX = :xDim
#    ,loopIterNameY = :yDim
#    ,loopIterNameZ = :zDim
#    ,loopXdim = 1
#    ,loopYdim= 1
#    ,loopZdim= 1
#   ,zOffset= :(($loopIterNameZ*gridDim().x))
#    ,zAdd =:(blockIdxX())
#   ,yOffset = :($loopIterNameY* blockDimY())
#   ,yAdd= :(threadIdxY())
#   ,xOffset= :($loopIterNameX * blockDimX())
#    ,xAdd= :(threadIdxX())
#   ,xCheck = :((x <= $arrDims[1]))
#   ,yCheck = :((y <= $arrDims[2]))
#   ,zCheck = :((z <= $arrDims[3]))
#   ,additionalActionAfterZ= :()
#    ,additionalActionAfterY= :()
#    ,additionalActionAfterX = :()
#    , is3d = true
#    ,ex
#    ,isFullBoundaryCheckX =false
#    , isFullBoundaryCheckY=false
#    , isFullBoundaryCheckZ=false
#    ,nobundaryCheckX=false
#    , nobundaryCheckY=false
#    , nobundaryCheckZ =false)
# #we will define expressions from deepest to most superficial

#   # xState = :(x= $xOffset +$xAdd)
#   # yState = :(y= $yOffset +$yAdd)
#   # zState = :(z= $zOffset +$zAdd)
#   # xState = :($xname= $xOffset +$xAdd)
#   # yState = :($yname= $yOffset +$yAdd)
#   # zState = :($zname= $zOffset +$zAdd)
#   xState = :($xname)
#   yState = :($yname)
#   zState = :($zname)

#   xCheck = :(true)
#   yCheck = :(true)
#   zCheck = :(true)

#   exp1 = getSubLoopPartialCheck(loopIterNameX,loopXdim,xState, additionalActionAfterX ,xCheck ,isFullBoundaryCheckX, nobundaryCheckX,ex)
#   exp2 = getSubLoopPartialCheck(loopIterNameY,loopYdim,yState, additionalActionAfterY ,yCheck ,isFullBoundaryCheckY, nobundaryCheckY,exp1)
#   exp3 = getSubLoopPartialCheck(loopIterNameZ,loopZdim,zState, additionalActionAfterZ ,zCheck ,isFullBoundaryCheckZ, nobundaryCheckZ,exp2)

#   #if it is 3d we will return x,y,z iterations if not only x and y
#   # if(is3d)
#   #   return exp3
#   # else
#   #   return exp2
#   # end

#   oneCheckEx = quote

#     CUDA.@cuprint("     aaaaaaaaaaa  \n")
#     CUDA.@cuprint("     aaaaaaaaaaa  \n")
#     CUDA.@cuprint("     aaaaaaaaaaa  \n")
#     CUDA.@cuprint("     aaaaaaaaaaa  \n")

#     @unroll for $loopIterNameX::UInt32 in  0:($loopXdim-1)  # we subtract one as we are intrsted only for those iterations that will not need to be bound checked
#       $xname= $xOffset +$xAdd
#         $ex
#         $additionalActionAfterX
#     end#for
#     if($xCheck)
#           $ex
#     end#if
#     $additionalActionAfterX
#   end
# return oneCheckEx

# end#generalizedIter3d


text:
using Statistics,Serialization,StatsPlots, DataFrames

# paras = ["no_pol","dist","dist_ho"]
badpol = true # do you have data for badpol and switches?
show_capacity = false # for infected plot

filename_prefix = "../figures/" # store the plots here

# set intial values as in coivd par ini
datat = 7 # interval of data collection
T = 744 #582  # time periods

# define period of plotting for economic plots
start_week = 2 # cut of first week
end_week = 102
# set periods for covid plots
inf_start_day = 1
inf_end_day = 600 # cut of after period 400
inf_start_week = 1
inf_end_week = 86 # 58*7 = 406 and 42*7 = 294

datapoint = fld(T,datat)
include("../data/vaccination_rollout/beta_l50_60//covid_par_ini.jl") # necessary for some basic parameter values (nhh)

plfull = plot()

pcts = ["100", "40", "25", "0"]

for p in ["beta_l5_", "beta_l50_", "beta_l5_alpha_l050_", "beta_l5_alpha_l050_alpha_o050_"]
    # scatterplot
    gdploss_means = []
    totcas_means = []
    pubacc_means = []

    for pct in pcts
        prefix = "../data/vaccination_rollout/$p$pct/"
        # store data scatterplot
        worker_results = deserialize(open("$(prefix)batchdata.dat"))
        gdploss = mean((results -> results[:togdploss]).(worker_results))
        totcas = mean(100*(results -> results[:totcas]).(worker_results))
        pubacch = (results -> results[:pubacctraj]).(worker_results)
        gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
        global nn = size(gdpmean)[1]
        pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
        push!(gdploss_means,gdploss)
        push!(totcas_means,totcas)
        push!(pubacc_means,pubacc)

        if "$p$pct" == "beta_l50_100"
            global point_A_x = gdploss
            global point_A_y = totcas
        end

        if "$p$pct" == "beta_l5_100"
            global point_B_x = gdploss
            global point_B_y = totcas
        end

        if "$p$pct" == "beta_l5_alpha_l050_100"
            global point_C_x = gdploss
            global point_C_y = totcas
        end

        if "$p$pct" == "beta_l5_alpha_l050_alpha_o050_100"
            global point_D_x = gdploss
            global point_D_y = totcas
        end

        if "$p$pct" == "beta_l5_alpha_l050_0"
            global point_C_0_x = gdploss
            global point_C_0_y = totcas
        end

        if "$p$pct" == "beta_l5_alpha_l050_25"
            global point_C_25_x = gdploss
            global point_C_25_y = totcas
        end

        if "$p$pct" == "beta_l5_alpha_l050_40"
            global point_C_40_x = gdploss
            global point_C_40_y = totcas
        end

        if "$p$pct" == "beta_l5_alpha_l050_100"
            global point_C_100_x = gdploss
            global point_C_100_y = totcas
        end

        #if "$p$pct" == "100_100"
        #    global point_E_x = gdploss
        #    global point_E_y = totcas
        #end
    end

    plot!(plfull, gdploss_means, totcas_means, linewidth = 2, linecolor = [:black], legend=false, xlabel = "average GDP loss [%]", ylabel = "mortality [%]") #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
    plot!(plfull, gdploss_means, totcas_means, linewidth = 2, seriestype = :scatter, color = [:black], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

end

### add specific points A B C
scatter!([point_A_x], [point_A_y-0.001], markersize = 0, markercolor = [:white], series_annotations = text.(["A"], :top))

scatter!([point_B_x], [point_B_y-0.001], markersize = 0, markercolor = [:white], series_annotations = text.(["B"], :top))

scatter!([point_C_x], [point_C_y-0.003], markersize = 0, markercolor = [:white], series_annotations = text.(["C"], :top))

scatter!([point_D_x], [point_D_y-0.001], markersize = 0, markercolor = [:white], series_annotations = text.(["D"], :top))

scatter!([point_C_0_x+0.08], [point_C_0_y+0.003], markersize = 0, markercolor = [:white], series_annotations = text.(["0%"], :top, 5))
scatter!([point_C_25_x+0.1], [point_C_25_y+0.004], markersize = 0, markercolor = [:white], series_annotations = text.(["25%"], :top, 5))
scatter!([point_C_40_x+0.08], [point_C_40_y+0.007], markersize = 0, markercolor = [:white], series_annotations = text.(["40%"], :top, 5))
scatter!([point_C_100_x-0.03], [point_C_100_y+0.007], markersize = 0, markercolor = [:white], series_annotations = text.(["100%"], :top, 5))

ylims!(plfull,0.005,0.25)
xlims!(plfull,0,4)

#ylims!(plfull,0.005,0.06)
#xlims!(plfull,0,8)

savefig(plfull,"$(filename_prefix)vaccination_rollout.pdf")

#ylims!(plarrows,0.005,0.13)
#xlims!(plarrows,0,12)


text:
using Signals
using Base.Test

include("signal.jl")
include("interp.jl")
include("window.jl")


text:
###
### **acb_poly.h** -- polynomials over the complex numbers
###

### Types, macros and constants

### Memory management
arbcall"void acb_poly_init(acb_poly_t poly)"
arbcall"void acb_poly_clear(acb_poly_t poly)"
arbcall"void acb_poly_fit_length(acb_poly_t poly, slong len)"
arbcall"void _acb_poly_set_length(acb_poly_t poly, slong len)"
arbcall"void _acb_poly_normalise(acb_poly_t poly)"
arbcall"void acb_poly_swap(acb_poly_t poly1, acb_poly_t poly2)"
arbcall"slong acb_poly_allocated_bytes(const acb_poly_t x)"

### Basic properties and manipulation
#mo arbcall"slong acb_poly_length(const acb_poly_t poly)" # clashes with Base.length
arbcall"slong acb_poly_degree(const acb_poly_t poly)"
arbcall"int acb_poly_is_zero(const acb_poly_t poly)"
arbcall"int acb_poly_is_one(const acb_poly_t poly)"
arbcall"int acb_poly_is_x(const acb_poly_t poly)"
arbcall"void acb_poly_zero(acb_poly_t poly)"
arbcall"void acb_poly_one(acb_poly_t poly)"
arbcall"void acb_poly_set(acb_poly_t dest, const acb_poly_t src)"
arbcall"void acb_poly_set_round(acb_poly_t dest, const acb_poly_t src, slong prec)"
arbcall"void acb_poly_set_trunc(acb_poly_t dest, const acb_poly_t src, slong n)"
arbcall"void acb_poly_set_trunc_round(acb_poly_t dest, const acb_poly_t src, slong n, slong prec)"
arbcall"void acb_poly_set_coeff_si(acb_poly_t poly, slong n, slong c)"
arbcall"void acb_poly_set_coeff_acb(acb_poly_t poly, slong n, const acb_t c)"
arbcall"void acb_poly_get_coeff_acb(acb_t v, const acb_poly_t poly, slong n)"
arbcall"void _acb_poly_shift_right(acb_ptr res, acb_srcptr poly, slong len, slong n)"
arbcall"void acb_poly_shift_right(acb_poly_t res, const acb_poly_t poly, slong n)"
arbcall"void _acb_poly_shift_left(acb_ptr res, acb_srcptr poly, slong len, slong n)"
arbcall"void acb_poly_shift_left(acb_poly_t res, const acb_poly_t poly, slong n)"
arbcall"void acb_poly_truncate(acb_poly_t poly, slong n)"
arbcall"slong acb_poly_valuation(const acb_poly_t poly)"

### Input and output
arbcall"void acb_poly_printd(const acb_poly_t poly, slong digits)"
#ns arbcall"void acb_poly_fprintd(FILE * file, const acb_poly_t poly, slong digits)"

### Random generation
#ns arbcall"void acb_poly_randtest(acb_poly_t poly, flint_rand_t state, slong len, slong prec, slong mag_bits)"

### Comparisons
arbcall"int acb_poly_equal(const acb_poly_t A, const acb_poly_t B)"
arbcall"int acb_poly_contains(const acb_poly_t poly1, const acb_poly_t poly2)"
#ni arbcall"int acb_poly_contains_fmpz_poly(const acb_poly_t poly1, const fmpz_poly_t poly2)"
#ni arbcall"int acb_poly_contains_fmpq_poly(const acb_poly_t poly1, const fmpq_poly_t poly2)"
arbcall"int _acb_poly_overlaps(acb_srcptr poly1, slong len1, acb_srcptr poly2, slong len2)"
arbcall"int acb_poly_overlaps(const acb_poly_t poly1, const acb_poly_t poly2)"
#ni arbcall"int acb_poly_get_unique_fmpz_poly(fmpz_poly_t z, const acb_poly_t x)"
arbcall"int acb_poly_is_real(const acb_poly_t poly)"

### Conversions
#ni arbcall"void acb_poly_set_fmpz_poly(acb_poly_t poly, const fmpz_poly_t re, slong prec)"
#ni arbcall"void acb_poly_set2_fmpz_poly(acb_poly_t poly, const fmpz_poly_t re, const fmpz_poly_t im, slong prec)"
arbcall"void acb_poly_set_arb_poly(acb_poly_t poly, const arb_poly_t re)"
arbcall"void acb_poly_set2_arb_poly(acb_poly_t poly, const arb_poly_t re, const arb_poly_t im)"
#ni arbcall"void acb_poly_set_fmpq_poly(acb_poly_t poly, const fmpq_poly_t re, slong prec)"
#ni arbcall"void acb_poly_set2_fmpq_poly(acb_poly_t poly, const fmpq_poly_t re, const fmpq_poly_t im, slong prec)"
arbcall"void acb_poly_set_acb(acb_poly_t poly, const acb_t src)"
arbcall"void acb_poly_set_si(acb_poly_t poly, slong src)"

### Bounds
arbcall"void _acb_poly_majorant(arb_ptr res, acb_srcptr poly, slong len, slong prec)"
arbcall"void acb_poly_majorant(arb_poly_t res, const acb_poly_t poly, slong prec)"

### Arithmetic
#mo arbcall"void _acb_poly_add(acb_ptr C, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong prec)" # clashes with _acb_vec_add
arbcall"void acb_poly_add(acb_poly_t C, const acb_poly_t A, const acb_poly_t B, slong prec)"
arbcall"void acb_poly_add_si(acb_poly_t C, const acb_poly_t A, slong B, slong prec)"
#mo arbcall"void _acb_poly_sub(acb_ptr C, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong prec)" # clashes with _acb_vec_sub
arbcall"void acb_poly_sub(acb_poly_t C, const acb_poly_t A, const acb_poly_t B, slong prec)"
arbcall"void acb_poly_add_series(acb_poly_t C, const acb_poly_t A, const acb_poly_t B, slong len, slong prec)"
arbcall"void acb_poly_sub_series(acb_poly_t C, const acb_poly_t A, const acb_poly_t B, slong len, slong prec)"
arbcall"void acb_poly_neg(acb_poly_t C, const acb_poly_t A)"
arbcall"void acb_poly_scalar_mul_2exp_si(acb_poly_t C, const acb_poly_t A, slong c)"
arbcall"void acb_poly_scalar_mul(acb_poly_t C, const acb_poly_t A, const acb_t c, slong prec)"
arbcall"void acb_poly_scalar_div(acb_poly_t C, const acb_poly_t A, const acb_t c, slong prec)"
arbcall"void _acb_poly_mullow_classical(acb_ptr C, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong n, slong prec)"
arbcall"void _acb_poly_mullow_transpose(acb_ptr C, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong n, slong prec)"
arbcall"void _acb_poly_mullow_transpose_gauss(acb_ptr C, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong n, slong prec)"
arbcall"void _acb_poly_mullow(acb_ptr C, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong n, slong prec)"
arbcall"void acb_poly_mullow_classical(acb_poly_t C, const acb_poly_t A, const acb_poly_t B, slong n, slong prec)"
arbcall"void acb_poly_mullow_transpose(acb_poly_t C, const acb_poly_t A, const acb_poly_t B, slong n, slong prec)"
arbcall"void acb_poly_mullow_transpose_gauss(acb_poly_t C, const acb_poly_t A, const acb_poly_t B, slong n, slong prec)"
arbcall"void acb_poly_mullow(acb_poly_t C, const acb_poly_t A, const acb_poly_t B, slong n, slong prec)"
arbcall"void _acb_poly_mul(acb_ptr C, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong prec)"
arbcall"void acb_poly_mul(acb_poly_t C, const acb_poly_t A1, const acb_poly_t B2, slong prec)"
arbcall"void _acb_poly_inv_series(acb_ptr Qinv, acb_srcptr Q, slong Qlen, slong len, slong prec)"
arbcall"void acb_poly_inv_series(acb_poly_t Qinv, const acb_poly_t Q, slong n, slong prec)"
arbcall"void _acb_poly_div_series(acb_ptr Q, acb_srcptr A, slong Alen, acb_srcptr B, slong Blen, slong n, slong prec)"
arbcall"void acb_poly_div_series(acb_poly_t Q, const acb_poly_t A, const acb_poly_t B, slong n, slong prec)"
arbcall"void _acb_poly_div(acb_ptr Q, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong prec)"
arbcall"void _acb_poly_rem(acb_ptr R, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong prec)"
arbcall"void _acb_poly_divrem(acb_ptr Q, acb_ptr R, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong prec)"
arbcall"int acb_poly_divrem(acb_poly_t Q, acb_poly_t R, const acb_poly_t A, const acb_poly_t B, slong prec)"
arbcall"void _acb_poly_div_root(acb_ptr Q, acb_t R, acb_srcptr A, slong len, const acb_t c, slong prec)"

### Composition
arbcall"void _acb_poly_taylor_shift_horner(acb_ptr g, const acb_t c, slong n, slong prec)"
arbcall"void acb_poly_taylor_shift_horner(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec)"
arbcall"void _acb_poly_taylor_shift_divconquer(acb_ptr g, const acb_t c, slong n, slong prec)"
arbcall"void acb_poly_taylor_shift_divconquer(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec)"
arbcall"void _acb_poly_taylor_shift_convolution(acb_ptr g, const acb_t c, slong n, slong prec)"
arbcall"void acb_poly_taylor_shift_convolution(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec)"
arbcall"void _acb_poly_taylor_shift(acb_ptr g, const acb_t c, slong n, slong prec)"
arbcall"void acb_poly_taylor_shift(acb_poly_t g, const acb_poly_t f, const acb_t c, slong prec)"
arbcall"void _acb_poly_compose_horner(acb_ptr res, acb_srcptr poly1, slong len1, acb_srcptr poly2, slong len2, slong prec)"
arbcall"void acb_poly_compose_horner(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong prec)"
arbcall"void _acb_poly_compose_divconquer(acb_ptr res, acb_srcptr poly1, slong len1, acb_srcptr poly2, slong len2, slong prec)"
arbcall"void acb_poly_compose_divconquer(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong prec)"
arbcall"void _acb_poly_compose(acb_ptr res, acb_srcptr poly1, slong len1, acb_srcptr poly2, slong len2, slong prec)"
arbcall"void acb_poly_compose(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong prec)"
arbcall"void _acb_poly_compose_series_horner(acb_ptr res, acb_srcptr poly1, slong len1, acb_srcptr poly2, slong len2, slong n, slong prec)"
arbcall"void acb_poly_compose_series_horner(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong n, slong prec)"
arbcall"void _acb_poly_compose_series_brent_kung(acb_ptr res, acb_srcptr poly1, slong len1, acb_srcptr poly2, slong len2, slong n, slong prec)"
arbcall"void acb_poly_compose_series_brent_kung(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong n, slong prec)"
arbcall"void _acb_poly_compose_series(acb_ptr res, acb_srcptr poly1, slong len1, acb_srcptr poly2, slong len2, slong n, slong prec)"
arbcall"void acb_poly_compose_series(acb_poly_t res, const acb_poly_t poly1, const acb_poly_t poly2, slong n, slong prec)"
arbcall"void _acb_poly_revert_series_lagrange(acb_ptr h, acb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void acb_poly_revert_series_lagrange(acb_poly_t h, const acb_poly_t f, slong n, slong prec)"
arbcall"void _acb_poly_revert_series_newton(acb_ptr h, acb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void acb_poly_revert_series_newton(acb_poly_t h, const acb_poly_t f, slong n, slong prec)"
arbcall"void _acb_poly_revert_series_lagrange_fast(acb_ptr h, acb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void acb_poly_revert_series_lagrange_fast(acb_poly_t h, const acb_poly_t f, slong n, slong prec)"
arbcall"void _acb_poly_revert_series(acb_ptr h, acb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void acb_poly_revert_series(acb_poly_t h, const acb_poly_t f, slong n, slong prec)"

### Evaluation
arbcall"void _acb_poly_evaluate_horner(acb_t y, acb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void acb_poly_evaluate_horner(acb_t y, const acb_poly_t f, const acb_t x, slong prec)"
arbcall"void _acb_poly_evaluate_rectangular(acb_t y, acb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void acb_poly_evaluate_rectangular(acb_t y, const acb_poly_t f, const acb_t x, slong prec)"
arbcall"void _acb_poly_evaluate(acb_t y, acb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void acb_poly_evaluate(acb_t y, const acb_poly_t f, const acb_t x, slong prec)"
arbcall"void _acb_poly_evaluate2_horner(acb_t y, acb_t z, acb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void acb_poly_evaluate2_horner(acb_t y, acb_t z, const acb_poly_t f, const acb_t x, slong prec)"
arbcall"void _acb_poly_evaluate2_rectangular(acb_t y, acb_t z, acb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void acb_poly_evaluate2_rectangular(acb_t y, acb_t z, const acb_poly_t f, const acb_t x, slong prec)"
arbcall"void _acb_poly_evaluate2(acb_t y, acb_t z, acb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void acb_poly_evaluate2(acb_t y, acb_t z, const acb_poly_t f, const acb_t x, slong prec)"

### Product trees
arbcall"void _acb_poly_product_roots(acb_ptr poly, acb_srcptr xs, slong n, slong prec)"
arbcall"void acb_poly_product_roots(acb_poly_t poly, acb_srcptr xs, slong n, slong prec)"
#ni arbcall"acb_ptr * _acb_poly_tree_alloc(slong len)"
#ni arbcall"void _acb_poly_tree_free(acb_ptr * tree, slong len)"
#ni arbcall"void _acb_poly_tree_build(acb_ptr * tree, acb_srcptr roots, slong len, slong prec)"

### Multipoint evaluation
arbcall"void _acb_poly_evaluate_vec_iter(acb_ptr ys, acb_srcptr poly, slong plen, acb_srcptr xs, slong n, slong prec)"
arbcall"void acb_poly_evaluate_vec_iter(acb_ptr ys, const acb_poly_t poly, acb_srcptr xs, slong n, slong prec)"
#ni arbcall"void _acb_poly_evaluate_vec_fast_precomp(acb_ptr vs, acb_srcptr poly, slong plen, acb_ptr * tree, slong len, slong prec)"
arbcall"void _acb_poly_evaluate_vec_fast(acb_ptr ys, acb_srcptr poly, slong plen, acb_srcptr xs, slong n, slong prec)"
arbcall"void acb_poly_evaluate_vec_fast(acb_ptr ys, const acb_poly_t poly, acb_srcptr xs, slong n, slong prec)"

### Interpolation
arbcall"void _acb_poly_interpolate_newton(acb_ptr poly, acb_srcptr xs, acb_srcptr ys, slong n, slong prec)"
arbcall"void acb_poly_interpolate_newton(acb_poly_t poly, acb_srcptr xs, acb_srcptr ys, slong n, slong prec)"
arbcall"void _acb_poly_interpolate_barycentric(acb_ptr poly, acb_srcptr xs, acb_srcptr ys, slong n, slong prec)"
arbcall"void acb_poly_interpolate_barycentric(acb_poly_t poly, acb_srcptr xs, acb_srcptr ys, slong n, slong prec)"
#ni arbcall"void _acb_poly_interpolation_weights(acb_ptr w, acb_ptr * tree, slong len, slong prec)"
#ni arbcall"void _acb_poly_interpolate_fast_precomp(acb_ptr poly, acb_srcptr ys, acb_ptr * tree, acb_srcptr weights, slong len, slong prec)"
arbcall"void _acb_poly_interpolate_fast(acb_ptr poly, acb_srcptr xs, acb_srcptr ys, slong len, slong prec)"
arbcall"void acb_poly_interpolate_fast(acb_poly_t poly, acb_srcptr xs, acb_srcptr ys, slong n, slong prec)"

### Differentiation
arbcall"void _acb_poly_derivative(acb_ptr res, acb_srcptr poly, slong len, slong prec)"
arbcall"void acb_poly_derivative(acb_poly_t res, const acb_poly_t poly, slong prec)"
arbcall"void _acb_poly_integral(acb_ptr res, acb_srcptr poly, slong len, slong prec)"
arbcall"void acb_poly_integral(acb_poly_t res, const acb_poly_t poly, slong prec)"

### Transforms
arbcall"void _acb_poly_borel_transform(acb_ptr res, acb_srcptr poly, slong len, slong prec)"
arbcall"void acb_poly_borel_transform(acb_poly_t res, const acb_poly_t poly, slong prec)"
arbcall"void _acb_poly_inv_borel_transform(acb_ptr res, acb_srcptr poly, slong len, slong prec)"
arbcall"void acb_poly_inv_borel_transform(acb_poly_t res, const acb_poly_t poly, slong prec)"
arbcall"void _acb_poly_binomial_transform_basecase(acb_ptr b, acb_srcptr a, slong alen, slong len, slong prec)"
arbcall"void acb_poly_binomial_transform_basecase(acb_poly_t b, const acb_poly_t a, slong len, slong prec)"
arbcall"void _acb_poly_binomial_transform_convolution(acb_ptr b, acb_srcptr a, slong alen, slong len, slong prec)"
arbcall"void acb_poly_binomial_transform_convolution(acb_poly_t b, const acb_poly_t a, slong len, slong prec)"
arbcall"void _acb_poly_binomial_transform(acb_ptr b, acb_srcptr a, slong alen, slong len, slong prec)"
arbcall"void acb_poly_binomial_transform(acb_poly_t b, const acb_poly_t a, slong len, slong prec)"

### Elementary functions
arbcall"void _acb_poly_pow_ui_trunc_binexp(acb_ptr res, acb_srcptr f, slong flen, ulong exp, slong len, slong prec)"
arbcall"void acb_poly_pow_ui_trunc_binexp(acb_poly_t res, const acb_poly_t poly, ulong exp, slong len, slong prec)"
arbcall"void _acb_poly_pow_ui(acb_ptr res, acb_srcptr f, slong flen, ulong exp, slong prec)"
arbcall"void acb_poly_pow_ui(acb_poly_t res, const acb_poly_t poly, ulong exp, slong prec)"
arbcall"void _acb_poly_pow_series(acb_ptr h, acb_srcptr f, slong flen, acb_srcptr g, slong glen, slong len, slong prec)"
arbcall"void acb_poly_pow_series(acb_poly_t h, const acb_poly_t f, const acb_poly_t g, slong len, slong prec)"
arbcall"void _acb_poly_pow_acb_series(acb_ptr h, acb_srcptr f, slong flen, const acb_t g, slong len, slong prec)"
arbcall"void acb_poly_pow_acb_series(acb_poly_t h, const acb_poly_t f, const acb_t g, slong len, slong prec)"
arbcall"void _acb_poly_sqrt_series(acb_ptr g, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sqrt_series(acb_poly_t g, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_rsqrt_series(acb_ptr g, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_rsqrt_series(acb_poly_t g, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_log_series(acb_ptr res, acb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void acb_poly_log_series(acb_poly_t res, const acb_poly_t f, slong n, slong prec)"
arbcall"void _acb_poly_log1p_series(acb_ptr res, acb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void acb_poly_log1p_series(acb_poly_t res, const acb_poly_t f, slong n, slong prec)"
arbcall"void _acb_poly_atan_series(acb_ptr res, acb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void acb_poly_atan_series(acb_poly_t res, const acb_poly_t f, slong n, slong prec)"
arbcall"void _acb_poly_exp_series_basecase(acb_ptr f, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_exp_series_basecase(acb_poly_t f, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_exp_series(acb_ptr f, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_exp_series(acb_poly_t f, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_exp_pi_i_series(acb_ptr f, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_exp_pi_i_series(acb_poly_t f, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_sin_cos_series_basecase(acb_ptr s, acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec, int times_pi)"
arbcall"void acb_poly_sin_cos_series_basecase(acb_poly_t s, acb_poly_t c, const acb_poly_t h, slong n, slong prec, int times_pi)"
arbcall"void _acb_poly_sin_cos_series_tangent(acb_ptr s, acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec, int times_pi)"
arbcall"void acb_poly_sin_cos_series_tangent(acb_poly_t s, acb_poly_t c, const acb_poly_t h, slong n, slong prec, int times_pi)"
arbcall"void _acb_poly_sin_cos_series(acb_ptr s, acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sin_cos_series(acb_poly_t s, acb_poly_t c, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_sin_series(acb_ptr s, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sin_series(acb_poly_t s, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_cos_series(acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_cos_series(acb_poly_t c, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_tan_series(acb_ptr g, acb_srcptr h, slong hlen, slong len, slong prec)"
arbcall"void acb_poly_tan_series(acb_poly_t g, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_sin_cos_pi_series(acb_ptr s, acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sin_cos_pi_series(acb_poly_t s, acb_poly_t c, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_sin_pi_series(acb_ptr s, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sin_pi_series(acb_poly_t s, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_cos_pi_series(acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_cos_pi_series(acb_poly_t c, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_cot_pi_series(acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_cot_pi_series(acb_poly_t c, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_sinh_cosh_series_basecase(acb_ptr s, acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sinh_cosh_series_basecase(acb_poly_t s, acb_poly_t c, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_sinh_cosh_series_exponential(acb_ptr s, acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sinh_cosh_series_exponential(acb_poly_t s, acb_poly_t c, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_sinh_cosh_series(acb_ptr s, acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sinh_cosh_series(acb_poly_t s, acb_poly_t c, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_sinh_series(acb_ptr s, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sinh_series(acb_poly_t s, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_cosh_series(acb_ptr c, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_cosh_series(acb_poly_t c, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_sinc_series(acb_ptr s, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_sinc_series(acb_poly_t s, const acb_poly_t h, slong n, slong prec)"

### Lambert W function
#ni arbcall"void _acb_poly_lambertw_series(acb_ptr res, acb_srcptr z, slong zlen, const fmpz_t k, int flags, slong len, slong prec)"
#ni arbcall"void acb_poly_lambertw_series(acb_poly_t res, const acb_poly_t z, const fmpz_t k, int flags, slong len, slong prec)"

### Gamma function
arbcall"void _acb_poly_gamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_gamma_series(acb_poly_t res, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_rgamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_rgamma_series(acb_poly_t res, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_lgamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_lgamma_series(acb_poly_t res, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_digamma_series(acb_ptr res, acb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void acb_poly_digamma_series(acb_poly_t res, const acb_poly_t h, slong n, slong prec)"
arbcall"void _acb_poly_rising_ui_series(acb_ptr res, acb_srcptr f, slong flen, ulong r, slong trunc, slong prec)"
arbcall"void acb_poly_rising_ui_series(acb_poly_t res, const acb_poly_t f, ulong r, slong trunc, slong prec)"

### Power sums
arbcall"void _acb_poly_powsum_series_naive(acb_ptr z, const acb_t s, const acb_t a, const acb_t q, slong n, slong len, slong prec)"
arbcall"void _acb_poly_powsum_series_naive_threaded(acb_ptr z, const acb_t s, const acb_t a, const acb_t q, slong n, slong len, slong prec)"
arbcall"void _acb_poly_powsum_one_series_sieved(acb_ptr z, const acb_t s, slong n, slong len, slong prec)"

### Zeta function
arbcall"void _acb_poly_zeta_em_choose_param(mag_t bound, ulong * N, ulong * M, const acb_t s, const acb_t a, slong d, slong target, slong prec)"
arbcall"void _acb_poly_zeta_em_bound1(mag_t bound, const acb_t s, const acb_t a, slong N, slong M, slong d, slong wp)"
arbcall"void _acb_poly_zeta_em_bound(arb_ptr vec, const acb_t s, const acb_t a, ulong N, ulong M, slong d, slong wp)"
arbcall"void _acb_poly_zeta_em_tail_naive(acb_ptr z, const acb_t s, const acb_t Na, acb_srcptr Nasx, slong M, slong len, slong prec)"
arbcall"void _acb_poly_zeta_em_tail_bsplit(acb_ptr z, const acb_t s, const acb_t Na, acb_srcptr Nasx, slong M, slong len, slong prec)"
arbcall"void _acb_poly_zeta_em_sum(acb_ptr z, const acb_t s, const acb_t a, int deflate, ulong N, ulong M, slong d, slong prec)"
arbcall"void _acb_poly_zeta_cpx_series(acb_ptr z, const acb_t s, const acb_t a, int deflate, slong d, slong prec)"
arbcall"void _acb_poly_zeta_series(acb_ptr res, acb_srcptr h, slong hlen, const acb_t a, int deflate, slong len, slong prec)"
arbcall"void acb_poly_zeta_series(acb_poly_t res, const acb_poly_t f, const acb_t a, int deflate, slong n, slong prec)"

### Other special functions
arbcall"void _acb_poly_polylog_cpx_small(acb_ptr w, const acb_t s, const acb_t z, slong len, slong prec)"
arbcall"void _acb_poly_polylog_cpx_zeta(acb_ptr w, const acb_t s, const acb_t z, slong len, slong prec)"
arbcall"void _acb_poly_polylog_cpx(acb_ptr w, const acb_t s, const acb_t z, slong len, slong prec)"
arbcall"void _acb_poly_polylog_series(acb_ptr w, acb_srcptr s, slong slen, const acb_t z, slong len, slong prec)"
arbcall"void acb_poly_polylog_series(acb_poly_t w, const acb_poly_t s, const acb_t z, slong len, slong prec)"
arbcall"void _acb_poly_erf_series(acb_ptr res, acb_srcptr z, slong zlen, slong n, slong prec)"
arbcall"void acb_poly_erf_series(acb_poly_t res, const acb_poly_t z, slong n, slong prec)"
arbcall"void _acb_poly_agm1_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)"
arbcall"void acb_poly_agm1_series(acb_poly_t res, const acb_poly_t z, slong n, slong prec)"
#mo arbcall"void _acb_poly_elliptic_k_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)" # alias to _acb_elliptic_k_series
#mo arbcall"void acb_poly_elliptic_k_series(acb_poly_t res, const acb_poly_t z, slong n, slong prec)" # alias to acb_elliptic_k_series
#mo arbcall"void _acb_poly_elliptic_p_series(acb_ptr res, acb_srcptr z, slong zlen, const acb_t tau, slong len, slong prec)" # alias to _acb_elliptic_p_series
#mo arbcall"void acb_poly_elliptic_p_series(acb_poly_t res, const acb_poly_t z, const acb_t tau, slong n, slong prec)" # alias to acb_elliptic_p_series

### Root-finding
arbcall"void _acb_poly_root_bound_fujiwara(mag_t bound, acb_srcptr poly, slong len)"
arbcall"void acb_poly_root_bound_fujiwara(mag_t bound, acb_poly_t poly)"
arbcall"void _acb_poly_root_inclusion(acb_t r, const acb_t m, acb_srcptr poly, acb_srcptr polyder, slong len, slong prec)"
arbcall"slong _acb_poly_validate_roots(acb_ptr roots, acb_srcptr poly, slong len, slong prec)"
arbcall"void _acb_poly_refine_roots_durand_kerner(acb_ptr roots, acb_srcptr poly, slong len, slong prec)"
arbcall"slong _acb_poly_find_roots(acb_ptr roots, acb_srcptr poly, acb_srcptr initial, slong len, slong maxiter, slong prec)"
arbcall"slong acb_poly_find_roots(acb_ptr roots, const acb_poly_t poly, acb_srcptr initial, slong maxiter, slong prec)"
arbcall"int _acb_poly_validate_real_roots(acb_srcptr roots, acb_srcptr poly, slong len, slong prec)"
arbcall"int acb_poly_validate_real_roots(acb_srcptr roots, const acb_poly_t poly, slong prec)"


text:
using Plots, LaTeXStrings, Statistics, JLD, BenchmarkTools, StatsBase, SpecialFunctions, ProgressMeter

Integral_answer=sqrt(π)/2 * erf(2)

function SSMC(f :: Function, a :: Real, b :: Real, N :: Integer )
    f₍ₓ₎= f.((b-a).*rand(N) .+ a)
    ∫f₍ₓ₎dx= (b-a)*(mean(f₍ₓ₎))
    error=std(f₍ₓ₎)/sqrt(N)
    return ∫f₍ₓ₎dx, error
end

function ISMC(f :: Function, g :: Function, y :: Function, a :: Real, b :: Real, N :: Integer)
    drivation= x-> f(x)/g(x)
    samples= drivation.(y.(rand(N)))
    ∫g₍ₓ₎dx= 1-ℯ^(-2)
    ∫f₍ₓ₎dx= ∫g₍ₓ₎dx * mean(samples)
    error= ∫g₍ₓ₎dx * std(samples)/sqrt(N)
    return ∫f₍ₓ₎dx, error
end

nList=Integer.([10^n for n in 1:8])
SSMC_res=[]
prog= Progress(length(nList))
n=0
for n in nList
    I, error=SSMC(x-> ℯ^(-x^2), 0, 2, n)
    bench = @benchmark SSMC(x->ℯ^(-x^2), 0, 2, n)
    push!(SSMC_res, (I, error, mean(bench).time, bench))
    next!(prog)
end
ISMC_res=[]
n=0
prog= Progress(length(nList))
for n in nList
    I, error=ISMC(x->ℯ^(-x^2), x -> ℯ^(-x), x -> -log(1-(1-ℯ^(-2))*x), 0, 2, n)
    bench = @benchmark ISMC(x->ℯ^(-x^2), x -> ℯ^(-x), x -> -log(1-(1-ℯ^(-2))*x), 0, 2, n)
    push!(ISMC_res, (I, error, mean(bench).time, bench))
    next!(prog)
end
save("../../chapter7/7.1/data.jld","ISMC", ISMC_res, "SSMC", SSMC_res, "samples", nList)
AnswerList=[Integral_answer for i in 1: length(nList)]
scatter(log.(nList), [getindex.(SSMC_res, 1) getindex.(ISMC_res, 1) AnswerList], markercolor=[:blue :purple :orange], label=nothing)
plot!(log.(nList), [getindex.(SSMC_res, 1) getindex.(ISMC_res, 1) AnswerList], linecolor=[:blue :purple :orange],
    dpi=600,labels=["SSMC" "ISMC" "Absolute Result"], ylabel="Integral", xlabel="Number of Samples",
        title=L"Monte\ Carlo\ Results\ for\ I= \int{}_{0}^{2}\ e^{-x^2}\ dx" )
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter7\\7.1\\Figs\\ResultsPlot.png")
scatter(log.(nList), [getindex.(SSMC_res, 2) getindex.(ISMC_res, 2)], markercolor=[:blue :purple], label= nothing)
plot!(log.(nList), [getindex.(SSMC_res, 2) getindex.(ISMC_res, 2)],
     labels=["SSMC" "ISMC"], xlabel="Number of Samples", ylabel=" Error", dpi=600, linecolor=[:blue :purple],
        title=L"Monte\ Carlo\ Errors\ of\ the\ Results\ for\ I= \int{}_{0}^{2}\ e^{-x^2}\ dx")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter7\\7.1\\Figs\\ErrorPlot.png")


text:
export VoxelGrid

"""
    VoxelGrid

Initialize VoxelGrid representation.

`voxels` should be Array of size `(N, N, N, B)` where `N` is the number of
voxels features and `B` is the batch size of VoxelGrid.

### Fields:

- `voxels`      - voxels features of VoxelGrid.

### Available Contructor:

- `VoxelGrid(voxels::AbstractArray)`
- `VoxelGrid(;voxelsAbstractArray)`
- `VoxelGrid(v::VoxelGrid)`
"""
mutable struct VoxelGrid{T<:Float32} <: AbstractObject
    voxels::AbstractArray{T,4}
end

# VoxelGrid(voxels::AbstractArray{Float32,4}) = VoxelGrid(voxels)

VoxelGrid(voxels::AbstractArray{T,4}) where {T} = VoxelGrid(Float32.(voxels))

function VoxelGrid(voxels::AbstractArray{T,3}) where {T}
    voxels = reshape(voxels, size(voxels)..., 1)
    return VoxelGrid(voxels)
end

VoxelGrid(; voxels) = VoxelGrid(voxels)
VoxelGrid(v::VoxelGrid) = VoxelGrid(v.voxels)

@functor VoxelGrid

Base.getindex(v::VoxelGrid, index::Number) = v.voxels[:, :, :, index]

function Base.show(io::IO, m::VoxelGrid{T}) where {T}
    print(
        io,
        "VoxelGrid{$(T)} Structure:\n    Batch size: ",
        size(m.voxels, 4),
        "\n    Voxels features: ",
        size(m.voxels, 1),
        "\n    Storage type: ",
        typeof(m.voxels),
    )
end

_assert_voxel(v::VoxelGrid) = all(0.0 .<= v.voxels .<= 1.0)


text:
## Time Intergral Methods designed for non-parallel programme

abstract TimeIntergral_onestep
abstract ModelParams
export ModelParams

function TimeForward!( p :: ModelParams, t :: Float64, dt :: Float64)
  error("No function for abstract types")
end

function fTrend( data:: Tuple, p:: ModelParams)
  error("No function for abstract types")
end

function eval!(tf :: TimeIntergral_onestep, p :: ModelParams,
  data :: Tuple, trends :: Tuple,
  t :: Float64, dt :: Float64)
  error("No function for abstract types")
end
export TimeForward

include("AdamsBashforth.jl")
include("RK4.jl")
include("DatasTrends.jl")


text:
using Distributions

function block_bootstrap_DGP(dfm::DynamicFactorModel, factor_innovations::Array{Float64, 2}, block_size=10)
end
block_bootstrap_DGP(dfm::DynamicFactorModel) = block_bootstrap_DGP(dfm, apply(randn, size(dfm.x)))

function parametric_bootstrap_DGP(dfm::DynamicFactorModel, factor_innovations::Array{Float64, 2})
    # draw x from a parametric bootstrap with the parameters (estimated factors and loadings) given in dfm
    # Note: this is a DGP for the factor equation (x = factors*loadings' + e where e are the innovations) only
    x = dfm.factors*dfm.loadings' + factor_innovations
end
parametric_bootstrap_DGP(dfm::DynamicFactorModel) = parametric_bootstrap_DGP(dfm, apply(randn, size(dfm.x)))  # if no arguments are applied all innovations are assumed standard normal

# static factor models:

function residual_bootstrap(dfm::DynamicFactorModel, B::Int, stat::Function)  # resample residuals TODO: what to do when we have a break?

    function resample_x_block(from_indices, to_indices)  # resample given indices, function can be used to resample by block of residuals
        dfm.factors[from_indices:to_indices, 1:dfm.number_of_factors]*dfm.loadings[:, 1:dfm.number_of_factors]' + dfm.factor_residuals[rand(Distributions.DiscreteUniform(from_indices, to_indices), length(from_indices:to_indices)), :]
    end
    function resample_x()
        break_points = [1 dfm.break_indices' length(dfm.y)+1]
        resampled_x = apply(vcat, [resample_x_block(break_points[i-1], break_points[i]-1) for i in 2:length(break_points)])
    end

    stats = Array(Float64, B)
    for b in 1:B
        resampled_x = resample_x()
        resampled_y = dfm.y  # TODO: resampling y is less interesting. This would be akin to checking a linear model for a break.
        resampled_w = dfm.w
        stats[b] = stat(DynamicFactorModel(resampled_y, resampled_w, resampled_x, dfm.number_of_factors, dfm.number_of_factors_criterion, dfm.factor_type, dfm.targeted_predictors, dfm.number_of_factor_lags, dfm.break_indices))
    end
    stats
end

function wild_bootstrap(dfm::DynamicFactorModel, B::Int, stat::Function)  # resample residuals and multiply by random Variable with mean 0 and variance 1, do that B times and calculate statistic stat
    stats = Array(Float64, B)
    for b in 1:B
        resampled_factor_residuals = dfm.factor_residuals[rand(Distributions.DiscreteUniform(1, length(dfm.y)), size(dfm.x, 1)), :]
        resampled_x = dfm.factors[:, 1:dfm.number_of_factors]*dfm.loadings[:, 1:dfm.number_of_factors]' + apply(vcat, [resampled_factor_residuals[t, :] .* randn() for t in 1:size(resampled_factor_residuals, 1)])
        resampled_y = dfm.y  # TODO: resampling y is less interesting. This would be akin to checking a linear model for a break.
        resampled_w = dfm.w  # TODO: maybe a function for how to resample w?
        stats[b] = stat(DynamicFactorModel(resampled_y, resampled_w, resampled_x, dfm.number_of_factors, dfm.number_of_factors_criterion, dfm.factor_type, dfm.targeted_predictors, dfm.number_of_factor_lags, dfm.break_indices))
    end
    stats
end

#y, x, f, lambdas, epsilon_x = factor_model_DGP(100, 60, 1; model="Breitung_Eickmeier_2011", b=0)
#w = reshape(ones(length(y)), (length(y), 1))
#dfm = DynamicFactorModel(y,w,x,1)


text:
"""
    dgtripcolor(x, y, z, t; zmin=nothing, zmax=nothing)

Draw a pseudocolor plot of discontinuous unstructured triangular data. `x`, and `y`, are
one-dimensional input arrays of the same length (the number of points). `t` is an integer array of
size `(3,nt)`, where `nt` is the number of triangles. `z` is an array of the same size as `t`
representing function values on each triangle (discontinuous across triangle edges). The coordinates
of the `i`th vertex of the `j`th triangle are given by `(x[t[i,j]], y[t[i,j]])`. The function value
at the `i`th vertex of the `j`th triangle is given by `z[i,j]`.

The colormap is specified by `cmap`. The colormap bounds can be set by specifying values for `zmin`
and `zmax`. If these parameters are set to `nothing`, then the minimum and maximum values of the `z`
array will be used.

For plotting continuous fields, see [`tripcolor`](@ref).
"""
function dgtripcolor(x, y, z, t, cmap::AbstractVector{T}; px, py, zmin=nothing, zmax=nothing, yflip=false, bg=nothing) where T
    nv = length(x)
    nt = size(t,2)
    @assert (size(x) == (nv,) && size(y) == (nv,) && size(z) == (3,nt) && size(t) == (3,nt)) "Bad input sizes"

    c = isnothing(bg) ? Matrix{T}(undef,px,py) : fill(bg,px,py)
    nc = length(cmap)

    xmin,xmax = extrema(x)
    ymin,ymax = extrema(y)
    w = xmax-xmin
    h = ymax-ymin
    isnothing(zmin) && (zmin = minimum(z))
    isnothing(zmax) && (zmax = maximum(z))

    for it=1:nt
        xt = @views x[t[:,it]]
        yt = @views y[t[:,it]]
        zt = @views z[:,it]

        detT = (yt[2]-yt[3])*(xt[1]-xt[3]) + (xt[3]-xt[2])*(yt[1]-yt[3])

        xmint,xmaxt = extrema(xt)
        ymint,ymaxt = extrema(yt)

        x0 = floor(Int,(xmint-xmin)/w*px)
        x1 = ceil(Int,(xmaxt-xmin)/w*px)
        y0 = floor(Int,(ymint-ymin)/h*py)
        y1 = ceil(Int,(ymaxt-ymin)/h*py)

        # Loop over all pixels in triangle's bounding box
        for i=x0:x1-1, j=y0:y1-1
            # Cartesian coordinates of pixel (i,j) (zero-indexed)
            cx = xmin + i/px*w
            cy = ymin + j/py*h
            # Convert from Cartesian to barycentric coordinates
            λ1 = ((yt[2] -yt[3])*(cx-xt[3]) + (xt[3]-xt[2])*(cy-yt[3]))/detT
            λ2 = ((yt[3] -yt[1])*(cx-xt[3]) + (xt[1]-xt[3])*(cy-yt[3]))/detT
            λ3 = 1 - λ1 - λ2
            tol = 1e-12
            # If pixel is outside of triangle, move on to next pixel
            (λ1>-tol && λ2<1+tol && λ2>-tol && λ2<1+tol && λ3>-tol && λ3<1+tol) || continue
            zval = zt[1]*λ1 + zt[2]*λ2 + zt[3]*λ3
            zval = clamp(zval,zmin,zmax)
            c[i+1,yflip ? py-j : j+1] = cmap[begin+getcolorind(zval,zmin,zmax,nc)]
        end
    end
    c
end

"""
    tripcolor(x, y, z, t; zmin=nothing, zmax=nothing)

Draw a pseudocolor plot of unstructured triangular data. `x`, `y`, and `z` are
one-dimensional input arrays of the same length (the number of points). `t` is an integer
array of size `(3,nt)`, where `nt` is the number of triangles. The coordinates of the `i`th
vertex of the `j`th triangle are given by `(x[t[i,j]], y[t[i,j]])`. The function value at the `j`th
vertex of the `i`th triangle is given by `z[i,j]`.

The colormap is specified by `cmap`. The colormap bounds can be set by specifying values for `zmin`
and `zmax`. If these parameters are set to `nothing`, then the minimum and maximum values of the `z`
array will be used.

For plotting discontinuous fields, see [`dgtripcolor`](@ref).
"""
function tripcolor(x, y, z, t, cmap::AbstractVector{T}; px, py, zmin=nothing, zmax=nothing, yflip=false, bg=nothing) where T
    nv = length(x)
    nt = size(t,2)
    @assert (size(x) == (nv,) && size(y) == (nv,) && size(z) == (nv,) && size(t) == (3,nt)) "Bad input sizes"

    dg_z = similar(t, eltype(z))
    for it=1:nt
        dg_z[:,it] .= @views z[t[:,it]]
    end
    dgtripcolor(x, y, dg_z, t, cmap; px, py, zmin, zmax, bg)
end

function getcolorind(level,lmin,lmax,nc)
    cfrac = (lmin == lmax) ? 0 : (level-lmin)/(lmax-lmin)
    round(Int,cfrac*(nc-1))
end


text:
"""
    PackageWatcher(modules) <: Rewriter

Rebuild when a source file in one of `modules` changes.

"""
struct PackageWatcher <: Rewriter
    modules
    handler
end

# FIXME: remove old constructor and handler field
PackageWatcher(modules) = PackageWatcher(modules, 1)

# TODO: update to work with event-based serving

function geteventsource(pkgwatcher::PackageWatcher, ch)
    watcher = LiveServer.SimpleWatcher(filename -> onsourcefilechanged(Path(filename), ch))
    for m in pkgwatcher.modules
        srcdir = joinpath(Path(pkgdir(m)), "src")
        for p in walkpath(srcdir)
            LiveServer.watch_file!(watcher, string(p))
        end
    end
    return watcher
end

struct SourceFileUpdated <: Event
    p::AbstractPath
end

function onsourcefilechanged(p, ch)
    @info "Source code file $p changed, revising and rebuilding..."
    Revise.revise()
    event = SourceFileUpdated(p)
    put!(ch, event)
end

function handle(server, ::ServerMode, ::SourceFileUpdated)
    # clear executecode cache
    for rewriter in server.project.rewriters
        if rewriter isa ExecuteCode
            reset!(rewriter)
        elseif rewriter isa Referencer
            reset!(rewriter)
        end
    end
    for p in keys(server.project.sources)
        # only REbuild
        if p in keys(server.project.outputs)
            addrewrite!(server, p)
            addbuild!(server, p)
        end
    end
end


text:
using DiffusionAnim
using Test

@testset "DiffusionAnim.jl" begin
    # Write your own tests here.
end


text:
function irf_ci_bootstrap(V::VAR, H::Int64, nrep::Int64, i::Bool=true)

  (T,K) = size(V.Y')
  IRFrmat = zeros(nrep, K^2*(H+1))
  u = V.ϵ*sqrt((T-V.p)/(T-V.p-K*V.p-1))          # rescaling residual (Stine, JASA 1987)
  for j in 1:nrep
    # Draw block of initial pre-sample values
    yr = zeros(K,T)                                       # bootstrap data
    pos = trunc(Integer, (rand()*(T-V.p+1)))+1            # position of initial draw
    yr[:,1:V.p] = V.Y'[pos:pos+V.p-1,:]'                   # drawing pre-sample obs
    # Draw innovations
    ur = zeros(K,T)                             # bootstrap innovations
    index = trunc(Integer, (rand(T-V.p)*(T-V.p))+1)         # index for innovation draws
    ur[:, V.p+1:T] = u[:, index]                  # drawing innovations
    # Recursively construct sample
    for i in (V.p+1):T
      V.i == true ? (yr[:,i] = V.β[:,1] + ur[:,i]) : (yr[:,i] = ur[:,i])
      for jj in 1:V.p
        V.i == true ? (yr[:,i] += V.β[:,(jj-1)*K + 2:jj*K+1]*yr[:,i-jj]) : (yr[:,i] += V.β[:,(jj-1)*K + 1:jj*K]*yr[:,i-jj])
      end
    end
    yr = yr'
    for i in 1:K                                   # demean yr bootstrap data
        yr[:,i] -= mean(yr[:,i])
    end
    #pr = V.p # also using lag length selection
    Vr = VAR(yr,V.p,true)
    # Bias correction: if the largest root of the companion matrix
    # is less than 1, do BIAS correction
    i = true && (any(abs(eigvals(get_VAR1_rep(Vr))).>=1) || bias_correction(Vr))
    #if ~ any(abs(eigvals(Ar)).>=1)
    #  Ar = asybc(Ar,SIGMAr,T,K,Vr.p) #############################CHECK
    #end
    #Ar = real(Ar)
    IRFr = irf(Vr, H, true) ##############PROBLEM HERE
    IRFrmat[j,:] = vec(IRFr')'
  end                              # end bootstrap
  # Calculate 95 perccent interval endpoints
  CILv = zeros(1,size(IRFrmat,2))
  CIHv = zeros(1,size(IRFrmat,2))
  for i = 1:size(IRFrmat,2)
    CILv[:,i] = quantile(vec(IRFrmat[:,i]),0.025)
    CIHv[:,i] = quantile(vec(IRFrmat[:,i]),0.975)
  end
  CIL  = reshape(CILv',H+1,K^2)';
  CIH  = reshape(CIHv',H+1,K^2)';
  return CIL, CIH
end


text:
module Adapt

export adapt

"""
    adapt(to, x)

Adapt a value `x` according to `to`. If no specific adaptions have been registered for `to`,
this call will be a no-op.

To alter the behavior, implement methods of `adapt_structure` and `adapt_storage` to
respectively define how to adapt structures, and the leaves of those structures.

For example, defining an adaptor for an environment where we can't have integers, and adding
a method to `adapt_storage` to appropriately convert those to floating-point numbers:

    julia> struct IntegerLessAdaptor end

    julia> Adapt.adapt_storage(::IntegerLessAdaptor, x::Int64) = Float64(x)

    julia> adapt(IntegerLessAdaptor(), 42)
    42.0

This will automatically work on known types too:

    julia> adapt(IntegerLessAdaptor(), tuple(1,2,3))
    (1.0, 2.0, 3.0)

If we want this to work with custom structures, we need to extend `adapt_structure`:

    julia> struct MyStructure
             x
           end

    julia> Adapt.adapt_structure(to, obj::MyStructure) = MyStructure(adapt(to, obj.x))

    julia> adapt(IntegerLessAdaptor(), MyStructure(42))
    MyStructure(42.0)
"""
adapt(to, x) = adapt_structure(to, x)

adapt_structure(to, x) = adapt_storage(to, x)
adapt_storage(to, x) = x

include("base.jl")
include("wrappers.jl")

end # module


text:
using MinkReduction
using Test

@testset "MinkReduction.jl" begin
    @test m = DeviousMat(10); a = m[1,:]; b = m[2,:]; c = m[3,:]; d,e,f = minkReduce(a,b,c); hcat(d,e,f) == [0.0 -1.0 0.0; 0.0 0.0 -1.0; 1.0 0.0 0.0]

end

max_stars_repo_name: glwhart/JuSym

max_issues_repo_name: glwhart/JuSym

max_forks_repo_name: glwhart/JuSym


text:
function limit_features(doc_terms::DocumentTermMatrix,
                        high::Int,
                        low::Int)
    doc_freqs = vec(sum(doc_terms.dtm, dims=2))

    # build mask to restrict terms
    mask = trues(length(doc_freqs))
    if high < 1
        mask .&= (doc_freqs .<= high)
    end
    if low > 0
        mask .&= (doc_freqs .>= low)
    end

    new_terms = doc_terms.terms[mask]

    return (doc_terms.dtm[mask, :], new_terms)
end

## Helper functions to build Corpus ##
_convert_bag_of_words(X::Dict{<:NGram, <:Integer}) =
    Dict(join(k, " ") => v for (k, v) in X)

build_corpus(X::Vector{<:Dict{<:NGram, <:Integer}}) =
    build_corpus(_convert_bag_of_words.(X))
build_corpus(X::Vector{<:Dict{S, <:Integer}}) where {S <: AbstractString} =
    Corpus(NGramDocument.(X))
build_corpus(X) = Corpus(TokenDocument.(X))

# based on https://github.com/zgornel/StringAnalysis.jl/blob/master/src/dtm.jl
# and https://github.com/JuliaText/TextAnalysis.jl/blob/master/src/dtm.jl
build_dtm(docs::Corpus) = build_dtm(docs, sort(collect(keys(lexicon(docs)))))
function build_dtm(docs::Corpus, terms::Vector{T}) where {T}
    # we are flipping the orientation of this matrix
    # so we get the `columnindices` from the TextAnalysis API
    row_indices = TextAnalysis.columnindices(terms)

    m = length(terms) # terms are rows
    n = length(docs)  # docs are columns

    rows = Vector{Int}(undef, 0) # terms
    columns = Vector{Int}(undef, 0) # docs
    values = Vector{Int}(undef, 0)
    for i in eachindex(docs.documents)
        doc = docs.documents[i]
        ngs = ngrams(doc)
        for ngram in keys(ngs)
            j = get(row_indices, ngram, 0)
            v = ngs[ngram]
            if j != 0
                push!(columns, i)
                push!(rows, j)
                push!(values, v)
            end
        end
    end
    if length(rows) > 0
        doc_term_mat = sparse(rows, columns, values, m, n)
    else
        doc_term_mat = spzeros(Int, m, n)
    end
    DocumentTermMatrix(doc_term_mat, terms, row_indices)
end

## General method to calculate IDF vector ##
function compute_idf(smooth_idf::Bool, doc_term_mat::SparseMatrixCSC{T}) where {T <: Real}
    n = size(doc_term_mat, 2) # docs are columns

    smooth_idf = Int(smooth_idf)
    documents_containing_term = vec(sum(doc_term_mat .> 0, dims=2)) .+ smooth_idf
    idf = log.((n + smooth_idf) ./ documents_containing_term) .+ 1

    return idf
end


text:
"""
    AAindex

An AAindex parser.

Provides simple functionality to read AAindex database files. See [Kawashima
S and Kanehisa M. (2000)](https://dx.doi.org/10.1093%2Fnar%2F28.1.374)
"""
module AAindex

using DataDeps
using JLD2
using StaticArrays

export
# Types
AbstractAAIndex, Metadata, Index, AMatrix,

# Functions
parse, aaindex_by_id, search

abstract type AbstractAAIndex end

include("./init.jl")
include("./data.jl")
include("./index.jl")
include("./parse.jl")
include("./search.jl")

end # module


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 71903690-5f0e-11eb-2435-a56a85d76fe0
using Plots

# ╔═╡ 7eef98bc-5f0e-11eb-0dc2-43a0e5715aa4
using PlutoUI

# ╔═╡ d5cced56-5f0c-11eb-36ce-37a624f02490
md"# Heat kernel

The heat kernel is defined by

$K(t,x,y) = \frac{1}{(4 ~ \pi ~ t)^{d/2}} ~ \exp\left(\frac{-\lvert x - y \rvert^2}{4 t}\right)$

with $d = 1,2,\cdots$ as the number of dimensions, time $t \geq 0$ and space $x \in \mathbb{R}^{d}$.
"

# ╔═╡ 4131c77a-5f0e-11eb-1aae-41015705d66e
function kernel(t, x, y; dim = 1)
    d = dim

    arg = -1 * (abs(x - y)^2)/(4*t)
    den = (4*pi*t)^(d/2)

    sol = (1/den) * exp(arg)

    return sol
end

# ╔═╡ 4118304e-5f0e-11eb-3844-5159c87d48b7
tspan = 0.001 : 0.001 : 1.0 # time

# ╔═╡ 40fa23ec-5f0e-11eb-2b02-73a75b7a6049
xspan = -3.0 : 0.1 : 3.0 # space

# ╔═╡ 401f2062-5f0e-11eb-2517-5900dd4b84ba
y = 0.0 # Offset

# ╔═╡ 679b9134-5f0e-11eb-159c-258752c5fcef
sol = zeros(length(xspan), length(tspan));

# ╔═╡ 6f5614c6-5f0e-11eb-3412-5381fbf426a0
for i = 1 : length(tspan)
    for k = 1 : length(xspan)
        sol[k,i] = kernel(tspan[i], xspan[k], y)
    end
end

# ╔═╡ ab00a9dc-5f0e-11eb-36e1-fd6e0d1acb68
@bind pos Slider(1:10:1000, default=1)

# ╔═╡ be7a3cf8-5f0e-11eb-3bb6-0d132e1b8e42
scatter(xspan, sol[:,pos], legend=false)

# ╔═╡ Cell order:
# ╟─d5cced56-5f0c-11eb-36ce-37a624f02490
# ╠═4131c77a-5f0e-11eb-1aae-41015705d66e
# ╠═4118304e-5f0e-11eb-3844-5159c87d48b7
# ╠═40fa23ec-5f0e-11eb-2b02-73a75b7a6049
# ╠═401f2062-5f0e-11eb-2517-5900dd4b84ba
# ╠═679b9134-5f0e-11eb-159c-258752c5fcef
# ╠═6f5614c6-5f0e-11eb-3412-5381fbf426a0
# ╠═71903690-5f0e-11eb-2435-a56a85d76fe0
# ╠═7eef98bc-5f0e-11eb-0dc2-43a0e5715aa4
# ╠═ab00a9dc-5f0e-11eb-36e1-fd6e0d1acb68
# ╠═be7a3cf8-5f0e-11eb-3bb6-0d132e1b8e42


text:
# Particle filtering an AR(1) model with measurement noise

μ0(θ) = MvNormal([θ[1]/(1-tanh(θ[2]))],fill(θ[3]^2/(1-tanh(θ[2])^2),1,1))
f(x,u,t,θ) = MvNormal(θ[1] .+ tanh(θ[2])*x,fill(θ[3]^2,1,1))
g(x,u,t,θ) = MvNormal(x,fill(θ[4]^2,1,1))

function q(x,y,u,t,θ)
    δ = θ[1]
    φ = tanh(θ[2])
    σ = θ[3]
    σε = θ[4]
    ν = δ .+ φ * x
    MvNormal((ν*σε^2 .+ y*σ^2)/(σε^2 + σ^2),fill((σ^2*σε^2)/(σ^2 + σε^2),1,1))
end

function ξ(x,y,u,t,θ)
    δ = θ[1]
    φ = tanh(θ[2])
    σ = θ[3]
    σε = θ[4]
    ν = δ .+ φ * x
    logpdf(MvNormal(ν,sqrt(σ^2 + σε^2)),y)
end

mar = ProposalStateSpaceModel(μ0,f,g,q,ξ)

δ = 1.0
φ = 0.95
ση = 1.0
σε = 1.0

θ0 = [δ;atanh(φ);ση;σε]

β  = [0.0;1.0]

T = 1000
N = 100

Random.seed!(1234)

X,Y = Filtering.simulate(mar,θ0,T)
X1 = X[1][1]

Random.seed!(1234)

@testset "Simulating from model" begin
    X,Y = Filtering.simulate(mar,θ0,T)

    # Test to ensure that the RNG seed is working
    @test X[1][1] == X1

    @test length(X) == T+1
    @test length(Y) == T

    μ = δ/(1-φ)
    σ = ση/sqrt(1-φ^2)

    # These are rough approximations, just to make sure we are
    # getting a time series that looks right. Perhaps they should
    # be replaced with proper hypothesis tests.

    @test abs(mean(hcat(X...)[1,:]) - μ) < 3.0
    @test abs(std(hcat(X...)[1,:]) - σ) < 2.0
    @test abs(autocor(hcat(X...)[1,:],[1])[1] - φ) < 0.2
end

@testset "Running Kalman filter" begin
    X,Y = Filtering.simulate(mar,θ0,T)

    @test X[1][1] == X1

    arh = ARH(β)
    kf = kalman_filter(arh,θ0,Y,map(x->[x],ones(T)))
    ks = kalman_smoother(arh,θ0,kf)
end

@testset "Running particle filter" begin
    X,Y = Filtering.simulate(mar,θ0,T)

    @test X[1][1] == X1

    p1 = Filtering.particle_filter(mar,θ0,Y,N)
    #p2 = Filtering.filtering_filter(mar,Y,θ0,N)
end

#=
@testset "Sampling from smoothing distribution" begin
    B = 1

    Xs = [Filtering.generate_realization(p2,mar,θ0) for i in 1:B]
end
=#


text:
mutable struct Configuration{T}
    solution::Bool
    value::Float64
    values::Dictionary{Int, T}
end

is_solution(c) = c.solution

get_value(c) = c.value
get_error(c) = is_solution(c) ? 0.0 : get_value(c)
get_values(c) = c.values
get_value(c, x) = get_values(c)[x]

set_value!(c, val) = c.value = val
# set_value!(c, x, val) = get_values(c)[x] = val
set_values!(c, values) = c.values = values
set_sat!(c, b) = c.solution = b

compute_cost(m, config::Configuration) = compute_cost(m, get_values(config))
compute_cost!(m, config::Configuration) = set_value!(config, compute_cost(m, config))

function Configuration(m::_Model, X)
    values = draw(m)
    val = compute_costs(m, values, X)
    sol = val ≈ 0.0
    opt = sol && !is_sat(m)
    return Configuration(sol, opt ? compute_objective(m, values) : val, values)
end


text:
module Geotherm

include("Geometry.jl")
include("Interpolate.jl")
include("Integrate.jl")

using BisectPy: bisect_right
import DimensionalData  # Need to load this to run code
using DimensionalData: AbstractDimMatrix, dims, @dim

using Geotherm.Geometry: Point2D, Point3D, Rectangle
using Geotherm.Interpolate: bilinear_interpolate
using Geotherm.Integrate: runge_kutta_iter

export Temp, Press, generate_trace

# This is different from `BisectPy.find_le`, which returns the value!
function find_le(a, x)
    i = bisect_right(a, x)
    return i > 1 ? i - 1 : nothing
end

# This is a helper function and should not be exported.
function inject_find_lower_bound(ps, ts, geothermal_gradient)
    return function (p, t)
        i, j = find_le(ps, p), find_le(ts, t)
        k, l = i + 1, j + 1
        interpolated_function = bilinear_interpolate(
            Point3D(ps[i], ts[j], geothermal_gradient[j, i]),  # Note the order of indices!
            Point3D(ps[i], ts[l], geothermal_gradient[l, i]),  # Note the order of indices!
            Point3D(ps[k], ts[j], geothermal_gradient[j, k]),  # Note the order of indices!
            Point3D(ps[k], ts[l], geothermal_gradient[l, k]),   # Note the order of indices!
        )
        return interpolated_function(p, t)
    end
end

@dim Temp "Temperature"
@dim Press "Pressure"

"""
    generate_trace(geothermal_gradient::AbstractDimMatrix, p0::Point2D, h=0.01, n=1000)

Generate a trace of `Point`s specifiying pressure and temperature.
"""
function generate_trace(
    geothermal_gradient::AbstractDimMatrix{T,<:Tuple{Temp,Press}},
    p0::Point2D,
    h = 0.01,
    n = 1000,
) where {T}
    ps, ts = dims(geothermal_gradient, (Press, Temp))
    trace = Point2D[p0]
    f = inject_find_lower_bound(ps, ts, geothermal_gradient)
    for i in 1:(n-1)
        push!(trace, runge_kutta_iter(trace[i], f, h))
    end
    return trace
end
generate_trace(
    geothermal_gradient::AbstractDimMatrix{T,<:Tuple{Press,Temp}},
    p0::Point2D,
    h = 0.01,
    n = 1000,
) where {T} = generate_trace(geothermal_gradient', p0, h, n)

end


text:
module FStrings
using Printf

export @f_str

# 0) Base version - has problems with `:` inside the Julia code section
# regex = r"(?<!\\)\{(.*?)(?::+?(.*?))?(?<!\\)\}"
# 1) Version that would require the Julia code section to escape `:` to `\:` inside {}
# regex = r"(?<!\\)\{(.*?)(?:(?<!\\):+?(.*?))?(?<!\\)\}"
# 2) Version that can keep normal `:` inside the Julia code, if there is one last `:` present as well
# regex = r"(?<!\\)\{(.*?)(?::+?([^:\n\r]*?))?(?<!\\)\}"
# 3) Version that can handle all cases by automatically matching only allowed format specifiers
# regex = r"(?<!\\)\{(.*?)(?::([^:\n\r^\(^\)^\[^\]]*?))?(?<!\\)\}"
# ...slightly improved matching:
regex = r"(?<!\\)\{(.*?)(?::([^:\n\r\(\)\[\]]*?[^:\n\r\(\)\[\]\d]))?(?<!\\)\}"

# Taken from: https://github.com/JuliaLang/julia/blob/3608c84e6093594fe86923339fc315231492484c/stdlib/Printf/src/Printf.jl
# printf specifiers:
#   %                       # start
#   (\d+\$)?                # arg (not supported)
#   [\-\+#0' ]*             # flags
#   (\d+)?                  # width
#   (\.\d*)?                # precision
#   (h|hh|l|ll|L|j|t|z|q)?  # modifier (ignored)
#   [diouxXeEfFgGaAcCsSp%]  # conversion
regex_valid_format_specifiers = r"^(\d+\$)?[\-\+#0' ]*(\d+)?(\.\d*)?(h|hh|l|ll|L|j|t|z|q)?[diouxXeEfFgGaAcCsSp%]$"

"""
    @f_str(string::AbstractString)

    Loose implementation of Python style `fstring` literal string interpolation
    based on `Printf.@sprintf`.

# Examples
```julia-repl
julia> using FStrings
julia> f"π ≈ {π:.2f}"
"π ≈ 3.14"
```

# Format Specifiers
Please refer to `Printf.@sprintf` for further details on the available
format specifiers. Also refer to the principle syntax of `fstring` via PEP 498.

# References
- [`Printf.@sprintf`]: https://docs.julialang.org/en/v1/stdlib/Printf/#Printf.@sprintf
- [PEP 498]: https://www.python.org/dev/peps/pep-0498/
"""
macro f_str(string::AbstractString)
    new_string = ""
    args = Any[]
    last_pos = 1
    for (idx, match_) = enumerate(eachmatch(regex, string))
        ext_string = string[last_pos:match_.offset-1]
        ext_string = replace(ext_string, "%"=>"%%")
        new_string *= ext_string
        last_pos = match_.offset + length(match_.match)

        expr_txt = match_.captures[1]
        if length(expr_txt) == 0
            @warn "Format string field without content present - skipping."
            continue
        end
        expr = Meta.parse(expr_txt)  # Needs to be valid Julia code

        frmt_txt = match_.captures[2]
        if frmt_txt === nothing
            frmt = "s"
        elseif frmt_txt == ""
            frmt = "s"
        #elseif frmt_txt[end] == 'b'
        # Potential extension for a binary format. Could use `bitstring`
        # to convert result of expr. Needs parsing of the leading `0`
        # and bitwidth count. I.e. f"{123:b}", f"{123:8b}", f"{123:08b}"
        else
            frmt = frmt_txt
        end
        # *) Handle invalid format strings already here
        if match(regex_valid_format_specifiers, frmt) === nothing
            exc_txt = "invalid printf format string: \"$frmt\""
            exc = :(throw(ArgumentError($exc_txt)))
            return exc
        end
        new_string *= "%$frmt"

        push!(args, expr)
    end
    if last_pos <= length(string)
        ext_string = string[last_pos:end]
        ext_string = replace(ext_string, "%"=>"%%")
        new_string *= ext_string
    end

    # `macro f_str` received escape characters (e.g. `"\n"`) with the backslash
    # escaped. Need to undo this before passing it on.
    new_string = unescape_string(new_string)

    # In case `new_string` is empty, the `@sprintf` throws `LoadError`.
    # Workaround for now:
    if length(new_string) == 0
        return new_string
    end

    # TODO: Didn't figure out how to escape this, such that
    # `ArgumentErrors` are thrown in the outer scope _without_
    # producing also a `LoadError` when there is a wrong format specifier.
    # If this would work, *) can be removed.
    # Keep in mind though, that `@sprintf` should be evaluated only once
    # in the outside code for performance reasons.
    ex = :(@sprintf($new_string, $(esc.(args)...)))
    return ex
end

end  # module


text:
using Documenter, BioStructures

makedocs(
    format = :html,
    sitename = "BioStructures.jl",
    pages = [
        "Home"         => "index.md",
        "Documentation"=> "documentation.md",
        "Contributing" => "contributing.md"
    ],
    authors = "Joe G Greener, The BioJulia Organisation and other contributors."
)

deploydocs(
    repo = "github.com/BioJulia/BioStructures.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
@enum WordClass begin
    Bad
    Questionable
    Overridden
end
struct BadWord
    class::WordClass
    word::String
end

const WordClassStrings = Dict(
    Bad => "Bad",
    Questionable => "Questionable",
    Overridden => "Overridden",
)

function Base.show(io::IO, w::BadWord)
    print(io, WordClassStrings[w.class], "(", w.word, ")")
end