text:
module GR3

import GR
import Libdl

macro triplet(t)
    :( Tuple{$t, $t, $t} )
end

macro ArrayToVector(ctype, data)
    return :( convert(Vector{$(esc(ctype))}, vec($(esc(data)))) )
end

mutable struct PNG s::Array{UInt8} end
mutable struct HTML s::AbstractString end

Base.show(io::IO, ::MIME"image/png", x::PNG) = write(io, x.s)
Base.show(io::IO, ::MIME"text/html", x::HTML) = print(io, x.s)

function _readfile(path)
    data = Array(UInt8, filesize(path))
    s = open(path, "r")
    read!(s, data)
end

mutable struct GR3Exception <: Exception
    msg::AbstractString
end
Base.showerror(io::IO, e::GR3Exception) = print(io, e.msg);

const msgs = [ "none", "invalid value", "invalid attribute", "init failed",
               "OpenGL error", "out of memory", "not initialized",
               "camera not initialized", "unknown file extension",
               "cannot open file", "export failed" ]

function _check_error()
  line = Cint[0]
  file = Ptr{UInt8}[0]
  error_code = ccall(GR.libGR3_ptr(:gr3_geterror), Int32, (Int32, Ptr{Cint}, Ptr{Ptr{UInt8}}), 1, line, file)
  if (error_code != 0)
    line = line[1]
    file = unsafe_string(file[1])
    if 0 <= error_code < length(msgs)
      msg = msgs[error_code + 1]
    else
      msg = "unknown error"
    end
    message = string("GR3 error (", file, ", l. ", line, "): ", msg)
    throw(GR3Exception(message))
  end
end

function init(attrib_list)
  ccall(GR.libGR3_ptr(:gr3_init), Int32, (Ptr{Int}, ), attrib_list)
  _check_error()
end
export init

function free(pointer)
  ccall(GR.libGR3_ptr(:gr3_free), Nothing, (Ptr{Nothing}, ), pointer)
  _check_error()
end
export free

function terminate()
  ccall(GR.libGR3_ptr(:gr3_terminate), Nothing, ())
  _check_error()
end
export terminate

function useframebuffer(framebuffer)
  ccall(GR.libGR3_ptr(:gr3_useframebuffer), Nothing, (UInt32, ), framebuffer)
  _check_error()
end
export useframebuffer

function usecurrentframebuffer()
  ccall(GR.libGR3_ptr(:gr3_usecurrentframebuffer), Nothing, ())
  _check_error()
end
export usecurrentframebuffer

function getimage(width, height, use_alpha=true)
  bpp = use_alpha ? 4 : 3
  bitmap = zeros(UInt8, width * height * bpp)
  ccall(GR.libGR3_ptr(:gr3_getimage),
        Int32,
        (Int32, Int32, Int32, Ptr{UInt8}),
        width, height, use_alpha, bitmap)
  _check_error()
  return bitmap
end
export getimage

function save(filename, width, height)
  ccall(GR.libGR3_ptr(:gr3_export),
        Int32,
        (Ptr{Cchar}, Int32, Int32),
        filename, width, height)
  _check_error()
  ext = splitext(filename)[end:end][1]
  if ext == ".png"
    content = PNG(_readfile(filename))
  elseif ext == ".html"
    content = HTML("<iframe src=\"files/$filename\" width=$width height=$height></iframe>")
  else
    content = nothing
  end
  return content
end
export save

function getrenderpathstring()
  val = ccall(GR.libGR3_ptr(:gr3_getrenderpathstring),
              Ptr{UInt8}, (), )
  _check_error()
  unsafe_string(val)
end
export getrenderpathstring

function drawimage(xmin, xmax, ymin, ymax, pixel_width, pixel_height, window)
  ccall(GR.libGR3_ptr(:gr3_drawimage),
        Int32,
        (Float32, Float32, Float32, Float32, Int32, Int32, Int32),
        xmin, xmax, ymin, ymax, pixel_width, pixel_height, window)
  _check_error()
end
export drawimage

function createmesh(n, vertices, normals, colors)
  mesh = Cint[0]
  _vertices = [ Float32(x) for x in vertices ]
  _normals = [ Float32(x) for x in normals ]
  _colors = [ Float32(x) for x in colors ]
  ccall(GR.libGR3_ptr(:gr3_createmesh),
        Int32,
        (Ptr{Cint}, Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}),
        mesh, n, @ArrayToVector(Float32, _vertices), @ArrayToVector(Float32, _normals), @ArrayToVector(Float32, _colors))
  _check_error()
  return mesh[1]
end
export createmesh

function createindexedmesh(num_vertices, vertices, normals, colors, num_indices, indices)
  mesh = Cint[0]
  _vertices = [ Float32(x) for x in vertices ]
  _normals = [ Float32(x) for x in normals ]
  _colors = [ Float32(x) for x in colors ]
  _indices = [ Float32(x) for x in indices ]
  ccall(GR.libGR3_ptr(:gr3_createindexedmesh),
        Int32,
        (Ptr{Cint}, Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Int32, Ptr{Int32}),
        mesh, num_vertices, @ArrayToVector(Float32, _vertices), @ArrayToVector(Float32, _normals), @ArrayToVector(Float32, _colors), num_indices, @ArrayToVector(Float32, _indices))
  _check_error()
  return mesh[1]
end
export createindexedmesh

function drawmesh(mesh::Int32, n, positions::@triplet(Real), directions::@triplet(Real), ups::@triplet(Real), colors::@triplet(Real), scales::@triplet(Real))
  _positions = [ Float32(x) for x in positions ]
  _directions = [ Float32(x) for x in directions ]
  _ups = [ Float32(x) for x in ups ]
  _colors = [ Float32(x) for x in colors ]
  _scales = [ Float32(x) for x in scales ]
  ccall(GR.libGR3_ptr(:gr3_drawmesh),
        Nothing,
        (Int32, Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}),
        mesh, n, @ArrayToVector(Float32, _positions), @ArrayToVector(Float32, _directions), @ArrayToVector(Float32, _ups), @ArrayToVector(Float32, _colors), @ArrayToVector(Float32, _scales))
  _check_error()
end
export drawmesh

function createheightmapmesh(heightmap, num_columns, num_rows)
  if num_columns * num_rows == length(heightmap)
    if ndims(heightmap) == 2
      heightmap = reshape(heightmap, num_columns * num_rows)
    end
    ccall(GR.libGR3_ptr(:gr3_createheightmapmesh),
          Nothing,
          (Ptr{Float32}, Int32, Int32),
          @ArrayToVector(Float32, heightmap), num_columns, num_rows)
    _check_error()
  else
    println("Array has incorrect length or dimension.")
  end
end
export createheightmapmesh

function drawheightmap(heightmap, num_columns, num_rows, positions, scales)
  if num_columns * num_rows == length(heightmap)
    if ndims(heightmap) == 2
      heightmap = reshape(heightmap, num_columns * num_rows)
    end
    _positions = [ Float32(x) for x in positions ]
    _scales = [ Float32(x) for x in scales ]
    ccall(GR.libGR3_ptr(:gr3_drawheightmap),
          Nothing,
          (Ptr{Float32}, Int32, Int32, Ptr{Float32}, Ptr{Float32}),
          @ArrayToVector(Float32, heightmap), num_columns, num_rows, @ArrayToVector(Float32, _positions), @ArrayToVector(Float32, _scales))
    _check_error()
  else
    println("Array has incorrect length or dimension.")
  end
end
export drawheightmap

function deletemesh(mesh)
  ccall(GR.libGR3_ptr(:gr3_deletemesh), Nothing, (Int32, ), mesh)
  _check_error()
end
export deletemesh

function setquality(quality)
  ccall(GR.libGR3_ptr(:gr3_setquality), Nothing, (Int32, ), quality)
  _check_error()
end
export setquality

function clear()
  ccall(GR.libGR3_ptr(:gr3_clear), Nothing, ())
  _check_error()
end
export clear

function cameralookat(camera_x, camera_y, camera_z,
                      center_x, center_y, center_z,
                      up_x, up_y, up_z)
  ccall(GR.libGR3_ptr(:gr3_cameralookat),
        Nothing,
        (Float32, Float32, Float32, Float32, Float32, Float32, Float32, Float32, Float32),
        camera_x, camera_y, camera_z, center_x, center_y, center_z, up_x, up_y, up_z)
  _check_error()
end
export cameralookat

function setcameraprojectionparameters(vertical_field_of_view, zNear, zFar)
  ccall(GR.libGR3_ptr(:gr3_setcameraprojectionparameters),
        Nothing,
        (Float32, Float32, Float32),
        vertical_field_of_view, zNear, zFar)
  _check_error()
end
export setcameraprojectionparameters

function setlightdirection(x, y, z)
  ccall(GR.libGR3_ptr(:gr3_setlightdirection),
        Nothing,
        (Float32, Float32, Float32),
        x, y, z)
  _check_error()
end
export setlightdirection

function drawcylindermesh(n, positions, directions, colors, radii, lengths)
  _positions = [ Float32(x) for x in positions ]
  _directions = [ Float32(x) for x in directions ]
  _colors = [ Float32(x) for x in colors ]
  _radii = [ Float32(x) for x in radii ]
  _lengths = [ Float32(x) for x in lengths ]
  ccall(GR.libGR3_ptr(:gr3_drawcylindermesh),
        Nothing,
        (Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}),
        n, @ArrayToVector(Float32, _positions), @ArrayToVector(Float32, _directions), @ArrayToVector(Float32, _colors), @ArrayToVector(Float32, _radii), @ArrayToVector(Float32, _lengths))
  _check_error()
end
export drawcylindermesh

function drawconemesh(n, positions, directions, colors, radii, lengths)
  _positions = [ Float32(x) for x in positions ]
  _directions = [ Float32(x) for x in directions ]
  _colors = [ Float32(x) for x in colors ]
  _radii = [ Float32(x) for x in radii ]
  _lengths = [ Float32(x) for x in lengths ]
  ccall(GR.libGR3_ptr(:gr3_drawconemesh),
        Nothing,
        (Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}),
        n, @ArrayToVector(Float32, _positions), @ArrayToVector(Float32, _directions), @ArrayToVector(Float32, _colors), @ArrayToVector(Float32, _radii), @ArrayToVector(Float32, _lengths))
  _check_error()
end
export drawconemesh

function drawspheremesh(n, positions, colors, radii)
  _positions = [ Float32(x) for x in positions ]
  _colors = [ Float32(x) for x in colors ]
  _radii = [ Float32(x) for x in radii ]
  ccall(GR.libGR3_ptr(:gr3_drawspheremesh),
        Nothing,
        (Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}),
        n, @ArrayToVector(Float32, _positions), @ArrayToVector(Float32, _colors), @ArrayToVector(Float32, _radii))
  _check_error()
end
export drawspheremesh

function drawcubemesh(n, positions, directions, ups, colors, scales)
  _positions = [ Float32(x) for x in positions ]
  _directions = [ Float32(x) for x in directions ]
  _ups = [ Float32(x) for x in ups ]
  _colors = [ Float32(x) for x in colors ]
  _scales = [ Float32(x) for x in scales ]
  ccall(GR.libGR3_ptr(:gr3_drawcubemesh),
        Nothing,
        (Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}),
        n, @ArrayToVector(Float32, _positions), @ArrayToVector(Float32, _directions), @ArrayToVector(Float32, _ups), @ArrayToVector(Float32, _colors), @ArrayToVector(Float32, _scales))
  _check_error()
end
export drawcubemesh

function setbackgroundcolor(red, green, blue, alpha)
  ccall(GR.libGR3_ptr(:gr3_setbackgroundcolor),
        Nothing,
        (Float32, Float32, Float32, Float32),
        red, green, blue, alpha)
  _check_error()
end
export setbackgroundcolor

function createisosurfacemesh(grid::Array{UInt16,3}, step::@triplet(Float64), offset::@triplet(Float64), isolevel::Int64)
  mesh = Cint[0]
  dim_x, dim_y, dim_z = size(grid)
  data = reshape(grid, dim_x * dim_y * dim_z)
  stride_x, stride_y, stride_z = strides(grid)
  step_x, step_y, step_z = [ float(x) for x in step ]
  offset_x, offset_y, offset_z = [ float(x) for x in offset ]
  err = ccall(GR.libGR3_ptr(:gr3_createisosurfacemesh),
              Int32,
              (Ptr{Cint}, Ptr{UInt16}, UInt16, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, Float64, Float64, Float64, Float64, Float64, Float64),
              mesh, @ArrayToVector(UInt16, data), UInt16(isolevel), dim_x, dim_y, dim_z, stride_x, stride_y, stride_z, step_x, step_y, step_z, offset_x, offset_y, offset_z)
  _check_error()
  return mesh[1]
end
export createisosurfacemesh

function surface(px, py, pz, option::Int)
  nx = length(px)
  ny = length(py)
  nz = length(pz)
  if ndims(pz) == 1
    out_of_bounds = nz != nx * ny
  elseif ndims(px) == ndims(py) == ndims(pz) == 2
    nx, ny = size(pz)
    out_of_bounds = size(px)[1] != ny || size(px)[2] != nx ||
                    size(py)[1] != ny || size(py)[2] != nx
  elseif ndims(pz) == 2
    out_of_bounds = size(pz)[1] != nx || size(pz)[2] != ny
  else
    out_of_bounds = true
  end
  if !out_of_bounds
    if option != GR.OPTION_3D_MESH
      if ndims(pz) == 2
        pz = reshape(pz, nx * ny)
      end
      _px = [ Float32(x) for x in px ]
      _py = [ Float32(y) for y in py ]
      _pz = [ Float32(z) for z in pz ]
    else
      _px = px'
      _py = py'
      _pz = pz
    end
    ccall(GR.libGR3_ptr(:gr3_surface),
          Nothing,
          (Int32, Int32, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Int32),
          nx, ny, @ArrayToVector(Float32, _px), @ArrayToVector(Float32, _py), @ArrayToVector(Float32, _pz), option)
    _check_error()
  else
    println("Arrays have incorrect length or dimension.")
  end
end

function volume(data::Array{Float64,3}, algorithm::Int64)
  dmin = Cdouble[-1]
  dmax = Cdouble[-1]
  nx, ny, nz = size(data)
  data = reshape(data, nx * ny * nz)
  ccall(GR.libGR3_ptr(:gr_volume),
        Nothing,
        (Cint, Cint, Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}),
        nx, ny, nz, data, algorithm, dmin, dmax)
  return dmin[1], dmax[1]
end

function createslicemeshes(grid; x::Union{Real, Nothing}=nothing, y::Union{Real, Nothing}=nothing, z::Union{Real, Nothing}=nothing, step::Union{Tuple{Real, Real, Real}, Nothing}=nothing, offset::Union{Tuple{Real, Real, Real}, Nothing}=nothing)
    if x === nothing && y === nothing && z === nothing
        x = 0.5
        y = 0.5
        z = 0.5
    end
    if typeof(grid[1,1,1]) <: Unsigned
        input_max = typemax(typeof(grid[1,1,1]))
    elseif typeof(grid[1,1,1]) <: Real
        input_max = convert(typeof(grid[1,1,1]), 1.0)
        grid = min.(grid, input_max)
    else
        println("grid must be three dimensional array of Real numbers")
        return(nothing)
    end
    scaling_factor = typemax(UInt16) / input_max
    grid = Array{UInt16, 3}(floor.(grid * scaling_factor))
    nx, ny, nz = size(grid)
    if step === nothing && offset === nothing
        step = (2.0/(nx-1), 2.0/(ny-1), 2.0/(nz-1))
        offset = (-1.0, -1.0, -1.0)
    elseif offset === nothing
        offset = (-step[1] * (nx-1) / 2.0,
                  -step[2] * (ny-1) / 2.0,
                  -step[3] * (nz-1) / 2.0)
    elseif step === nothing
        step = (-offset[1] * 2.0 / (nx-1),
                -offset[2] * 2.0 / (ny-1),
                -offset[3] * 2.0 / (nz-1))
    end
    stride_x, stride_y, stride_z = 1, nx, nx*ny
    dim_x, dim_y, dim_z = convert(Tuple{UInt32, UInt32, UInt32}, size(grid))
    step_x, step_y, step_z = convert(Tuple{Float64, Float64, Float64}, step)
    offset_x, offset_y, offset_z = convert(Tuple{Float64, Float64, Float64}, offset)

    if x !== nothing
        x = convert(UInt32, floor(clamp(x, 0, 1) * nx))
        mesh = Cint[0]
        ccall(GR.libGR3_ptr(:gr3_createxslicemesh),
            Nothing,
            (Ptr{UInt32}, Ptr{UInt16}, UInt32,
            UInt32, UInt32, UInt32,
            UInt32, UInt32, UInt32,
            Float64, Float64, Float64,
            Float64, Float64, Float64),
            mesh, grid, x,
            dim_x, dim_y, dim_z,
            stride_x, stride_y, stride_z,
            step_x, step_y, step_z,
            offset_x, offset_y, offset_z)
        _check_error()
        _mesh_x = mesh[1]
    else
        _mesh_x = nothing
    end

    if y !== nothing
        y = convert(UInt32, floor(clamp(y, 0, 1) * ny))
        mesh = Cint[0]
        ccall(GR.libGR3_ptr(:gr3_createyslicemesh),
            Nothing,
            (Ptr{UInt32}, Ptr{UInt16}, UInt32,
            UInt32, UInt32, UInt32,
            UInt32, UInt32, UInt32,
            Float64, Float64, Float64,
            Float64, Float64, Float64),
            mesh, grid, y,
            dim_x, dim_y, dim_z,
            stride_x, stride_y, stride_z,
            step_x, step_y, step_z,
            offset_x, offset_y, offset_z)
        _check_error()
        _mesh_y = mesh[1]
    else
        _mesh_y = nothing
    end

    if z !== nothing
        z = convert(UInt32, floor(clamp(z, 0, 1) * nz))
        mesh = Cint[0]
        ccall(GR.libGR3_ptr(:gr3_createzslicemesh),
            Nothing,
            (Ptr{UInt32}, Ptr{UInt16}, UInt32,
            UInt32, UInt32, UInt32,
            UInt32, UInt32, UInt32,
            Float64, Float64, Float64,
            Float64, Float64, Float64),
            mesh, grid, z,
            dim_x, dim_y, dim_z,
            stride_x, stride_y, stride_z,
            step_x, step_y, step_z,
            offset_x, offset_y, offset_z)
        _check_error()
        _mesh_z = mesh[1]
    else
        _mesh_z = nothing
    end

    return(_mesh_x, _mesh_y, _mesh_z)
end
export createslicemeshes

function createxslicemesh(grid, x::Real=0.5; step::Union{Tuple{Real, Real, Real}, Nothing}=nothing, offset::Union{Tuple{Real, Real, Real}, Nothing}=nothing)
    return createslicemeshes(grid, x=x, step=step, offset=offset)[1]
end
export createxslicemesh

function createyslicemesh(grid, y::Real=0.5; step::Union{Tuple{Real, Real, Real}, Nothing}=nothing, offset::Union{Tuple{Real, Real, Real}, Nothing}=nothing)
    return createslicemeshes(grid, y=y, step=step, offset=offset)[2]
end
export createyslicemesh

function createzslicemesh(grid, z::Real=0.5; step::Union{Tuple{Real, Real, Real}, Nothing}=nothing, offset::Union{Tuple{Real, Real, Real}, Nothing}=nothing)
    return createslicemeshes(grid, z=z, step=step, offset=offset)[3]
end
export createzslicemesh

function drawxslicemesh(grid, x::Real=0.5; step::Union{Tuple{Real, Real, Real}, Nothing}=nothing, offset::Union{Tuple{Real, Real, Real}, Nothing}=nothing, position::Tuple{Real, Real, Real}=(0, 0, 0), direction::Tuple{Real, Real, Real}=(0, 0, 1), up::Tuple{Real, Real, Real}=(0, 1, 0), color::Tuple{Real, Real, Real}=(1, 1, 1), scale::Tuple{Real, Real, Real}=(1, 1, 1))
    mesh = createxslicemesh(grid, x, step=step, offset=offset)
    drawmesh(mesh, 1, position, direction, up, color, scale)
    deletemesh(mesh)
end
export drawxslicemesh

function drawyslicemesh(grid, y::Real=0.5; step::Union{Tuple{Real, Real, Real}, Nothing}=nothing, offset::Union{Tuple{Real, Real, Real}, Nothing}=nothing, position::Tuple{Real, Real, Real}=(0, 0, 0), direction::Tuple{Real, Real, Real}=(0, 0, 1), up::Tuple{Real, Real, Real}=(0, 1, 0), color::Tuple{Real, Real, Real}=(1, 1, 1), scale::Tuple{Real, Real, Real}=(1, 1, 1))
    mesh = createxslicemesh(grid, y, step=step, offset=offset)
    drawmesh(mesh, 1, position, direction, up, color, scale)
    deletemesh(mesh)
end
export drawyslicemesh

function drawzslicemesh(grid, z::Real=0.5; step::Union{Tuple{Real, Real, Real}, Nothing}=nothing, offset::Union{Tuple{Real, Real, Real}, Nothing}=nothing, position::Tuple{Real, Real, Real}=(0, 0, 0), direction::Tuple{Real, Real, Real}=(0, 0, 1), up::Tuple{Real, Real, Real}=(0, 1, 0), color::Tuple{Real, Real, Real}=(1, 1, 1), scale::Tuple{Real, Real, Real}=(1, 1, 1))
    mesh = createxslicemesh(grid, z, step=step, offset=offset)
    drawmesh(mesh, 1, position, direction, up, color, scale)
    deletemesh(mesh)
end
export drawzslicemesh

function drawslicemeshes(data; x::Union{Real, Nothing}=nothing, y::Union{Real, Nothing}=nothing, z::Union{Real, Nothing}=nothing, step::Union{Tuple{Real, Real, Real}, Nothing}=nothing, offset::Union{Tuple{Real, Real, Real}, Nothing}=nothing, position::Tuple{Real, Real, Real}=(0, 0, 0), direction::Tuple{Real, Real, Real}=(0, 0, 1), up::Tuple{Real, Real, Real}=(0, 1, 0), color::Tuple{Real, Real, Real}=(1, 1, 1), scale::Tuple{Real, Real, Real}=(1, 1, 1))
    meshes = createslicemeshes(data, x=x, y=y, z=z, step=step, offset=offset)
    for mesh in meshes
        if mesh !== nothing
            drawmesh(mesh, 1, position, direction, up, color, scale)
            deletemesh(mesh)
        end
    end
end
export drawslicemeshes

function setlightsources(num_lights, directions, colors)
  @assert length(directions) == length(colors)
  _directions = [ Float32(x) for x in directions ]
  _colors = [ Float32(x) for x in colors ]
  ccall(GR.libGR3_ptr(:gr3_setlightsources),
        Nothing,
        (Int32, Ptr{Float32}, Ptr{Float32}),
        num_lights, @ArrayToVector(Float32, _directions), @ArrayToVector(Float32, _colors))
end
export setlightsources

function getlightsources(num_lights)
  directions = Vector{Cfloat}(undef, num_lights * 3)
  colors = Vector{Cfloat}(undef, num_lights * 3)
  ccall(GR.libGR3_ptr(:gr3_getlightsources),
        Nothing,
        (Int32, Ptr{Cfloat}, Ptr{Cfloat}),
        num_lights, directions, colors)
  return directions, colors
end
export getlightsources

const IA_END_OF_LIST = 0
const IA_FRAMEBUFFER_WIDTH = 1
const IA_FRAMEBUFFER_HEIGHT = 2

const ERROR_NONE = 0
const ERROR_INVALID_VALUE = 1
const ERROR_INVALID_ATTRIBUTE = 2
const ERROR_INIT_FAILED = 3
const ERROR_OPENGL_ERR = 4
const ERROR_OUT_OF_MEM = 5
const ERROR_NOT_INITIALIZED = 6
const ERROR_CAMERA_NOT_INITIALIZED = 7
const ERROR_UNKNOWN_FILE_EXTENSION = 8
const ERROR_CANNOT_OPEN_FILE = 9
const ERROR_EXPORT = 10

const QUALITY_OPENGL_NO_SSAA  = 0
const QUALITY_OPENGL_2X_SSAA  = 2
const QUALITY_OPENGL_4X_SSAA  = 4
const QUALITY_OPENGL_8X_SSAA  = 8
const QUALITY_OPENGL_16X_SSAA = 16
const QUALITY_POVRAY_NO_SSAA  = 0+1
const QUALITY_POVRAY_2X_SSAA  = 2+1
const QUALITY_POVRAY_4X_SSAA  = 4+1
const QUALITY_POVRAY_8X_SSAA  = 8+1
const QUALITY_POVRAY_16X_SSAA = 16+1

const DRAWABLE_OPENGL = 1
const DRAWABLE_GKS = 2

const SURFACE_DEFAULT     =  0
const SURFACE_NORMALS     =  1
const SURFACE_FLAT        =  2
const SURFACE_GRTRANSFORM =  4
const SURFACE_GRCOLOR     =  8
const SURFACE_GRZSHADED   = 16

end # module


text:
using DataFrames, StatsBase
"""
    function comp_refs(dat...; nmask = 100, nthread=7)
---
This is to compare several references for imputation on a same target.
The last file in `dat` is the target.
All the files before this are references.

## Data required
- Genotypes in `final report` format of the 3 platforms.
- Map files for 3 the 3 platforms

## Methods
- Select shared loci
- Random mask 100, by defaut, loci of the target platform
- Impute with two references in turn.
- Compare the call rates of the two imputation.
- Make a report.
"""
function comp_refs(dat...; nmask = 100, window=300)
    length(dat) < 2 && (@error "Specify at least one ref, one target")
    tgt = dat[end]
    refs = dat[1:end-1]
    tbim = read_bim("$tgt.bim")

    lstsnp = begin
        shared = Set(tbim.snp)
        for r in refs
            tt = read_bim("$r.bim")
            shared = intersect(shared, Set(tt.snp))
        end
        sample(collect(shared), nmask, replace=false)
    end
    smpl = Set(lstsnp)

    dir, _ = splitdir(tgt)
    # sample SNP to be masked
    open("$dir/mask.snp", "w") do io
        for snp in tbim.snp
            snp ∈ smpl && println(io, snp)
        end
    end

    #- plink extract for later comparison
    @info "Extract $nmask SNP test set"
    run(pipeline(`$plink
			--sheep
			--bfile $tgt
			--extract $dir/mask.snp
			--recode vcf-iid
			--out $dir/mask`, devnull))

    @info "Prepare $tgt.vcf.gz with $nmask SNP removed"
    run(pipeline(`$plink
			--sheep
			--bfile $tgt
			--exclude $dir/mask.snp
			--recode vcf-iid bgz
			--out $tgt`, devnull))

    @info "Impute the references about their few missing genotypes"
    for r in refs
        @info "Prepare reference $r"
        run(pipeline(`$plink
			--sheep
			--bfile $r
			--recode vcf-iid bgz
			--out $r`, devnull))
        run(pipeline(`java -jar $beagle
                        window=$window
			ne=$nsNe
			gt=$r.vcf.gz
			out=$r-p`, `grep Window`))
    end                         # p for phased

    @info "Impute $tgt with each of the references"
    for r in refs
        @info "Impute with $r-p"
        run(pipeline(`java -jar $beagle
                        window=$window
			ne=$nsNe
			ref=$r-p.vcf.gz
			gt=$tgt.vcf.gz
			out=$r-i`, `grep Window`))
        run(pipeline(`$plink --sheep
			     --vcf $r-i.vcf.gz
			     --extract $dir/mask.snp
			     --make-bed
			     --recode vcf-iid
			     --out $r-m`, devnull))
    end
    for r in refs
        @info "Imputation results with $r"
        comp_vcf("$r-m.vcf", "$dir/mask.vcf")
    end
end

"""
    function comp_vcf(va, vb)
---
Compare two plain (not zipped) VCF files of their alleles line by line.
Return proportion of agreed allles and genotypes.
"""
function comp_vcf(va, vb)
    @info join(["Comparing", "- $va", "- $vb"], "\n")
    fa = open(va, "r")
    fb = open(vb, "r")
    mg = ma = 0       # m:misimputed; g:genotype; a:allele

    # skip headers
    n1 = n2 = 0
    for line in eachline(fa)
        if line[2] ≠ '#'
            n1 = length(split(line))
            break
        end
    end
    for line in eachline(fb)
        if line[2] ≠ '#'
            n2 = length(split(line))
            break
        end
    end
    n1 ≠ n2 && @error "Not same ID set"
    nid = n1 - 9
    nms = nlc = 0

    for la in eachline(fa)
        nlc += 1
        ta = split(la)
        tb = split(readline(fb))
        snpa, snpb = (ta[3], tb[3])
        snpa ≠ snpb && @error "Not same SNP set"
        a1, a2 = (ta[4][1], ta[5][1])
        b1, b2 = (tb[4][1], tb[5][1])
        flip = a1 ≠ b1 || a2 ≠ b2
        flip && @warn "$snpa minor allele flipped"
        for i in 10:n1
            a, b, x, y = Int.((ta[i][1], ta[i][3], tb[i][1], tb[i][3]))
            if x == 46 || y == 46 || a == 46 || b == 46
                nms += 1         # Int('.') = 46
            else
                if flip
                    (a ≠ y) && (ma += 1)
                    (b ≠ x) && (ma += 1)
                else
                    (a ≠ x) && (ma += 1)
                    (b ≠ y) && (ma += 1)
                end
                ((a + b) ≠ (x + y)) && (mg += 1)
            end
        end
    end

    # report
    println(lpad("            Total ID: ", 40), nid)
    println(lpad("           Total SNP: ", 40), nlc)
    println(lpad("   Genotypes missing: ", 40), nms/nid/nlc)
    println(lpad("Genotypes misimputed: ", 40), mg/nid/nlc)
    println(lpad("  Alleles misimputed: ", 40), ma/nid/nlc/2)

    close(fa)
    close(fb)
end

"""
    function read_map(mp)
---
## Description

This read a map of 3 columns:

- SNP name
- Chromosome number
- Base pair position

into a dataframe.

## Note

Only the first appeared SNP will be included in the map.
Its other repeat(s) will be omitted.
"""
function read_map(mp)
    SNP = Set{String}()
    df = DataFrame(snp=String[], chr=Int8[], bp=Int32[])
    for line in eachline(mp)
        snp, chr, bp = split(line)
        chr = parse(Int8, chr)
        bp = parse(Int32, bp)
        if snp ∉ SNP            # remove duplicate
            push!(df, [snp chr bp])
            push!(SNP, snp)
        end
    end
    df
end

"""
    function read_bim(bim)
---
To read a plink bim file, assuming all chromosomes are integer.
The 6 columns:

1. Chromosome
2. SNP name
3. Morgen
4. BP position
5. Allele 1
6. Allele 2
"""
function read_bim(bim)
    df = DataFrame(chr=Int8[],
                   snp=String[],
                   mgn=Float64[], # morgen
                   bp=Int32[],
                   a1=Char[],
                   a2=Char[])
    for line in eachline(bim)
        chr, snp, mgn, bp, a1, a2 = split(line)
        chr = parse(Int8, chr)
        mgn = parse(Float64, mgn)
        bp  = parse(Int32, bp)
        push!(df, [chr snp mgn bp a1[1] a2[1]])
    end
    df
end

"""
    make_bed(dir, lmap, target)
---
## Description
This function will merge the genotypes files of GSGT ver-2 format in `dir` into one
`vcf` file according linkage map `lmap`.

## dir
In the future, put all the files from a same platform in the same folder.  This
folder should contains the genotype files of GSGT ver-2 format.  Nothing else.

## lmap
This file contains 3 columns:
1. SNP name
2. Chromsome number.  Here they 1-26 autosomes only.
3. Base pair position.

The SNP are ordered on their chromosome numbers and bp positions.

## Example
    NSG.make_bed("dat/raw/a17k", "dat/maps/a17k.map", "dat/plink/a17k")

Above will merge all the files in `dat/raw/a17k` into `dat/plink/a17k.{bed,bim,fam}`.
"""
function make_bed(dir, lmap, target)
    title("Merge files in $dir into $target.bed")
    item("Check parameters")
    if !isdir(dir)
        warning("$dir doesn't exist")
        return
    end

    dir = abspath(dir)
    files = readdir(dir, join=true)
    if length(files) == 0
        warning("No file in $dir")
    end

    if !isfile(lmap)
        warning("Map $lmap doesn't exist")
        return
    end
    lmap = abspath(lmap)

    target =  begin             # check files to be saved
        pt, target = splitdir(target)
        if length(pt) == 0
            warning("You are saving the bed files in the current dir")
            pt = pwd()
        else
            pt = abspath(pt)
            isdir(pt) || mkpath(pt)
        end
        joinpath(pt, target)
    end

    done("OK")

    item("Create beagle files")
    tmap = read_map(lmap)
    ID, GT = read_fr(files...)

    temp = mktempdir(".")
    write_bgl(tmap, ID, GT, temp)

    item("Create VCF files")
    for chr in 1:26
        run(pipeline(`java -jar $bin_dir/beagle2vcf.jar $chr $temp/$chr.mrk $temp/$chr.bgl -`,
                     "$temp/$chr.vcf"))
    end

    cp("$temp/1.vcf", "$target.vcf", force = true)
    open("$target.vcf", "a") do io
        for chr in 2:26
            for line in eachline("$temp/$chr.vcf")
                line[1] ≠ '#' && println(io, line)
            end
        end
    end
    done()

    item("Make $target.bed")
    _ = read(run(`$plink --sheep --vcf $target.vcf --const-fid --out $target`), String)
    rm(temp, recursive=true, force=true)
    done()
end

"""
    function write_bgl(lmp, ID, GT, dir)
---
Write linkage map `lmp`, ID vector `ID`, and genotype dictionary `GT`
into Beagle files in `dir.
"""
function write_bgl(lmp, ID, GT, dir)
    CHR = lmp.chr[1]
    BP = -1
    mrk = open("$dir/$CHR.mrk", "w")
    bgl = open("$dir/$CHR.bgl", "w")
    print(bgl, "I\tid\t")
    println(bgl, join(repeat(1:length(ID), inner=2), '\t'))
    for (snp, chr, bp) in eachrow(lmp)
        bp == BP && continue
        BP = bp
        if chr != CHR
            CHR = chr
            close(mrk)
            close(bgl)
            mrk = open("$dir/$chr.mrk", "w")
            bgl = open("$dir/$chr.bgl", "w")
            print(bgl, "I\tid\t")
            println(bgl, join(repeat(ID, inner=2), '\t'))
        end
        if haskey(GT, snp)
            alleles = Set(collect(GT[snp]))
            setdiff!(alleles, '-') # remove missing values
            length(alleles) !=2 && continue
            print(bgl, "M $snp ")
            println(bgl, join(collect(GT[snp]), ' '))
            println(mrk, "$snp $bp ", join(alleles, ' '))
        end
    end
    close(mrk)
    close(bgl)
end

"""
    function read_fr(files...)
---
read `ID` vector  and genotypes dictionary `GT`(`SNP => alleles`) for the `files...`,
which are of final report format.
"""
function read_fr(files...)
    ID = String[]
    GT = Dict{String, String}()
    for f in files
        open(f, "r") do io
            for line in eachline(io) # skip header
                (line == "[Data]") && break
            end
            line = readline(io) # ID line
            append!(ID, split(line))
            for line in eachline(io)
                t = split(line)
                if haskey(GT, t[1])
                    GT[t[1]] *= join(t[2:end])
                else
                    GT[t[1]]  = join(t[2:end])
                end
            end
        end
    end
    if length(unique(length.(values(GT)))) > 1
        @error "The final report files are of different format"
    end
    ID, GT
end

ext: jl


text:
using ChaosTools, PkgBenchmark

@benchgroup "lyapunov_discrete" begin
        ds = Systems.henon()
        @bench "henon" lyapunov($ds, 10000,
        d0 = 1e-9, threshold = 1e-6, Ttr = 100)

        ds = Systems.towel()
        @bench "towel" lyapunov($ds, 10000,
        d0 = 1e-9, threshold = 1e-6, Ttr = 100)

        ds = Systems.coupledstandardmaps(5)
        @bench "CSM" lyapunov($ds, 10000,
        d0 = 1e-9, threshold = 1e-6, Ttr = 100)
    end

@benchgroup "lyapunov_continuous" begin
    ds = Systems.lorenz()
    @bench "lorenz" lyapunov($ds, 1000.0, d0 = 1e-9, threshold = 1e-6)
    @bench "lorenz2" lyapunov($ds, 1000.0, d0 = 1e-9, threshold = 1e-6)

    ds = Systems.duffing(β = -1, ω = 1, f = 0.3)
    @bench "duffing" lyapunov($ds, 1000.0, d0 = 1e-9, threshold = 1e-6)
    @bench "duffing2" lyapunov($ds, 1000.0, d0 = 1e-9, threshold = 1e-6)
end


text:
## PARTE 1
cd(@__DIR__)
function play_game(input, stop)
    turns = zeros(Int, stop)
    turns[1:length(input)] .= input
    foreach(length(input)+1:stop) do i
        j = findprev(==(turns[i-1]), turns, i-2)
        isnothing(j) && return
        turns[i] = i - j - 1
    end
    turns
end

##
n = 200000
#input = [18, 11, 9, 0, 5, 1]
input = [0, 3, 6]
@time turns = play_game(input, n)
n => last(turns)


text:
using Luxor
using Test

# testing that functions that accept eg `:fill` also can be used with keyword `action=:fill`

expressions = [
:(arc(O, 70, 0, π/2, :stroke)),
:(arc(O, 70, 0, π/2, action=:stroke)),

:(arc(O, 30, 0, π/2, :fill)),
:(arc(O, 30, 0, π/2, action=:fill)),

:(arc2r(O, Point(-30, 50), Point(20, 40), :stroke)),
:(arc2r(O, Point(-30, 50), Point(20, 40), action=:stroke)),

:(arc2sagitta(Point(30, 50), Point(20, 50), 40, :stroke)),
:(arc2sagitta(Point(30, 50), Point(20, 50), 40, action=:stroke)),

:(box(-30, -30, 50, 50, :fill)),
:(box(-30, -30, 50, 50, action=:fill)),

:(box([O, Point(24, 20), Point(20, 30)], :fill)),
:(box([O, Point(24, 20), Point(20, 30)], action=:fill)),

:(box(BoundingBox(box([O, Point(24, 20), Point(20, 30)])), :fill)),
:(box(BoundingBox(box([O, Point(24, 20), Point(20, 30)])), action=:fill)),

:(begin
    t = Table(3, 3)
    pts = box(t, 2, 1, action=:stroke, vertices=false, reversepath=true)
    pts = box(t, 2, 1, :stroke, vertices=false, reversepath=true)
    t = Tiler(50, 50, 3, 3)
    pts = box(t, 5, action=:stroke, vertices=false, reversepath=true)
    pts = box(t, 5, :stroke, vertices=false, reversepath=true)
end),

:(box(Point(-30, -50), Point(20, 50), :fill)),
:(box(Point(-25, -20), Point(23, 20), action=:fill)),

:(box(O, 30, 60, 10, :stroke)),
:(box(O, 35, 60, 10, action = :stroke)),

:(box(O, 30, 80, [10, 10, 20, 20], :stroke)),
:(box(O, 35, 80, [10, 10, 20, 20], action = :stroke)),

:(pts = box(O, 20, 20, :stroke; reversepath=true, vertices=true); poly(pts, :fill)),
:(pts = box(O, 20, 20, action=:stroke; reversepath=true, vertices=true); poly(pts, action=:fill)),

:(carc(O, 60, π, π/3, :stroke)),
:(carc(O, 60, π, π/3, action=:stroke)),

:(carc(O, 50, π, π/3, :stroke)),
:(carc(O, 50, π, π/3, action=:stroke)),

:(carc2r(O, Point(30, 30), Point(-20, -20), :stroke)),
:(carc2r(O, Point(30, 30), Point(-20, -20), action=:stroke)),

:(carc2sagitta(Point(-20, -20), Point(30, 30), 30, :stroke)),
:(carc2sagitta(Point(-20, -20), Point(30, 30), 30, action=:stroke)),
:(circle(O, 30, :stroke)),
:(circle(O, 30, action=:stroke)),
:(circle(20, 10, 30, :stroke)),
:(circle(20, 10, 30, action=:stroke)),
:(circle(O, Point(20, -20), :stroke)),
:(circle(O, Point(20, -20), action=:stroke)),
:(circle(O, Point(20, -20), Point(-10, 30), :stroke)),
:(circle(O, Point(20, -20), Point(-10, 30), action=:stroke)),
:(circlepath(O, 20, :stroke; reversepath=true)),
:(circlepath(O, 20, action=:stroke; reversepath=true)),
:(crescent(Point(10, 10), 18, 25, :stroke; vertices=false, reversepath=true, steps=20)),
:(crescent(Point(10, 10), 18, 25, action=:stroke; vertices=false, reversepath=true, steps=20)),
:(crescent(O, 20, Point(30, 30), 30, :stroke; vertices=false, reversepath=true)),
:(crescent(O, 20, Point(30, 30), 30, action=:stroke; vertices=false, reversepath=true)),

:(ellipse(10, 20, 30, 40, :stroke)),
:(ellipse(10, 20, 30, 40, action=:stroke)),

:(ellipse(O, 40, 50, :stroke)),
:(ellipse(O, 30, 30, action=:stroke)),

:(ellipse(O - (10, 10), O - (10, 10), 50, :stroke; stepvalue=π/30, vertices=false, reversepath=true)),
:(ellipse(O - (10, 10), O - (10, 10), 50, action=:stroke; stepvalue=π/30, vertices=false, reversepath=true)),

:(ellipse(O - (10, 10), O - (10, 10), O, :stroke; stepvalue=π/13, vertices=false, reversepath=true)),
:(ellipse(O - (10, 10), O - (10, 10), O, action=:stroke; stepvalue=π/13, vertices=false, reversepath=true)),

:(epitrochoid(41, 7, 23, :stroke; close=false, stepby=π/60, vertices=false)),
:(epitrochoid(41, 7, 23, action=:stroke; close=false, stepby=π/60, vertices=false)),

:(hypotrochoid(44, 7, 27, :stroke; close=false, stepby=π/39, vertices=false)),
:(hypotrochoid(44, 7, 27, action=:stroke; close=false, stepby=π/39, vertices=false)),

:(line(O, Point(30, -40), :stroke)),
:(line(O, Point(30, -40), action=:stroke)),

:(ngon(20, 10, 20, 5, :stroke)),
:(ngon(20, 10, 20, 5, action=:stroke)),
:(ngon(20, 20, 20, 6, π/3, :stroke)),
:(ngon(20, 20, 20, 6, π/3, action=:stroke)),
:(ngon(O, 23, 7, π/3, :fill, vertices=false, reversepath=true)),
:(ngon(O, 23, 7, π/3, action=:fill, vertices=false, reversepath=true)),

:(ngonside(O, 50, 7, π/3, :stroke)),
:(ngonside(O, 50, 7, π/3, action=:stroke)),
:(ngonside(O, 30, 7, :stroke)),
:(ngonside(O, 30, 7, action=:stroke)),

:(pie(O, 23, 0, π, :fill)),
:(pie(O, 23, 0, π, action=:fill)),
:(pie(20, 10, 50, π/3, 0, :stroke)),
:(pie(20, 10, 50, π/3, 0, action=:stroke)),
:(pie(50, 0, π/2, :stroke)),
:(pie(50, 0, π/2, action=:stroke)),

:(poly(BoundingBox(box(O, 20, 20)), :stroke)),
:(poly(BoundingBox(box(O, 20, 20)), action=:stroke)),

:(poly(ngon(O, 50, 7), :stroke, close=true, reversepath=false)),
:(poly(ngon(O, 50, 7), action=:stroke, close=true, reversepath=false)),

:(polycross(O, 23, 3, 0.7, π/3, :fill, vertices=false, reversepath=true)),
:(polycross(O, 23, 3, 0.7, π/3, action=:fill, vertices=false, reversepath=true)),

:(prettypoly(box(O, 20, 20), close=true, :stroke)),
:(prettypoly(box(O, 20, 20), action=:stroke, close=true)),

:(prettypoly(ngon(O, 30, 5, vertices=true), close=true, :none, () -> (star(O, 16, 7, 0.4, 0.0, action=:stroke)))),
:(prettypoly(ngon(O, 30, 5, vertices=true), close=true, action=:none, () -> (star(O, 16, 7, 0.4, 0.0, action=:stroke)))),

:(prettypoly(BoundingBox() * 0.05, :fill)),
:(prettypoly(BoundingBox() * 0.05, action=:fill)),

:(rect(-20, -30, 15, 29, :stroke)),
:(rect(-20, -30, 15, 29, action=:stroke)),

:(rect(O, 26, 31, :stroke; reversepath=true, vertices=false)),
:(rect(O, 26, 31, action=:stroke; reversepath=true, vertices=false)),

:(sector(O, 28, 45, 0, 3π/2, :stroke)),
:(sector(O, 28, 45, 0, 3π/2, action=:stroke)),

:(sector(18, 40, 0, 3π/2, :stroke)),
:(sector(18, 40, 0, 3π/2, action=:stroke)),

:(sector(O, 18, 45, 0, 3π/2, 10, :stroke)),
:(sector(O, 18, 45, 0, 3π/2, 10, action=:stroke)),

:(sector(18, 25, 0, 3π/2, 10, :stroke)),
:(sector(18, 25, 0, 3π/2, 10, action=:stroke)),

:(spiral(20, -2, :stroke; stepby=π/30, period=12π, vertices=false, log=false)),
:(spiral(20, -2, action=:stroke; stepby=π/30, period=12π, vertices=false, log=false)),

:(squircle(O, 23, 23, :stroke; rt=0.6, vertices=false, stepby=π/30, reversepath=true)),
:(squircle(O, 23, 23, action=:stroke; rt=0.6, vertices=false, stepby=π/30, reversepath=true)),

:(star(20, 10, 23, 6, 0.5, 0.0, :fill; vertices=false, reversepath=true)),
:(star(20, 10, 23, 6, 0.5, 0.0, action=:fill; vertices=false, reversepath=true)),

:(star(O, 33, 7, 0.3, π, :stroke; vertices=false, reversepath=true)),
:(star(O, 33, 7, 0.3, π, action=:stroke; vertices=false, reversepath=true)),

:(textoutlines("text", Point(0, 0), action=:path; halign=:left, valign=:middle, startnewpath=true); fillpath()),
:(textoutlines("text", Point(0, 0), :stroke, halign=:left, valign=:middle, startnewpath=true)),
]

fname = "action-dispatch-tests.pdf"
Drawing(2000, 2000, fname)
origin()
background("black")
fontsize(10)
tiles = Tiler(2000, 2000, 12, 12)
for (pos, n) in tiles
    @layer begin
        translate(pos)
        if n <= length(expressions)
            sethue("white")
            textwrap(string(expressions[n]), tiles.tilewidth * 0.7, O - (tiles.tilewidth/2, tiles.tileheight/2), leading=10)
            sethue("gold")
            eval(expressions[n])
        end
    end
end

@test finish() == true
println("...finished test: output in $(fname)")


text:
#TODO
# compute partial derivatives by chain-ruling built-in T-Dist pdf and cdf with derivatives of sub-functions

include("kernel.jl")
include("transforms.jl")
include("examples.jl")
include("buffers.jl")
include("priors.jl")
include("settings.jl")

using Distributions
using Printf
using SpecialFunctions
using Plots
using Polynomials
using LinearAlgebra

#export prob, partial_theta, partial_lambda, partial_z0, posterior_theta, posterior_lambda, checkDerivative

"""
Exemplifies how slow regular K can be (when for loop are used)
"""
function time_func_artifact(θ::Float64, setting::setting{Array{Float64, 2}, Array{Float64, 1}})
    s = setting.s
    s0 = setting.s0
    X = setting.X
    X0 = setting.X0
    z = setting.z
    reset_timer!()
    @timeit "Eθ" Eθ = K(s0, s0, θ, rbf)
    @timeit "Σθ" Σθ =  K(s, s, θ, rbf)
    @timeit "fasterΣθ" Σθ =  fastK(s, s, θ, rbf_single)
    @timeit "Bθ" Bθ =  K(s0, s, θ, rbf)
    @timeit "choleskyΣθ" choleskyΣθ = cholesky(Σθ)
    @timeit "choleskyXΣX" choleskyXΣX = cholesky(Hermitian(X'*(choleskyΣθ\X)))
    @timeit "Dθ" Dθ = Eθ - Bθ*(choleskyΣθ\Bθ')
    @timeit "Hθ" Hθ = X0 - Bθ*(choleskyΣθ\X)
    @timeit "Cθ" Cθ = Dθ + Hθ*(choleskyXΣX\Hθ')
    @timeit "Eθ_prime" Eθ_prime = K(s0, s0, θ, rbf_prime)
    @timeit "Eθ_prime2" Eθ_prime2 = K(s0, s0, θ, rbf_prime2)
    @timeit "Σθ_prime" Σθ_prime = K(s, s, θ, rbf_prime)
    @timeit "Σθ_prime_faster" Σθ_prime = fastK(s, s, θ, rbf_prime_single)
    @timeit "Σθ_prime2" Σθ_prime2 = K(s, s, θ, rbf_prime2)
    @timeit "Σθ_prime2_faster" Σθ_prime2 = fastK(s, s, θ, rbf_prime2_single)
    @timeit "Bθ_prime" Bθ_prime = K(s0, s, θ, rbf_prime)
    @timeit "Bθ_prime2" Bθ_prime2 = K(s0, s, θ, rbf_prime2)
    print_timer()
    return θ_params(Eθ, Σθ, Bθ, Dθ, Hθ, Cθ, Eθ_prime, Eθ_prime2, Σθ_prime, Σθ_prime2, Bθ_prime, Bθ_prime2, choleskyΣθ, choleskyXΣX)
end

"""
Compute derivative of beta hat with respect to theta
"""
function compute_betahat_prime_theta(choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, choleskyXΣX::Cholesky{Float64,Array{Float64, 2}}, expr_mid, Σθ_prime, X::Array{Float64,2}, gλz, Σθ_inv_X)
    AA = choleskyXΣX\(expr_mid)*(choleskyXΣX\(X'*(choleskyΣθ\gλz)))
    BB = - (choleskyXΣX\(X'*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\gλz)))))
    βhat_prime_theta = AA + BB
    βhat_prime_theta = reshape(βhat_prime_theta, size(βhat_prime_theta, 1), size(βhat_prime_theta, 2)) #turn 1D array into 2D array
end

"""
Compute second derivative of beta hat
Q: - derivative of Sigma_theta^-1 (function handle)
dQ: derivative of Q (function handle)
gλz: g_lambda(z)
"""
function compute_βhat_prime2_theta(choleskyXΣX::Cholesky{Float64,Array{Float64, 2}}, choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, expr_mid, X::Array{Float64,2}, Q, dQ, gλz)
    XQX = X'*Q(X) #precompute
    Qgλz = Q(gλz) #precompute
    βhatEPXR1 = - (choleskyXΣX\(-expr_mid * (choleskyXΣX\((XQX * (choleskyXΣX\(X'*(choleskyΣθ\gλz))))))))
    βhatEPXR2 = choleskyXΣX\((X'*dQ(X*(choleskyXΣX\(X'*(choleskyΣθ\gλz))))))
    βhatEPXR3  = - (choleskyXΣX\((XQX*(choleskyXΣX\(-expr_mid * (choleskyXΣX\(X'*(choleskyΣθ\(gλz)))))))))
    βhatEPXR4 = - (choleskyXΣX\(XQX*(choleskyXΣX\(X'*Qgλz))))
    βhatEPXR5 = (choleskyXΣX\(-expr_mid * (choleskyXΣX\(X'*Qgλz))))
    βhatEPXR6 = - (choleskyXΣX\(X'*dQ(gλz)))
    βhat_prime2_theta = βhatEPXR1 + βhatEPXR2 + βhatEPXR3 + βhatEPXR4 + βhatEPXR5 + βhatEPXR6
end

"""
First derivative of qtilde with respect to theta
"""
function compute_qtilde_prime_theta(gλz, X::Array{Float64,2}, βhat, βhat_prime_theta, choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, Σθ_prime::Array{Float64,2})
    meanvv = gλz - X*βhat
    rr = X*βhat_prime_theta
    AA = (-rr)' * (choleskyΣθ \ meanvv)
    BB = - meanvv' * (choleskyΣθ \ (Σθ_prime * (choleskyΣθ \ meanvv)))
    CC =  meanvv' * (choleskyΣθ \ (-rr))
    qtilde_prime_theta = AA .+ BB .+ CC
end

"""
Second derivative of qtilde with respect to theta
"""
function compute_qtilde_prime2_theta(choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, X::Array{Float64,2}, meanvv, Q, dQ, βhat_prime_theta, βhat_prime2_theta)
    Xdβ = X*βhat_prime_theta
    qtildeEXPR1 = -meanvv'*dQ(meanvv)
    qtildeEXPR2 = 2*(X*βhat_prime_theta)'*Q(meanvv)
    qtildeEXPR3 = 2*Xdβ'*(choleskyΣθ\(Xdβ))
    qtildeEXPR4 = 2*meanvv'*Q(Xdβ)
    qtildeEXPR5 = -2*meanvv'*(choleskyΣθ\(X*βhat_prime2_theta))
    qtilde_prime2_theta = qtildeEXPR1 .+ qtildeEXPR2 .+ qtildeEXPR3 .+ qtildeEXPR4 .+ qtildeEXPR5
end

"""
First derivative of Htheta with respect to theta
"""
function compute_Hθ_prime(Bθ_prime::Array{Float64,2}, Σθ_inv_X::Array{Float64,2}, choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, Σθ_prime::Array{Float64,2}, Bθ::Array{Float64,2}, X::Array{Float64,2})
    #compute Hθ_prime
    AA = -Bθ_prime*Σθ_inv_X
    #BB = Bθ*Σθ\(Σθ_prime*(Σθ\X)) displaying the bug
    BB = Bθ*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\X)))
    Hθ_prime = AA + BB
end

"""
Second derivative of Htheta with respect to theta
"""
function compute_Hθ_prime2(Bθ::Array{Float64,2}, Bθ_prime::Array{Float64,2}, Bθ_prime2::Array{Float64,2}, Σθ_inv_X, choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, X, Q, dQ)
    AA = - Bθ_prime2*(choleskyΣθ\X)
    BB = 2*Bθ_prime*Q(X)
    CC = Bθ*dQ(X)
    AA+BB+CC
end

"""
First derivative of m_theta with respect to theta
"""
function compute_m_prime_theta(Bθ::Array{Float64,2}, Bθ_prime::Array{Float64,2}, choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, Σθ_prime::Array{Float64,2}, gλz, βhat, βhat_prime_theta,Hθ::Array{Float64,2}, Hθ_prime::Array{Float64,2})
    AA = Bθ_prime*(choleskyΣθ\gλz)
    BB = - Bθ*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\gλz)))
    CC = Hθ_prime*βhat
    DD = Hθ*βhat_prime_theta
    m_prime_theta = AA + BB + CC + DD
end

"""
First derivative of m_theta with respect to theta
"""
function compute_m_prime2_theta(Bθ::Array{Float64,2}, Bθ_prime::Array{Float64,2}, Bθ_prime2::Array{Float64,2}, choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, Σθ_prime::Array{Float64,2}, gλz, βhat, βhat_prime_theta,
    βhat_prime2_theta, Hθ, Hθ_prime, Hθ_prime2, Q, dQ)
    AA = Bθ_prime2*(choleskyΣθ\gλz) - 2*Bθ_prime*(Q(gλz))
    BB = -Bθ*(dQ(gλz))
    CC = Hθ_prime2*βhat+2*Hθ_prime*βhat_prime_theta
    DD = Hθ*βhat_prime2_theta

    m_prime_theta2 = AA + BB + CC + DD
end

"""
First derivative of D_theta with respect to theta
"""
function compute_Dθ_prime(choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, Bθ::Array{Float64,2}, Eθ_prime::Array{Float64,2}, Σθ_prime::Array{Float64,2}, Bθ_prime::Array{Float64,2})
    sigma_inv_B = choleskyΣθ \ Bθ' #precomputation
    AA = Eθ_prime - Bθ_prime * sigma_inv_B
    BB = sigma_inv_B' * Σθ_prime * sigma_inv_B
    CC = - sigma_inv_B' * Bθ_prime'
    Dθ_prime = AA + BB + CC
end

"""
Second derivative of D_theta with respect to theta
"""
function compute_Dθ_prime2(choleskyΣθ::Cholesky{Float64,Array{Float64, 2}}, Bθ::Array{Float64,2}, Bθ_prime::Array{Float64,2}, Bθ_prime2::Array{Float64,2}, Eθ_prime2::Array{Float64,2}, Q, dQ)
    AA = Eθ_prime2 - Bθ_prime2*(choleskyΣθ\Bθ') - Bθ*(choleskyΣθ\Bθ_prime2')
    BB = Bθ*dQ(Bθ') + 2*Bθ*Q(Bθ_prime') + 2*Bθ_prime*Q(Bθ')
    CC = -2*Bθ_prime*(choleskyΣθ\Bθ_prime')
    AA + BB + CC
end

"""
First derivative of C_theta with respect to theta
"""
function compute_Cθ_prime(Dθ_prime::Array{Float64,2}, Hθ::Array{Float64,2},  Hθ_prime::Array{Float64,2}, choleskyXΣX::Cholesky{Float64,Array{Float64, 2}}, Σθ_inv_X::Array{Float64,2}, Σθ_prime::Array{Float64,2})
    AA = Dθ_prime + Hθ_prime*(choleskyXΣX\Hθ')
    BB = Hθ*(choleskyXΣX\(Σθ_inv_X'*Σθ_prime*Σθ_inv_X))*(choleskyXΣX\Hθ')
    CC = Hθ*(choleskyXΣX\(Hθ_prime'))
    C_theta_prime = AA + BB + CC
end

"""
Second derivative of C_theta with respect to theta
"""
function compute_Cθ_prime2(Dθ_prime2::Array{Float64,2}, Hθ::Array{Float64,2}, Hθ_prime::Array{Float64,2}, Hθ_prime2::Array{Float64,2}, choleskyXΣX::Cholesky{Float64,Array{Float64, 2}}, dPinv, d2Pinv)
    AA = Dθ_prime2 + Hθ_prime2 * (choleskyXΣX\Hθ')
    BB = Hθ * (choleskyXΣX\Hθ_prime2')
    CC = Hθ*(d2Pinv(Hθ'))
    DD = 2*(Hθ_prime * (choleskyXΣX\Hθ_prime') + Hθ_prime * dPinv(Hθ') + Hθ*(dPinv(Hθ_prime'))  )
    AA+BB+CC+DD
end

"""
p(z0| theta, lambda, z). Computes T-Distribution explicitly (and unstably).
"""
function prob_artifact(θ, λ, setting, type = "Gaussian")
    s = setting.s; s0 = setting.s0; X = setting.X; X0 = setting.X0; z = setting.z; n = size(X, 1); p = size(X, 2); k = size(X0, 1)  #unpack setting
    theta_params = funcθ(θ, train, test, type)
    g = boxCox #boxCox by default
    Eθ = theta_params.Eθ
    Σθ = theta_params.Σθ
    Bθ = theta_params.Bθ
    Dθ = theta_params.Dθ
    Hθ = theta_params.Hθ
    Cθ = theta_params.Cθ
    Σθ_inv_X = theta_params.Σθ_inv_X
    choleskyΣθ = theta_params.choleskyΣθ
    choleskyXΣX = theta_params.choleskyXΣX

    βhat = (X'*(choleskyΣθ\X))\(X'*(choleskyΣθ\g(z, λ)))
    qtilde = (expr = g(z, λ)-X*βhat; expr'*(choleskyΣθ\expr))
    m = Bθ*(choleskyΣθ\g(z, λ)) + Hθ*βhat
    #t = LocationScale(m[1], sqrt(qtilde[1]*Cθ[1]/(n-p)), TDist(n-p))
    #p = z0 -> Distributions.pdf(t, g(z0, λ))
    #return p
    cc = gamma((n-p+k)/2)/gamma((n-p)/2)/pi^(k/2) #constant term
    expr = z0 -> g(z0, λ) .- m
    return z0 -> cc*(det(qtilde*Cθ)^(-1/2))*(1+expr(z0)'*((qtilde*Cθ)\expr(z0)))^(-(n-p+k)/2)
end

"""
Compute derivative of p(z0|theta, lambda, z) w.r.t theta. This function should compute things that depend on BOTH theta and lambda. It
takes in all pertinent theta-dependent quantities as inputs.

WARNING: not tested, because want to get rid of cc (the unstably computed constant term), and redo some computations to use built-in cdf and pdf
"""
function partial_theta(θ::Float64, λ::Float64, train::trainingData{A, B}, test::testingData{A}, transforms, theta_params::Union{θ_params{A},
    Cholesky{Float64, A}, θ_param_derivs{A, Cholesky{Float64, A}}, Nothing} = nothing, type::String = "Gaussian") where A<:Array{Float64, 2} where B<:Array{Float64, 1}
    if theta_params === nothing
        #println("WARNING: recompute theta_params in partial_theta")
        theta_params = funcθ(θ, train, test, type)
        #(Eθ, Σθ, Bθ, Dθ, Hθ, Cθ, Eθ_prime, Eθ_prime2, Σθ_prime, Σθ_prime2, Bθ_prime, Bθ_prime2, choleskyΣθ, choleskyXΣX) = getθ_Params(theta_params)
        #(s, s0, X, X0, z,  n, p, k) = getSettingParams(setting)
    end
    if type == "Gaussian" || type == "Turan"
        Eθ = theta_params.Eθ
        Σθ = theta_params.Σθ
        Bθ = theta_params.Bθ
        Dθ = theta_params.Dθ
        Hθ = theta_params.Hθ
        Cθ = theta_params.Cθ
        Σθ_inv_X = theta_params.Σθ_inv_X
        choleskyΣθ = theta_params.choleskyΣθ
        choleskyXΣX = theta_params.choleskyXΣX
    else
        throw(ArgumentError("Quadrature type undefined. Please enter \"Gaussian\" or \"Turan\" for last arg."))
    end
    if type == "Turan" #load higher derivatives
        Eθ_prime = theta_params.Eθ_prime
        Eθ_prime2 = theta_params.Eθ_prime2
        Σθ_prime = theta_params.Σθ_prime
        Σθ_prime2 = theta_params.Σθ_prime2
        Bθ_prime = theta_params.Bθ_prime
        Bθ_prime2 = theta_params.Bθ_prime2
        Dθ_prime = theta_params.Dθ_prime
        Dθ_prime2 = theta_params.Dθ_prime2
        Hθ_prime = theta_params.Hθ_prime
        Hθ_prime2 = theta_params.Hθ_prime2
        Cθ_prime = theta_params.Cθ_prime
        Cθ_prime2 = theta_params.Cθ_prime2
        tripleΣ = theta_params.tripleΣ
    end

    s = train.s; s0 = test.s0; X = train.X; X0 = test.X0; z = train.z; n = size(X, 1); p = size(X, 2); k = size(X0, 1)  #unpack

    g = boxCox #boxCox by default
    gλz = g(z, λ)

    @timeit "βhat" βhat = (X'*(choleskyΣθ\X))\(X'*(choleskyΣθ\g(z, λ)))
    @timeit "qtilde" qtilde = (expr = g(z, λ)-X*βhat; expr'*(choleskyΣθ\expr))
    @timeit "m" m = Bθ*(choleskyΣθ\g(z, λ)) + Hθ*βhat
    meanvv = gλz - X*βhat

    dg = transforms.df
    jac = z0 -> abs(reduce(*, map(x -> dg(x, λ), z0)))

    cc = gamma((n-p+k)/2)/gamma((n-p)/2)/pi^(k/2) #constant term
    t = LocationScale(m[1], sqrt(qtilde[1]*Cθ[1]/(n-p)), TDist(n-p))
    vanillat = LocationScale(0, 1, TDist(n-p))
    vanillaT_pdf = z0 -> Distributions.pdf.(vanillat, z0)
    T_pdf = z0 -> Distributions.pdf.(t, g(z0, λ)) #no Jacobian
    main_pdf = z0 -> Distributions.pdf.(t, g(z0, λ)) * jac(z0)
    main_cdf = z0 -> Distributions.cdf.(t, g(z0, λ))
    #This block of code computes first and second derivatives.
    if type == "Turan"
        @timeit "expr_mid" expr_mid = X'*(choleskyΣθ\(Σθ_prime * Σθ_inv_X))#precompute
        #first derivatives
        @timeit "some first derivs, which depend on θ AND λ" begin
        βhat_prime_theta = compute_betahat_prime_theta(choleskyΣθ, choleskyXΣX, expr_mid, Σθ_prime, X, gλz, Σθ_inv_X)
        qtilde_prime_theta = compute_qtilde_prime_theta(gλz, X, βhat, βhat_prime_theta, choleskyΣθ, Σθ_prime)
        m_prime_theta = compute_m_prime_theta(Bθ, Bθ_prime, choleskyΣθ,Σθ_prime, gλz, βhat, βhat_prime_theta,Hθ, Hθ_prime)
        end
        #copied from func(theta, examples) in structures.jl
        dQ = Y -> (choleskyΣθ\((Σθ_prime2*(choleskyΣθ\Y)))) - 2*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\(Y))))))
        Q = Y -> (choleskyΣθ\(Σθ_prime*(choleskyΣθ\(Y))))
        #second derivatives
        #@timeit "second derivs" begin
        @timeit "βhat_prime2_theta" βhat_prime2_theta  = compute_βhat_prime2_theta(choleskyXΣX, choleskyΣθ, expr_mid, X, Q, dQ, gλz)
        @timeit "qtilde_prime2_theta" qtilde_prime2_theta = compute_qtilde_prime2_theta(choleskyΣθ, X, meanvv, Q, dQ, βhat_prime_theta, βhat_prime2_theta)
        #@timeit "Hθ_prime2" Hθ_prime2 = compute_Hθ_prime2(Bθ, Bθ_prime, Bθ_prime2, Σθ_inv_X, choleskyΣθ, X, Q, dQ)
        @timeit "m_prime2_theta" m_prime2_theta = compute_m_prime2_theta(Bθ, Bθ_prime, Bθ_prime2, choleskyΣθ, Σθ_prime, gλz, βhat, βhat_prime_theta,
        βhat_prime2_theta, Hθ, Hθ_prime, Hθ_prime2, Q, dQ)
        #@timeit "Dθ_prime2" Dθ_prime2 = compute_Dθ_prime2(choleskyΣθ, Bθ, Bθ_prime, Bθ_prime2, Eθ_prime2, Q, dQ)
        #@timeit "Cθ_prime2" Cθ_prime2 = compute_Cθ_prime2(Dθ_prime2, Hθ, Hθ_prime, Hθ_prime2, choleskyXΣX, dPinv, d2Pinv)
        #end
        @timeit "main deriv partial_theta" begin
        #compute derivative of main expression

        expr = z0 -> g(z0, λ) .- m
        qC = qtilde*Cθ
        bilinearform = z0 -> 1 .+ expr(z0)'*(qC\(expr(z0)))
        detqC = det(qC)
        qC_prime_theta =  qtilde_prime_theta .* Cθ + qtilde .* Cθ_prime
        AA = -0.5 * detqC^(-1/2) * tr(qC\(qC_prime_theta))
        qC_inv_prime_theta = Y -> - qC\(qC_prime_theta * (qC\Y))

        dbilinearform = z0 -> -m_prime_theta'*(qC\(expr(z0))) .+ expr(z0)' * qC_inv_prime_theta(expr(z0)) .- expr(z0)'*(qC\m_prime_theta)
        EE = z0 -> bilinearform(z0)^(-(n-p+k)/2)
        FF = z0 -> detqC^(-1/2) * (-(n-p+k)/2) * (bilinearform(z0))^(-(n-p+k+2)/2)
        dT_pdf = z0 -> (cc*(AA*EE(z0) .+ FF(z0)*(dbilinearform(z0))))[1]
        dmain_pdf = z0 -> (cc*(AA*EE(z0) .+ FF(z0)*(dbilinearform(z0))))[1] * jac(z0)
        end
        @timeit "main second deriv partial_theta" begin
        #compute second derivative of main expression
        qC_prime2_theta = qtilde_prime2_theta .* Cθ + qtilde .* Cθ_prime2 + 2* Cθ_prime .* qtilde_prime_theta
        qC_inv_prime2_theta = Y -> 2*(qC\(qC_prime_theta*(qC\(qC_prime_theta*(qC\Y))))) - qC\(qC_prime2_theta*(qC\Y))

        PP = 0.25 * (detqC)^(-1/2) * (tr(qC\qC_prime_theta))^2
        QQ = -0.5*detqC^(-1/2) * tr(qC_inv_prime_theta(qC_prime_theta) + qC\(qC_prime2_theta))
        dAA = PP+QQ #second derivative of det(qC)^-1/2

        d2bilinearform = z0 -> (- m_prime2_theta'*(qC\expr(z0)) .- expr(z0)'*(qC\m_prime2_theta) .+ expr(z0)'*qC_inv_prime2_theta(expr(z0))
                                .+ 2*(m_prime_theta'*(qC\(m_prime_theta)) .- m_prime_theta'*(qC_inv_prime_theta(expr(z0))) .- expr(z0)'*(qC_inv_prime_theta(m_prime_theta)))
        )
        bformpower = z0 -> bilinearform(z0)^(-(n-p+k)/2)
        dbformpower = z0 -> -((n-p+k)/2)*bilinearform(z0)^(-(n-p+k+2)/2)*dbilinearform(z0)
        d2bformpower = z0 -> (n-p+k)/2 * (n-p+k+2)/2 * bilinearform(z0)^(-(n-p+k+4)/2)*dbilinearform(z0)^2 - (n-p+k)/2*bilinearform(z0)^(-(n-p+k+2)/2)*d2bilinearform(z0)
        d2main_pdf = z0 -> (cc* (dAA*bformpower(z0) .+ 2*AA*dbformpower(z0) .+ detqC^(-1/2)*d2bformpower(z0)))[1] * jac(z0)
        end

        #derivatives of g(z0)-m/sqrt(qC/(n-p)) w.r.t theta for chain rule application
        locationscale = z0 -> expr(z0) ./ sqrt(qC/(n-p))
        dlocationscale = z0 -> - m_prime_theta ./ sqrt(qC ./ (n-p)) - 0.5 .* expr(z0) ./ (qC ./ (n-p))^1.5 .* qC_prime_theta/(n-p)
        d2locationscale = z0 ->  - m_prime2_theta ./ sqrt(qC/(n-p)) .+ m_prime_theta ./(qC/(n-p))^1.5 .* qC_prime_theta/(n-p) .+ expr(z0) * sqrt(n-p) .* (0.75*(qC)^(-5/2) .* (qC_prime_theta)^2 - 0.5 * qC^(-3/2) .* qC_prime2_theta)
        dmain_cdf = z0 -> vanillaT_pdf(expr(z0)/(sqrt(qC/(n-p)))) * dlocationscale(z0) #T_pdf(z0) .* dlocationscale(z0) #+ main_cdf(z0) .* (-0.5 ./ qC .* qC_prime_theta)

        #Below we write down some helper functions, which will allow us to compute derivatives using the built-in TDist CDF and PDF.
        #Observe that d4' = d5 and d6' = d7
        d4 = z0 -> vanillaT_pdf(expr(z0)/(sqrt(qC/(n-p))))/sqrt(qC/(n-p))
        d5 = z0 -> dT_pdf(z0)
        d6 = z0 -> vanillaT_pdf(expr(z0)/(sqrt(qC/(n-p))))
        d7 = z0 -> sqrt(qC)*(dT_pdf(z0)/sqrt(n-p) .+ 0.5*vanillaT_pdf(expr(z0)/(sqrt(qC/(n-p)))) ./ (qC)^(3/2) .* qC_prime_theta)

        d2main_cdf = z0 -> d7(z0) .* dlocationscale(z0) .+ vanillaT_pdf(expr(z0)/(sqrt(qC/(n-p))))*d2locationscale(z0)

        return [[main_pdf, dmain_pdf, d2main_pdf], [main_cdf, dmain_cdf, d2main_cdf]]
        #return ((main_pdf, dmain_pdf, d2main_pdf), (main_cdf, dmain_cdf, d2main_cdf, d4, d5, d6, d7))
        #return ((main_pdf, dmain_pdf, d2main_pdf), (main_cdf, dmain_cdf, d2main_cdf), (locationscale, dlocationscale, d2locationscale))

        #currently pdf and cdf are decoupled, so that calling one followe by the other results in
        #repeated computation of certain quantities. However coupling them would result in extra
        #computation if only one is needed.
    end
    #return (vec(Bθ_prime), vec(Bθ_prime2))
    return [[main_pdf], [main_cdf]]
    #return (dmain, d2main)
    #return (dbformpower, d2bformpower)
    #return (dbilinearform, d2bilinearform)
    #return (qC_prime_theta, qC_prime2_theta)
end
    #return (vec(qC_inv_prime_theta(I)), vec(qC_inv_prime2_theta(I)))
    #return (qtilde_prime_theta, qtilde_prime2_theta)
    #return (Cθ_prime, Cθ_prime2)
    #return (vec(Bθ_prime), vec(Bθ_prime2))
    #return (m_prime_theta, m_prime2_theta)
    #return (vec(Eθ_prime), vec(Eθ_prime2))
    #return (vec(Bθ), vec(Bθ_prime))
    #return (vec(Bθ_prime), vec(Bθ_prime2))
    #return (vec(Q(X)), vec(dQ(X)))
    #return (vec(Dθ), vec(Dθ_prime))
    #return (vec(Dθ_prime), vec(Dθ_prime2))
    #return (vec(dP(X')), vec(d2P(X')))
    #return (vec(dPinv(X')), vec(d2Pinv(X')))
    #main_deriv = (AA*EE(z0) + FF(z0)*(BB(z0) + CC(z0) + DD(z0)))
    #return (vec(Σθ), vec(Σθ_prime))
    #return (βhat, βhat_prime_theta)
    #return ([qtilde], [qtilde_prime_theta])
    #return (vec(Eθ), vec(Eθ_prime))
    #return (vec(Bθ), vec(Bθ_prime))
    #return (vec(Σθ), vec(Σθ_prime))
    #return (vec(Hθ), vec(Hθ_prime))
    #return (m, m_prime_theta)
    #return (vec(Dθ), vec(Dθ_prime))
    #return (vec(Cθ), vec(C_theta_prime))
    #return (vec(qC_inv), vec(qC_inv_prime_theta))
    #return ([detqC^(-1/2)], [AA])
    #return (bilinearform, z0 -> BB(z0) + CC(z0) + DD(z0))

"""
Compute derivative of p(z0|theta, lambda, z) w.r.t lambda
"""
function partial_lambda(θ, λ, setting)
    (s, s0, X, X0, z, n, p, k, Eθ, Σθ, Bθ, choleskyΣθ, choleskyXΣX, Dθ, Hθ, Cθ, Eθ_prime,Eθ_prime2, Σθ_prime, Σθ_prime2, Bθ_prime, Bθ_prime2) = func(θ, setting)

    #transformation and its derivatives
    g = boxCox #boxCox by default
    dg = boxCoxPrime
    dgλ = boxCoxPrime_lambda
    dgλ2 = boxCoxPrime_lambda2
    dgλx = boxCoxMixed_lambda_z

    βhat = choleskyXΣX\(X'*(choleskyΣθ\g(z, λ)))
    qtilde = (expr = g(z, λ)-X*βhat; expr'*(choleskyΣθ\expr))
    m = Bθ*(choleskyΣθ\g(z, λ)) + Hθ*βhat

    #compute βhat_prime_lambda
    dgλz = dgλ(z, λ)
    gλz = g(z, λ)
    βhat_prime_lambda = choleskyXΣX\(X'*(choleskyΣθ\dgλz))

    #compute qtilde_prime_lambda
    AA = dgλz - X*βhat_prime_lambda
    BB = gλz - X*βhat
    qtilde_prime_lambda = 2*AA'*(choleskyΣθ\BB)
    #compute m_prime_lambda
    m_prime_lambda = Bθ*(choleskyΣθ\dgλz)+Hθ*βhat_prime_lambda

    #compute main expression
    jac = z0 -> abs(reduce(*, map(x -> dg(x, λ), z0)))
    qC = qtilde*Cθ
    detqC = det(qC)
    expr = z0 -> (g(z0, λ) .- m)

    bilinearform = z0 -> 1 .+ expr(z0)'*(qC\(expr(z0)))

    #compute derivative of Jacobian
    function djac(z0)
        dgλx_vec = z0 -> map(zi -> dgλx(zi, λ), z0)
        dg_vec = z0 -> map(zi -> dg(zi, λ), z0)
        local prod = z0 -> reduce(*, dg_vec(z0))
        local p = prod(z0)
        p != 0 ? (AAA = dgλx_vec(z0);
        BBB = dg_vec(z0);
        ff = i -> p * AAA[i]/BBB[i] ;
        ee = reduce(+, map(ff, collect(1:length(z0))));
        dg(2, 2) < 0 ? (-1)^(length(z0))*ee : ee) : 0
    end

    #return jac, djac

    #compute derivative of main expression
    EXPR1 = z0 -> jac(z0)
    EXPR2 = z0 -> detqC^(-1/2)
    EXPR3 = z0 -> bilinearform(z0)^(-(n-p+k)/2)
    dEXPR1 = z0 -> djac(z0)
    dEXPR2 = z0 -> -0.5 * detqC^(-1/2)*tr(qC\(qtilde_prime_lambda*Cθ))

    qCinv = qC\I
    qC_inv_prime_lambda = - qC \ ((qtilde_prime_lambda*Cθ) * qCinv)

    dbilinearform = z0 -> (local AA = dgλ(z0, λ) - m_prime_lambda;
                    local BB = (g(z0, λ) - m);
                    2*AA'*(qC \ BB) + BB' * qC_inv_prime_lambda * BB
                    )
    dEXPR3 = z0 -> -(n-p+k)/2 * bilinearform(z0)^(-(n-p+k+2)/2) * dbilinearform(z0)

    #product rule
    main = z0 -> (EXPR1(z0) * EXPR2(z0) * EXPR3(z0))
    main_deriv = z0 -> (dEXPR1(z0)*EXPR2(z0)*EXPR3(z0) + dEXPR2(z0)*EXPR1(z0)*EXPR3(z0)+dEXPR3(z0)*EXPR1(z0)*EXPR2(z0))
    #sub-functions (dependence on z0)
    #return (EXPR1, dEXPR1)
    #return (EXPR2, dEXPR2)
    #return (EXPR3, dEXPR3)
    return (main, main_deriv)

    #sub-constants
    #return (qtilde, qtilde_prime_lambda)
    #return (βhat, βhat_prime_lambda)
    #return (gλz, dgλz)
    #return (qCinv, qC_inv_prime_lambda)
    #return (m, m_prime_lambda)
end
"""
Compute deriative of p(z0|theta, lambda, z) with respect to z0
"""
function partial_z0(θ, λ, setting, g = boxCox, dg = boxCoxPrime, dg2 = boxCoxPrime2)
    (s, s0, X, X0, z, n, p, k, Eθ, Σθ, Bθ, choleskyΣθ, choleskyXΣX, Dθ, Hθ, Cθ, Eθ_prime,Eθ_prime2, Σθ_prime, Σθ_prime2, Bθ_prime, Bθ_prime2) = func(θ, setting)

    βhat = choleskyXΣX\(X'*(choleskyΣθ\g(z, λ)))
    qtilde = (expr = g(z, λ)-X*βhat; expr'*(choleskyΣθ\expr))
    m = Bθ*(choleskyΣθ\g(z, λ)) + Hθ*βhat

    jac = z0 -> abs(reduce(*, map(x -> dg(x, λ), z0)))

    djac = z0 ->(local n = length(z0);
                reshape(
                (local dd = map(x->dg(x, λ), z0);
                local p = abs(reduce(*, dd));
                (f = i -> (dd[i] !=0 ? p/dd[i]*sign(dd[i])*(-1)^(((sign(dd[mod(i+1, n)+1])+1)/2)+1)*
                dg2(z0[i], λ) : 0) ;
                map(f, collect(1:n)))), n, 1))

    qC = qtilde*Cθ
    detqC = det(qC)

    ccc = gamma((n-p+k)/2)/gamma((n-p)/2)/pi^(k/2)/detqC^(1/2) #constant

    expr = z0 -> (g(z0, λ) .- m)'*(qC\(g(z0, λ) .- m))

    dexpr = z0 -> (local n = length(z0);
                    2 * Diagonal(dg(z0, λ)) * (qC\(g(z0, λ) .- m));
                    )

    #expr2 = z0 -> (1 .+ (g(z0, λ) .- m)'*(qC\(g(z0, λ) .- m)))^(-(n-p+k)/2)
    #dexpr2 = z0 -> (-(n-p+k)/2) * (1 .+ (g(z0, λ) .- m)'*(qC\(g(z0, λ) .- m)))^(-(n-p+k+2)/2) * 2 * Diagonal(dg(z0, λ)) * (qC\(g(z0, λ) .- m))

    expr3 = z0 -> ccc*jac(z0)*(1 .+ (g(z0, λ) .- m)'*(qC\(g(z0, λ) .- m)))^(-(n-p+k)/2)
    dexpr3 = z0 -> ccc*(jac(z0)*(-(n-p+k)/2) * (1 .+ (g(z0, λ) .- m)'*(qC\(g(z0, λ) .- m)))^(-(n-p+k+2)/2) * 2 * Diagonal(dg(z0, λ)) * (qC\(g(z0, λ) .- m))
                    + djac(z0)*(1 .+ (g(z0, λ) .- m)'*(qC\(g(z0, λ) .- m)))^(-(n-p+k)/2))

    #return (z0 -> g(z0, λ), z0 -> [dg(z0, λ)[1] 0 ; 0 dg(z0, λ)[2]])
    #return (expr, dexpr)
    #return (jac, djac)
    return (expr3, dexpr3)
end

"""
Compute p(theta, lambda| z), possibly in addition to derivatives with respect to theta. Differs from paper in that Jacobian term computation is delayed
until the matrix of weights is formed.
"""
#function posterior_theta(θ::Float64, λ::Float64, pθ, dpθ, dpθ2, pλ, setting::setting{Array{Float64,2}, Array{Float64, 1}}, theta_params::Union{θ_params{Array{Float64, 2},
function posterior_theta(θ::Float64, λ::Float64, priorθ, priorλ, train::trainingData{A, B}, test::testingData{A}, transforms, theta_params::Union{θ_params{A,
    Cholesky{Float64, A}}, θ_param_derivs{A, Cholesky{Float64,A}}, Nothing} = nothing, type::String="Gaussian") where A<:Array{Float64, 2} where B<:Array{Float64, 1}
    if theta_params === nothing
        #println("WARNING: recompute theta_params in partial_theta")
        theta_params = funcθ(θ, train, test, type)
        #(Eθ, Σθ, Bθ, Dθ, Hθ, Cθ, Eθ_prime, Eθ_prime2, Σθ_prime, Σθ_prime2, Bθ_prime, Bθ_prime2, choleskyΣθ, choleskyXΣX) = getθ_Params(theta_params)
        #(s, s0, X, X0, z,  n, p, k) = getSettingParams(setting)
    end
    if type == "Gaussian" || type == "Turan"
        Eθ = theta_params.Eθ
        Σθ = theta_params.Σθ
        Bθ = theta_params.Bθ
        Dθ = theta_params.Dθ
        Hθ = theta_params.Hθ
        Cθ = theta_params.Cθ
        Σθ_inv_X = theta_params.Σθ_inv_X
        choleskyΣθ = theta_params.choleskyΣθ
        choleskyXΣX = theta_params.choleskyXΣX
    else
        throw(ArgumentError("Quadrature type undefined. Please enter \"Gaussian\" or \"Turan\" for last arg."))
    end
    if type == "Turan" #load higher derivatives
        Eθ_prime = theta_params.Eθ_prime
        Eθ_prime2 = theta_params.Eθ_prime2
        Σθ_prime = theta_params.Σθ_prime
        Σθ_prime2 = theta_params.Σθ_prime2
        Bθ_prime = theta_params.Bθ_prime
        Bθ_prime2 = theta_params.Bθ_prime2
        Dθ_prime = theta_params.Dθ_prime
        Dθ_prime2 = theta_params.Dθ_prime2
        Hθ_prime = theta_params.Hθ_prime
        Hθ_prime2 = theta_params.Hθ_prime2
        Cθ_prime = theta_params.Cθ_prime
        Cθ_prime2 = theta_params.Cθ_prime2
        tripleΣ = theta_params.tripleΣ
    end
    #unpack args
    pθ = priorθ.f; dpθ = priorθ.df; dpθ2 = priorθ.d2f; pλ = priorλ.f
    s = train.s; s0 = test.s0; X = train.X; X0 = test.X0; z = train.z; n = size(X, 1); p = size(X, 2); k = size(X0, 1)
    g = transforms.f; dg = transforms.df
    #compute auxiliary quantities
    gλz = g(z, λ)
    βhat = (X'*(choleskyΣθ\X))\(X'*(choleskyΣθ\gλz))
    qtilde = (expr = gλz-X*βhat; expr'*(choleskyΣθ\expr))
    m = Bθ*(choleskyΣθ\gλz) + Hθ*βhat
    Σθ_inv_X = Σθ\X
    meanvv = gλz - X*βhat

    EXPR1 = det(choleskyΣθ)^(-1/2)
    #println("det Sigma theta ^-1/2: ", EXPR1)
    EXPR2 = det(choleskyXΣX)^(-1/2)
    EXPR3 = qtilde^(-(n-p)/2)
    EXPR4 = pθ(θ)*pλ(λ)
    main = EXPR1*EXPR2*EXPR3*EXPR4

    if type=="Turan"
        #compute betahat_prime_theta
        expr_mid = X'*(choleskyΣθ\Σθ_prime)*(Σθ_inv_X)
        βhat_prime_theta = compute_betahat_prime_theta(choleskyΣθ, choleskyXΣX, expr_mid, Σθ_prime, X, gλz, Σθ_inv_X)

        #compute qtilde_prime_theta
        qtilde_prime_theta  = compute_qtilde_prime_theta(gλz, X, βhat, βhat_prime_theta, choleskyΣθ, Σθ_prime)

        trΣθqΣθ_prime = tr(choleskyΣθ\Σθ_prime) #precompute
        dEXPR1 = -0.5 * det(choleskyΣθ)^(-1/2)*trΣθqΣθ_prime
        XΣdΣΣX = Σθ_inv_X' * Σθ_prime * Σθ_inv_X  #precompute
        trexpr = tr(choleskyXΣX\(-XΣdΣΣX))#precompute
        dEXPR2 = -0.5*det(choleskyXΣX)^(-1/2)*trexpr
        dEXPR3 = -((n-p)/2)*qtilde^(-(n-p+2)/2)*qtilde_prime_theta
        dEXPR4 = dpθ(θ)*pλ(λ)

        dmain = (dEXPR1 * EXPR2 * EXPR3 * EXPR4 .+ dEXPR2*EXPR1*EXPR3*EXPR4 .+ dEXPR3*EXPR1*EXPR2*EXPR4 .+ dEXPR4*EXPR1*EXPR2*EXPR3)[1]

        #-====================== dmain2 (second derivative)===============================
        d2EXPR1 = 0.25 * EXPR1 * trΣθqΣθ_prime^2 -0.5 * EXPR1* tr(choleskyΣθ\Σθ_prime2 - choleskyΣθ\(Σθ_prime*(choleskyΣθ\Σθ_prime)))

        dQ = Y -> (choleskyΣθ\((Σθ_prime2*(choleskyΣθ\Y)))) - 2*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\(Y))))))
        dPinv = Y -> choleskyXΣX\(X'*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\(X*(choleskyXΣX\Y))))))

        d2EXPR2 = 0.25 * EXPR2 * trexpr^2 .+ 0.5*EXPR2*tr(dPinv(XΣdΣΣX) .+ choleskyXΣX\(X'*dQ(X)))

        Q = Y -> (choleskyΣθ\(Σθ_prime*(Σθ\(Y))))
        βhat_prime2_theta = compute_βhat_prime2_theta(choleskyXΣX, choleskyΣθ, expr_mid, X, Q, dQ, gλz)
        qtilde_prime2_theta = compute_qtilde_prime2_theta(choleskyΣθ, X, meanvv, Q, dQ, βhat_prime_theta, βhat_prime2_theta)
        d2EXPR3 = ((n-p)/2)*((n-p+2)/2)*qtilde^(-(n-p+4)/2)*qtilde_prime_theta^2 - (n-p)/2 * qtilde^(-(n-p+2)/2)*qtilde_prime2_theta
        d2EXPR4 = dpθ2(θ)*pλ(λ)
        d2main = (d2EXPR1*EXPR2*EXPR3*EXPR4 .+ d2EXPR2*EXPR1*EXPR3*EXPR4 .+ d2EXPR3*EXPR1*EXPR2*EXPR4 .+ d2EXPR4*EXPR1*EXPR2*EXPR3
                    .+ 2*(dEXPR1*dEXPR2*EXPR3*EXPR4 .+ dEXPR1*EXPR2*dEXPR3*EXPR4 .+ dEXPR1*EXPR2*EXPR3*dEXPR4 .+ EXPR1*dEXPR2*dEXPR3*EXPR4
                    .+ EXPR1*dEXPR2*EXPR3*dEXPR4 .+ EXPR1*EXPR2*dEXPR3*dEXPR4))[1]
        if false #debug
            println("===============================================================")
            println("theta");println(θ);println("betahat");println(βhat);println("glambdaz");println(gλz);println("qtilde")
            println(qtilde);println("m");println(m);println("expr_mid");println(expr_mid);println("βhat_prime2_theta")
            println(βhat_prime2_theta);println(" qtilde_prime2_theta");println(qtilde_prime2_theta)
            println("pθ(θ)"); println(pθ(θ)); println("pλ(λ)");println(pλ(λ));println("EXPR1");println(EXPR1);println("EXPR2");println(EXPR2)
            println("EXPR3"); println(EXPR3);println("EXPR4");println(EXPR4);println("dEXPR1"); println(dEXPR1)
            println("dEXPR2"); println(dEXPR2);println("dEXPR3");println(dEXPR3);println("dEXPR4")
            println(dEXPR4); println("d2EXPR1");println(d2EXPR1);println("d2EXPR2");println(d2EXPR2)
            println("d2EXPR3");println(d2EXPR3);println("d2EXPR4");println(d2EXPR4);println("main")
            println(main);println("dmain");println(dmain);println("d2main");println(d2main)
            println("===============================================================")
        end
        return (main, dmain, d2main)
    elseif type == "Gaussian"
        return main
    else
        throw(ArgumentError("Quadrature type undefined. Please use \"Gaussian\" or \"Turan\" for last arg"))
    end

    #return (qtilde_prime_theta, qtilde_prime2_theta)
    #return (choleskyXΣX\(expr_mid)*(choleskyXΣX\(X'*(Σθ\gλz))), βhatEPXR1+βhatEPXR2+βhatEPXR3+βhatEPXR4)
    #return (- (choleskyXΣX\(X'*(Σθ\(Σθ_prime*(Σθ\gλz))))), βhatEPXR5 + βhatEPXR6)
    #return (Q(gλz), dQ(gλz))
    #return (βhat_prime_theta, βhat_prime2_theta)
    #return (vec(choleskyXΣX\X), vec(dPinv(X)))
    #return (Σθ\(Σθ_prime*(Σθ\X)), dQ(X))
    #return (dEXPR2, d2EXPR2)
    #return (main, dmain)

end

"""
Compute derivative of p(theta, lambda| z) with respect to theta
"""
function posterior_lambda(θ, λ, pθ, pλ, dpλ, setting)
    (s, s0, X, X0, z, n, p, k, Eθ, Σθ, Bθ, choleskyΣθ, choleskyXΣX, Dθ, Hθ, Cθ, Eθ_prime,Eθ_prime2, Σθ_prime, Σθ_prime2, Bθ_prime, Bθ_prime2) = func(θ, setting)

    g = boxCox
    dg = boxCoxPrime
    dgλ = boxCoxPrime_lambda
    dgλx = boxCoxMixed_lambda_z

    βhat = (X'*(choleskyΣθ\X))\(X'*(choleskyΣθ\g(z, λ)))
    qtilde = (expr = g(z, λ)-X*βhat; expr'*(choleskyΣθ\expr))

    #compute βhat_prime_lambda
    dgλz = dgλ(z, λ)
    gλz = g(z, λ)
    βhat_prime_lambda = choleskyXΣX\(X'*(choleskyΣθ\dgλz))

    #compute qtilde_prime_lambda
    AA = dgλz - X*βhat_prime_lambda
    BB = gλz - X*βhat
    qtilde_prime_lambda = 2*AA'*(choleskyΣθ\BB)

    jac = z0 -> (abs(reduce( *, map(x -> dg(x, λ), z0))))
    jacz = jac(z)

    #compute derivative of Jacobian
    function djac(z0)
        dgλx_vec = z0 -> map(zi -> dgλx(zi, λ), z0)
        dg_vec = z0 -> map(zi -> dg(zi, λ), z0)
        local prod = z0 -> reduce(*, dg_vec(z0))
        local p = prod(z0)
        p != 0 ? (AAA = dgλx_vec(z0);
        BBB = dg_vec(z0);
        ff = i -> p * AAA[i]/BBB[i] ;
        ee = reduce(+, map(ff, collect(1:length(z0))));
        dg(2, 2) < 0 ? (-1)^(length(z0))*ee : ee) : 0
    end

    ccc = det(choleskyΣθ)^(-1/2)*det(choleskyXΣX)^(-1/2)*pθ(θ)
    EXPR1 = qtilde^(-(n-p)/2)
    EXPR2 = jacz^(1-p/n)
    EXPR3 = pλ(λ)
    dEXPR1 = -(n-p)/2*qtilde^(-(n-p+2)/2)*qtilde_prime_lambda
    dEXPR2 = (1-p/n)*jacz^(-p/n)*djac(z)
    dEXPR3 = dpλ(λ)
    main = ccc * EXPR1*EXPR2*EXPR3
    dmain = ccc*(dEXPR1*EXPR2*EXPR3 + dEXPR2*EXPR1*EXPR3 + dEXPR3*EXPR1*EXPR2)
    return (main, dmain)
end

"""
Use Taylor's Theorem with Remainder to check
validity of computed derivative. More specifically, check
that error of linear approximation decays like O(h^2)

last: smallest power of exp^(-x) we wish to compute in scheme
"""
function checkDerivative(f, df, x0, first = 3, last = 12, num = 10)
    f0 = f(x0)
    df0 = df(x0)
    if size(x0, 2)>1
        dx = rand(size(x0, 1), size(x0, 2))
    else
        dx = rand(size(x0, 1))
    end
    h = collect(first:(last-first)/num:last)
    for i=1:length(h)
        h[i] = exp(-h[i])
    end
    A = zeros(length(h))
    for i = 1:length(h)
        #println(x0)
        #println(h[i]*dx)
        fi = f(x0 .+ h[i]*dx)

        if false #debug
            println("dx: ", dx)
            println("fi: ", fi)
            println("f0: ", f0)
            println("df0: ", df0)
            println("increment", h[i]*dx)
        end
        try
            A[i] = norm((fi .- f0) .- df0' * (h[i] * dx))
        catch DimensionMismatch #we catch the case when f: R^1 -> R^n, in which case  df0'*dx will yield an error
            #println("caught in check deriv")
            A[i] = norm((fi .- f0) .- df0 .* (h[i] * dx))
        end
    end
    #println(A)
    #println(h)

    return (log.(h), log.(A))
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg