text:
module Cairo

import Base.Sys

# For libcairo
using Cairo_jll
# For libpangocairo
using Pango_jll
# For libgobject
using Glib_jll

# Deprecate old library variables
Base.@deprecate_binding _jl_libcairo      Cairo.libcairo   false
Base.@deprecate_binding _jl_libgobject    Cairo.libgobject false
Base.@deprecate_binding _jl_libpango      Cairo.libpango   false
Base.@deprecate_binding _jl_libpangocairo Cairo.libpango   false

using Colors

import Graphics
using Graphics: BoundingBox, GraphicsContext, GraphicsDevice
import Graphics: arc, clip, clip_preserve, close_path, creategc, device_to_user!, device_to_user_distance!, fill_preserve, height, line_to, move_to, new_path, new_sub_path, paint, rectangle, rel_line_to, rel_move_to, reset_clip, restore, rotate, save, scale, set_dash, set_line_width, set_source, set_source_rgb, set_source_rgba, stroke, stroke_preserve, stroke_transformed, stroke_transformed_preserve, textwidth, translate, user_to_device!, user_to_device_distance!, width, circle, reset_transform
import Base: copy, fill

libcairo_version = VersionNumber(unsafe_string(
      ccall((:cairo_version_string,Cairo.libcairo),Cstring,()) ))
libpango_version = VersionNumber(unsafe_string(
      ccall((:pango_version_string,Cairo.libpango),Cstring,()) ))
if !Sys.iswindows()
    libpangocairo_version = VersionNumber(unsafe_string(
          ccall((:pango_version_string,Cairo.libpangocairo),Cstring,()) ))
    libgobject_version = VersionNumber(
          unsafe_load(cglobal((:glib_major_version, Cairo.libgobject), Cuint)),
          unsafe_load(cglobal((:glib_minor_version, Cairo.libgobject), Cuint)),
          unsafe_load(cglobal((:glib_micro_version, Cairo.libgobject), Cuint)))
end

import Base.show

include("constants.jl")

export
    # drawing surface and context types
    CairoSurface, CairoContext, CairoPattern,

    # surface constructors
    CairoRGBSurface, CairoPDFSurface, CairoEPSSurface, CairoXlibSurface,
    CairoARGBSurface, CairoSVGSurface, CairoImageSurface, CairoQuartzSurface,
    CairoWin32Surface, CairoScriptSurface, CairoRecordingSurface,
    CairoPSSurface, surface_create_similar,

    # surface and context management
    finish, destroy, status, get_source,
    creategc, save, restore, show_page, width, height,

    # pattern
    pattern_create_radial, pattern_create_linear,
    pattern_add_color_stop_rgb, pattern_add_color_stop_rgba,
    pattern_set_filter, pattern_set_extend, pattern_get_surface,

    # mesh patterns (version > 1.12)
    CairoPatternMesh,
    mesh_pattern_begin_patch, mesh_pattern_end_patch,
    mesh_pattern_move_to, mesh_pattern_line_to,
    mesh_pattern_curve_to,
    mesh_pattern_set_corner_color_rgb,
    mesh_pattern_set_corner_color_rgba,

    # drawing attribute manipulation
    set_antialias, get_antialias,
    set_fill_type, set_line_width, set_dash,
    set_source_rgb, set_source_rgba, set_source_surface, set_line_type,
    set_line_cap, set_line_join,
    set_operator, get_operator, set_source,
    CairoMatrix,

    # coordinate systems
    reset_transform, rotate, scale, translate, user_to_device!,
    device_to_user!, user_to_device_distance!, device_to_user_distance!,
    get_matrix, set_matrix,

    # clipping
    clip, clip_preserve, reset_clip,

    # fill, stroke, path, and shape commands
    fill, fill_preserve, new_path, new_sub_path, close_path, paint, paint_with_alpha, stroke,
    stroke_preserve, stroke_transformed, stroke_transformed_preserve,
    move_to, line_to, rel_line_to, rel_move_to,
    rectangle, circle, arc, arc_negative,
    curve_to, rel_curve_to,
    path_extents,

    # path copy
    copy_path, copy_path_flat, convert_cairo_path_data,

    # other path operations
    get_current_point, has_current_point,

    # text
    text,
    update_layout, show_layout, get_layout_size, layout_text,
    set_text, set_latex,
    set_font_face, set_font_size, select_font_face,
    textwidth, textheight, text_extents,
    TeXLexer, tex2pango, show_text, text_path,

    # images
    write_to_png, image, read_from_png,

    # push+pop group
    push_group, pop_group

## head docstring following

"
Surfaces, the canvas you are painting on

    CairoSurface, CairoRGBSurface, CairoPDFSurface, CairoEPSSurface,
    CairoXlibSurface, CairoARGBSurface, CairoSVGSurface,
    CairoImageSurface, CairoQuartzSurface,
    CairoWin32Surface, CairoScriptSurface, CairoRecordingSurface,
    CairoPSSurface,

Context, the handle to coordinate transformation, paint+Color

    CairoContext, finish, destroy, status, get_source,
    creategc, save, restore, show_page, width, height

Path creation API

    new_path, new_sub_path, close_path,
    move_to, line_to, rel_line_to, rel_move_to,
    rectangle, circle, arc, arc_negative,
    curve_to, rel_curve_to,
    path_extents,

Stroking and painting API

    fill, fill_preserve, paint, paint_with_alpha, stroke,
    stroke_preserve, stroke_transformed, stroke_transformed_preserve

    CairoContext, finish, destroy, status, get_source,
    creategc, save, restore, show_page, width, height

"
Cairo

function write_to_stream_callback(s::IO, buf::Ptr{UInt8}, len::UInt32)
    n = unsafe_write(s,buf,len)
    Int32((n == len) ? 0 : 11)
end

get_stream_callback(::Type{T}) where T = @cfunction write_to_stream_callback  Int32 (Ref{T}, Ptr{UInt8}, UInt32)

function read_from_stream_callback(s::IO, buf::Ptr{UInt8}, len::UInt32)
    # wrap the provided buf into a julia Array
    b1 = unsafe_wrap(Array,buf,len)

    # read from stream
    nb = readbytes!(s,b1,len)

    # provide a return status
    (nb == len) ? STATUS_SUCCESS : STATUS_READ_ERROR
end

get_readstream_callback(::Type{T}) where T = @cfunction read_from_stream_callback Int32 (Ref{T}, Ptr{UInt8}, UInt32)

abstract type CairoSurface{T<:Union{UInt32,RGB24,ARGB32}} <: GraphicsDevice end

mutable struct CairoSurfaceBase{T<:Union{UInt32,RGB24,ARGB32}} <: CairoSurface{T}
    ptr::Ptr{Nothing}
    width::Float64
    height::Float64

    function CairoSurface{T}(ptr::Ptr{Nothing}, w, h) where {T}
        self = new{T}(ptr, w, h)
        finalizer(destroy, self)
        self
    end
    function CairoSurface{T}(ptr::Ptr{Nothing}) where {T}
        ccall(
          (:cairo_surface_reference,libcairo),
          Ptr{Nothing}, (Ptr{Nothing}, ), ptr)
        self = new{T}(ptr)
        finalizer(destroy, self)
        self
    end
end

mutable struct CairoSurfaceImage{T<:Union{UInt32,RGB24,ARGB32}} <: CairoSurface{T}
    ptr::Ptr{Nothing}
    width::Float64
    height::Float64
    data::Matrix{T}

    function CairoSurface{T}(ptr::Ptr{Nothing}, w, h, data::Matrix{T}) where {T}
        self = new{T}(ptr, w, h, data)
        finalizer(destroy, self)
        self
    end
end

mutable struct CairoSurfaceIOStream{T<:Union{UInt32,RGB24,ARGB32}} <: CairoSurface{T}
    ptr::Ptr{Nothing}
    width::Float64
    height::Float64
    stream::IO

    function CairoSurface{T}(ptr::Ptr{Nothing}, w, h, stream::IO) where {T}
        self = new{T}(ptr, w, h, stream)
        finalizer(destroy, self)
        self
    end
end
function get_stream_ptr(surf::T) where {T<:CairoSurfaceIOStream}
    GC.@preserve surf begin
        return unsafe_load(Ptr{Ptr{Nothing}}(pointer_from_objref(surf) + fieldoffset(T, 4)))
    end
end

CairoSurface(ptr, w, h) = CairoSurface{UInt32}(ptr, w, h)
CairoSurface(ptr, w, h, data) = CairoSurface{eltype(data)}(ptr, w, h, data)
CairoSurface(ptr, w, h, stream::IO) = CairoSurface{UInt32}(ptr, w, h, stream)
CairoSurface(ptr) = CairoSurface{UInt32}(ptr)

width(surface::CairoSurface) = surface.width
height(surface::CairoSurface) = surface.height

function destroy(surface::CairoSurface)
    if surface.ptr == C_NULL
        return
    end
    ccall((:cairo_surface_destroy,libcairo), Nothing, (Ptr{Nothing},), surface.ptr)
    surface.ptr = C_NULL
    nothing
end

# function resize(surface::CairoSurface, w, h)
#     if OS_NAME == :Linux
#         CairoXlibSurfaceSetSize(surface.ptr, w, h)
#     elseif OS_NAME == :Darwin
#     elseif OS_NAME == :Windows
#     else
#         error("Unsupported operating system")
#     end
#     surface.width = w
#     surface.height = h
# end

for name in (:finish,:flush,:mark_dirty)
    @eval begin
        $name(surface::CairoSurface) =
            ccall(($(string("cairo_surface_",name)),libcairo),
                  Nothing, (Ptr{Nothing},), surface.ptr)
    end
end

function status(surface::CairoSurface)
    ccall((:cairo_surface_status,libcairo),
          Int32, (Ptr{Nothing},), surface.ptr)
end

function CairoImageSurface(w::Real, h::Real, format::Integer)
    ptr = ccall((:cairo_image_surface_create,libcairo),
                Ptr{Nothing}, (Int32,Int32,Int32), format, w, h)
    CairoSurface(ptr, w, h)
end

CairoRGBSurface(w::Real, h::Real) = CairoImageSurface(w, h, FORMAT_RGB24)
CairoARGBSurface(w::Real, h::Real) = CairoImageSurface(w, h, FORMAT_ARGB32)
CairoARGBSurface(img) = CairoImageSurface(img, FORMAT_ARGB32)
CairoRGBSurface(img) = CairoImageSurface(img, FORMAT_RGB24)

function CairoImageSurface(img::Array{UInt32,2}, format::Integer; flipxy::Bool = true)
    if flipxy
        img = permutedims(img, (2,1))
    end
    w,h = size(img)
    stride = format_stride_for_width(format, w)
    @assert stride == 4w
    ptr = ccall((:cairo_image_surface_create_for_data,libcairo),
                Ptr{Nothing}, (Ptr{Nothing},Int32,Int32,Int32,Int32),
                img, format, w, h, stride)
    CairoSurface(ptr, w, h, img)
end

function CairoImageSurface(img::Matrix{T}) where {T<:Union{RGB24,ARGB32}}
    w,h = size(img)
    stride = format_stride_for_width(format(T), w)
    @assert stride == 4w
    ptr = ccall((:cairo_image_surface_create_for_data,libcairo),
                Ptr{Nothing}, (Ptr{Nothing},Int32,Int32,Int32,Int32),
                img, format(T), w, h, stride)
    CairoSurface(ptr, w, h, img)
end

format(::Type{RGB24}) = FORMAT_RGB24
format(::Type{ARGB32}) = FORMAT_ARGB32
format(surf::CairoSurface{T}) where {T<:Union{RGB24,ARGB32}} = T

## PDF ##

function CairoPDFSurface(stream::T, w::Real, h::Real) where {T<:IO}
    callback = get_stream_callback(T)
    surf = CairoSurface(C_NULL, w, h, stream)
    GC.@preserve surf begin
        surf.ptr = ccall((:cairo_pdf_surface_create_for_stream,libcairo), Ptr{Nothing},
                         (Ptr{Nothing}, Ptr{Nothing}, Float64, Float64), callback,
                         get_stream_ptr(surf), w, h)
    end
    return surf
end

function CairoPDFSurface(filename::AbstractString, w_pts::Real, h_pts::Real)
    ptr = ccall((:cairo_pdf_surface_create,libcairo), Ptr{Nothing},
                (Ptr{UInt8},Float64,Float64), String(filename), w_pts, h_pts)
    CairoSurface(ptr, w_pts, h_pts)
end

## EPS ##

function CairoEPSSurface(stream::T, w::Real, h::Real) where {T<:IO}
    callback = get_stream_callback(T)
    surf = CairoSurface(C_NULL, w, h, stream)
    GC.@preserve surf begin
        surf.ptr = ccall((:cairo_ps_surface_create_for_stream,libcairo), Ptr{Nothing},
                         (Ptr{Nothing}, Ptr{Nothing}, Float64, Float64),
                         callback, get_stream_ptr(surf), w, h)
        ccall((:cairo_ps_surface_set_eps,libcairo), Nothing,
              (Ptr{Nothing},Int32), surf.ptr, 1)
    end
    return surf
end

function CairoEPSSurface(filename::AbstractString, w_pts::Real, h_pts::Real)
    ptr = ccall((:cairo_ps_surface_create,libcairo), Ptr{Nothing},
                (Ptr{UInt8},Float64,Float64), String(filename), w_pts, h_pts)
    ccall((:cairo_ps_surface_set_eps,libcairo), Nothing,
          (Ptr{Nothing},Int32), ptr, 1)
    CairoSurface(ptr, w_pts, h_pts)
end

## PS ##

function CairoPSSurface(stream::T, w::Real, h::Real) where {T<:IO}
    callback = get_stream_callback(T)
    surf = CairoSurface(C_NULL, w, h, stream)
    GC.@preserve surf begin
        surf.ptr = ccall((:cairo_ps_surface_create_for_stream,libcairo), Ptr{Nothing},
                         (Ptr{Nothing}, Ptr{Nothing}, Float64, Float64),
                         callback, get_stream_ptr(surf), w, h)
        ccall((:cairo_ps_surface_set_eps,libcairo), Nothing,
              (Ptr{Nothing},Int32), surf.ptr, 0)
    end
    return surf
end

function CairoPSSurface(filename::AbstractString, w_pts::Real, h_pts::Real)
    ptr = ccall((:cairo_ps_surface_create,libcairo), Ptr{Nothing},
                (Ptr{UInt8},Float64,Float64), String(filename), w_pts, h_pts)
    ccall((:cairo_ps_surface_set_eps,libcairo), Nothing,
          (Ptr{Nothing},Int32), ptr, 0)
    CairoSurface(ptr, w_pts, h_pts)
end

## Xlib ##

function CairoXlibSurface(display, drawable, visual, w, h)
    ptr = ccall((:cairo_xlib_surface_create,libcairo), Ptr{Nothing},
                (Ptr{Nothing}, Int, Ptr{Nothing}, Int32, Int32),
                display, drawable, visual, w, h)
    CairoSurface(ptr, w, h)
end

CairoXlibSurfaceSetSize(surface, w, h) =
    ccall((:cairo_xlib_surface_set_size,libcairo), Nothing,
          (Ptr{Nothing}, Int32, Int32),
          surface, w, h)

## Quartz ##
function CairoQuartzSurface(context, w, h)
    ptr = ccall((:cairo_quartz_surface_create_for_cg_context,libcairo),
          Ptr{Nothing}, (Ptr{Nothing}, UInt32, UInt32), context, w, h)
    CairoSurface(ptr, w, h)
end

## Win32 ##

function CairoWin32Surface(hdc,w,h)
    ptr = ccall((:cairo_win32_surface_create, libcairo),
                Ptr{Nothing}, (Ptr{Nothing},), hdc)
    CairoSurface(ptr,w,h)
end

## SVG ##

function CairoSVGSurface(stream::T, w::Real, h::Real) where {T<:IO}
    callback = get_stream_callback(T)
    surf = CairoSurface(C_NULL, w, h, stream)
    GC.@preserve surf begin
        surf.ptr = ccall((:cairo_svg_surface_create_for_stream,libcairo), Ptr{Nothing},
                         (Ptr{Nothing}, Ptr{Nothing}, Float64, Float64), callback,
                         get_stream_ptr(surf), w, h)
    end
    return surf
end

function CairoSVGSurface(filename::AbstractString, w::Real, h::Real)
    ptr = ccall((:cairo_svg_surface_create,libcairo), Ptr{Nothing},
                (Ptr{UInt8},Float64,Float64), String(filename), w, h)
    CairoSurface(ptr, w, h)
end

## PNG ##

function read_from_png(filename::AbstractString)
    ptr = ccall((:cairo_image_surface_create_from_png,libcairo),
                Ptr{Nothing}, (Ptr{UInt8},), String(filename))
    w = ccall((:cairo_image_surface_get_width,libcairo),
              Int32, (Ptr{Nothing},), ptr)
    h = ccall((:cairo_image_surface_get_height,libcairo),
              Int32, (Ptr{Nothing},), ptr)
    CairoSurface(ptr, w, h)
end

function write_to_png(surface::CairoSurface, stream::T) where {T<:IO}
    callback = get_stream_callback(T)
    ccall((:cairo_surface_write_to_png_stream,libcairo), status_t,
          (Ptr{UInt8},Ptr{Nothing},Any), surface.ptr, callback, stream)
end

function write_to_png(surface::CairoSurface, filename::AbstractString)
    ccall((:cairo_surface_write_to_png,libcairo), status_t,
          (Ptr{UInt8},Ptr{UInt8}), surface.ptr, String(filename))
end

show(io::IO, ::MIME"image/png", surface::CairoSurface) =
   write_to_png(surface, io)

function read_from_png(stream::T) where {T<:IO}
    callback = get_readstream_callback(T)
    ptr = ccall((:cairo_image_surface_create_from_png_stream, Cairo.libcairo),
                Ptr{Nothing}, (Ptr{Nothing},Ref{IO}), callback, stream)
    w = ccall((:cairo_image_surface_get_width,Cairo.libcairo),
              Int32, (Ptr{Nothing},), ptr)
    h = ccall((:cairo_image_surface_get_height,Cairo.libcairo),
              Int32, (Ptr{Nothing},), ptr)
    Cairo.CairoSurface(ptr, w, h)
end

## Generic ##

function surface_create_similar(s::CairoSurface, w = width(s), h = height(s))
    ptr = ccall((:cairo_surface_create_similar,libcairo), Ptr{Nothing},
                (Ptr{Nothing}, Int32, Int32, Int32),
                s.ptr, CONTENT_COLOR_ALPHA, w, h)
    CairoSurface(ptr, w, h)
end

# Utilities

function format_stride_for_width(format::Integer, width::Integer)
    ccall((:cairo_format_stride_for_width,libcairo), Int32,
          (Int32,Int32), format, width)
end

## Scripting (only target IO stream)

mutable struct CairoScript <: GraphicsDevice
    ptr::Ptr{Nothing}
    stream::IO

    function CairoScript(stream::T) where {T<:IO}
        callback = get_stream_callback(T)
        self = new(C_NULL, stream)
        GC.@preserve self begin
            self.ptr = ccall((:cairo_script_create_for_stream,libcairo), Ptr{Nothing},
                             (Ptr{Nothing}, Ptr{Nothing}), callback, get_stream_ptr(self))
        end
        finalizer(destroy, self)
        self
    end
end
function get_stream_ptr(s::T) where {T<:CairoScript}
    GC.@preserve s begin
        return unsafe_load(Ptr{Ptr{Nothing}}(pointer_from_objref(s) + fieldoffset(T, 2)))
    end
end

function destroy(s::CairoScript)
    if s.ptr == C_NULL
        return
    end
    ccall((:cairo_device_destroy,libcairo), Nothing, (Ptr{Nothing},), s.ptr)
    s.ptr = C_NULL
    nothing
end

function CairoScriptSurface(stream::IO, w::Real, h::Real)
    s = CairoScript(stream)
    ptr = ccall((:cairo_script_surface_create,libcairo), Ptr{Nothing},
                (Ptr{Nothing},Int32,Float64,Float64),s.ptr ,CONTENT_COLOR_ALPHA, w, h)
    CairoSurface(ptr, w, h, stream)
end

function CairoScriptSurface(stream::IO, sc::CairoSurface)
    s = CairoScript(stream)
    ptr = ccall((:cairo_script_surface_create_for_target,libcairo), Ptr{Nothing},
                (Ptr{Nothing},Ptr{Nothing}),s.ptr, sc.ptr)
    CairoSurface(ptr, sc.width, sc.height)
end

mutable struct CairoRectangle
    x0::Float64
    y0::Float64
    x1::Float64
    y1::Float64
end

CairoRectangle() = CairoRectangle(0.0, 0.0, 0.0, 0.0)

function CairoRecordingSurface(content::Int32,extents::CairoRectangle)
    ptr = ccall((:cairo_recording_surface_create,libcairo), Ptr{Nothing},
                (Int32,Ptr{Nothing}),content, Ref(extents))
    CairoSurface(ptr)
end
function CairoRecordingSurface(content::Int32)
    ptr = ccall((:cairo_recording_surface_create,libcairo), Ptr{Nothing},
                (Int32,Ptr{Nothing}),content, C_NULL)
    CairoSurface(ptr)
end

CairoRecordingSurface() = CairoRecordingSurface(CONTENT_COLOR_ALPHA)

function script_from_recording_surface(s::CairoScript,r::CairoSurface)
    ccall((:cairo_script_from_recording_surface,libcairo), Int32,
                (Ptr{Nothing},Ptr{Nothing}),s.ptr, r.ptr)
end
# -----------------------------------------------------------------------------

mutable struct CairoContext <: GraphicsContext
    ptr::Ptr{Nothing}
    surface::CairoSurface
    layout::Ptr{Nothing} # cache PangoLayout

    function CairoContext(surface::CairoSurface)
        ptr = ccall((:cairo_create,libcairo),
                    Ptr{Nothing}, (Ptr{Nothing},), surface.ptr)
        layout = ccall((:pango_cairo_create_layout,libpangocairo),
                       Ptr{Nothing}, (Ptr{Nothing},), ptr)
        self = new(ptr, surface, layout)
        finalizer(destroy, self)
        self
    end
    function CairoContext(ptr::Ptr{Nothing})
        ccall((:cairo_reference,libcairo),
                   Ptr{Nothing}, (Ptr{Nothing},), ptr)
        surface_p = ccall((:cairo_get_target,libcairo),
                   Ptr{Nothing}, (Ptr{Nothing},), ptr)
        surface = CairoSurface(surface_p)
        layout = ccall((:pango_cairo_create_layout,libpangocairo),
                  Ptr{Nothing}, (Ptr{Nothing},), ptr)
        self = new(ptr,surface,layout)
        finalizer(destroy, self)
        self
    end

end

creategc(s::CairoSurface) = CairoContext(s)

function destroy(ctx::CairoContext)
    if ctx.ptr == C_NULL
        return
    end
    ccall((:g_object_unref,libgobject), Nothing, (Ptr{Nothing},), ctx.layout)
    _destroy(ctx)
    ctx.ptr = C_NULL
    nothing
end

 width(ctx::CairoContext) =  width(ctx.surface)
height(ctx::CairoContext) = height(ctx.surface)

function copy(ctx::CairoContext)
    surf = surface_create_similar(ctx.surface)
    c = creategc(surf)
    set_source_surface(c, ctx.surface)
    paint(c)
    set_matrix(c, get_matrix(ctx))
    c
end

# Copy a rectangular region
function copy(ctx::CairoContext, bb::BoundingBox)
    w = width(bb)
    h = height(bb)
    surf = surface_create_similar(ctx.surface, ceil(Int,w), ceil(Int,h))
    c = creategc(surf)
    set_source_surface(c, ctx.surface, -bb.xmin, -bb.ymin)
    rectangle(c, 0, 0, w, h)
    fill(c)
    set_matrix(c, get_matrix(ctx))
    c
end

for (NAME, FUNCTION) in Any[(:_destroy, :cairo_destroy),
                         (:save, :cairo_save),
                         (:restore, :cairo_restore),
                         (:show_page, :cairo_show_page),
                         (:clip, :cairo_clip),
                         (:clip_preserve, :cairo_clip_preserve),
                         (:reset_clip, :cairo_reset_clip),
                         (:reset_transform, :cairo_identity_matrix),
                         (:fill, :cairo_fill),
                         (:fill_preserve, :cairo_fill_preserve),
                         (:new_path, :cairo_new_path),
                         (:new_sub_path, :cairo_new_sub_path),
                         (:close_path, :cairo_close_path),
                         (:paint, :cairo_paint),
                         (:stroke_transformed, :cairo_stroke),
                         (:stroke_transformed_preserve, :cairo_stroke_preserve)]
    @eval begin
        $NAME(ctx::CairoContext) =
            ccall(($(Expr(:quote,FUNCTION)),libcairo),
                  Nothing, (Ptr{Nothing},), ctx.ptr)
    end
end

function stroke(ctx::CairoContext)
    save(ctx)
    # use uniform scale for stroking
    reset_transform(ctx)
    ccall((:cairo_stroke, libcairo), Nothing, (Ptr{Nothing},), ctx.ptr)
    restore(ctx)
end

function stroke_preserve(ctx::CairoContext)
    save(ctx)
    reset_transform(ctx)
    ccall((:cairo_stroke_preserve, libcairo), Nothing, (Ptr{Nothing},), ctx.ptr)
    restore(ctx)
end

function paint_with_alpha(ctx::CairoContext, a)
    ccall((:cairo_paint_with_alpha, libcairo),
          Nothing, (Ptr{Nothing}, Float64), ctx.ptr, a)
end

function get_operator(ctx::CairoContext)
    Int(ccall((:cairo_get_operator,libcairo), Int32, (Ptr{Nothing},), ctx.ptr))
end

for (NAME, FUNCTION) in Any[(:set_fill_type, :cairo_set_fill_rule),
                         (:set_operator, :cairo_set_operator),
                         (:set_line_cap, :cairo_set_line_cap),
                         (:set_line_join, :cairo_set_line_join)]
    @eval begin
        $NAME(ctx::CairoContext, i0::Integer) =
            ccall(($(Expr(:quote,FUNCTION)),libcairo),
                  Nothing, (Ptr{Nothing},Int32), ctx.ptr, i0)
    end
end

for (NAME, FUNCTION) in Any[(:set_line_width, :cairo_set_line_width),
                         (:rotate, :cairo_rotate),
                         (:set_font_size, :cairo_set_font_size)]
    @eval begin
        $NAME(ctx::CairoContext, d0::Real) =
            ccall(($(Expr(:quote,FUNCTION)),libcairo),
                  Nothing, (Ptr{Nothing},Float64), ctx.ptr, d0)
    end
end

for (NAME, FUNCTION) in Any[(:line_to, :cairo_line_to),
                         (:move_to, :cairo_move_to),
                         (:rel_line_to, :cairo_rel_line_to),
                         (:rel_move_to, :cairo_rel_move_to),
                         (:scale, :cairo_scale),
                         (:translate, :cairo_translate)]
    @eval begin
        $NAME(ctx::CairoContext, d0::Real, d1::Real) =
            ccall(($(Expr(:quote,FUNCTION)),libcairo),
                  Nothing, (Ptr{Nothing},Float64,Float64), ctx.ptr, d0, d1)
    end
end

for (NAME, FUNCTION) in Any[(:curve_to, :cairo_curve_to),
                         (:rel_curve_to, :cairo_rel_curve_to)]
    @eval begin
        $NAME(ctx::CairoContext, d0::Real, d1::Real, d2::Real, d3::Real, d4::Real, d5::Real) =
            ccall(($(Expr(:quote,FUNCTION)),libcairo),
                  Nothing, (Ptr{Nothing},Float64,Float64,Float64,Float64,Float64,Float64), ctx.ptr, d0, d1, d2, d3, d4, d5)
    end
end

for (NAME, FUNCTION) in Any[(:arc, :cairo_arc),
                         (:arc_negative, :cairo_arc_negative)]
    @eval begin
        $NAME(ctx::CairoContext, xc::Real, yc::Real, radius::Real, angle1::Real, angle2::Real) =
            ccall(($(Expr(:quote,FUNCTION)),libcairo),
                  Nothing, (Ptr{Nothing},Float64,Float64,Float64,Float64,Float64),
                  ctx.ptr, xc, yc, radius, angle1, angle2)
    end
end

set_source_rgb(ctx::CairoContext, r::Real, g::Real, b::Real) =
    ccall((:cairo_set_source_rgb,libcairo),
          Nothing, (Ptr{Nothing},Float64,Float64,Float64), ctx.ptr, r, g, b)

set_source_rgba(ctx::CairoContext, r::Real, g::Real, b::Real, a::Real) =
    ccall((:cairo_set_source_rgba,libcairo), Nothing,
          (Ptr{Nothing},Float64,Float64,Float64,Float64),
          ctx.ptr, r, g, b, a)

function set_source(ctx::CairoContext, c::Color)
    rgb = convert(RGB, c)
    set_source_rgb(ctx, rgb.r, rgb.g, rgb.b)
end

function set_source(ctx::CairoContext, ac::TransparentColor)
    rgba = convert(RGBA, ac)
    set_source_rgba(ctx, rgba.r, rgba.g, rgba.b, rgba.alpha)
end

set_source(dest::CairoContext, src::CairoContext) = set_source_surface(dest, src.surface)

set_source(dest::CairoContext, src::CairoSurface) = set_source_surface(dest, src)

rectangle(ctx::CairoContext, x::Real, y::Real, w::Real, h::Real) =
    ccall((:cairo_rectangle,libcairo), Nothing,
          (Ptr{Nothing},Float64,Float64,Float64,Float64),
          ctx.ptr, x, y, w, h)

function set_dash(ctx::CairoContext, dashes::Vector{Float64}, offset::Real = 0.0)
    ccall((:cairo_set_dash,libcairo), Nothing,
          (Ptr{Nothing},Ptr{Float64},Int32,Float64), ctx.ptr, dashes, length(dashes), offset)
end

function set_source_surface(ctx::CairoContext, s::CairoSurface, x::Real = 0.0, y::Real = 0.0)
    ccall((:cairo_set_source_surface,libcairo), Nothing,
          (Ptr{Nothing},Ptr{Nothing},Float64,Float64), ctx.ptr, s.ptr, x, y)
end

function set_source(ctx::CairoContext, s::CairoSurface, x::Real, y::Real)
    set_source_surface(ctx, s, x, y)
end

# cairo_path data and functions

mutable struct CairoPath_t
    status::Cairo.status_t
    data::Ptr{Float64}
    num_data::UInt32
end

mutable struct CairoPath <: GraphicsDevice
    ptr::Ptr{CairoPath_t}

    function CairoPath(ptr::Ptr{Nothing})
        self = new(ptr)
        finalizer(destroy, self)
        self
    end
end

# Abstract, contains type (moveto,lineto,curveto,closepath) and points
mutable struct CairoPathEntry
    element_type::UInt32
    points::Array{Float64,1}
end

function destroy(path::CairoPath)
    if path.ptr == C_NULL
        return
    end
    ccall((:cairo_path_destroy,libcairo), Nothing, (Ptr{Nothing},), path.ptr)
    path.ptr = C_NULL
    nothing
end

function copy_path(ctx::CairoContext)
    ptr = ccall((:cairo_copy_path, libcairo),
                    Ptr{Nothing}, (Ptr{Nothing},),ctx.ptr)
    path = CairoPath(ptr)
    finalizer(destroy, path)
    path
end

function copy_path_flat(ctx::CairoContext)
    ptr = ccall((:cairo_copy_path_flat, libcairo),
                    Ptr{Nothing}, (Ptr{Nothing},),ctx.ptr)
    path = CairoPath(ptr)
    finalizer(destroy, path)
    path
end

function convert_cairo_path_data(p::CairoPath)
    c = unsafe_load(p.ptr)

    # The original data (pointed by c.data) is an array of Unions. We
    # define here by Float64 (most data is) and reinterpret in the header.

    path_data = CairoPathEntry[]
    @static if VERSION >= v"0.7"
        c_data = unsafe_wrap(Array, c.data, (Int(c.num_data*2), 1), own=false)
    else
        c_data = unsafe_wrap(Array, c.data, (Int(c.num_data*2), 1), false)
    end

    data_index = 1
    while data_index <= ((c.num_data)*2)

        # read header (reinterpret a Float64 to UInt64 and split to UInt32 x 2)
        element_length = reinterpret(UInt64,c_data[data_index]) >> 32
        element_type = reinterpret(UInt64,c_data[data_index]) & 0xffffffff

        # copy points x,y
        points = Vector{Float64}(undef, (element_length - 1) * 2)
        for i=1:(element_length-1)*2
            points[i] = c_data[data_index+i+1]
        end

        g = CairoPathEntry(element_type,points)
        push!(path_data,g)

        # goto next element
        data_index += (element_length*2)

    end
    path_data
end

# other path operations

function get_current_point(ctx::CairoContext)

    x = Ref{Cdouble}(0)
    y = Ref{Cdouble}(0)
    ccall((:cairo_get_current_point, libcairo),
            Nothing, (Ptr{Nothing},Ref{Cdouble},Ref{Cdouble}),ctx.ptr,x,y)

    x[],y[]
end

function has_current_point(ctx::CairoContext)
    Bool(ccall((:cairo_has_current_point, libcairo),
            Cint, (Ptr{Nothing},),ctx.ptr))
end

# user<->device coordinate translation

for (fname,cname) in ((:user_to_device!,:cairo_user_to_device),
                      (:device_to_user!,:cairo_device_to_user),
                      (:user_to_device_distance!,:cairo_user_to_device_distance),
                      (:device_to_user_distance!,:cairo_device_to_user_distance))
    @eval begin
        function ($fname)(ctx::CairoContext, p::Vector{Float64})
            ccall(($(Expr(:quote,cname)),libcairo),
                  Nothing, (Ptr{Nothing}, Ptr{Float64}, Ptr{Float64}),
                  ctx.ptr, pointer(p,1), pointer(p,2))
            p
        end
    end
end

function image(ctx::CairoContext, s::CairoSurface, x, y, w, h)
    rectangle(ctx, x, y, w, h)
    save(ctx)
    translate(ctx, x, y)
    scale(ctx, w/s.width, h/s.height)
    set_source_surface(ctx, s, 0, 0)
    if abs(w) > s.width && abs(h) > s.height
        # use NEAREST filter when stretching an image
        # it's usually better to see pixels than a blurry mess when viewing
        # a small image
        p = get_source(ctx)
        pattern_set_filter(p, FILTER_NEAREST)
    end
    fill(ctx)
    restore(ctx)
end

image(ctx::CairoContext, img::Array{UInt32,2}, x, y, w, h) =
    image(ctx, CairoRGBSurface(img), x, y, w, h)

function push_group(ctx::CairoContext)
    if ctx.ptr == C_NULL
        return
    end
    ccall((:cairo_push_group, libcairo), Nothing, (Ptr{Nothing},),ctx.ptr)
    nothing
end

function pop_group(ctx::CairoContext)
    if ctx.ptr == C_NULL
        return
    end
    ptr = ccall((:cairo_pop_group, libcairo), Ptr{Nothing}, (Ptr{Nothing},),ctx.ptr)
    pattern = CairoPattern(ptr)
    finalizer(destroy, pattern)
    pattern
end

# -----------------------------------------------------------------------------

mutable struct CairoPattern
    ptr::Ptr{Nothing}
end

function CairoPattern(s::CairoSurface)
    ptr = ccall((:cairo_pattern_create_for_surface, libcairo),
                    Ptr{Nothing}, (Ptr{Nothing},), s.ptr)
    # Ideally we'd check the status, but at least for certain releases of the library
    # the return value seems not to be set properly (random values are returned)
#     status = ccall((:cairo_pattern_status, libcairo),
#                     Cint, (Ptr{Nothing},), s.ptr)
#     if status != 0
#         error("Error creating Cairo pattern: ", bytestring(
#               ccall((:cairo_status_to_string, libcairo),
#                     Ptr{UInt8}, (Cint,), status)))
#     end
    pattern = CairoPattern(ptr)
    finalizer(destroy, pattern)
    pattern
end

set_source(dest::CairoContext, src::CairoPattern) =
    ccall((:cairo_set_source, libcairo),
          Nothing, (Ptr{Nothing}, Ptr{Nothing}), dest.ptr, src.ptr)

function get_source(ctx::CairoContext)
    CairoPattern(ccall((:cairo_get_source,libcairo),
                       Ptr{Nothing}, (Ptr{Nothing},), ctx.ptr))
end

function pattern_set_filter(p::CairoPattern, f)
    ccall((:cairo_pattern_set_filter,libcairo), Nothing,
          (Ptr{Nothing},Int32), p.ptr, f)
end

function pattern_set_extend(p::CairoPattern, val)
    ccall((:cairo_pattern_set_extend,libcairo), Nothing,
          (Ptr{Nothing},Int32), p.ptr, val)
end

function pattern_create_radial(cx0::Real, cy0::Real, radius0::Real, cx1::Real, cy1::Real, radius1::Real)
    ptr = ccall((:cairo_pattern_create_radial, libcairo),
                    Ptr{Nothing}, (Float64,Float64,Float64,Float64,Float64,Float64),cx0,cy0,radius0,cx1,cy1,radius1)
    pattern = CairoPattern(ptr)
    finalizer(destroy, pattern)
    pattern
end

function pattern_create_linear(x0::Real, y0::Real, x1::Real, y1::Real)
    ptr = ccall((:cairo_pattern_create_linear, libcairo),
                    Ptr{Nothing}, (Float64,Float64,Float64,Float64),x0,y0,x1,y1)
    pattern = CairoPattern(ptr)
    finalizer(destroy, pattern)
    pattern
end

function pattern_add_color_stop_rgb(pat::CairoPattern, offset::Real, red::Real, green::Real, blue::Real)
    ccall((:cairo_pattern_add_color_stop_rgb, libcairo),
                    Nothing, (Ptr{Nothing},Float64,Float64,Float64,Float64),pat.ptr,offset,red,green,blue)
end

function pattern_add_color_stop_rgba(pat::CairoPattern, offset::Real, red::Real, green::Real, blue::Real, alpha::Real)
    ccall((:cairo_pattern_add_color_stop_rgba, libcairo),
                    Nothing, (Ptr{Nothing},Float64,Float64,Float64,Float64,Float64),pat.ptr,offset,red,green,blue,alpha)
end

function pattern_get_surface(pat::CairoPattern)
    ptrref = Ref{Ptr{Nothing}}()
    status = ccall((:cairo_pattern_get_surface, libcairo), Cint,
                   (Ptr{Nothing}, Ref{Ptr{Nothing}}), pat.ptr, ptrref)
    if status == STATUS_PATTERN_TYPE_MISMATCH
        error("Cannot get surface from a non-surface pattern.")
    end
    ptr = ptrref.x

    ccall((:cairo_surface_reference, libcairo), Ptr{Nothing}, (Ptr{Nothing},), ptr)
    typ = ccall((:cairo_surface_get_type, libcairo), Cint, (Ptr{Nothing},), ptr)

    w = 0.0
    h = 0.0
    if typ == CAIRO_SURFACE_TYPE_IMAGE
        w = ccall((:cairo_image_surface_get_width, libcairo),
                  Int32, (Ptr{Nothing},), ptr)
        h = ccall((:cairo_image_surface_get_height, libcairo),
                  Int32, (Ptr{Nothing},), ptr)
    end
    return CairoSurface(ptr, w, h)
end

function destroy(pat::CairoPattern)
    if pat.ptr == C_NULL
        return
    end
    ccall((:cairo_pattern_destroy,libcairo), Nothing, (Ptr{Nothing},), pat.ptr)
    pat.ptr = C_NULL
    nothing
end

# mesh pattern

# create mesh pattern
function CairoPatternMesh()
    ptr = ccall((:cairo_pattern_create_mesh, libcairo),
                    Ptr{Nothing}, ())
    pattern = CairoPattern(ptr)
    #status = ccall((:cairo_pattern_status, libcairo),
    #                Cint, (Ptr{Nothing},), pattern.ptr)
    #if status != 0
    #    error("Error creating Cairo pattern: ", bytestring(
    #          ccall((:cairo_status_to_string, libcairo),
    #                Ptr{Uint8}, (Cint,), status)))
    #end
    finalizer(destroy, pattern)
    pattern
end

#for (NAME, FUNCTION) in Any[(:set_line_width, :cairo_set_line_width),
#                         (:rotate, :cairo_rotate),
#                         (:set_font_size, :cairo_set_font_size)]
#    @eval begin
#        $NAME(ctx::CairoContext, d0::Real) =
#            ccall(($(Expr(:quote,FUNCTION)),libcairo),
#                  Nothing, (Ptr{Nothing},Float64), ctx.ptr, d0)
#    end
#end

for (NAME, FUNCTION) in Any[(:mesh_pattern_begin_patch, :cairo_mesh_pattern_begin_patch),
                         (:mesh_pattern_end_patch, :cairo_mesh_pattern_end_patch)]
    @eval begin
        $NAME(pattern::CairoPattern) =
            ccall(($(Expr(:quote,FUNCTION)),libcairo),
                  Nothing, (Ptr{Nothing},), pattern.ptr)
    end
end

for (NAME, FUNCTION) in Any[(:mesh_pattern_line_to, :cairo_mesh_pattern_line_to),
                         (:mesh_pattern_move_to, :cairo_mesh_pattern_move_to)]
    @eval begin
        $NAME(pattern::CairoPattern, d0::Real, d1::Real) =
            ccall(($(Expr(:quote,FUNCTION)),libcairo),
                  Nothing, (Ptr{Nothing},Float64,Float64), pattern.ptr, d0, d1)
    end
end

for (NAME, FUNCTION) in Any[(:mesh_pattern_curve_to, :cairo_mesh_pattern_curve_to)]

    @eval begin
        $NAME(pattern::CairoPattern, d0::Real, d1::Real, d2::Real, d3::Real, d4::Real, d5::Real) =
            ccall(($(Expr(:quote,FUNCTION)),libcairo),
                  Nothing, (Ptr{Nothing},Float64,Float64,Float64,Float64,Float64,Float64), pattern.ptr, d0, d1, d2, d3, d4, d5)
    end
end

function mesh_pattern_set_corner_color_rgb(pat::CairoPattern, corner_num::Real, red::Real, green::Real, blue::Real)
    ccall((:cairo_mesh_pattern_set_corner_color_rgb, libcairo),
                    Nothing, (Ptr{Nothing},Int32,Float64,Float64,Float64),pat.ptr,corner_num,red,green,blue)
end

function mesh_pattern_set_corner_color_rgba(pat::CairoPattern, corner_num::Real, red::Real, green::Real, blue::Real, alpha::Real)
    ccall((:cairo_mesh_pattern_set_corner_color_rgba, libcairo),
                    Nothing, (Ptr{Nothing},Int32,Float64,Float64,Float64,Float64),pat.ptr,corner_num,red,green,blue,alpha)
end

# ----

set_antialias(ctx::CairoContext, a) =
    ccall((:cairo_set_antialias,libcairo), Nothing,
          (Ptr{Nothing},Cint), ctx.ptr, a)

get_antialias(ctx::CairoContext) =
    ccall((:cairo_get_antialias,libcairo), Cint,
          (Ptr{Nothing},), ctx.ptr)

# -----------------------------------------------------------------------------

struct CairoMatrix
    xx::Float64
    yx::Float64
    xy::Float64
    yy::Float64
    x0::Float64
    y0::Float64
end

CairoMatrix() = CairoMatrix(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

function get_matrix(ctx::CairoContext)
    m = [CairoMatrix()]
    ccall((:cairo_get_matrix, libcairo), Nothing, (Ptr{Nothing}, Ptr{Nothing}), ctx.ptr, m)
    m[1]
end

function set_matrix(ctx::CairoContext, m::CairoMatrix)
    ccall((:cairo_set_matrix, libcairo), Nothing, (Ptr{Nothing}, Ptr{Nothing}), ctx.ptr, [m])
end

function set_matrix(p::CairoPattern, m::CairoMatrix)
    ccall((:cairo_pattern_set_matrix, libcairo), Nothing, (Ptr{Nothing}, Ptr{Nothing}), p.ptr, [m])
end

# -----------------------------------------------------------------------------
function set_line_type(ctx::CairoContext, nick::AbstractString)
    if nick == "solid"
        dash = Float64[]
    elseif nick == "dotted" || nick == "dot"
        dash = [1.,3.]
    elseif nick == "dotdashed"
        dash = [1.,3.,4.,4.]
    elseif nick == "longdashed"
        dash = [6.,6.]
    elseif nick == "shortdashed" || nick == "dash" || nick == "dashed"
        dash = [4.,4.]
    elseif nick == "dotdotdashed"
        dash = [1.,3.,1.,3.,4.,4.]
    elseif nick == "dotdotdotdashed"
        dash = [1.,3.,1.,3.,1.,3.,4.,4.]
    else
        error("unknown line type ", nick)
    end
    set_dash(ctx, dash)
end

# -----------------------------------------------------------------------------
# text commands

function set_font_face(ctx::CairoContext, str::AbstractString)
    fontdesc = ccall((:pango_font_description_from_string,libpango),
                     Ptr{Nothing}, (Ptr{UInt8},), String(str))
    ccall((:pango_layout_set_font_description,libpango), Nothing,
          (Ptr{Nothing},Ptr{Nothing}), ctx.layout, fontdesc)
    ccall((:pango_font_description_free,libpango), Nothing,
          (Ptr{Nothing},), fontdesc)
end

function set_text(ctx::CairoContext, text::AbstractString, markup::Bool = false)
    if markup
        ccall((:pango_layout_set_markup,libpango), Nothing,
            (Ptr{Nothing},Ptr{UInt8},Int32), ctx.layout, String(text), -1)
    else
        # clear attributes (from previous set markups)
        ccall((:pango_layout_set_attributes,libpango), Nothing,
            (Ptr{Nothing},Ptr{Nothing}), ctx.layout, C_NULL)

        ccall((:pango_layout_set_text,libpango), Nothing,
            (Ptr{Nothing},Ptr{UInt8},Int32), ctx.layout, String(text), -1)
    end
    text
end

function get_layout_size(ctx::CairoContext)
    w = Vector{Int32}(undef, 2)
    ccall((:pango_layout_get_pixel_size,libpango), Nothing,
          (Ptr{Nothing},Ptr{Int32},Ptr{Int32}), ctx.layout, pointer(w,1), pointer(w,2))
    w
end

function update_layout(ctx::CairoContext)
    ccall((:pango_cairo_update_layout,libpangocairo), Nothing,
          (Ptr{Nothing},Ptr{Nothing}), ctx.ptr, ctx.layout)
end

function show_layout(ctx::CairoContext)
    ccall((:pango_cairo_show_layout,libpangocairo), Nothing,
          (Ptr{Nothing},Ptr{Nothing}), ctx.ptr, ctx.layout)
end

text_extents(ctx::CairoContext,value::AbstractString) =
    text_extents!(ctx,value, Matrix{Float64}(undef, 6, 1))

function text_extents!(ctx::CairoContext,value::AbstractString,extents)
    ccall((:cairo_text_extents, libcairo),
          Nothing, (Ptr{Nothing}, Ptr{UInt8}, Ptr{Float64}),
          ctx.ptr, String(value), extents)
    extents
end

function path_extents(ctx::CairoContext)
    dx1 = Cdouble[0]
    dx2 = Cdouble[0]
    dy1 = Cdouble[0]
    dy2 = Cdouble[0]

    ccall((:cairo_path_extents, libcairo),
          Nothing, (Ptr{Nothing}, Ptr{Cdouble}, Ptr{Cdouble},
          Ptr{Cdouble}, Ptr{Cdouble}),
          ctx.ptr, dx1, dy1, dx2, dy2)

    return(dx1[1],dy1[1],dx2[1],dy2[1])
end

function show_text(ctx::CairoContext,value::AbstractString)
    ccall((:cairo_show_text, libcairo),
          Nothing, (Ptr{Nothing}, Ptr{UInt8}),
          ctx.ptr, String(value))
end

function text_path(ctx::CairoContext,value::AbstractString)
    ccall((:cairo_text_path, libcairo),
          Nothing, (Ptr{Nothing}, Ptr{UInt8}),
          ctx.ptr, String(value))
end

function select_font_face(ctx::CairoContext,family::AbstractString,slant,weight)
    ccall((:cairo_select_font_face, libcairo),
          Nothing, (Ptr{Nothing}, Ptr{UInt8},
                 font_slant_t, font_weight_t),
          ctx.ptr, String(family),
          slant, weight)
end

function align2offset(a::AbstractString)
    if     a == "center" return 0.5
    elseif a == "left"   return 0.0
    elseif a == "right"  return 1.0
    elseif a == "top"    return 0.0
    elseif a == "bottom" return 1.0
    end
    @assert false
end

function text(ctx::CairoContext, x::Real, y::Real, str::AbstractString;
              halign::AbstractString = "left", valign::AbstractString = "bottom", angle::Real = 0, markup::Bool=false)
    move_to(ctx, x, y)
    save(ctx)
    reset_transform(ctx)
    rotate(ctx, -angle*pi/180.)

    set_text(ctx, str, markup)
    update_layout(ctx)

    extents = get_layout_size(ctx)
    dxrel = -align2offset(halign)
    dyrel = align2offset(valign)
    rel_move_to(ctx, dxrel*extents[1], -dyrel*extents[2])

    show_layout(ctx)
    restore(ctx)
    w, h = Graphics.device_to_user(ctx, extents[1], extents[2])
    BoundingBox(x+dxrel*w, x+(dxrel+1)*w, y-dyrel*h, y+(1-dyrel)*h)
end

function textwidth(ctx::CairoContext, str::AbstractString, markup::Bool = false)
    set_text(ctx, str, markup)
    extents = get_layout_size(ctx)
    extents[1]
end

function textheight(ctx::CairoContext, str::AbstractString, markup::Bool = false)
    set_text(ctx, str, markup)
    extents = get_layout_size(ctx)
    extents[2]
end

set_latex(ctx::CairoContext, str::AbstractString, fontsize::Real) = set_text(ctx, tex2pango(str, fontsize), true)

mutable struct TeXLexer
    str::String
    len::Int
    pos::Int
    token_stack::Array{String,1}

    function TeXLexer(str::AbstractString)
        s = String(str)
        new(s, lastindex(s), 1, String[])
    end
end

function get_token(self::TeXLexer)
    if self.pos > self.len
        return nothing
    end

    if length(self.token_stack) > 0
        return pop!(self.token_stack)
    end

    str = self.str[self.pos:end]
    re_control_sequence = r"^\\[a-zA-Z]+[ ]?|^\\[^a-zA-Z][ ]?"
    m = match(re_control_sequence, str)
    if m !== nothing
        token = m.match
        self.pos = self.pos + sizeof(token)
        # consume trailing space
        if length(token) > 2 && token[end] == ' '
            token = token[1:end-1]
        end
    else
        @static if VERSION >= v"0.7"
            token, self.pos = iterate(self.str, self.pos)
        else
            token, self.pos = next(self.str, self.pos)
        end
        token = string(token)
    end

    return token
end

function put_token(self::TeXLexer, token)
    push!(self.token_stack, token)
end

function peek(self::TeXLexer)
    token = get_token(self)
    put_token(self, token)
    return token
end

function map_text_token(token::AbstractString)
    if haskey(_text_token_dict, token)
        return _text_token_dict[token]
    else
        return get(_common_token_dict, token, token)
    end
end

function map_math_token(token::AbstractString)
    if haskey(_math_token_dict, token)
        return _math_token_dict[token]
    else
        return get(_common_token_dict, token, token)
    end
end

function math_group(lexer::TeXLexer)
    output = ""
    bracketmode = false
    while true
        token = get_token(lexer)
        if token == nothing
            break
        end

        if token == "{"
            bracketmode = true
        elseif token == "}"
            break
        else
            output = string(output, map_math_token(token))
            if !bracketmode
                break
            end
        end
    end
    return output
end

#font_code = [ "\\f0", "\\f1", "\\f2", "\\f3" ]

function tex2pango(str::AbstractString, fontsize::Real)
    output = ""
    mathmode = true
    font_stack = Any[]
    font = 1
    script_size = fontsize/1.618034

    lexer = TeXLexer(str)
    while true
        token = get_token(lexer)
        if token == nothing
            break
        end

        more_output = ""

        if token == "\$"
#            mathmode = !mathmode
            more_output = "\$"
        elseif token == "{"
            push!(font_stack, font)
        elseif token == "}"
            old_font = pop!(font_stack)
            if old_font != font
                font = old_font
#                more_output = font_code[font]
            end
        elseif token == "\\rm"
            font = 1
#            more_output = font_code[font]
        elseif token == "\\it"
            font = 2
#            more_output = font_code[font]
        elseif token == "\\bf"
            font = 3
#            more_output = font_code[font]
        elseif !mathmode
            more_output = map_text_token(token)
        elseif token == "_"
            more_output = string("<sub><span font=\"$script_size\">", math_group(lexer), "</span></sub>")
            #if peek(lexer) == "^"
            #    more_output = string("\\mk", more_output, "\\rt")
            #end
        elseif token == "^"
            more_output = string("<sup><span font=\"$script_size\">", math_group(lexer), "</span></sup>")
            #if peek(lexer) == "_"
            #    more_output = string("\\mk", more_output, "\\rt")
            #end
        else
            more_output = map_math_token(token)
        end

        output = string(output, more_output)
    end
    return output
end

@deprecate text(ctx::CairoContext,x::Real,y::Real,str::AbstractString,fontsize::Real,halign::AbstractString,valign,angle)    text(ctx,x,y,set_latex(ctx,str,fontsize),halign=halign,valign=valign,angle=angle,markup=true)
@deprecate layout_text(ctx::CairoContext, str::AbstractString, fontsize::Real)       set_latex(ctx, str, fontsize)
@deprecate textwidth(ctx::CairoContext, str::AbstractString, fontsize::Real)         textwidth(ctx, tex2pango(str, fontsize), true)
@deprecate textheight(ctx::CairoContext, str::AbstractString, fontsize::Real)        textheight(ctx, tex2pango(str, fontsize), true)
@deprecate cairo_write_to_ios_callback(s::Ptr{Nothing}, buf::Ptr{UInt8}, len::UInt32)   write_to_ios_callback(s, buf, len)
@deprecate cairo_write_to_stream_callback(s::IO, buf::Ptr{UInt8}, len::UInt32)       write_to_stream_callback(s, buf, len)
@deprecate text_extents(ctx::CairoContext,value::AbstractString,extents) text_extents!(ctx,value,extents)

end  # module


text:
# Copyright 2017 Stephen L. Smith and Frank Imeson
# Modified copyright 2020 A. Kudryavtsev and M. Khachay
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Three different parameter settings for the GTSP solver
default -- moderate runtime.  ~ 10-20 seconds for 100 sets
fast -- very quick, lower quality solution
slow -- slower, but potentially higher quality solution
"""
function parameter_settings(num_vertices, num_sets, sets, problem_instance, args)
    args = Dict(args)
    mode = get(args, :mode, "default")
    ############ default setting #####################
    if mode == "default"
        num_iterations = get(args, :num_iterations, 60) * num_sets
        param = Dict(
            :cold_trials => get(args, :trials, 5),
            :warm_trials => get(args, :restarts, 3),
            :max_time => get(args, :max_time, 360),
            :init_tour => get(args, :init_tour, "rand"),
            :prob_reopt => get(args, :reopt, 1.0),
            :accept_percentage => 0.05,
            :prob_accept => 10.0 / num_iterations,
            :num_iterations => num_iterations,
            :latest_improvement => num_iterations / 2,
            :first_improvement => num_iterations / 4,
            :max_removals => min(100, max(round(Int64, 0.3 * num_sets), 1)),
            :insertions => ["randpdf", "cheapest"],
        )

        ################## very_fast  ##########################
    elseif mode == "fast"

        num_iterations = get(args, :num_iterations, 60) * num_sets
        param = Dict(
            :cold_trials => get(args, :trials, 3),
            :warm_trials => get(args, :restarts, 2),
            :max_time => get(args, :max_time, 300),
            :init_tour => get(args, :init_tour, "insertion"),
            :prob_reopt => get(args, :reopt, 0.2),
            :accept_percentage => 0.05,
            :prob_accept => 10.0 / num_iterations,
            :num_iterations => num_iterations,
            :latest_improvement => num_iterations / 4,
            :first_improvement => num_iterations / 6,
            :max_removals => min(20, max(round(Int64, 0.1 * num_sets), 1)),
            :insertions => ["randpdf"],
        )

        ################## attempt slow search  ##########################
    elseif mode == "slow"
        num_iterations = get(args, :num_iterations, 150) * num_sets
        param = Dict(
            :cold_trials => get(args, :trials, 10),
            :warm_trials => get(args, :restarts, 5),
            :max_time => get(args, :max_time, 1200),
            :init_tour => get(args, :init_tour, "rand"),
            :prob_reopt => get(args, :reopt, 1.0),
            :accept_percentage => 0.05,
            :prob_accept => 10.0 / num_iterations,
            :num_iterations => num_iterations,
            :latest_improvement => num_iterations / 3,
            :first_improvement => num_iterations / 6,
            :max_removals => max(round(Int64, 0.4 * num_sets), 1),
            :insertions => ["randpdf", "cheapest"],
        )

    else
        error("mode not recognized.  Use default, fast, or slow")
    end

    # insertion algs
    if get(args, :insertion_algs, "default") == "cheapest"
        param[:insertions] = ["randpdf", "cheapest"]
    elseif get(args, :insertion_algs, "default") == "randpdf"
        param[:insertions] = ["randpdf"]
    end

    # insertion powers
    if get(args, :insertion_algs, "default") == "classic"
        param[:insertion_powers] = [-10.0, 0.0, 10.0]
    else
        param[:insertion_powers] = [-10.0, -1.0, -0.5, 0.0, 0.5, 1.0, 10.0]
    end

    # removal algs and powers
    if get(args, :removal_algs, "default") == "classic"
        param[:removals] = ["worst"]
        param[:removal_powers] = [0.0, 10.0]
    else
        param[:removals] = ["distance", "worst", "segment"]
        param[:removal_powers] = [-0.5, 0.0, 0.5, 1.0, 10.0]
    end

    # parameters that are common for all modes
    param[:mode] = mode
    param[:problem_instance] = split(problem_instance, "/")[end]
    param[:num_sets] = num_sets
    param[:num_vertices] = num_vertices
    param[:output_file] = get(args, :output, "None")
    param[:print_output] = get(args, :verbose, 3)
    param[:epsilon] = get(args, :epsilon, 0.5)
    param[:noise] = get(args, :noise, "Add")
    param[:adaptive_iter] = 1
    param[:print_time] = 5
    param[:budget] = get(args, :budget, typemin(Int64))
    param[:timeout] = false
    param[:budget_met] = false
    param[:min_set] = min_set(sets)
    param[:min_removals] = (param[:max_removals] > 1 ? 2 : 1)
    param[:seed] = get(args, :seed, -1)
    print_params(param)

    return param
end

ext: jl

max_stars_repo_name: AndreiKud/PCGLNS

max_issues_repo_name: AndreiKud/PCGLNS

max_forks_repo_name: AndreiKud/PCGLNS


text:
# For use in geodesic routines

abstract _geodesic

type null_geodesic <: _geodesic
end

# The following functions are generally named after the associated C API
# functions, but without the pj prefix.

immutable ProjUV
    u::Cdouble
    v::Cdouble
end

"forward projection from Lat/Lon to X/Y (only supports 2 dimensions)"
function _fwd!(lonlat::Vector{Cdouble}, proj_ptr::Ptr{Void})
    xy = ccall((:pj_fwd, libproj), ProjUV, (ProjUV, Ptr{Void}), ProjUV(lonlat[1], lonlat[2]), proj_ptr)
    _errno() == 0 || error("forward projection error: $(_strerrno())")
    lonlat[1] = xy.u; lonlat[2] = xy.v
    lonlat
end

"Row-wise projection from Lat/Lon to X/Y (only supports 2 dimensions)"
function _fwd!(lonlat::Array{Cdouble,2}, proj_ptr::Ptr{Void})
    for i=1:size(lonlat,1)
        xy = ccall((:pj_fwd, libproj), ProjUV, (ProjUV, Ptr{Void}),
                   ProjUV(lonlat[i,1], lonlat[i,2]), proj_ptr)
        lonlat[i,1] = xy.u; lonlat[i,2] = xy.v
    end
    _errno() == 0 || error("forward projection error: $(_strerrno())")
    lonlat
end

"inverse projection from X/Y to Lat/Lon (only supports 2 dimensions)"
function _inv!(xy::Vector{Cdouble}, proj_ptr::Ptr{Void})
    lonlat = ccall((:pj_inv, libproj), ProjUV, (ProjUV, Ptr{Void}),
                   ProjUV(xy[1], xy[2]), proj_ptr)
    _errno() == 0 || error("inverse projection error: $(_strerrno())")
    xy[1] = lonlat.u; xy[2] = lonlat.v
    xy
end

"Row-wise projection from X/Y to Lat/Lon (only supports 2 dimensions)"
function _inv!(xy::Array{Cdouble,2}, proj_ptr::Ptr{Void})
    for i=1:size(xy,1)
        lonlat = ccall((:pj_inv, libproj), ProjUV, (ProjUV, Ptr{Void}),
                       ProjUV(xy[i,1], xy[i,2]), proj_ptr)
        xy[i,1] = lonlat.u; xy[i,2] = lonlat.v
    end
    _errno() == 0 || error("inverse projection error: $(_strerrno())")
    xy
end

function _init_plus(proj_string::String)
    proj_ptr = ccall((:pj_init_plus, libproj), Ptr{Void}, (Cstring,), proj_string)
    if proj_ptr == C_NULL
        # TODO: use context?
        error("Could not parse projection: \"$proj_string\": $(_strerrno())")
    end
    proj_ptr
end

"Free C datastructure associated with a projection. For internal use!"
function _free(proj_ptr::Ptr{Void})
    @assert proj_ptr != C_NULL
    ccall((:pj_free, libproj), Void, (Ptr{Void},), proj_ptr)
end

"Get human readable error string from proj.4 error code"
function _strerrno(code::Cint)
    unsafe_string(ccall((:pj_strerrno, libproj), Cstring, (Cint,), code))
end

"Get global errno string in human readable form"
function _strerrno()
    _strerrno(_errno())
end

"Get error number"
function _errno()
    unsafe_load(ccall((:pj_get_errno_ref, libproj), Ptr{Cint}, ()))
end

"Get projection definition string in the proj.4 plus format"
function _get_def(proj_ptr::Ptr{Void})
    @assert proj_ptr != C_NULL
    opts = 0 # Apparently obsolete argument, not used in current proj source
    unsafe_string(ccall((:pj_get_def, libproj), Cstring, (Ptr{Void}, Cint), proj_ptr, opts))
end

"Low level interface to libproj transform, C_NULL can be passed in for z, if it's 2-dimensional"
function _transform!(src_ptr::Ptr{Void}, dest_ptr::Ptr{Void}, point_count::Integer, point_stride::Integer,
                     x::Ptr{Cdouble}, y::Ptr{Cdouble}, z::Ptr{Cdouble})
    @assert src_ptr != C_NULL && dest_ptr != C_NULL
    err = ccall((:pj_transform, libproj), Cint, (Ptr{Void}, Ptr{Void}, Clong, Cint, Ptr{Cdouble}, Ptr{Cdouble},
                Ptr{Cdouble}), src_ptr, dest_ptr, point_count, point_stride, x, y, z)
    err != 0 && error("transform error: $(_strerrno(err))")
end

function _transform!(src_ptr::Ptr{Void}, dest_ptr::Ptr{Void}, position::Vector{Cdouble})
    @assert src_ptr != C_NULL && dest_ptr != C_NULL
    ndim = length(position)
    @assert ndim >= 2

    x = pointer(position)
    y = x + sizeof(Cdouble)
    z = (ndim == 2) ? Ptr{Cdouble}(C_NULL) : x + 2*sizeof(Cdouble)

    _transform!(src_ptr, dest_ptr, 1, 1, x, y, z)
    position
end
function _transform!(src_ptr::Ptr{Void}, dest_ptr::Ptr{Void}, position::Array{Cdouble,2})
    @assert src_ptr != C_NULL && dest_ptr != C_NULL
    npoints, ndim = size(position)
    @assert ndim >= 2

    x = pointer(position)
    y = x + sizeof(Cdouble)*npoints
    z = (ndim == 2) ? Ptr{Cdouble}(C_NULL) : x + 2*sizeof(Cdouble)*npoints

    _transform!(src_ptr, dest_ptr, npoints, 1, x, y, z)
    position
end

function _is_latlong(proj_ptr::Ptr{Void})
    @assert proj_ptr != C_NULL
    ccall((:pj_is_latlong, libproj), Cint, (Ptr{Void},), proj_ptr) != 0
end

function _is_geocent(proj_ptr::Ptr{Void})
    @assert proj_ptr != C_NULL
    ccall((:pj_is_geocent, libproj), Cint, (Ptr{Void},), proj_ptr) != 0
end

"Get a string describing the underlying version of libproj in use"
function _get_release()
    unsafe_string(ccall((:pj_get_release, libproj), Cstring, ()))
end

"""
Fetch the internal definition of the spheroid as a tuple (a, es), where

    a = major_axis
    es = eccentricity squared

"""
function _get_spheroid_defn(proj_ptr::Ptr{Void})
    major_axis = Ref{Cdouble}()
    eccentricity_squared = Ref{Cdouble}()
    ccall((:pj_get_spheroid_defn, libproj), Void, (Ptr{Void}, Ptr{Cdouble}, Ptr{Cdouble}),
          proj_ptr, major_axis, eccentricity_squared)
    major_axis[], eccentricity_squared[]
end

"Returns true if the two datums are identical, otherwise false."
function _compare_datums(p1_ptr::Ptr{Void}, p2_ptr::Ptr{Void})
    Bool(ccall((:pj_compare_datums, libproj), Cint, (Ptr{Void}, Ptr{Void}), p1_ptr, p2_ptr))
end

# Unused/untested

# """
# Return the lat/long coordinate system on which a projection is based.
# If the coordinate system passed in is latlong, a clone of the same will be returned.
# """
# function _latlong_from_proj(proj_ptr::Ptr{Void})
#     ccall((:pj_latlong_from_proj, libproj), Ptr{Void}, (Ptr{Void},), proj_ptr)
# end


text:
import Logging
import LoggingExtras
using LinearAlgebra

"""
    Sets a two-logger TeeLogger as the global_logger:
    - the current logger
    - a logger to a file with debug minimum level

# Arguments
- `filepath_p` : file to which the logs will be created (overwritten)
"""
function init_logging(filepath_p)
    """
    formats output to only print source_info for level > Warn
    """
    function logger_metafmt(level::Logging.LogLevel, _module, group, id, file, line)
        @nospecialize
        color = Logging.default_logcolor(level)
        prefix = string(level == Logging.Warn ? "Warning" : string(level), ':')
        suffix::String = ""
        level < Logging.Warn && return color, prefix, suffix
        _module !== nothing && (suffix *= "$(_module)")
        if file !== nothing
            _module !== nothing && (suffix *= " ")
            suffix *= Base.contractuser(file)::String
            if line !== nothing
                suffix *= ":$(isa(line, UnitRange) ? "$(first(line))-$(last(line))" : line)"
            end
        end
        !isempty(suffix) && (suffix = "@ " * suffix)
        return color, prefix, suffix
    end

    log_file = endswith(filepath_p, ".log") ? filepath_p : filepath_p * ".log"
    file_logger = Logging.ConsoleLogger(open(log_file, "w"), Logging.Debug, meta_formatter=logger_metafmt)
    #file_logger = LoggingExtras.MinLevelLogger(LoggingExtras.FileLogger(log_file), Logging.Debug) #print souce_info at each log line
    logger = LoggingExtras.TeeLogger(Logging.current_logger(), file_logger)

    Base.global_logger(logger)
end

"""
    redirect_to_file(f::Function, file_p::String, mode_p="w")

Execute function `f` while redirecting C and Julia level stdout to the file file_p.
Note that `file_p` is open in write mode by default.

# Arguments
- `f::Function` : the function to execute
- `file_p::String` : name of the file to print to
- `mode_p` : open mode of the file (defaults to "w")
"""
function redirect_to_file(f::Function, file_p::String, mode_p="w")
    open(file_p, mode_p) do file_l
        redirect_to_file(f, file_l)
    end
end
function redirect_to_file(f::Function, io::IO)
    Base.Libc.flush_cstdio()
    redirect_stdout(io) do
        f()
        Base.Libc.flush_cstdio()
    end
end

function pretty_print(io::IO, d::AbstractDict; spacing=1, sort_p::Bool=false)
    d_l = d
    if sort_p
        d_l = sort(d)
    end
    for (k,v) in d_l
        if typeof(v) <: AbstractDict
            str_k = "$(repr(k)) => "
            print(io, join(fill(" ", spacing)), str_k, "\n")
            pretty_print(io, v, spacing=spacing+1+length(str_k), sort_p=sort_p)
        else
            print(io, join(fill(" ", spacing)), k, " => ", v, "\n")
        end
    end
    nothing
end

function rm_files(dir::String, f)
    if isdir(dir)
        matches = filter(f, readdir(dir))
        paths = joinpath.(dir, matches)
        files_to_rm = filter(isfile, paths)
        foreach(rm, files_to_rm)
    end
end

function rm_non_prefixed(dir::String, prefix::String)
    rm_files(dir, !startswith(prefix))
end

function rm_prefixed(dir::String, prefix::String)
    rm_files(dir, startswith(prefix))
end

"""
Like Base.split() but returns a vector of String instead of a vector of SubString{String}
"""
function split_str(str::AbstractString, dlm; keepempty::Bool=true)
    return String.(split(str, dlm, keepempty=keepempty))
end

function split_with_space(str::String)
    result = String[];
    if length(str) > 0
        start_with_quote = startswith(str, "\"");
        buffer_quote = split(str, keepempty=false, "\"");
        i = 1;
        while i <= length(buffer_quote)
            if i > 1 || !start_with_quote
                str2 = buffer_quote[i];
                buffer_space = split(str2, keepempty=false);
                for str3 in buffer_space
                    push!(result, str3);
                end
                i += 1;
            end
            if i <= length(buffer_quote)
                push!(result, buffer_quote[i]);
                i += 1;
            end
        end
    end
    return result;
end

function is_different(a, b)::Bool
    return a != b
end
function is_different(a::Number, b::Number, eps=1e-09)::Bool
    return abs(a-b) >= eps
end
function in_bounds(value, lower, upper, eps=1e-09)::Bool
    return (lower - eps) <= value <= (upper + eps)
end

"""
normalize_zero : If True, [0 0 0] will be normalized to [1/3 1/3 1/3] else, to [Nan Nan Nan]
Note : [-1 1] is normalized to [-0.5 0.5]
"""
function normalize_values(d::AbstractDict, normalize_zero::Bool=true)
    values_l = collect(values(d))
    if normalize_zero && !is_different(norm(values_l), 0.)
        normalized_values = fill(1. / length(values_l), length(values_l))
    else
        normalized_values = normalize(collect(values(d)), 1)
    end

    return Dict(zip( keys(d), normalized_values ))
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# run the algorithm to test 2 layer boostNet

dataset = "set5" # "set2", "set5"
train_data, train_label, validation_data, validation_label, test_data, test_label =
  share(load_data(dataset))

model = boostNet(train_data, train_label', m=2, T=1000, B=10)

say()
p = predict(model, train_data)
say("train error: ", mean(p .!= train_label'))

p = predict(model, test_data)
say("test error: ", mean(p .!= test_label'))

say("love huifang")

max_stars_repo_name: innerlee/Halfspace

max_issues_repo_name: innerlee/Halfspace

max_forks_repo_name: innerlee/Halfspace


text:
using Random

struct RandomPartitioner <: AbstractPartitioner
    width::Int
    n::Int
    masks::Vector{Vector{Int}}

    function RandomPartitioner(width::Int, n::Int; seed::Union{Nothing,Int}=nothing)
        !isnothing(seed) && seed < 0 && throw(DomainError(seed, "Seed must be non-negative"))
        width < 0 && throw(DomainError(width, "Width must be non-negative"))
        n < 1 && throw(DomainError(n, "Partition may not be less then one"))

        rng = isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)

        ordering = randperm(rng, width)
        nmasks = ceil(Int, width / n)

        # Could I avoid allocating the masks?
        masks = collect(Iterators.partition(ordering, n))

        new(width, n, masks)
    end
end

function Base.iterate(itr::PartitionerIterator{RandomPartitioner,<:AbstractVector}, state=1)
    state > length(itr.partitioner.masks) && return nothing

    return itr.X[itr.partitioner.masks[state]], state + 1
end

function Base.iterate(itr::PartitionerIterator{RandomPartitioner,<:AbstractMatrix}, state=1)
    state > length(itr.partitioner.masks) && return nothing

    return itr.X[:, itr.partitioner.masks[state]], state + 1
end

function random_partitioning(X::T, n::Int; seed::Union{Nothing,Int}=nothing) where {T <: AbstractVector}
    return partition(RandomPartitioner(length(X), n; seed), X)
end

function random_partitioning(X::T, n::Int; seed::Union{Nothing,Int}=nothing) where {T <: AbstractMatrix}
    return partition(RandomPartitioner(size(X, 2), n; seed), X)
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
# Shared Memory (part of B.2)

export @LocalMemory

lmem_id = 0

macro LocalMemory(T, dims)
    # FIXME: generating a unique id in the macro is incorrect, as multiple parametrically typed
    #        functions will alias the id (and the size might be a parameter). but incrementing in
    #        the @generated function doesn't work, as it is supposed to be pure and identical
    #        invocations will erroneously share (and even cause multiple lmem globals).
    id = gensym("static_lmem")

    quote
        len = prod($(esc(dims)))
        ptr = emit_localmemory(Val($(QuoteNode(id))), $(esc(T)), Val(len))
        oneDeviceArray($(esc(dims)), ptr)
    end
end

# get a pointer to local memory, with known (static) or zero length (dynamic shared memory)
@generated function emit_localmemory(::Val{id}, ::Type{T}, ::Val{len}=Val(0)) where {id,T,len}
    Context() do ctx
        eltyp = convert(LLVMType, T; ctx)
        T_ptr = convert(LLVMType, LLVMPtr{T,AS.Local}; ctx)

        # create a function
        llvm_f, _ = create_function(T_ptr)

        # create the global variable
        mod = LLVM.parent(llvm_f)
        gv_typ = LLVM.ArrayType(eltyp, len)
        gv = GlobalVariable(mod, gv_typ, GPUCompiler.safe_name(string(id)), AS.Local)
        if len > 0
            linkage!(gv, LLVM.API.LLVMInternalLinkage)
            initializer!(gv, null(gv_typ))
        end
        # TODO: Make the alignment configurable
        alignment!(gv, Base.datatype_alignment(T))

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            ptr = gep!(builder, gv, [ConstantInt(0; ctx), ConstantInt(0; ctx)])

            untyped_ptr = bitcast!(builder, ptr, T_ptr)

            ret!(builder, untyped_ptr)
        end

        call_function(llvm_f, LLVMPtr{T,AS.Local})
    end
end


text:
@testset "Parallel.Distance" begin
    g4 = PathDiGraph(5)
    adjmx1 = [0 1 0; 1 0 1; 0 1 0] # graph
    adjmx2 = [0 1 0; 1 0 1; 1 1 0] # digraph
    a1 = SimpleGraph(adjmx1)
    a2 = SimpleDiGraph(adjmx2)
    distmx1 = [Inf 2.0 Inf; 2.0 Inf 4.2; Inf 4.2 Inf]
    distmx2 = [Inf 2.0 Inf; 3.2 Inf 4.2; 5.5 6.1 Inf]

    for g in testgraphs(a1)
        z = @inferred(LightGraphs.eccentricity(g, distmx1))
        y = @inferred(Parallel.eccentricity(g, distmx1))
        @test isapprox(y, z)
        @test @inferred(LightGraphs.diameter(y)) == @inferred(Parallel.diameter(g, distmx1)) == 6.2
        @test @inferred(LightGraphs.periphery(y)) == @inferred(Parallel.periphery(g, distmx1)) == [1, 3]
        @test @inferred(LightGraphs.radius(y)) == @inferred(Parallel.radius(g, distmx1)) == 4.2
        @test @inferred(LightGraphs.center(y)) == @inferred(Parallel.center(g, distmx1)) == [2]
    end

    for g in testdigraphs(a2)
        z = @inferred(LightGraphs.eccentricity(g, distmx2))
        y = @inferred(Parallel.eccentricity(g, distmx2))
        @test isapprox(y, z)
        @test @inferred(LightGraphs.diameter(y)) == @inferred(Parallel.diameter(g, distmx2)) == 6.2
        @test @inferred(LightGraphs.periphery(y)) == @inferred(Parallel.periphery(g, distmx2)) == [1]
        @test @inferred(LightGraphs.radius(y)) == @inferred(Parallel.radius(g, distmx2)) == 4.2
        @test @inferred(LightGraphs.center(y)) == @inferred(Parallel.center(g, distmx2)) == [2]
    end
end


text:
@testset "functions/summary.jl" begin
    # Write your tests here.
    # Test Array
    #x = [1.0, 2.0, -3.0]
    x=[-5.,-4.,-3.,-2.,-1.,0,1.,2.,3.,4.,5.]
    y=[1.,2.,3.,4.,5.,6.,7.,8.,9.,-10.,11.]
    z=[1.,1.,1.]

    @test abs_energy(x) == 110.
    @test abs_energy(y) == 506.

    @test absolute_maximum(x) == 5.
    @test absolute_maximum(y) == 11.

    @test absolute_sum_of_changes(x) == 10
    @test absolute_sum_of_changes(y) == 48

    @test agg_autocorrelation(x,mean,2) == 0.6212121212121213#0.6213
    @test agg_autocorrelation(y,mean,2) == -0.18300395256916996#-0.1825

    @test autocorrelation(x, 2) == 0.5151515151515152#0.5151
    @test autocorrelation(y,2) == 0.015546772068511216#0.0155

    @test benford_correlation(x) == 0.64055170771251
    @test benford_correlation(y) == 0.8641230464899805

    @test binned_entropy(y,3) == 0.3821755938143392
    @test binned_entropy(x,3) == 0.7167794185129729

    @test c3(x,5) == 0.
    @test c3(x,4) == 0.
    @test c3(y,5) == 66.
    @test c3(x,6) == 0.

    @test change_quantiles(x, 0., 1., mean) == 0
    @test change_quantiles(y, 0., 1., var) == 10.266666666666667

    @test cid_ce(x, true) == 0.9534625892455924
    @test cid_ce(y, false) == 28.460498941515414
    @test cid_ce(z, true) == 0.

    @test count_above(x, 4.) == 0.09090909090909091
    @test count_above(y, 3.) == 0.6363636363636364

    @test count_above_mean(x) == 5
    @test count_above_mean(y) == 6

    @test count_below(x, 0.) == 0.45454545454545453
    @test count_below(y, 1.) == 0.09090909090909091

    @test count_below_mean(x) == 5
    @test count_below_mean(y) == 5

    @test first_location_of_maximum(x) == 1
    @test first_location_of_maximum(y) == 1

    @test first_location_of_minimum(x) == 0.09090909090909091
    @test first_location_of_minimum(y) == 0.9090909090909091

    @test longest_strike_above_mean(x) == 5
    @test longest_strike_above_mean(y) == 5

    @test longest_strike_below_mean(x) == 5
    @test longest_strike_below_mean(y) == 4

    @test has_duplicate(x) == false
    @test has_duplicate(y) == false

    @test has_duplicate_max(x) == false
    @test has_duplicate_max(y) == false

    @test has_duplicate_min(x) == false
    @test has_duplicate_min(y) == false

    @test large_standard_deviation(x, 2.) == false
    @test large_standard_deviation(y, 2.) == false

    @test last_location_of_maximum(x) == 0.09090909090909091
    @test last_location_of_maximum(y) == 0.09090909090909091

    @test last_location_of_minimum(x) == 1
    @test last_location_of_minimum(y) == 0.18181818181818182

    @test mean_abs_change(x) == 1.0
    @test mean_abs_change(y) == 4.8

    @test mean_change(x) == 1.
    @test mean_change(y) == 1.

    @test mean_n_absolute_max(x, 3) == 4.666666666666667
    @test mean_n_absolute_max(y, 2) == 10.5

    @test mean_second_derivative_central(x) == 0.5555555555555556
    @test mean_second_derivative_central(y) == 1.0

    @test number_crossing_m(x, 2.0) == 1
    @test number_crossing_m(y, 7.0) == 3

    # @test percentage_of_reoccurring_datapoints_to_all_datapoints(x) == 0.0
    # @test percentage_of_reoccurring_datapoints_to_all_datapoints(y) == 0.0
    #
    # @test percentage_of_reoccurring_values_to_all_values(x) == 0.0
    # @test percentage_of_reoccurring_values_to_all_values(y) == 0.0

    @test range_count(x, 0., 2.) == 2
    @test range_count(y, 6., 9.) == 3

    @test ratio_beyond_r_sigma(x, 1.) == 1.6363636363636365
    @test ratio_beyond_r_sigma(y, 2.) == 1.2892561983471074

    # @test ratio_value_number_to_time_series_length(x) == 1.
    # @test ratio_value_number_to_time_series_length(y) == 1.

    @test root_mean_square(x) == 3.1622776601683795
    @test root_mean_square(y) == 6.782329983125268

#   @test sample_entropy(x) == NaN != NaN ...
    @test sample_entropy(y) == 0.11778303565638351

    # @test sum_of_reoccurring_data_points(x) == 0
    # @test sum_of_reoccurring_data_points(y) == 0

end


text:
"""
Model of 2-state synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""

function mdl_machine_ode!(device_states,
                          output_ode,
                          current_r,
                          current_i,
                          Sbase::Float64,
                          f0::Float64,
                          device::DynGenerator{BaseMachine, S, A, TG, P})  where {S <: Shaft,
                                                                                A <: AVR,
                                                                                TG <: TurbineGov,
                                                                                P <: PSS}

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.machine]
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = device.inner_vars[VR_gen_var]
    V_tI = device.inner_vars[VI_gen_var]

    #Get parameters
    R = device.machine.R
    Xd_p = device.machine.Xd_p
    eq_p = device.machine.eq_p
    BaseMVA = device.machine.MVABase

    #RI to dq transformation
    V_dq = ri_dq(δ)*[V_tR; V_tI]

    #Obtain electric variables
    i_d = (1.0/(R^2 + Xd_p^2)) * ( Xd_p*(eq_p - V_dq[2]) - R*V_dq[1])  #15.36
    i_q = (1.0/(R^2 + Xd_p^2)) * ( Xd_p*V_dq[1] + R*(eq_p - V_dq[2]) ) #15.36
    Pe = (V_dq[1]+ R*i_d)*i_d + (V_dq[2] + R*i_q)*i_q      #15.35

    #Update inner_vars
    device.inner_vars[τe_var] = Pe #Model assume ω approx 1.0

    #Compute current from the generator to the grid
    I_RI = (BaseMVA/Sbase)*dq_ri(δ)*[i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 4-state (One d- and One q-) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""

function mdl_machine_ode!(device_states,
                          output_ode,
                          current_r,
                          current_i,
                          Sbase::Float64,
                          f0::Float64,
                          device::DynGenerator{OneDOneQMachine, S, A, TG, P})  where {S <: Shaft,
                                                                                A <: AVR,
                                                                                TG <: TurbineGov,
                                                                                P <: PSS}

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.machine]

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ed_p = internal_states[2]

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.machine]
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = device.inner_vars[VR_gen_var]
    V_tI = device.inner_vars[VI_gen_var]
    Vf = device.inner_vars[Vf_var]

    #Get parameters
    R = device.machine.R
    Xd = device.machine.Xd
    Xq = device.machine.Xq
    Xd_p = device.machine.Xd_p
    Xq_p = device.machine.Xq_p
    Td0_p = device.machine.Td0_p
    Tq0_p = device.machine.Tq0_p
    BaseMVA = device.machine.MVABase

    #RI to dq transformation
    V_dq = ri_dq(δ)*[V_tR; V_tI]

    #Obtain electric variables
    i_d = (1.0/(R^2 + Xd_p*Xq_p)) * ( Xq_p*(eq_p - V_dq[2]) + R*(ed_p - V_dq[1]) )  #15.32
    i_q = (1.0/(R^2 + Xd_p*Xq_p)) * ( -Xd_p*(ed_p - V_dq[1]) + R*(eq_p - V_dq[2]) )  #15.32
    Pe = (V_dq[1]+ R*i_d)*i_d + (V_dq[2] + R*i_q)*i_q      #15.35

    #Compute ODEs
    output_ode[local_ix[1]] = (1.0/Td0_p)*(-eq_p - (Xd - Xd_p)*i_d + Vf)     #15.29 eq_p
    output_ode[local_ix[2]] = (1.0/Tq0_p)*(-ed_p + (Xq - Xq_p)*i_q)          #15.30 ed_p

    #Update inner_vars
    device.inner_vars[τe_var] = Pe #We should divide by ω, so we should receive ω

    #Compute current from the generator to the grid
    I_RI = (BaseMVA/Sbase)*dq_ri(δ)*[i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 6-state (MarconatoMachine) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""

function mdl_machine_ode!(device_states,
                          output_ode,
                          current_r,
                          current_i,
                          Sbase::Float64,
                          f0::Float64,
                          device::DynGenerator{MarconatoMachine, S, A, TG, P})  where {S <: Shaft,
                                                                                A <: AVR,
                                                                                TG <: TurbineGov,
                                                                                P <: PSS}

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.machine]

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    ψq = internal_states[1]
    ψd = internal_states[2]
    eq_p = internal_states[3]
    ed_p = internal_states[4]
    eq_pp = internal_states[5]
    ed_pp = internal_states[6]

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.machine]
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = device.inner_vars[VR_gen_var]
    V_tI = device.inner_vars[VI_gen_var]
    Vf = device.inner_vars[Vf_var]

    #Get parameters
    R = device.machine.R
    Xd = device.machine.Xd
    Xq = device.machine.Xq
    Xd_p = device.machine.Xd_p
    Xq_p = device.machine.Xq_p
    Xd_pp = device.machine.Xd_pp
    Xq_pp = device.machine.Xq_pp
    Td0_p = device.machine.Td0_p
    Tq0_p = device.machine.Tq0_p
    Td0_pp = device.machine.Td0_pp
    Tq0_pp = device.machine.Tq0_pp
    T_AA = device.machine.T_AA
    γd = device.machine.γd
    γq = device.machine.γq
    BaseMVA = device.machine.MVABase

    #RI to dq transformation
    V_dq = ri_dq(δ)*[V_tR; V_tI]

    #Obtain electric variables
    i_d = (1.0/Xd_pp)*(eq_pp - ψd)      #15.18
    i_q = (1.0/Xq_pp)*(-ed_pp - ψq)     #15.18
    τ_e = ψd*i_q - ψq*i_d               #15.6

    #Compute ODEs
    output_ode[local_ix[1]] = 2*π*f0*(R*i_q - ω*ψd + V_dq[2])                        #15.9 ψq
    output_ode[local_ix[2]] = 2*π*f0*(R*i_d + ω*ψq + V_dq[1])                        #15.9 ψd
    output_ode[local_ix[3]] = ( (1.0/Td0_p)*(-eq_p - (Xd - Xd_p - γd)*i_d
                                + (1-(T_AA/Td0_p))*Vf) )                             #15.16 eq_p
    output_ode[local_ix[4]] = (1.0/Tq0_p)*(-ed_p + (Xq - Xq_p - γq)*i_q)             #15.16 ed_p
    output_ode[local_ix[5]] = ( (1.0/Td0_pp)*(-eq_pp + eq_p
                                - (Xd_p - Xd_pp + γd)*i_d + (T_AA/Td0_p)*Vf) )       #15.16 eq_pp
    output_ode[local_ix[6]] = (1.0/Tq0_pp)*(-ed_pp + ed_p + (Xq_p - Xq_pp + γq)*i_q) #15.16 ed_pp

    #Update inner_vars
    device.inner_vars[τe_var] = τ_e

    #Compute current from the generator to the grid
    I_RI = (BaseMVA/Sbase)*dq_ri(δ)*[i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 4-state (SimpleMarconatoMachine) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""

function mdl_machine_ode!(device_states,
                          output_ode,
                          current_r,
                          current_i,
                          Sbase::Float64,
                          f0::Float64,
                          device::DynGenerator{SimpleMarconatoMachine, S, A, TG, P})  where {S <: Shaft,
                                                                                A <: AVR,
                                                                                TG <: TurbineGov,
                                                                                P <: PSS}

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.machine]

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ed_p = internal_states[2]
    eq_pp = internal_states[3]
    ed_pp = internal_states[4]

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.machine]
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = device.inner_vars[VR_gen_var]
    V_tI = device.inner_vars[VI_gen_var]
    Vf = device.inner_vars[Vf_var]

    #Get parameters
    R = device.machine.R
    Xd = device.machine.Xd
    Xq = device.machine.Xq
    Xd_p = device.machine.Xd_p
    Xq_p = device.machine.Xq_p
    Xd_pp = device.machine.Xd_pp
    Xq_pp = device.machine.Xq_pp
    Td0_p = device.machine.Td0_p
    Tq0_p = device.machine.Tq0_p
    Td0_pp = device.machine.Td0_pp
    Tq0_pp = device.machine.Tq0_pp
    T_AA = device.machine.T_AA
    γd = device.machine.γd
    γq = device.machine.γq
    BaseMVA = device.machine.MVABase

    #RI to dq transformation
    V_dq = ri_dq(δ)*[V_tR; V_tI]

    #Obtain electric variables
    #i_dq = inv([[R -Xq_pp]; [Xq_pp R]])*[ed_pp - V_dq[1]; eq_pp - V_dq[2]]
    i_d = (1.0/(R^2 + Xd_pp*Xq_pp)) * (Xq_pp*(eq_pp - V_dq[2])  +  R*(ed_pp - V_dq[1]))      #15.25
    i_q = (1.0/(R^2 + Xd_pp*Xq_pp)) * (-Xd_pp*(ed_pp - V_dq[1]) +  R*(eq_pp - V_dq[2]))      #15.25
    τ_e = (V_dq[1]+ R*i_d)*i_d + (V_dq[2] + R*i_q)*i_q               #15.6

    #Compute ODEs
    output_ode[local_ix[1]] = ( (1.0/Td0_p)*(-eq_p - (Xd - Xd_p - γd)*i_d
                                + (1-(T_AA/Td0_p))*Vf) )                             #15.16 eq_p
    output_ode[local_ix[2]] = (1.0/Tq0_p)*(-ed_p + (Xq - Xq_p - γq)*i_q)             #15.16 ed_p
    output_ode[local_ix[3]] = ( (1.0/Td0_pp)*(-eq_pp + eq_p
                                - (Xd_p - Xd_pp + γd)*i_d + (T_AA/Td0_p)*Vf) )       #15.16 eq_pp
    output_ode[local_ix[4]] = (1.0/Tq0_pp)*(-ed_pp + ed_p + (Xq_p - Xq_pp + γq)*i_q) #15.16 ed_pp

    #Update inner_vars
    device.inner_vars[τe_var] = τ_e

    #Compute current from the generator to the grid
    I_RI = (BaseMVA/Sbase)*dq_ri(δ)*[i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 6-state (AndersonFouadMachine) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""

function mdl_machine_ode!(device_states,
                          output_ode,
                          current_r,
                          current_i,
                          Sbase::Float64,
                          f0::Float64,
                          device::DynGenerator{AndersonFouadMachine, S, A, TG, P})  where {S <: Shaft,
                                                                                A <: AVR,
                                                                                TG <: TurbineGov,
                                                                                P <: PSS}

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.machine]

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    ψq = internal_states[1]
    ψd = internal_states[2]
    eq_p = internal_states[3]
    ed_p = internal_states[4]
    eq_pp = internal_states[5]
    ed_pp = internal_states[6]

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.machine]
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = device.inner_vars[VR_gen_var]
    V_tI = device.inner_vars[VI_gen_var]
    Vf = device.inner_vars[Vf_var]

    #Get parameters
    R = device.machine.R
    Xd = device.machine.Xd
    Xq = device.machine.Xq
    Xd_p = device.machine.Xd_p
    Xq_p = device.machine.Xq_p
    Xd_pp = device.machine.Xd_pp
    Xq_pp = device.machine.Xq_pp
    Td0_p = device.machine.Td0_p
    Tq0_p = device.machine.Tq0_p
    Td0_pp = device.machine.Td0_pp
    Tq0_pp = device.machine.Tq0_pp
    BaseMVA = device.machine.MVABase

    #RI to dq transformation
    V_dq = ri_dq(δ)*[V_tR; V_tI]

    #Obtain electric variables
    i_d = (1.0/Xd_pp)*(eq_pp - ψd)      #15.18
    i_q = (1.0/Xq_pp)*(-ed_pp - ψq)     #15.18
    τ_e = ψd*i_q - ψq*i_d               #15.6

    #Compute ODEs
    output_ode[local_ix[1]] = 2*π*f0*(R*i_q - ω*ψd + V_dq[2])                        #15.9 ψq
    output_ode[local_ix[2]] = 2*π*f0*(R*i_d + ω*ψq + V_dq[1])                        #15.9 ψd
    output_ode[local_ix[3]] = (1.0/Td0_p)*(-eq_p + (Xd - Xd_p)*i_d + Vf)             #15.19 eq_p
    output_ode[local_ix[4]] = (1.0/Tq0_p)*(-ed_p + (Xq - Xq_p)*i_q)                  #15.19 ed_p
    output_ode[local_ix[5]] = (1.0/Td0_pp)*(-eq_pp + eq_p - (Xd_p - Xd_pp)*i_d)      #15.19 eq_pp
    output_ode[local_ix[6]] = (1.0/Tq0_pp)*(-ed_pp + ed_p + (Xq_p - Xq_pp)*i_q)      #15.19 ed_pp

    #Update inner_vars
    device.inner_vars[τe_var] = τ_e

    #Compute current from the generator to the grid
    I_RI = (BaseMVA/Sbase)*dq_ri(δ)*[i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 4-state (SimpleAFMachine) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""

function mdl_machine_ode!(device_states,
                          output_ode,
                          current_r,
                          current_i,
                          Sbase::Float64,
                          f0::Float64,
                          device::DynGenerator{SimpleAFMachine, S, A, TG, P})  where {S <: Shaft,
                                                                                A <: AVR,
                                                                                TG <: TurbineGov,
                                                                                P <: PSS}

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.machine]

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ed_p = internal_states[2]
    eq_pp = internal_states[3]
    ed_pp = internal_states[4]

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.machine]
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = device.inner_vars[VR_gen_var]
    V_tI = device.inner_vars[VI_gen_var]
    Vf = device.inner_vars[Vf_var]

    #Get parameters
    R = device.machine.R
    Xd = device.machine.Xd
    Xq = device.machine.Xq
    Xd_p = device.machine.Xd_p
    Xq_p = device.machine.Xq_p
    Xd_pp = device.machine.Xd_pp
    Xq_pp = device.machine.Xq_pp
    Td0_p = device.machine.Td0_p
    Tq0_p = device.machine.Tq0_p
    Td0_pp = device.machine.Td0_pp
    Tq0_pp = device.machine.Tq0_pp
    BaseMVA = device.machine.MVABase

    #RI to dq transformation
    V_dq = ri_dq(δ)*[V_tR; V_tI]

    #Obtain electric variables
    i_d = (1.0/(R^2 + Xd_pp*Xq_pp)) * (Xq_pp*(eq_pp - V_dq[2])  +  R*(ed_pp - V_dq[1]))      #15.25
    i_q = (1.0/(R^2 + Xd_pp*Xq_pp)) * (-Xd_pp*(ed_pp - V_dq[1]) +  R*(eq_pp - V_dq[2]))      #15.25
    τ_e = (V_dq[1]+ R*i_d)*i_d + (V_dq[2] + R*i_q)*i_q               #15.6

    #Compute ODEs
    output_ode[local_ix[1]] = (1.0/Td0_p)*(-eq_p + (Xd - Xd_p)*i_d + Vf)            #15.19 eq_p
    output_ode[local_ix[2]] = (1.0/Tq0_p)*(-ed_p + (Xq - Xq_p)*i_q)                 #15.19 ed_p
    output_ode[local_ix[3]] = (1.0/Td0_pp)*(-eq_pp + eq_p - (Xd_p - Xd_pp)*i_d)     #15.19 eq_pp
    output_ode[local_ix[4]] = (1.0/Tq0_pp)*(-ed_pp + ed_p + (Xq_p - Xq_pp)*i_q)     #15.19 ed_pp

    #Update inner_vars
    device.inner_vars[τe_var] = τ_e

    #Compute current from the generator to the grid
    I_RI = (BaseMVA/Sbase)*dq_ri(δ)*[i_d; i_q]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 5-state (FullOrderMachine) synchronous machine in Julia.
Refer to Power System Dynamics: Stability and Control, by J. Machowski, J. Bialek and J. Bumby,
or Power System Stability and Control by P. Kundur, for the equations
"""

function mdl_machine_ode!(device_states,
                          output_ode,
                          current_r,
                          current_i,
                          Sbase::Float64,
                          f0::Float64,
                          device::DynGenerator{FullMachine, S, A, TG, P})  where {S <: Shaft,
                                                                                A <: AVR,
                                                                                TG <: TurbineGov,
                                                                                P <: PSS}

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.machine]

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    ψd = internal_states[1]
    ψq = internal_states[2]
    ψf = internal_states[1]
    ψ1d = internal_states[2]
    ψ1q = internal_states[3]

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.machine]
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = device.inner_vars[VR_gen_var]
    V_tI = device.inner_vars[VI_gen_var]
    Vf = device.inner_vars[Vf_var]
    #ψd = device.inner_vars[ψd_var]
    #ψq = device.inner_vars[ψq_var]

    #Get parameters
    R = device.machine.R
    R_f = device.machine.R_f
    R_1d= device.machine.R_1d
    R_1q = device.machine.R_1q
    inv_d_fluxlink = device.machine.inv_d_fluxlink
    inv_q_fluxlink = device.machine.inv_q_fluxlink
    BaseMVA = device.machine.MVABase

    #RI to dq transformation
    V_dq = ri_dq(δ)*[V_tR; V_tI]

    #Obtain electric currents] using Flux Linkage equations
    i_df1d = inv_d_fluxlink*[ψd; ψf; ψ1d]  #11.18 in Machowski (3.127, 3.130 and 3.131 in Kundur)
    i_q1q = inv_q_fluxlink*[ψq; ψ1q]        #11.19 in Machowski (3.128 and 3.132 in Kundur)

    #Compute ODEs
    output_ode[local_ix[1]] = 2*pi*f0*(ψq*ω + R*i_df1d[1] + V_dq[1])  #11.32 in Machowski or 3.120 in Kundur
    output_ode[local_ix[2]] = 2*pi*f0*(-ψd*ω + R*i_q1q[1] + V_dq[2])
    output_ode[local_ix[3]] = -R_f*i_df1d[2] + Vf      #11.33 in Machowski or 3.123 in Kundur
    output_ode[local_ix[4]] = -R_1d*i_df1d[3]           #11.33 in Machowski or 3.124 in Kundur
    output_ode[local_ix[5]] = -R_1q*i_q1q[2]            #11.33 in Machowski or 3.125 in Kundur

    #Update stator fluxes and torque
    #ψd = (1.0/ω)*(R*i_q1q[1] + V_dq[2])         #11.32 in Machowski or 3.121 in Kundur
    #ψq = (1.0/ω)*(-R*i_df1d[1] - V_dq[1])       #11.32 in Machowski or 3.120 in Kundur
    τ_e = ψd*i_q1q[1] - ψq*i_df1d[1]            #15.6 in Milano or 3.117 in Kundur

    #Update inner_vars
    device.inner_vars[τe_var] = τ_e
    device.inner_vars[ψd_var] = ψd
    device.inner_vars[ψq_var] = ψq

    #Compute current from the generator to the grid
    I_RI = (BaseMVA/Sbase)*dq_ri(δ)*[i_df1d[1]; i_q1q[1]]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end

"""
Model of 3-state (FullOrderMachine) synchronous machine in Julia.
Refer to Power System Dynamics: Stability and Control, by J. Machowski, J. Bialek and J. Bumby,
or Power System Stability and Control by P. Kundur, for the equations
"""

function mdl_machine_ode!(device_states,
                          output_ode,
                          current_r,
                          current_i,
                          Sbase::Float64,
                          f0::Float64,
                          device::DynGenerator{SimpleFullMachine, S, A, TG, P})  where {S <: Shaft,
                                                                                A <: AVR,
                                                                                TG <: TurbineGov,
                                                                                P <: PSS}

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.machine]

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    ψf = internal_states[1]
    ψ1d = internal_states[2]
    ψ1q = internal_states[3]

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.machine]
    δ = device_states[external_ix[1]]
    ω = device_states[external_ix[2]]

    #Obtain inner variables for component
    V_tR = device.inner_vars[VR_gen_var]
    V_tI = device.inner_vars[VI_gen_var]
    Vf = device.inner_vars[Vf_var]
    ψd = device.inner_vars[ψd_var]
    ψq = device.inner_vars[ψq_var]

    #Get parameters
    R = device.machine.R
    R_f = device.machine.R_f
    R_1d= device.machine.R_1d
    R_1q = device.machine.R_1q
    inv_d_fluxlink = device.machine.inv_d_fluxlink
    inv_q_fluxlink = device.machine.inv_q_fluxlink
    BaseMVA = device.machine.MVABase

    #RI to dq transformation
    V_dq = ri_dq(δ)*[V_tR; V_tI]

    #Obtain electric currents] using Flux Linkage equations
    i_df1d = inv_d_fluxlink*[ψd; ψf; ψ1d]  #11.18 in Machowski (3.127, 3.130 and 3.131 in Kundur)
    i_q1q = inv_q_fluxlink*[ψq; ψ1q]        #11.19 in Machowski (3.128 and 3.132 in Kundur)

    #Compute ODEs
    output_ode[local_ix[1]] = -R_f*i_df1d[2] + Vf      #11.33 in Machowski or 3.123 in Kundur
    output_ode[local_ix[2]] = -R_1d*i_df1d[3]           #11.33 in Machowski or 3.124 in Kundur
    output_ode[local_ix[3]] = -R_1q*i_q1q[2]            #11.33 in Machowski or 3.125 in Kundur

    #Update stator fluxes and torque
    ψd = (1.0/ω)*(R*i_q1q[1] + V_dq[2])         #11.32 in Machowski or 3.121 in Kundur
    ψq = (1.0/ω)*(-R*i_df1d[1] - V_dq[1])       #11.32 in Machowski or 3.120 in Kundur
    τ_e = ψd*i_q1q[1] - ψq*i_df1d[1]            #15.6 in Milano or 3.117 in Kundur

    #Update inner_vars
    device.inner_vars[τe_var] = τ_e
    device.inner_vars[ψd_var] = ψd
    device.inner_vars[ψq_var] = ψq

    #Compute current from the generator to the grid
    I_RI = (BaseMVA/Sbase)*dq_ri(δ)*[i_df1d[1]; i_q1q[1]]

    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]

    return
end