text:
module GLUT

import GetC.@getCFun

#Handle platform invariances...
const libGlut = @windows? "freeglut" : @linux? "libglut" : @osx? "/System/Library/Frameworks/GLUT.framework/GLUT" : error("Platform not recognized. Try setting GLUT library name yourself")

@getCFun libGlut glutInit glutInit(pargc::Ptr{Int32}, argv::Ptr{Ptr{Uint8}})::Void
glutInit() = glutInit(convert(Ptr{Int32},pointer([1])), ["a"])
export glutInit

@getCFun libGlut glutInitWindowPosition glutInitWindowPosition(x::Int32, y::Int32)::Void
export glutInitWindowPosition
@getCFun libGlut glutInitWindowSize glutInitWindowSize(width::Int32, height::Int32)::Void
export glutInitWindowSize
@getCFun libGlut glutInitDisplayMode glutInitDisplayMode(displayMode::Uint32)::Void
export glutInitDisplayMode
@getCFun libGlut glutInitDisplayString glutInitDisplayString(displayMode::Ptr{Uint8})::Void
export glutInitDisplayString

@getCFun libGlut glutMainLoop glutMainLoop()::Void
export glutMainLoop

@getCFun libGlut glutCreateWindow glutCreateWindow(title::Ptr{Uint8})::Int32
export glutCreateWindow
@getCFun libGlut glutCreateSubWindow glutCreateSubWindow(window::Int32, x::Int32, y::Int32, width::Int32, height::Int32)::Int32
export glutCreateSubWindow
@getCFun libGlut glutDestroyWindow glutDestroyWindow(window::Int32)::Void
export glutDestroyWindow
@getCFun libGlut glutSetWindow glutSetWindow(window::Int32)::Void
export glutSetWindow
@getCFun libGlut glutGetWindow glutGetWindow()::Int32
export glutGetWindow
@getCFun libGlut glutSetWindowTitle glutSetWindowTitle(title::Ptr{Uint8})::Void
export glutSetWindowTitle
@getCFun libGlut glutSetIconTitle glutSetIconTitle(title::Ptr{Uint8})::Void
export glutSetIconTitle
@getCFun libGlut glutReshapeWindow glutReshapeWindow(width::Int32, height::Int32)::Void
export glutReshapeWindow
@getCFun libGlut glutPositionWindow glutPositionWindow(x::Int32, y::Int32)::Void
export glutPositionWindow
@getCFun libGlut glutShowWindow glutShowWindow()::Void
export glutShowWindow
@getCFun libGlut glutHideWindow glutHideWindow()::Void
export glutHideWindow
@getCFun libGlut glutIconifyWindow glutIconifyWindow()::Void
export glutIconifyWindow
@getCFun libGlut glutPushWindow glutPushWindow()::Void
export glutPushWindow
@getCFun libGlut glutPopWindow glutPopWindow()::Void
export glutPopWindow
@getCFun libGlut glutFullScreen glutFullScreen()::Void
export glutFullScreen

@getCFun libGlut glutPostWindowRedisplay glutPostWindowRedisplay(window::Int32)::Void
export glutPostWindowRedisplay
@getCFun libGlut glutPostRedisplay glutPostRedisplay()::Void
export glutPostRedisplay
@getCFun libGlut glutSwapBuffers glutSwapBuffers()::Void
export glutSwapBuffers

@getCFun libGlut glutTimerFunc glutTimerFunc(time::Uint32, callback::Ptr{None}, value::Int32)::Void
export glutTimerFunc
@getCFun libGlut glutIdleFunc glutIdleFunc(callback::Ptr{None})::Void
export glutIdleFunc

@getCFun libGlut glutKeyboardFunc glutKeyboardFunc(callback::Ptr{None})::Void
export glutKeyboardFunc
@getCFun libGlut glutKeyboardUpFunc glutKeyboardUpFunc(callback::Ptr{None})::Void
export glutKeyboardUpFunc
@getCFun libGlut glutSpecialFunc glutSpecialFunc(callback::Ptr{None})::Void
export glutSpecialFunc
@getCFun libGlut glutSpecialUpFunc glutSpecialUpFunc(callback::Ptr{None})::Void
export glutSpecialUpFunc
@getCFun libGlut glutSetKeyRepeat glutSetKeyRepeat(repeatMode::Cint)::Cint
export glutSetKeyRepeat

@getCFun libGlut glutReshapeFunc glutReshapeFunc(callback::Ptr{None})::Void
export glutReshapeFunc
@getCFun libGlut glutVisibilityFunc glutVisibilityFunc(callback::Ptr{None})::Void
export glutVisibilityFunc
@getCFun libGlut glutDisplayFunc glutDisplayFunc(callback::Ptr{None})::Void
export glutDisplayFunc
@getCFun libGlut glutMouseFunc glutMouseFunc(callback::Ptr{None})::Void
export glutMouseFunc
@getCFun libGlut glutMotionFunc glutMotionFunc(callback::Ptr{None})::Void
export glutMotionFunc
@getCFun libGlut glutPassiveMotionFunc glutPassiveMotionFunc(callback::Ptr{None})::Void
export glutPassiveMotionFunc
@getCFun libGlut glutEntryFunc glutEntryFunc(callback::Ptr{None})::Void
export glutEntryFunc
@getCFun libGlut glutCloseFunc glutCloseFunc(callback::Ptr{None})::Void
export glutCloseFunc

@getCFun libGlut glutGetProcAddress glutGetProcAddress(name::Ptr{Uint8})::Ptr{Void}
export glutGetProcAddress
@getCFun libGlut glutGet glutGet(state::Uint16)::Cint
export glutGet

const  FREEGLUT                        = 1
const  GLUT_API_VERSION                = 4
const  FREEGLUT_VERSION_2_0            = 1
const  GLUT_XLIB_IMPLEMENTATION        = 13

const  GLUT_KEY_F1                     = 0x0001
const  GLUT_KEY_F2                     = 0x0002
const  GLUT_KEY_F3                     = 0x0003
const  GLUT_KEY_F4                     = 0x0004
const  GLUT_KEY_F5                     = 0x0005
const  GLUT_KEY_F6                     = 0x0006
const  GLUT_KEY_F7                     = 0x0007
const  GLUT_KEY_F8                     = 0x0008
const  GLUT_KEY_F9                     = 0x0009
const  GLUT_KEY_F10                    = 0x000A
const  GLUT_KEY_F11                    = 0x000B
const  GLUT_KEY_F12                    = 0x000C
const  GLUT_KEY_LEFT                   = 0x0064
const  GLUT_KEY_UP                     = 0x0065
const  GLUT_KEY_RIGHT                  = 0x0066
const  GLUT_KEY_DOWN                   = 0x0067
const  GLUT_KEY_PAGE_UP                = 0x0068
const  GLUT_KEY_PAGE_DOWN              = 0x0069
const  GLUT_KEY_HOME                   = 0x006A
const  GLUT_KEY_END                    = 0x006B
const  GLUT_KEY_INSERT                 = 0x006C

const  GLUT_LEFT_BUTTON                = 0x0000
const  GLUT_MIDDLE_BUTTON              = 0x0001
const  GLUT_RIGHT_BUTTON               = 0x0002
const  GLUT_DOWN                       = 0x0000
const  GLUT_UP                         = 0x0001
const  GLUT_LEFT                       = 0x0000
const  GLUT_ENTERED                    = 0x0001

const  GLUT_RGB                        = 0x0000
const  GLUT_RGBA                       = 0x0000
const  GLUT_INDEX                      = 0x0001
const  GLUT_SINGLE                     = 0x0000
const  GLUT_DOUBLE                     = 0x0002
const  GLUT_ACCUM                      = 0x0004
const  GLUT_ALPHA                      = 0x0008
const  GLUT_DEPTH                      = 0x0010
const  GLUT_STENCIL                    = 0x0020
const  GLUT_MULTISAMPLE                = 0x0080
const  GLUT_STEREO                     = 0x0100
const  GLUT_LUMINANCE                  = 0x0200

const  GLUT_MENU_NOT_IN_USE            = 0x0000
const  GLUT_MENU_IN_USE                = 0x0001
const  GLUT_NOT_VISIBLE                = 0x0000
const  GLUT_VISIBLE                    = 0x0001
const  GLUT_HIDDEN                     = 0x0000
const  GLUT_FULLY_RETAINED             = 0x0001
const  GLUT_PARTIALLY_RETAINED         = 0x0002
const  GLUT_FULLY_COVERED              = 0x0003

const  GLUT_WINDOW_X                   = 0x0064
const  GLUT_WINDOW_Y                   = 0x0065
const  GLUT_WINDOW_WIDTH               = 0x0066
const  GLUT_WINDOW_HEIGHT              = 0x0067
const  GLUT_WINDOW_BUFFER_SIZE         = 0x0068
const  GLUT_WINDOW_STENCIL_SIZE        = 0x0069
const  GLUT_WINDOW_DEPTH_SIZE          = 0x006A
const  GLUT_WINDOW_RED_SIZE            = 0x006B
const  GLUT_WINDOW_GREEN_SIZE          = 0x006C
const  GLUT_WINDOW_BLUE_SIZE           = 0x006D
const  GLUT_WINDOW_ALPHA_SIZE          = 0x006E
const  GLUT_WINDOW_ACCUM_RED_SIZE      = 0x006F
const  GLUT_WINDOW_ACCUM_GREEN_SIZE    = 0x0070
const  GLUT_WINDOW_ACCUM_BLUE_SIZE     = 0x0071
const  GLUT_WINDOW_ACCUM_ALPHA_SIZE    = 0x0072
const  GLUT_WINDOW_DOUBLEBUFFER        = 0x0073
const  GLUT_WINDOW_RGBA                = 0x0074
const  GLUT_WINDOW_PARENT              = 0x0075
const  GLUT_WINDOW_NUM_CHILDREN        = 0x0076
const  GLUT_WINDOW_COLORMAP_SIZE       = 0x0077
const  GLUT_WINDOW_NUM_SAMPLES         = 0x0078
const  GLUT_WINDOW_STEREO              = 0x0079
const  GLUT_WINDOW_CURSOR              = 0x007A

const  GLUT_SCREEN_WIDTH               = 0x00C8
const  GLUT_SCREEN_HEIGHT              = 0x00C9
const  GLUT_SCREEN_WIDTH_MM            = 0x00CA
const  GLUT_SCREEN_HEIGHT_MM           = 0x00CB
const  GLUT_MENU_NUM_ITEMS             = 0x012C
const  GLUT_DISPLAY_MODE_POSSIBLE      = 0x0190
const  GLUT_INIT_WINDOW_X              = 0x01F4
const  GLUT_INIT_WINDOW_Y              = 0x01F5
const  GLUT_INIT_WINDOW_WIDTH          = 0x01F6
const  GLUT_INIT_WINDOW_HEIGHT         = 0x01F7
const  GLUT_INIT_DISPLAY_MODE          = 0x01F8
const  GLUT_ELAPSED_TIME               = 0x02BC
const  GLUT_WINDOW_FORMAT_ID           = 0x007B

const  GLUT_HAS_KEYBOARD               = 0x0258
const  GLUT_HAS_MOUSE                  = 0x0259
const  GLUT_HAS_SPACEBALL              = 0x025A
const  GLUT_HAS_DIAL_AND_BUTTON_BOX    = 0x025B
const  GLUT_HAS_TABLET                 = 0x025C
const  GLUT_NUM_MOUSE_BUTTONS          = 0x025D
const  GLUT_NUM_SPACEBALL_BUTTONS      = 0x025E
const  GLUT_NUM_BUTTON_BOX_BUTTONS     = 0x025F
const  GLUT_NUM_DIALS                  = 0x0260
const  GLUT_NUM_TABLET_BUTTONS         = 0x0261
const  GLUT_DEVICE_IGNORE_KEY_REPEAT   = 0x0262
const  GLUT_DEVICE_KEY_REPEAT          = 0x0263
const  GLUT_HAS_JOYSTICK               = 0x0264
const  GLUT_OWNS_JOYSTICK              = 0x0265
const  GLUT_JOYSTICK_BUTTONS           = 0x0266
const  GLUT_JOYSTICK_AXES              = 0x0267
const  GLUT_JOYSTICK_POLL_RATE         = 0x0268

const  GLUT_OVERLAY_POSSIBLE           = 0x0320
const  GLUT_LAYER_IN_USE               = 0x0321
const  GLUT_HAS_OVERLAY                = 0x0322
const  GLUT_TRANSPARENT_INDEX          = 0x0323
const  GLUT_NORMAL_DAMAGED             = 0x0324
const  GLUT_OVERLAY_DAMAGED            = 0x0325

const  GLUT_VIDEO_RESIZE_POSSIBLE      = 0x0384
const  GLUT_VIDEO_RESIZE_IN_USE        = 0x0385
const  GLUT_VIDEO_RESIZE_X_DELTA       = 0x0386
const  GLUT_VIDEO_RESIZE_Y_DELTA       = 0x0387
const  GLUT_VIDEO_RESIZE_WIDTH_DELTA   = 0x0388
const  GLUT_VIDEO_RESIZE_HEIGHT_DELTA  = 0x0389
const  GLUT_VIDEO_RESIZE_X             = 0x038A
const  GLUT_VIDEO_RESIZE_Y             = 0x038B
const  GLUT_VIDEO_RESIZE_WIDTH         = 0x038C
const  GLUT_VIDEO_RESIZE_HEIGHT        = 0x038D

const  GLUT_NORMAL                     = 0x0000
const  GLUT_OVERLAY                    = 0x0001

const  GLUT_ACTIVE_SHIFT               = 0x0001
const  GLUT_ACTIVE_CTRL                = 0x0002
const  GLUT_ACTIVE_ALT                 = 0x0004

const  GLUT_CURSOR_RIGHT_ARROW         = 0x0000
const  GLUT_CURSOR_LEFT_ARROW          = 0x0001
const  GLUT_CURSOR_INFO                = 0x0002
const  GLUT_CURSOR_DESTROY             = 0x0003
const  GLUT_CURSOR_HELP                = 0x0004
const  GLUT_CURSOR_CYCLE               = 0x0005
const  GLUT_CURSOR_SPRAY               = 0x0006
const  GLUT_CURSOR_WAIT                = 0x0007
const  GLUT_CURSOR_TEXT                = 0x0008
const  GLUT_CURSOR_CROSSHAIR           = 0x0009
const  GLUT_CURSOR_UP_DOWN             = 0x000A
const  GLUT_CURSOR_LEFT_RIGHT          = 0x000B
const  GLUT_CURSOR_TOP_SIDE            = 0x000C
const  GLUT_CURSOR_BOTTOM_SIDE         = 0x000D
const  GLUT_CURSOR_LEFT_SIDE           = 0x000E
const  GLUT_CURSOR_RIGHT_SIDE          = 0x000F
const  GLUT_CURSOR_TOP_LEFT_CORNER     = 0x0010
const  GLUT_CURSOR_TOP_RIGHT_CORNER    = 0x0011
const  GLUT_CURSOR_BOTTOM_RIGHT_CORNER = 0x0012
const  GLUT_CURSOR_BOTTOM_LEFT_CORNER  = 0x0013
const  GLUT_CURSOR_INHERIT             = 0x0064
const  GLUT_CURSOR_NONE                = 0x0065
const  GLUT_CURSOR_FULL_CROSSHAIR      = 0x0066

const  GLUT_RED                        = 0x0000
const  GLUT_GREEN                      = 0x0001
const  GLUT_BLUE                       = 0x0002

const  GLUT_KEY_REPEAT_OFF             = 0x0000
const  GLUT_KEY_REPEAT_ON              = 0x0001
const  GLUT_KEY_REPEAT_DEFAULT         = 0x0002

const  GLUT_JOYSTICK_BUTTON_A          = 0x0001
const  GLUT_JOYSTICK_BUTTON_B          = 0x0002
const  GLUT_JOYSTICK_BUTTON_C          = 0x0004
const  GLUT_JOYSTICK_BUTTON_D          = 0x0008

const  GLUT_GAME_MODE_ACTIVE           = 0x0000
const  GLUT_GAME_MODE_POSSIBLE         = 0x0001
const  GLUT_GAME_MODE_WIDTH            = 0x0002
const  GLUT_GAME_MODE_HEIGHT           = 0x0003
const  GLUT_GAME_MODE_PIXEL_DEPTH      = 0x0004
const  GLUT_GAME_MODE_REFRESH_RATE     = 0x0005
const  GLUT_GAME_MODE_DISPLAY_CHANGED  = 0x0006

export  FREEGLUT
export  GLUT_API_VERSION
export  FREEGLUT_VERSION_2_0
export  GLUT_XLIB_IMPLEMENTATION

export  GLUT_KEY_F1
export  GLUT_KEY_F2
export  GLUT_KEY_F3
export  GLUT_KEY_F4
export  GLUT_KEY_F5
export  GLUT_KEY_F6
export  GLUT_KEY_F7
export  GLUT_KEY_F8
export  GLUT_KEY_F9
export  GLUT_KEY_F10
export  GLUT_KEY_F11
export  GLUT_KEY_F12
export  GLUT_KEY_LEFT
export  GLUT_KEY_UP
export  GLUT_KEY_RIGHT
export  GLUT_KEY_DOWN
export  GLUT_KEY_PAGE_UP
export  GLUT_KEY_PAGE_DOWN
export  GLUT_KEY_HOME
export  GLUT_KEY_END
export  GLUT_KEY_INSERT

export  GLUT_LEFT_BUTTON
export  GLUT_MIDDLE_BUTTON
export  GLUT_RIGHT_BUTTON
export  GLUT_DOWN
export  GLUT_UP
export  GLUT_LEFT
export  GLUT_ENTERED

export  GLUT_RGB
export  GLUT_RGBA
export  GLUT_INDEX
export  GLUT_SINGLE
export  GLUT_DOUBLE
export  GLUT_ACCUM
export  GLUT_ALPHA
export  GLUT_DEPTH
export  GLUT_STENCIL
export  GLUT_MULTISAMPLE
export  GLUT_STEREO
export  GLUT_LUMINANCE

export  GLUT_MENU_NOT_IN_USE
export  GLUT_MENU_IN_USE
export  GLUT_NOT_VISIBLE
export  GLUT_VISIBLE
export  GLUT_HIDDEN
export  GLUT_FULLY_RETAINED
export  GLUT_PARTIALLY_RETAINED
export  GLUT_FULLY_COVERED

export  GLUT_WINDOW_X
export  GLUT_WINDOW_Y
export  GLUT_WINDOW_WIDTH
export  GLUT_WINDOW_HEIGHT
export  GLUT_WINDOW_BUFFER_SIZE
export  GLUT_WINDOW_STENCIL_SIZE
export  GLUT_WINDOW_DEPTH_SIZE
export  GLUT_WINDOW_RED_SIZE
export  GLUT_WINDOW_GREEN_SIZE
export  GLUT_WINDOW_BLUE_SIZE
export  GLUT_WINDOW_ALPHA_SIZE
export  GLUT_WINDOW_ACCUM_RED_SIZE
export  GLUT_WINDOW_ACCUM_GREEN_SIZE
export  GLUT_WINDOW_ACCUM_BLUE_SIZE
export  GLUT_WINDOW_ACCUM_ALPHA_SIZE
export  GLUT_WINDOW_DOUBLEBUFFER
export  GLUT_WINDOW_RGBA
export  GLUT_WINDOW_PARENT
export  GLUT_WINDOW_NUM_CHILDREN
export  GLUT_WINDOW_COLORMAP_SIZE
export  GLUT_WINDOW_NUM_SAMPLES
export  GLUT_WINDOW_STEREO
export  GLUT_WINDOW_CURSOR

export  GLUT_SCREEN_WIDTH
export  GLUT_SCREEN_HEIGHT
export  GLUT_SCREEN_WIDTH_MM
export  GLUT_SCREEN_HEIGHT_MM
export  GLUT_MENU_NUM_ITEMS
export  GLUT_DISPLAY_MODE_POSSIBLE
export  GLUT_INIT_WINDOW_X
export  GLUT_INIT_WINDOW_Y
export  GLUT_INIT_WINDOW_WIDTH
export  GLUT_INIT_WINDOW_HEIGHT
export  GLUT_INIT_DISPLAY_MODE
export  GLUT_ELAPSED_TIME
export  GLUT_WINDOW_FORMAT_ID

export  GLUT_HAS_KEYBOARD
export  GLUT_HAS_MOUSE
export  GLUT_HAS_SPACEBALL
export  GLUT_HAS_DIAL_AND_BUTTON_BOX
export  GLUT_HAS_TABLET
export  GLUT_NUM_MOUSE_BUTTONS
export  GLUT_NUM_SPACEBALL_BUTTONS
export  GLUT_NUM_BUTTON_BOX_BUTTONS
export  GLUT_NUM_DIALS
export  GLUT_NUM_TABLET_BUTTONS
export  GLUT_DEVICE_IGNORE_KEY_REPEAT
export  GLUT_DEVICE_KEY_REPEAT
export  GLUT_HAS_JOYSTICK
export  GLUT_OWNS_JOYSTICK
export  GLUT_JOYSTICK_BUTTONS
export  GLUT_JOYSTICK_AXES
export  GLUT_JOYSTICK_POLL_RATE

export  GLUT_OVERLAY_POSSIBLE
export  GLUT_LAYER_IN_USE
export  GLUT_HAS_OVERLAY
export  GLUT_TRANSPARENT_INDEX
export  GLUT_NORMAL_DAMAGED
export  GLUT_OVERLAY_DAMAGED

export  GLUT_VIDEO_RESIZE_POSSIBLE
export  GLUT_VIDEO_RESIZE_IN_USE
export  GLUT_VIDEO_RESIZE_X_DELTA
export  GLUT_VIDEO_RESIZE_Y_DELTA
export  GLUT_VIDEO_RESIZE_WIDTH_DELTA
export  GLUT_VIDEO_RESIZE_HEIGHT_DELTA
export  GLUT_VIDEO_RESIZE_X
export  GLUT_VIDEO_RESIZE_Y
export  GLUT_VIDEO_RESIZE_WIDTH
export  GLUT_VIDEO_RESIZE_HEIGHT

export  GLUT_NORMAL
export  GLUT_OVERLAY

export  GLUT_ACTIVE_SHIFT
export  GLUT_ACTIVE_CTRL
export  GLUT_ACTIVE_ALT

export  GLUT_CURSOR_RIGHT_ARROW
export  GLUT_CURSOR_LEFT_ARROW
export  GLUT_CURSOR_INFO
export  GLUT_CURSOR_DESTROY
export  GLUT_CURSOR_HELP
export  GLUT_CURSOR_CYCLE
export  GLUT_CURSOR_SPRAY
export  GLUT_CURSOR_WAIT
export  GLUT_CURSOR_TEXT
export  GLUT_CURSOR_CROSSHAIR
export  GLUT_CURSOR_UP_DOWN
export  GLUT_CURSOR_LEFT_RIGHT
export  GLUT_CURSOR_TOP_SIDE
export  GLUT_CURSOR_BOTTOM_SIDE
export  GLUT_CURSOR_LEFT_SIDE
export  GLUT_CURSOR_RIGHT_SIDE
export  GLUT_CURSOR_TOP_LEFT_CORNER
export  GLUT_CURSOR_TOP_RIGHT_CORNER
export  GLUT_CURSOR_BOTTOM_RIGHT_CORNER
export  GLUT_CURSOR_BOTTOM_LEFT_CORNER
export  GLUT_CURSOR_INHERIT
export  GLUT_CURSOR_NONE
export  GLUT_CURSOR_FULL_CROSSHAIR

export  GLUT_RED
export  GLUT_GREEN
export  GLUT_BLUE

export  GLUT_KEY_REPEAT_OFF
export  GLUT_KEY_REPEAT_ON
export  GLUT_KEY_REPEAT_DEFAULT

export  GLUT_JOYSTICK_BUTTON_A
export  GLUT_JOYSTICK_BUTTON_B
export  GLUT_JOYSTICK_BUTTON_C
export  GLUT_JOYSTICK_BUTTON_D

export  GLUT_GAME_MODE_ACTIVE
export  GLUT_GAME_MODE_POSSIBLE
export  GLUT_GAME_MODE_WIDTH
export  GLUT_GAME_MODE_HEIGHT
export  GLUT_GAME_MODE_PIXEL_DEPTH
export  GLUT_GAME_MODE_REFRESH_RATE
export  GLUT_GAME_MODE_DISPLAY_CHANGED

end


text:
using JuLIP: Atoms,
             neighbourlist,
             site_energies,
             sites,
             JVec

using JuLIP.Potentials: site_virial,
                        evaluate,
                        evaluate_d,
                        evaluate_d!

using NeighbourLists: max_neigs

using KahanSummation

import JuLIP: site_energies,
              energy,
              forces,
              virial,
              cutoff

import NBodyIPs: evaluate_many!,
                 evaluate_many_d!,
                 evaluate_many_ed!,
                 eval_site_nbody!,
                 _grad_len2pos!,
                 TempEdV

function cutoff(V::AbstractEnvIP)
   @assert cutoff(Vn(V)) <= cutoff(Vr(V))
   return cutoff(Vr(V::EnvIP))
end

n_fun(V::EnvIP, n) = (V.t == 0) ? 1.0 : n^V.t

function n_fun_d(V::EnvIP, n)
   if V.t == 0
      return 0.0
   elseif V.t == 1
      return 1.0
   else
      return V.t * n^(V.t-1)
   end
end

site_ns(V::EnvIP, at) = n_fun.(Ref(V), site_energies(Vn(V), at))

function site_ns_ed(V::EnvIP, at)
   Vns = site_energies(Vn(V), at)
   return n_fun.(Ref(V), Vns), n_fun_d.(Ref(V), Vns)
end

function site_n_d!(dVn, V::EnvIP, r, R, Ni, dNi)
   for n = 1:length(r)
      dVn[n] = 0.5 * dNi * evaluate_d(Vn(V), r[n]) * R[n] / r[n]
   end
   return dVn
end

site_energies(V::AbstractEnvIP, at::Atoms) =
      site_ns(V, at) .* site_energies(Vr(V), at)

energy(V::AbstractEnvIP, at::Atoms) = sum_kbn(site_energies(V, at))

function forces(V::AbstractEnvIP{N}, at::Atoms{T}) where {N, T}

   # compute the n values
   Ns, dNs = site_ns_ed(V, at)
   # compute the inner v values
   Vs = site_energies(Vr(V), at)

   # now assemble site forces and use those to create
   # total forces by mixing N and V
   cutoff_n = cutoff(Vn(V))
   nlist = neighbourlist(at, cutoff(V))  # this checks that cutoff(Vn) <= cutoff(Vr)
   maxneigs = max_neigs(nlist)
   F = zeros(JVec{T}, length(at))
   dVsite = zeros(JVec{T}, maxneigs)
   dVn = zeros(JVec{T}, maxneigs)

   for (i, j, r, R) in sites(nlist)
      # compute the site energy gradients
      fill!(dVsite, zero(JVec{T}))
      eval_site_nbody!(
            Val(N), i, j, R, cutoff(V), false,
            (out, R, i, J, temp) -> evaluate_d!(out, Vr(V), R, i, J),
            dVsite, nothing )   # dVsite == out, nothing == temp
      # compute the neighbour count gradients
      site_n_d!(dVn, V, r, R, Ns[i], dNs[i])

      # write site energy gradient into forces
      for n = 1:length(j)
         f = Ns[i] * dVsite[n] + Vs[i] * dVn[n]
         F[j[n]] -= f
         F[i] += f
      end
   end
   return F
end

function virial(V::AbstractEnvIP{N}, at::Atoms{T}) where {N, T}
   # compute the n values
   Ns, dNs = site_ns_ed(V, at)
   # compute the inner v values
   Vs = site_energies(Vr(V), at)

   # now assemble site forces and use those to create
   # total forces by mixing N and V
   cutoff_n = cutoff(Vn(V))
   nlist = neighbourlist(at, cutoff(V))  # this checks that cutoff(Vn) <= cutoff(Vr)
   maxneigs = max_neigs(nlist)
   dVsite = zeros(JVec{T}, maxneigs)
   dVn = zeros(JVec{T}, maxneigs)

   S = @SMatrix zeros(3,3)

   for (i, j, r, R) in sites(nlist)
      # compute the site energy gradients
      fill!(dVsite, zero(JVec{T}))
      eval_site_nbody!(
            Val(N), i, j, R, cutoff(V), false,
            (out, R, ii, J, temp) -> evaluate_d!(out, Vr(V), R, ii, J),
            dVsite, nothing )   # dVsite == out, nothing == temp
      # compute the neighbour count gradients
      site_n_d!(dVn, V, r, R, Ns[i], dNs[i])

      # convert the two into a single site energy gradient
      # so that we can call site_virial on it
      for n = 1:length(j)
         dVsite[n] = Ns[i] * dVsite[n] + Vs[i] * dVn[n]
      end

      S += site_virial(dVsite, R)
   end
   return S
end

function energy(B::Vector{TB}, at::Atoms{T}, typewarn=true
                ) where {TB <: EnvIP{N}, T} where {N}
   if typewarn
      !isconcretetype(TB) && @warn("TB is not a leaf type")
   end

   Br = [Vr(b) for b in B]

   rcut = cutoff(B[1])
   nlist = neighbourlist(at, rcut)
   temp = zeros(T, length(B))
   E = zeros(T, length(B))
   Etemp = zeros(T, length(B))

   # all the site energies => should be trivial in terms of cost
   Ns = [ site_ns(V, at) for V in B ]

   for (i, j, r, R) in sites(nlist)
      # evaluate all the site energies at the same time
      # for each simples, write the nB energies into temp
      # then add them to E, which is just passed through all the
      # various loops, so no need to update it here again
      fill!(Etemp, zero(T))
      eval_site_nbody!(Val(N), i, j, R, rcut, false,
                       (out, R, ii, J, temp) -> evaluate_many!(out, Br, R, ii, J),
                       Etemp, nothing)
      #
      for nb = 1:length(B)
         E[nb] += Etemp[nb] * Ns[nb][i]
      end
   end
   return E
end

function forces(B::AbstractVector{TB}, at::Atoms{T}, typewarn=true
              )where {TB <: EnvIP{N}, T} where {N}
   if typewarn
      !isconcretetype(TB) && @warn("TB is not a leaf type")
   end

   Br = [Vr(b) for b in B]

   rcut = cutoff(B[1])
   nlist = neighbourlist(at, rcut)
   maxneigs = max_neigs(nlist)
   nedges = (N*(N-1))÷2
   nB = length(B)
   # forces
   F =      [ zeros(JVec{T}, length(at)) for n = 1:nB ]
   # site gradient
   dVsite = [ zeros(JVec{T}, maxneigs)   for n = 1:nB ]

   # extras dfor Env
   dVn = zeros(JVec{T}, maxneigs)
   Etemp = zeros(T, length(B))

   # compute the N-components
   Ns = [ site_ns(V, at) for V in B ]
   dNs = [ site_ns_ed(V, at)[2] for V in B ]

   for (i, j, r, R) in sites(nlist)
      # clear dVsite and Etemp
      for n = 1:nB; fill!(dVsite[n], zero(JVec{T})); end
      fill!(Etemp, zero(T))
      # fill site energy and dVsite
      eval_site_nbody!(Val(N), i, j, R, rcut, false,
                       (out, R, ii, J, temp) -> evaluate_many!(out, Br, R, ii, J),
                       Etemp, nothing)
      eval_site_nbody!(Val(N), i, j, R, rcut, false,
                       (out, R, ii, J, temp) -> evaluate_many_d!(out, Br, R, ii, J),
                       dVsite, nothing)

      # write it into the force vectors
      for ib = 1:nB, n = 1:length(j)
         site_n_d!(dVn, B[ib], r, R, Ns[ib][i], dNs[ib][i])
         f = (Ns[ib][i] * dVsite[ib][n] + dVn[n] * Etemp[ib])
         F[ib][j[n]] -= f
         F[ib][i] += f
      end
   end
   return F
end

function virial(B::AbstractVector{TB}, at::Atoms{T}, typewarn=true
              )where {TB <: EnvIP{N}, T} where {N}
   if typewarn
      !isconcretetype(TB) && @warn("TB is not a leaf type")
   end

   Br = [Vr(b) for b in B]

   rcut = cutoff(B[1])
   nlist = neighbourlist(at, rcut)
   maxneigs = max_neigs(nlist)
   nedges = (N*(N-1))÷2
   nB = length(B)
   # virials (main output)
   S = fill((@SMatrix zeros(3,3)), length(B))
   # site gradient
   dVsite = [ zeros(JVec{T}, maxneigs)   for n = 1:length(B) ]

   # extras dfor Env
   dVn = zeros(JVec{T}, maxneigs)
   Etemp = zeros(T, length(B))

   # compute the N-components
   Ns = [ site_ns(V, at) for V in B ]
   dNs = [ site_ns_ed(V, at)[2] for V in B ]

   for (i, j, r, R) in sites(nlist)
      # clear dVsite and Etemp
      for n = 1:nB; fill!(dVsite[n], zero(JVec{T})); end
      fill!(Etemp, zero(T))
      # fill site energy and dVsite
      eval_site_nbody!(Val(N), i, j, R, rcut, false,
                       (out, R, ii, J, temp) -> evaluate_many!(out, Br, R, ii, J),
                       Etemp, nothing)
      eval_site_nbody!(Val(N), i, j, R, rcut, false,
                       (out, R, ii, J, temp) -> evaluate_many_d!(out, Br, R, ii, J),
                       dVsite, nothing)

      # use Etemp (site energies), the Ns and dVn to convert dVsite in
      # proper site energies
      for ib = 1:length(B)
         site_n_d!(dVn, B[ib], r, R, Ns[ib][i], dNs[ib][i])
         for n = 1:length(j)
            dVsite[ib][n] = Ns[ib][i] * dVsite[ib][n] + Etemp[ib] * dVn[n]
         end
      end

      # update the virials
      for iB = 1:length(B)
         S[iB] += site_virial(dVsite[iB], R)
      end
   end
   return S
end

# ============================================================
#     Re-Implement from Scratch for Faster `EnvPolys`
# ============================================================

function site_energies(V::EnvPoly{N}, at::Atoms{T}) where {N, T <: AbstractFloat}
   Es = zeros(T, length(at))

   # all the site energies => should be trivial in terms of cost
   Ns = site_energies(V.Vn, at)

   # storage
   Etemp = zeros(T, length(V.Vr))

   rcut = cutoff(V.Vr[1])
   for (i, j, r, R) in sites(at, rcut)
      fill!(Etemp, zero(T))
      eval_site_nbody!(Val(N), i, j, R, rcut, false,
                       (out, R, ii, J, temp) -> evaluate_many!(out, V.Vr, R, ii, J),
                       Etemp, nothing)
      for P = 1:length(V.Vr)
         Es[i] += Ns[i]^(P-1) * Etemp[P]
      end
   end

   return Es
end

function forces(V::EnvPoly{N}, at::Atoms{T}) where {N, T}
   rcut = cutoff(V.Vr[1])
   nlist = neighbourlist(at, rcut)
   maxneigs = max_neigs(nlist)
   nedges = (N*(N-1))÷2
   nVr = length(V.Vr)

   # forces
   F = zeros(JVec{T}, length(at))
   # site gradient
   dVsite = [ zeros(JVec{T}, maxneigs) for _ = 1:nVr ]

   # extras for Env
   Etemp = zeros(T, nVr)
   out = TempEdV(Etemp, dVsite)

   # compute the N-components
   Ns = site_energies(V.Vn, at)

   for (i, j, r, R) in sites(nlist)
      # clear dVsite and Etemp
      for n = 1:nVr; fill!(dVsite[n], zero(JVec{T})); end
      fill!(Etemp, zero(T))
      # fill site energy and dVsite
      eval_site_nbody!(Val(N), i, j, R, rcut, false,
                       (out, R, ii, J, temp) -> evaluate_many_ed!(out, V.Vr, R, ii, J),
                       out, nothing)

      # write it into the force vectors
      # first for P = 0
      for n = 1:length(j)
         F[j[n]] -= dVsite[1][n]
         F[i] += dVsite[1][n]
      end

      for P = 2:nVr, n = 1:length(j)
         # energy contribution:
         #   Es[i] += Ns[i]^(P-1) * Etemp[P]
         dVn = 0.5 * evaluate_d(V.Vn, r[n]) * R[n] / r[n]
         f = (  Ns[i]^(P-1) * dVsite[P][n]
              + (P-1) * Ns[i]^(P-2) * dVn * Etemp[P] )
         F[j[n]] -= f
         F[i] += f
      end
   end
   return F
end

# =================================================================
#     Re-Implement from Scratch again for Even Faster `EnvPolys`
#     with ClusterBLDesc
# =================================================================

# [!!!EXPERIMENTAL!!!]

using NBodyIPs: ClusterBLDesc
using NBodyIPs.Polys: NBPoly, StNBPoly

function many_to_sites!(Es, i, j, J, Eloc, Ns)
   for P = 1:length(Eloc)
      Es[i] += Ns[i]^(P-1) * Eloc[P]
      for a in J
         Es[j[a]] += Ns[j[a]]^(P-1) * Eloc[P]
      end
   end
   return Eloc
end

function site_energies(V::EnvPoly{N, TVR}, at::Atoms{T}
         ) where {TVR <: StNBPoly{N,TD}} where {N, T <: AbstractFloat, TD <: ClusterBLDesc}

   Es = zeros(T, length(at))
   # precompute the environment terms
   Ns = site_energies(V.Vn, at)
   # storage
   Etemp = zeros(T, length(V.Vr))

   rcut = cutoff(V.Vr[1])
   for (i, j, r, R) in sites(at, rcut)
      fill!(Etemp, zero(T))
      eval_site_nbody!(Val(N), i, j, R, rcut, true,
           (out, R, ii, J, temp) -> many_to_sites!(Es, ii, j, J, evaluate_many!(out, V.Vr, R, ii, J), Ns),
           Etemp, nothing)
   end

   return Es
end

# function forces(V::EnvPoly{N}, at::Atoms{T}) where {N, T}
#    rcut = cutoff(V.Vr[1])
#    nlist = neighbourlist(at, rcut)
#    maxneigs = max_neigs(nlist)
#    nedges = (N*(N-1))÷2
#    nVr = length(V.Vr)
#
#    # forces
#    F = zeros(JVec{T}, length(at))
#    # site gradient
#    dVsite = [ zeros(JVec{T}, maxneigs) for _ = 1:nVr ]
#
#    # extras for Env
#    Etemp = zeros(T, nVr)
#    out = TempEdV(Etemp, dVsite)
#
#    # compute the N-components
#    Ns = site_energies(V.Vn, at)
#
#    for (i, j, r, R) in sites(nlist)
#       # clear dVsite and Etemp
#       for n = 1:nVr; fill!(dVsite[n], zero(JVec{T})); end
#       fill!(Etemp, zero(T))
#       # fill site energy and dVsite
#       eval_site_nbody!(Val(N), i, j, R, rcut, false,
#                        (out, R, ii, J, temp) -> evaluate_many_ed!(out, V.Vr, R, ii, J),
#                        out, nothing)
#
#       # write it into the force vectors
#       # first for P = 0
#       for n = 1:length(j)
#          F[j[n]] -= dVsite[1][n]
#          F[i] += dVsite[1][n]
#       end
#
#       for P = 2:nVr, n = 1:length(j)
#          # energy contribution:
#          #   Es[i] += Ns[i]^(P-1) * Etemp[P]
#          dVn = 0.5 * evaluate_d(V.Vn, r[n]) * R[n] / r[n]
#          f = (  Ns[i]^(P-1) * dVsite[P][n]
#               + (P-1) * Ns[i]^(P-2) * dVn * Etemp[P] )
#          F[j[n]] -= f
#          F[i] += f
#       end
#    end
#    return F
# end

ext: jl

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
################################################################################
#
#  Ringmorphisms of orders onto finite fields
#
################################################################################

# NfOrd -> FqNmod
mutable struct NfOrdToFqNmodMor <: Map{NfOrd, FqNmodFiniteField, HeckeMap, NfOrdToFqNmodMor}
  header::MapHeader{NfOrd, FqNmodFiniteField}
  poly_of_the_field::gfp_poly
  P::NfOrdIdl
  powers::Vector{nf_elem}

  function NfOrdToFqNmodMor()
    r = new()
    r.header = MapHeader{NfOrd, FqNmodFiniteField}()
    return r
  end

  function NfOrdToFqNmodMor(O::NfOrd, F::FqNmodFiniteField, g::gfp_poly)
    # assume that F = F_p[X]/(g) and g is a factor of f mod p

    z = new()
    d = degree(nf(O))
    p = characteristic(F)
    a = gen(nf(O))
    tmp_gfp_poly = parent(g)()
    z.poly_of_the_field = g

    powers = Vector{nf_elem}(undef, d)

    powers[1] = a
    for i in 2:d
      powers[i] = powers[i - 1] * a
    end
    z.powers = powers

    function _image(x::NfOrdElem)
      u = F()
      gg = parent(nf(O).pol)(elem_in_nf(x))::fmpq_poly
      fmpq_poly_to_gfp_poly_raw!(tmp_gfp_poly, gg)
      ccall((:nmod_poly_rem, libflint), Nothing,
            (Ref{gfp_poly}, Ref{gfp_poly}, Ref{gfp_poly}, Ptr{Nothing}),
            tmp_gfp_poly, tmp_gfp_poly, g, pointer_from_objref(F)+sizeof(fmpz))
      ccall((:fq_nmod_set, libflint), Nothing,
            (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}),
            u, tmp_gfp_poly, F)
      ccall((:fq_nmod_reduce, libflint), Nothing,
            (Ref{fq_nmod}, Ref{FqNmodFiniteField}), u, F)
      return u
    end

    # The lift is even simpler!
    function _preimage(y::fq_nmod)
      zz = O()
      zz.elem_in_nf = nf(O)(coeff(y, 0))
      for i in 2:d
        add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
      end
      zz.elem_in_nf = mod(zz.elem_in_nf, p)
      return zz
    end

    z.header = MapHeader{NfOrd, FqNmodFiniteField}(O, F, _image, _preimage)
    return z
  end

  function NfOrdToFqNmodMor(O::NfOrd, P::NfOrdIdl)
    z = NfOrdToFqNmodMor()
    z.P = P
    p = minimum(P)
    a, g, b = get_residue_field_data(P)
    psmall = Int(p)
    R = GF(psmall, cached = false)
    Rx, x = PolynomialRing(R, "_\$", cached = false)
    F = FqNmodFiniteField(Rx(g), Symbol("_\$"), false)
    d = degree(g)
    n = degree(O)
    imageofbasis = Vector{fq_nmod}(undef, n)
    powers = Vector{nf_elem}(undef, d)
    c = Rx()
    for i in 1:n
      ib = F()
      @assert d == ncols(b[i])
      for j in 1:d
        setcoeff!(c, j - 1, b[i][1, j])
      end
      ccall((:fq_nmod_set, libflint), Nothing, (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}), ib, c, F)
      imageofbasis[i] = ib
    end

    powers[1] = a.elem_in_nf
    for i in 2:d
      powers[i] = powers[i - 1] * a.elem_in_nf
    end
    z.powers = powers

    tempF = F()

    function _image(x::NfOrdElem)
      v = coordinates(x, copy = false)
      zz = zero(F)
      for i in 1:n
        ccall((:fq_nmod_mul_fmpz, libflint), Nothing,
              (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}),
              tempF, imageofbasis[i], v[i], F)
        add!(zz, zz, tempF)
      end
      return zz
    end

    function _preimage(y::fq_nmod)
      zz = O()
      zz.elem_in_nf = nf(O)(coeff(y, 0))
      for i in 2:d
        add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
      end
      zz.elem_in_nf = mod(zz.elem_in_nf, p)
      return zz
    end

    z.header = MapHeader{NfOrd, FqNmodFiniteField}(O, F, _image, _preimage)
    return z
  end
end

function preimage(f::NfOrdToFqNmodMor, y::fq_nmod)
  O = domain(f)
  p = minimum(f.P)
  powers = f.powers
  d = length(powers)
  zz = O()
  zz.elem_in_nf = nf(O)(coeff(y, 0))
  for i in 2:d
    add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
  end
  zz.elem_in_nf = mod(zz.elem_in_nf, p)
  return zz
end

# S is the type of the order, T the type of the ideal and U the elem_type of the order, which define the quotient ring
mutable struct AbsOrdQuoMap{S, T, U} <: Map{S, AbsOrdQuoRing{S, T}, HeckeMap, AbsOrdQuoMap}
  header::MapHeader{S, AbsOrdQuoRing{S, T}}

  function AbsOrdQuoMap{S, T, U}(O::S, Q::AbsOrdQuoRing{S, T}) where {S, T, U}
    z = new()

    function _image(x::U)
      return Q(x)
    end

    function _image(x::FacElem)
      return mod(x, Q)
    end

    function _preimage(x::AbsOrdQuoRingElem{S, T, U})
      x.elem = mod(x.elem, parent(x))
      return x.elem
    end

    z.header = MapHeader(O, Q, _image, _preimage)
    return z
  end
end

function AbsOrdQuoMap(O::S, Q::AbsOrdQuoRing{S, T}) where {S, T}
  U = elem_type(O)
  return AbsOrdQuoMap{S, T, U}(O, Q)
end

const NfOrdQuoMap = AbsOrdQuoMap{NfOrd, NfOrdIdl, NfOrdElem}

function Mor(O::NfOrd, F::FqNmodFiniteField, y::fq_nmod)
  return NfOrdToFqNmodMor(O, F, y)
end

function Mor(O::NfOrd, F::FqNmodFiniteField, h::gfp_poly)
  return NfOrdToFqNmodMor(O, F, h)
end

function evaluate(f::fmpz_poly, r::fq_nmod)
  #Horner - stolen from Claus

  if length(f) == 0
    return parent(r)()
  end

  l = f.length-1
  s = parent(r)(coeff(f, l))
  for i =l-1:-1:0
    s = s*r + parent(r)(coeff(f, i))
  end
  return s
end

function evaluate!(z::fq_nmod, f::fmpz_poly, r::fq_nmod)
  #Horner - stolen from Claus

  zero!(z)

  if length(f) == 0
    return z
  end

  l = f.length-1
  set!(z, parent(r)(coeff(f, l)))
  #s = parent(r)(coeff(f, l))
  for i =l-1:-1:0
    mul!(z, z, r)
    add!(z, z, parent(r)(coeff(f, i)))
    #s = s*r + parent(r)(coeff(f, i))
  end
  return z
end

function _get_coeff_raw(x::fq_nmod, i::Int)
  u = ccall((:nmod_poly_get_coeff_ui, libflint), UInt, (Ref{fq_nmod}, Int), x, i)
  return u
end

function _get_coeff_raw(x::fq, i::Int)
  t = fmpz()
  ccall((:fmpz_poly_get_coeff_fmpz, libflint), Nothing, (Ref{fmpz}, Ref{fq}, Int), t, x, i)
  return t
end

# this is super slow
# improve
function (f::NfOrdQuoMap)(I::NfOrdIdl)
  O = domain(f)
  Q = codomain(f)
  B = Q.ideal + I
  nB = norm(B, copy = false)
  if isone(nB)
    return one(Q)
  end
  _assure_weakly_normal_presentation(B)
  nQ = norm(Q.ideal, copy = false)
  if _normmod(nQ, B.gen_two) == nB
    return Q(B.gen_two)
  end
  b = basis(B, copy = false)
  range1nQ2 = fmpz(1):nQ^2
  z = O()
  while true
    z = rand!(z, b, range1nQ2)
    if _normmod(nQ, z) == nB
      break
    end
  end
  return Q(z)
end

function (f::NfOrdQuoMap)(p::PolyElem{NfOrdElem})
  F = codomain(f)
  Fx,_ = PolynomialRing(F, "_\$", cached = false)

  ar = NfOrdElem[ coeff(p, i) for i in 0:degree(p) ]

  z = Fx(map(f, ar))

  return z
end

mutable struct NfOrdToFqMor <: Map{NfOrd, FqFiniteField, HeckeMap, NfOrdToFqMor}
  header::MapHeader{NfOrd, FqFiniteField}
  poly_of_the_field::gfp_fmpz_poly
  P::NfOrdIdl
  fastpath::Bool
  # Some temporary variables
  tmp_gfp_fmpz_poly::gfp_fmpz_poly
  t_fmpz_poly::fmpz_poly
  t_fmpz::fmpz
  a::nf_elem

  function NfOrdToFqMor()
    z = new()
    return z
  end

  function NfOrdToFqMor(O::NfOrd, F::FqFiniteField, g::gfp_fmpz_poly)
    # assume that F = F_p[X]/(g) and g is a factor of f mod p

    z = new()
    z.fastpath = true
    p = characteristic(F)
    z.tmp_gfp_fmpz_poly = parent(g)()
    z.t_fmpz_poly = fmpz_poly()
    z.t_fmpz = fmpz()

    z.a = gen(nf(O))
    z.poly_of_the_field = g

    z.header = MapHeader{NfOrd, FqFiniteField}(O, F)# _image, _preimage)

    return z
  end
end

function NfOrdToFqMor(O::NfOrd, P::NfOrdIdl)#, g::fmpz_poly, a::NfOrdElem, b::Vector{fmpz_mat})
  z = NfOrdToFqMor()
  z.fastpath = false
  z.P = P
  a, g, b = get_residue_field_data(P)
  p = minimum(P)
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "_\$", cached = false)
  F = FqFiniteField(Rx(g), Symbol("_\$"), false)
  d = degree(g)
  n = degree(O)
  imageofbasis = Vector{fq}(undef, n)
  powers = Vector{nf_elem}(undef, d)
  c = Rx()

  for i in 1:n
    ib = F()
    @assert d == ncols(b[i])
    for j in 1:d
      setcoeff!(c, j - 1, b[i][1, j])
    end
    ccall((:fq_set, libflint), Nothing, (Ref{fq}, Ref{gfp_fmpz_poly}, Ref{FqFiniteField}), ib, c, F)
    imageofbasis[i] = ib
  end

  powers[1] = a.elem_in_nf
  for i in 2:d
    powers[i] = powers[i - 1] * a.elem_in_nf
  end

  tempF = F()

  function _image(x::NfOrdElem)
    v = coordinates(x, copy = false)
    zz = zero(F)
    for i in 1:n
      ccall((:fq_mul_fmpz, libflint), Nothing,
            (Ref{fq}, Ref{fq}, Ref{fmpz}, Ref{FqFiniteField}),
            tempF, imageofbasis[i], v[i], F)
      add!(zz, zz, tempF)
    end
    return zz
  end

  function _preimage(y::fq)
    zz = O()
    zz.elem_in_nf = nf(O)(coeff(y, 0))
    for i in 2:d
      add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
    end
    zz.elem_in_nf = mod(zz.elem_in_nf, p)
    return zz
  end

  z.header = MapHeader{NfOrd, FqFiniteField}(O, F, _image, _preimage)
  return z
end

function image(f::NfOrdToFqMor, x::NfOrdElem)
  if f.fastpath
    F = codomain(f)
    O = domain(f)
    u = F()
    gg = parent(nf(O).pol)(elem_in_nf(x, copy = false))::fmpq_poly
    fmpq_poly_to_gfp_fmpz_poly_raw!(f.tmp_gfp_fmpz_poly, gg, f.t_fmpz_poly, f.t_fmpz)
    ccall((:fmpz_mod_poly_rem, libflint), Nothing, (Ref{gfp_fmpz_poly}, Ref{gfp_fmpz_poly}, Ref{gfp_fmpz_poly}, Ref{Nemo.fmpz_mod_ctx_struct}), f.tmp_gfp_fmpz_poly, f.tmp_gfp_fmpz_poly, f.poly_of_the_field, f.tmp_gfp_fmpz_poly.parent.base_ring.ninv)
    ccall((:fq_set, libflint), Nothing, (Ref{fq}, Ref{gfp_fmpz_poly}, Ref{FqFiniteField}), u, f.tmp_gfp_fmpz_poly, F)
    ccall((:fq_reduce, libflint), Nothing, (Ref{fq}, Ref{FqFiniteField}), u, F)
    return u
  else
    return f.header.image(x)::fq
  end
end

function preimage(f::NfOrdToFqMor, x::fq)
  if f.fastpath
    O = domain(f)
    F = codomain(f)
    zz = nf(O)()

    a = f.a
    # TODO: Do something more clever here
    for i in 0:degree(F)-1
      zz = zz + coeff(x, i)*a^i
    end

    return O(zz, false)::NfOrdElem
  else
    @assert isdefined(f.header, :preimage)
    return f.header.preimage(x)::NfOrdElem
  end
end

function Mor(O::NfOrd, F::FqFiniteField, h::gfp_fmpz_poly)
  return NfOrdToFqMor(O, F, h)
end

################################################################################
#
#  ResidueField degree 1 primes
#
################################################################################

mutable struct NfOrdToGFMor <: Map{NfOrd, GaloisField, HeckeMap, NfOrdToFqNmodMor}
  header::MapHeader{NfOrd, GaloisField}
  poly_of_the_field::gfp_poly
  P::NfOrdIdl

  function NfOrdToGFMor()
    r = new()
    r.header = MapHeader{NfOrd, GaloisField}()
    return r
  end

  function NfOrdToGFMor(O::NfOrd, F::GaloisField, g::gfp_poly)
    # assume that F = F_p[X]/(g) and g is a factor of f mod p of degree 1

    z = new()
    tmp_gfp_poly = parent(g)()
    z.poly_of_the_field = g
		local _image
    let g = g, tmp_gfp_poly = tmp_gfp_poly, O = O, F = F
    	function _image(x::NfOrdElem)
      	gg = parent(nf(O).pol)(elem_in_nf(x))::fmpq_poly
      	fmpq_poly_to_gfp_poly_raw!(tmp_gfp_poly, gg)
      	ccall((:nmod_poly_rem, libflint), Nothing,
        	    (Ref{gfp_poly}, Ref{gfp_poly}, Ref{gfp_poly}, Ptr{Nothing}),
          	  tmp_gfp_poly, tmp_gfp_poly, g, pointer_from_objref(F)+sizeof(fmpz))
      	return coeff(tmp_gfp_poly, 0)
			end
    end

    z.header = MapHeader{NfOrd, GaloisField}(O, F, _image)
    return z
  end

  function NfOrdToGFMor(O::NfOrd, P::NfOrdIdl)
    z = NfOrdToGFMor()
    z.P = P
    p = minimum(P)
    a, g, b = get_residue_field_data(P)
    psmall = Int(p)
		n = degree(O)
    F = GF(psmall, cached = false)
    imageofbasis = Vector{gfp_elem}(undef, n)
    for i in 1:n
      imageofbasis[i] = F(b[i][1, 1])
    end

		local _image
		let imageofbasis = imageofbasis, F = F, n = n
   		function _image(x::NfOrdElem)
      	v = coordinates(x, copy = false)
      	tempF = zero(UInt)
      	for i in 1:n
          tempF += mul_mod(imageofbasis[i].data, v[i], F)
       	end
      	return F(tempF)
			end
    end

    z.header = MapHeader{NfOrd, GaloisField}(O, F, _image)
    return z
  end
end

function preimage(f::NfOrdToGFMor, a::gfp_elem)
  return domain(f)(a.data)
end

Mor(O::NfOrd, F::GaloisField, g::gfp_poly) = NfOrdToGFMor(O, F, g)

mutable struct NfOrdToGFFmpzMor <: Map{NfOrd, Nemo.GaloisFmpzField, HeckeMap, NfOrdToGFFmpzMor}
  header::MapHeader{NfOrd, Nemo.GaloisFmpzField}
  poly_of_the_field::gfp_fmpz_poly
  P::NfOrdIdl

  function NfOrdToGFFmpzMor()
    r = new()
    return r
  end

  function NfOrdToGFFmpzMor(O::NfOrd, F::Nemo.GaloisFmpzField, g::gfp_fmpz_poly)
    # assume that F = F_p[X]/(g) and g is a factor of f mod p of degree 1

    z = new()
    tmp_gfp_poly = parent(g)()
    z.poly_of_the_field = g
		local _image
    let g = g, tmp_gfp_poly = tmp_gfp_poly, O = O, F = F
    	function _image(x::NfOrdElem)
      	gg = parent(nf(O).pol)(elem_in_nf(x))::fmpq_poly
      	fmpq_poly_to_gfp_fmpz_poly_raw!(tmp_gfp_poly, gg)
				rem!(tmp_gfp_poly, tmp_gfp_poly, g)
      	return coeff(tmp_gfp_poly, 0)
			end
    end

    z.header = MapHeader{NfOrd, Nemo.GaloisFmpzField}(O, F, _image)
    return z
  end

  function NfOrdToGFFmpzMor(O::NfOrd, P::NfOrdIdl)
    z = NfOrdToGFFmpzMor()
    z.P = P
    p = minimum(P)
    a, g, b = get_residue_field_data(P)
		n = degree(O)
    F = GF(p, cached = false)
    imageofbasis = Vector{Nemo.gfp_fmpz_elem}(undef, n)
    for i in 1:n
      imageofbasis[i] = F(b[i][1, 1])
    end

    tempF = F()
		local _image
		let tempF = tempF, imageofbasis = imageofbasis, F = F, n = n
   		function _image(x::NfOrdElem)
      	v = coordinates(x, copy = false)
      	zz = zero(F)
      	for i in 1:n
          mul!(tempF, imageofbasis[i], v[i])
        	add!(zz, zz, tempF)
      	end
      	return zz
			end
    end

    z.header = MapHeader{NfOrd, Nemo.GaloisFmpzField}(O, F, _image)
    return z
  end
end

function preimage(f::NfOrdToGFFmpzMor, a::Nemo.gfp_fmpz_elem)
  return domain(f)(lift(a))
end

Mor(O::NfOrd, F::Nemo.GaloisFmpzField, h::gfp_fmpz_poly) = NfOrdToGFFmpzMor(O, F, h)

################################################################################
#
#  Extend to number field
#
################################################################################

mutable struct NfToFinFldMor{T} <: Map{AnticNumberField, T, HeckeMap, NfToFinFldMor{T}}
  header::MapHeader{AnticNumberField, T}

  function NfToFinFldMor{T}() where T
    r = new{T}()
    r.header = MapHeader{AnticNumberField, T}()
    return r
  end
end

function extend(f::T, K::AnticNumberField) where T <: Union{NfOrdToFqNmodMor, NfOrdToFqMor, NfOrdToGFMor, NfOrdToGFFmpzMor}
  nf(domain(f)) != K && error("Number field is not the number field of the order")

  z = NfToFinFldMor{typeof(codomain(f))}()

  z.header.domain = K
  z.header.codomain = f.header.codomain

  pia = anti_uniformizer(f.P)
  O = domain(f)
  P = f.P

  _image = function(x::nf_elem)
    !iszero(x) && valuation(x, P) < 0 && error("Element not p-integral")
    m = denominator(x, domain(f))
    l = valuation(m, P)
    if l == 0
      return f(O(m*x, false))//f(O(m))
    else
      pial = pia^l
      pialm = pial * m
      w = pialm * x
      return f(O(w, false))//f(O(pialm, false))
    end
  end

  _preimage = function(x)
    return elem_in_nf(preimage(f, x))
  end

  z.header.image = _image
  z.header.preimage = _preimage

  return z
end

#=
function (f::Union{NfOrdToFqNmodMor, NfOrdToFqMor, NfOrdToGFMor, NfOrdToGFFmpzMor})(p::PolyElem{NfOrdElem})
  return map_coefficients(f, p)
end
=#
@doc Markdown.doc"""
    extend_easy(f::Hecke.NfOrdToFqNmodMor, K::AnticNumberField) -> NfToFqNmodMor

For a residue field map from a prime ideal, extend the domain of the map
to the entire field.
Requires the prime ideal to be coprime to the index, unramified and
over a small integer. The resulting map can very efficiently be
evaluated using `image(map, elem)`.
The resulting map can be applied to
  * `nf_elem`
  * `FacElem{nf_elem}`
Will throw a `BadPrime` exception if applied to an element in the
field with a $p$ in the denominator. In the case of `FacElem`, zero
is also not permitted (and will produce a `BadPrime` error).
"""
function extend_easy(f::Hecke.NfOrdToFqNmodMor, K::AnticNumberField)
  nf(domain(f)) != K && error("Number field is not the number field of the order")
  return NfToFqNmodMor_easy(f, K)
end

#a stop-gap, mainly for non-monic polynomials
function extend_easy(f::Hecke.NfOrdToFqMor, K::AnticNumberField)
  return NfToFqMor_easy(f, K)
end

function extend_easy(f::Hecke.NfOrdToGFMor, K::AnticNumberField)
  return NfToGFMor_easy(f, K)
end

function extend_easy(f::Hecke.NfOrdToGFFmpzMor, K::AnticNumberField)
  return NfToGFFmpzMor_easy(f, K)
end

mutable struct NfToFqMor_easy <: Map{AnticNumberField, FqFiniteField, HeckeMap, NfToFqMor_easy}
  header::MapHeader{AnticNumberField, FqFiniteField}
  Fq::FqFiniteField
  s::fq
  t::gfp_fmpz_poly
  function NfToFqMor_easy(a::Map, k::AnticNumberField)
    r = new()
    r.Fq = codomain(a)
    r.header = MapHeader(k, r.Fq)
    r.s = r.Fq()
    r.t = PolynomialRing(GF(characteristic(r.Fq), cached = false), cached = false)[1]()
    return r
  end
end

function image(mF::NfToFqMor_easy, a::FacElem{nf_elem, AnticNumberField}, quo::Int = 0)
  Fq = mF.Fq
  q = one(Fq)
  t = mF.t
  s = mF.s
  for (k, v) = a.fac
    vv = v
    if quo != 0
      vv = v % quo
      if vv < 0
        vv += quo
      end
    end
    @assert vv < order(Fq)  #please complain if this is triggered
    if !iszero(vv)
      if denominator(k) % characteristic(Fq) == 0
        throw(BadPrime(characteristic(Fq)))
      end
      _nf_to_fq!(s, k, Fq, t)
      if iszero(s)
        throw(BadPrime(1))
      end
      if vv < 0
        ccall((:fq_inv, libflint), Nothing, (Ref{fq}, Ref{fq}, Ref{FqFiniteField}), s, s, Fq)
        vv = -vv
      end
      ccall((:fq_pow_ui, libflint), Nothing, (Ref{fq}, Ref{fq}, Int, Ref{FqFiniteField}), s, s, vv, Fq)
      mul!(q, q, s)
    end
  end
  return q
end

function image(mF::NfToFqMor_easy, a::nf_elem, n_quo::Int = 0)
  Fq = mF.Fq
  q = Fq()
  if denominator(a) % characteristic(Fq) == 0
    throw(BadPrime(characteristic(Fq)))
  end
  _nf_to_fq!(q, a, Fq, mF.t)
  return q
end

mutable struct NfToFqNmodMor_easy <: Map{AnticNumberField, FqNmodFiniteField, HeckeMap, NfToFqNmodMor_easy}
  header::MapHeader{AnticNumberField, FqNmodFiniteField}
  Fq::FqNmodFiniteField
  s::fq_nmod
  t::gfp_poly
  function NfToFqNmodMor_easy(a::Map, k::AnticNumberField)
    r = new()
    r.Fq = codomain(a)
    r.header = MapHeader(k, r.Fq)
    r.s = r.Fq()
    r.t = PolynomialRing(GF(UInt(characteristic(r.Fq)), cached=false), cached=false)[1]()
    return r
  end
end

function image(mF::NfToFqNmodMor_easy, a::FacElem{nf_elem, AnticNumberField}, quo::Int = 0)
  Fq = mF.Fq
  q = one(Fq)
  t = mF.t
  s = mF.s
  oFq = order(Fq) # fmpz
  small_mod = UInt(0)
  char_Fq = characteristic(Fq)

  if quo != 0
    small_mod = UInt(quo)
  end

  for (k, v) = a.fac
    # I want to map k^v to F. I can reduce mod q (reduction modulo q - 1 is
    # done by the power function itself.

    inver = false

    if v < 0
      v = -v
      inver = true
    end

    local vv::UInt

    if quo != 0
      if v > small_mod
        vv = fmpz_mod_ui(v, small_mod)
      else
        vv = UInt(v)
      end
    end

    if (quo != 0 && vv != 0) || !iszero(v)
      # We have something to do
      if iszero(denominator(k) % char_Fq)
        throw(BadPrime(characteristic(Fq)))
      end
      _nf_to_fq!(s, k, Fq, t)
      if iszero(s)
        throw(BadPrime(1))
      end
      if inver
        ccall((:fq_nmod_inv, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), s, s, Fq)
      end
      if quo != 0
        ccall((:fq_nmod_pow_ui, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Int, Ref{FqNmodFiniteField}), s, s, vv, Fq)
      else
        ccall((:fq_nmod_pow, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}), s, s, v, Fq)
      end
      mul!(q, q, s)
    end
  end
  return q
end

function image(mF::NfToFqNmodMor_easy, a::FacElem{nf_elem, AnticNumberField}, D::Vector, cached::Bool, quo::Int = 0)
  Fq = mF.Fq
  q = one(Fq)
  t = mF.t
  s = mF.s
  i = 0
  char = UInt(Fq.n)
  small_mod = UInt(0)

  if quo != 0
    small_mod = UInt(quo)
  end

  for (k, v) in a.fac
    i += 1

    inver = false

    if v < 0
      v = -v
      inver = true
    end

    local vv::UInt

    if quo != 0
      if v > small_mod
        vv = fmpz_mod_ui(v, small_mod)
      else
        vv = UInt(v)
      end
    end

    # We always reduce, so do the test first
    if !cached && iszero(fmpz_mod_ui(denominator(k), char))
      throw(BadPrime(characteristic(Fq)))
    end

    if !((quo != 0 && vv != 0) || !iszero(v))
      if !cached
        nf_elem_to_gfp_poly!(t, k)
        D[i] = zero(parent(t))
        set!(D[i], t)
      end
    end

    if (quo != 0 && vv != 0) || !iszero(v)
      if cached
        s = zero(Fq)
        ccall((:fq_nmod_set, libflint), Nothing,
          (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}), s, D[i], Fq)
        _reduce(s)
      else
        nf_elem_to_gfp_poly!(t, k)
        #tt = deepcopy(t)
        if isassigned(D, i)
          y = D[i]
          if y.mod_n == t.mod_n
            y.parent = t.parent
            set!(y, t)
          else
            y.mod_n = t.mod_n
            y.mod_ninv = t.mod_ninv
            y.mod_norm = t.mod_norm
            y.parent = t.parent
            set!(y, t)
          end
        else
          D[i] = zero(parent(t))
          set!(D[i], t)
        end
        s = zero(Fq)
        ccall((:fq_nmod_set, libflint), Nothing,
          (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}), s, D[i], Fq)
        _reduce(s)
      end
      if iszero(s)
        throw(BadPrime(1))
      end

      if inver
        ccall((:fq_nmod_inv, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), s, s, Fq)
      end
      if quo != 0
        ccall((:fq_nmod_pow_ui, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Int, Ref{FqNmodFiniteField}), s, s, vv, Fq)
      else
        ccall((:fq_nmod_pow, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}), s, s, v, Fq)
      end
      mul!(q, q, s)

      #if vv < 0
      #  s = inv(s)
      #  vv = -vv
      #end
      #s = s^Int(vv)
      #q = mul!(q, q, s)
    end
  end
  return q
end

function image(mF::NfToFqNmodMor_easy, a::nf_elem, n_quo::Int = 0)
  Fq = mF.Fq
  q = Fq()
  if denominator(a) % characteristic(Fq) == 0
    throw(BadPrime(characteristic(Fq)))
  end
  _nf_to_fq!(q, a, Fq, mF.t)
  return q
end

function _nf_to_gfp_elem(b::nf_elem, a_tmp::gfp_poly, def_pol::gfp_poly)
  nf_elem_to_gfp_poly!(a_tmp, b)
  rem!(a_tmp, a_tmp, def_pol)
  return coeff(a_tmp, 0)
end

function _nf_to_gfp_elem(b::nf_elem, a_tmp::gfp_fmpz_poly, def_pol::gfp_fmpz_poly)
  nf_elem_to_gfp_fmpz_poly!(a_tmp, b)
  rem!(a_tmp, a_tmp, def_pol)
  return coeff(a_tmp, 0)
end

mutable struct NfToGFMor_easy <: Map{AnticNumberField, GaloisField, HeckeMap, NfToGFMor_easy}
  header::MapHeader{AnticNumberField, GaloisField}
  Fq::GaloisField
  defining_pol::gfp_poly
  s::gfp_elem
  t::gfp_poly
  function NfToGFMor_easy(a::NfOrdToGFMor, k::AnticNumberField)
    r = new()
    r.Fq = codomain(a)
    r.defining_pol = a.poly_of_the_field
    r.header = MapHeader(k, r.Fq)
    r.s = r.Fq()
    r.t = zero(parent(r.defining_pol))
    return r
  end
end

function image(mF::NfToGFMor_easy, a::FacElem{nf_elem, AnticNumberField}, quo::Int = 0)
  Fq = mF.Fq
  p = mF.defining_pol
  q = one(Fq)
  t = mF.t
  for (k, v) = a.fac
    vv = v
    if quo != 0
      vv = v %quo
      if vv < 0
        vv += quo
      end
    end
    @assert vv < order(Fq)  #please complain if this is triggered
    if !iszero(vv)
      if denominator(k) % characteristic(Fq) == 0
        throw(BadPrime(characteristic(Fq)))
      end
      s = _nf_to_gfp_elem(k, t, p)
      if iszero(s)
        throw(BadPrime(1))
      end
      if vv < 0
        s = inv(s)
        vv = -vv
      end
      s = s^vv
      q = mul!(q, q, s)
    end
  end
  return q
end

function image(mF::NfToGFMor_easy, a::FacElem{nf_elem, AnticNumberField}, D::Vector, cached::Bool, quo::Int = 0)
  # cached == true also implies that all the denominators are coprime
  Fq = mF.Fq
  p = mF.defining_pol
  q = one(Fq)
  t = mF.t
  i = 0
  pminusone = Fq.n - 1
  @assert ismonic(p)
  evaluateat = -coeff(p, 0)
  for (k, v) in a.fac
    i += 1
    if v > 0 && v < pminusone
      vv = UInt(v)
    else
      vv = fmpz_mod_ui(v, pminusone)
    end
    if quo != 0
      vv = vv % quo # vv will always be positive
    end
    @assert vv < Fq.n  #please complain if this is triggered

    # We always have to reduce, so check first
    if !cached && (fmpz_mod_ui(denominator(k), Fq.n) == 0)
      throw(BadPrime(characteristic(Fq)))
    end

    if iszero(vv) && !cached
      D[i] = zero(parent(t))
      nf_elem_to_gfp_poly!(t, k)
      set!(D[i], t)
    end

    if !iszero(vv)
      if cached
        s = evaluate_raw(D[i], evaluateat)
      else
        nf_elem_to_gfp_poly!(t, k)
        #tt = deepcopy(t)
        if isassigned(D, i)
          y = D[i]
          if y.mod_n == t.mod_n
            y.parent = t.parent
            set!(y, t)
          else
            y.mod_n = t.mod_n
            y.mod_ninv = t.mod_ninv
            y.mod_norm = t.mod_norm
            y.parent = t.parent
            set!(y, t)
          end
        else
          D[i] = zero(parent(t))
          set!(D[i], t)
        end
        s = evaluate(t, evaluateat)
      end
      #s = _nf_to_gfp_elem(k, t, p)
      if iszero(s)
        throw(BadPrime(1))
      end
      if vv < 0
        s = inv(s)
        vv = -vv
      end
      s = s^Int(vv)
      q = mul!(q, q, s)
    end
  end
  return q
end

function image(mF::NfToGFMor_easy, a::nf_elem, n_quo::Int = 0)
  Fq = mF.Fq
  p = mF.defining_pol
  t = mF.t
  if denominator(a) % characteristic(Fq) == 0
    throw(BadPrime(characteristic(Fq)))
  end
  return _nf_to_gfp_elem(a, t, p)
end

function image(mF::NfToGFMor_easy, a::nf_elem, D::Vector, cached::Bool, n_quo::Int = 0)
  Fq = mF.Fq
  p = mF.defining_pol
  t = mF.t

  @assert ismonic(p)
  evaluateat = -coeff(p, 0)

  if denominator(a) % characteristic(Fq) == 0
    throw(BadPrime(characteristic(Fq)))
  end
  if cached
    @assert length(D) == 1
    s = evaluate_raw(D[1], evaluateat)
    #rem!(t, D[1], p)
    #s = coeff(t, 0)
  else
    nf_elem_to_gfp_poly!(t, a)
    D[1] = deepcopy(t)
    #rem!(t, t, p)
    #s = coeff(t, 0)
    s = evaluate_raw(t, evaluateat)
  end
  return s
end

mutable struct NfToGFFmpzMor_easy <: Map{AnticNumberField, Nemo.GaloisFmpzField, HeckeMap, NfToGFFmpzMor_easy}
  header::MapHeader{AnticNumberField, Nemo.GaloisFmpzField}
  Fq::Nemo.GaloisFmpzField
  defining_pol::Nemo.gfp_fmpz_poly
  s::Nemo.gfp_fmpz_elem
  t::Nemo.gfp_fmpz_poly
  function NfToGFFmpzMor_easy(a::NfOrdToGFFmpzMor, k::AnticNumberField)
    r = new()
    r.Fq = codomain(a)
    r.header = MapHeader(k, r.Fq)
    r.s = r.Fq()
    r.defining_pol = a.poly_of_the_field
    r.t = zero(parent(a.poly_of_the_field))
    return r
  end
end

function image(mF::NfToGFFmpzMor_easy, a::FacElem{nf_elem, AnticNumberField}, quo::Int = 0)
  Fq = mF.Fq
  p = mF.defining_pol
  q = one(Fq)
  t = mF.t
  for (k, v) = a.fac
    vv = v
    if quo != 0
      vv = v %quo
      if vv < 0
        vv += quo
      end
    end
    @assert vv < order(Fq)  #please complain if this is triggered
    if !iszero(vv)
      if denominator(k) % characteristic(Fq) == 0
        throw(BadPrime(characteristic(Fq)))
      end
      s = _nf_to_gfp_fmpz_elem(k, t, p)
      if iszero(s)
        throw(BadPrime(1))
      end
      if vv < 0
        s = inv(s)
        vv = -vv
      end
      s = s^vv
      mul!(q, q, s)
    end
  end
  return q
end

function image(mF::NfToGFFmpzMor_easy, a::nf_elem, n_quo::Int = 0)
  Fq = mF.Fq
  p = mF.defining_pol
  t = mF.t
  if denominator(a) % characteristic(Fq) == 0
    throw(BadPrime(characteristic(Fq)))
  end
  return _nf_to_gfp_fmpz_elem(a, t, p)
end

################################################################################
#
#  AbsOrdToAlgAssMor type
#
################################################################################

mutable struct AbsOrdToAlgAssMor{S, T} <: Map{S, AlgAss{T}, HeckeMap, AbsOrdToAlgAssMor}
  header::MapHeader

  function AbsOrdToAlgAssMor{S, T}(O::S, A::AlgAss{T}, _image::Function, _preimage::Function) where {S <: Union{ NfAbsOrd, AlgAssAbsOrd }, T}
    z = new{S, T}()
    z.header = MapHeader(O, A, _image, _preimage)
    return z
  end
end

function AbsOrdToAlgAssMor(O::Union{ NfAbsOrd, AlgAssAbsOrd }, A::AlgAss{T}, _image, _preimage) where {T}
  return AbsOrdToAlgAssMor{typeof(O), T}(O, A, _image, _preimage)
end

# Helper

function mul!(z::gfp_elem, x::gfp_elem, y::fmpz)
  R = parent(x)
  d = ccall((:fmpz_fdiv_ui, libflint), UInt, (Ref{fmpz}, UInt), y, R.n)
  r = ccall((:n_mulmod2_preinv, libflint), UInt, (UInt, UInt, UInt, UInt),
             x.data, d, R.n, R.ninv)
  z.data = r
  return z
end

function mul_mod(x::UInt, y::fmpz, R)
  d = ccall((:fmpz_fdiv_ui, libflint), UInt, (Ref{fmpz}, UInt), y, R.n)
  r = ccall((:n_mulmod2_preinv, libflint), UInt, (UInt, UInt, UInt, UInt),
             x, d, R.n, R.ninv)
  return r
end

function mul!(z::Nemo.gfp_fmpz_elem, x::Nemo.gfp_fmpz_elem, y::fmpz)
  R = parent(x)
  ccall((:fmpz_mod, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}),
        z.data, y, R.n)

  ccall((:fmpz_mod_mul, libflint), Nothing,
        (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{Nemo.fmpz_mod_ctx_struct}),
        z.data, x.data, z.data, R.ninv)
  return z
end

function rem!(z::gfp_poly, a::gfp_poly, b::gfp_poly)
  ccall((:nmod_poly_rem, libflint), Nothing,
        	    (Ref{gfp_poly}, Ref{gfp_poly}, Ref{gfp_poly}, Ptr{Nothing}),
          	  z, a, b, pointer_from_objref(base_ring(z))+sizeof(fmpz))
  return z
end

function evaluate_raw(x::gfp_poly, y::gfp_elem)
  z = ccall((:nmod_poly_evaluate_nmod, libflint), UInt,
              (Ref{gfp_poly}, UInt), x, y.data)
  return parent(y)(z)
end

function fmpz_mod_ui(x::fmpz, y::UInt)
  return ccall((:fmpz_fdiv_ui, libflint), UInt, (Ref{fmpz}, UInt), x, y)
end


text:
using Distributed

num_feeds = 50
num_procs =  2
proc_delay = 0.01

for i in 1:num_procs
	println("started proc $(i)")
	sshflags=`-i /Users/robert/src/julia/juliaSnippets/distributed/ms-keys`
	addprocs(["robert@mbp2018.local"]; tunnel=false, sshflags=sshflags, max_parallel=10, dir="/Users/robert/src/julia/juliaSnippets//distributed", topology=:all_to_all)
end

println("defining struct everywhere")
@everywhere struct Ctrl
	n::Int64
	quit::Bool
end

println("loading echo() everywhere")
@everywhere function echo()
	println("echo() started on $(myid())")
	c = Channel{Ctrl}(10)
	@async while true
		d = take!(c)
		if d.quit == true
			println("$(myid()) received quit, exitting, closing channel")
			close(c)
			return
		end
		println("$(myid()) received $(d.n)")
	end
	println("echo() on $(myid()) returning RemoteChannel")
	return c
end

# create an array of remote channels
x = Array{RemoteChannel,1}(undef, num_procs + 1)

# create a remote Int channel in pid 2
for i in 1:num_procs + 1
	x[i]= RemoteChannel(echo, i)
end

println("there are $(procs())")
for i in procs()
	println("starting echo on $(i)")
	remotecall(echo, i)
end

println("feeding...")
while true
	for i=1:num_feeds
		chan = 1 + (i % size(x)[1])
		println("sending value $(i) to $(chan)")
		try
			if (i % 13 == 0)
				put!(x[chan], Ctrl(i, true))
			else
				put!(x[chan], Ctrl(i, false))
			end
		catch e
			if isa(e, InvalidStateException)
				println("channel $(i) is closed, nothing sent")
				deleteat!(x,chan)
			end
		end
		sleep(proc_delay)
	end
	break
end

println("done, cleaning up.")
println("removing $(workers())")
t = rmprocs(workers()...)
wait(t)
println("done.")

max_stars_repo_name: rabarar/juliaSnippets

max_issues_repo_name: rabarar/juliaSnippets

max_forks_repo_name: rabarar/juliaSnippets


text:
module Signals

  using JSON, HTTP, Dates, TimeZones

  include("Utils.jl")
  include("BaseAPI.jl")

  export SignalsAPI,
         get_leaderboard,
         upload_predictions,
         submission_status,
         public_user_profile,
         daily_user_performances,
         daily_submissions_performances,
         ticker_universe,
         download_validation_data,
         stake_get

  const SIGNALS_DOM = "https://numerai-signals-public-data.s3-us-west-2.amazonaws.com"
  const TOURNAMENT = 11

  struct SignalsAPI <: BaseAPI
    public_id::Union{String,Nothing}
    secret_key::Union{String,Nothing}
    tournament::Int
    signals_dom::String
    ticker_universe_url::String
    historical_data_url::String
  end

  function SignalsAPI(public_id::Union{String,Nothing},
                      secret_key::Union{String,Nothing};
                      tournament=TOURNAMENT)
    return SignalsAPI(public_id, secret_key, tournament, SIGNALS_DOM,
                                    "$(SIGNALS_DOM)/latest_universe.csv",
                                    "$(SIGNALS_DOM)/signals_train_val_bbg.csv")
  end

  """
    get_account(api::SignalsAPI)::Dict{String,Any}

  Get all information about your account
  """
  get_account(api::SignalsAPI)::Dict{String,Any} = _get_account(api)

  """
    get_models(api::SignalsAPI)::Dict{String,String}

  Get mapping of account model names to model ids for convenience
  """
  get_models(api::SignalsAPI)::Dict{String,String} = _get_models(api)

  """
    get_current_round(api::SignalsAPI)::Union{Real,Nothing}

  Get number of the current active round
  """
  get_current_round(api::SignalsAPI)::Union{Real,Nothing} = _get_current_round(api)

  """
    get_account_transactions(api::SignalsAPI)::Dict{String,Vector}

  Get all your account deposits and withdrawals
  """
  get_account_transactions(api::SignalsAPI)::Dict{String,Vector} = _get_account_transactions(api)

  """
      get_leaderboard(api::SignalsAPI; limit::Int=50, offset::Int=0)::Vector{Dict}

  Get the current Numerai Signals leaderboard

  # Arguments
  - `limit::Int=50`: number of items to return (optional, defaults to 50)
  - `offset::Int=0`: number of items to skip (optional, defaults to 0)

  # Example
  ```julia-repl
  julia> get_leaderboard(signal_api)
  ```
  """
  function get_leaderboard(api::SignalsAPI; limit::Int=50, offset::Int=0)::Vector{Dict}
    query = """
      query(\$limit: Int! \$offset: Int!) {
        signalsLeaderboard(limit: \$limit offset: \$offset) {
          prevRank
          rank
          sharpe
          today
          username
          mmc
          mmcRank
          nmrStaked
        }
      }
    """
    variables = Dict( "limit" => limit, "offset" => offset )
    data = raw_query(api, query, variables=variables)["data"]["signalsLeaderboard"]
    return data
  end

  """
      upload_predictions(api::SignalsAPI, file_path::String;
                          model_id::Union{String,Nothing}=nothing)::String

  Upload predictions from file
  """
  function upload_predictions(api::SignalsAPI, file_path::String;
                                model_id::Union{String,Nothing}=nothing)::String
    @info "uploading predictions..."
    auth_query = """
      query(\$filename: String! \$modelId: String) {
        submissionUploadSignalsAuth(filename: \$filename modelId: \$modelId) {
              filename
              url
          }
      }
    """
    variables = Dict( "filename" => basename(file_path), "modelId" => model_id )
    submission_resp = raw_query(api, auth_query, variables=variables, authorization=true)
    submission_auth = submission_resp["data"]["submissionUploadSignalsAuth"]
    headers = haskey(ENV, "NUMERAI_COMPUTE_ID") ? ["x_compute_id"=>ENV["NUMERAI_COMPUTE_ID"],] : []

    open(file_path, "r") do io
      HTTP.request("PUT", submission_auth["url"], headers, read(io))
    end

    create_query = """
      mutation(\$filename: String! \$modelId: String) {
        createSignalsSubmission(filename: \$filename modelId: \$modelId) {
          id
          firstEffectiveDate
        }
      }
    """
    variables = Dict( "filename" => submission_auth["filename"], "modelId" => model_id )
    create = raw_query(api, create_query, variables=variables, authorization=true)
    return create["data"]["createSignalsSubmission"]["id"]
  end

  """
      submission_status(api::SignalsAPI; model_id::Union{String,Nothing}=nothing)::Dict

  Submission status of the last submission associated with the account
  """
  function submission_status(api::SignalsAPI; model_id::Union{String,Nothing}=nothing)::Vector{Dict}
    query = """
      query(\$modelId: String) {
          model(modelId: \$modelId) {
            latestSignalsSubmission {
              id
              filename
              firstEffectiveDate
              userId
              submissionIp
              submittedCount
              filteredCount
              invalidTickers
              hasHistoric
              historicMean
              historicStd
              historicSharpe
              historicMaxDrawdown
            }
          }
        }
    """
    variables = Dict( "modelId" =>  model_id )
    data = raw_query(api, query, variables=variables, authorization=true)
    return data["data"]["model"]["latestSignalsSubmission"]
  end

  """
      public_user_profile(api::SignalsAPI, username::String)::Dict

  Fetch the public Numerai Signals profile of a user
  """
  function public_user_profile(api::SignalsAPI, username::String)::Dict
    query = """
      query(\$username: String!) {
        signalsUserProfile(username: \$username) {
          rank
          id
          startDate
          username
          bio
          sharpe
          totalStake
        }
      }
    """
    variables = Dict( "username" => username )
    profile = raw_query(api, query, variables=variables)["data"]["signalsUserProfile"]
    profile["startDate"] = parse_datetime_string(profile["startDate"])
    profile["totalStake"] = parse_float_string(profile["totalStake"])
    return profile
  end

  """
      daily_user_performances(api::SignalsAPI, username::String)::Vector{Dict}

  Fetch daily Numerai Signals performance of a user
  """
  function daily_user_performances(api::SignalsAPI, username::String)::Vector{Dict}
    query = """
      query(\$username: String!) {
        signalsUserProfile(username: \$username) {
          dailyUserPerformances {
            rank
            date
            sharpe
            mmcRep
            reputation
          }
        }
      }
    """
    variables = Dict( "username" => username )
    data = raw_query(api, query, variables=variables)["data"]["signalsUserProfile"]
    performances = data["dailyUserPerformances"]
    for perf in performances
      perf["date"] = parse_datetime_string(perf["date"])
    end
    return performances
  end

  """
      daily_submissions_performances(api::SignalsAPI, username::String)::Vector{Dict}

  Fetch daily Numerai Signals performance of a user's submissions
  """
  function daily_submissions_performances(api::SignalsAPI, username::String)::Vector{Dict}
    query = """
      query(\$username: String!) {
        signalsUserProfile(username: \$username) {
          dailySubmissionPerformances {
            date
            returns
            submissionTime
            correlation
            mmc
            roundNumber
            corrRep
            mmcRep
          }
        }
      }
    """
    variables = Dict( "username" => username )
    data = raw_query(api, query, variables=variables)["data"]["signalsUserProfile"]
    performances = data["dailySubmissionPerformances"]
    for perf in performances
      perf["date"] = parse_datetime_string(perf["date"])
      perf["submissionTime"] = parse_datetime_string(perf["submissionTime"])
    end

    performances = filter(perf -> !isnothing(perf["date"]), performances)

    sort!(performances, by=x -> (x["roundNumber"], x["date"]), rev=true)

    return performances
  end

  """
      ticker_universe(api::SignalsAPI)::Vector{String}

  Fetch universe of accepted tickers
  """
  function ticker_universe(api::SignalsAPI)::Vector{String}
    result = HTTP.request("GET", api.ticker_universe_url)
    tickers = String(result.body)
    tickers = split(tickers, "\n")
    tickers = [strip(ticker) for ticker in tickers
                if (strip(ticker) ≠ "bloomberg_ticker") & (strip(ticker) ≠ "")]
    return tickers
  end

  """
      download_validation_data(api::SignalsAPI;
                                dest_path::String=".",
                                dest_filename::Union{String,Nothing}=nothing,
                                show_progress_bar::Bool=true)::String

  Download CSV file with historical targets and ticker universe
  """
  function download_validation_data(api::SignalsAPI;
                                    dest_path::String=".",
                                    dest_filename::Union{String,Nothing}=nothing,
                                    show_progress_bar::Bool=true)::String
    if isnothing(dest_filename)
      dest_filename = "numerai_signals_historical.csv"
    end
    dataset_path = joinpath(dest_path, dest_filename)

    ensure_directory_exists(dest_path)
    download_file(api.historical_data_url, dataset_path, show_progress_bar=show_progress_bar)

    return dataset_path
  end

  """
      stake_get(api::SignalsAPI, username::String)::Real

  Get current stake for a given users
  """
  function stake_get(api::SignalsAPI, username::String)::Real
    data = public_user_profile(api, username)
    return data["totalStake"]
  end

end # module


text:
module Day11
export solvePart1, solvePart2Slow, solvePart2

serialNumber = 6548

function findCellPower(x::Int, y::Int, serialNumber::Int)
    rackId = x + 10
    power = rackId * y
    power += serialNumber
    power *= rackId
    power = div(power, 100) % 10
    power -= 5
    return power
end

function findAllCellPowers(serialNumber::Int)
    map(xy -> findCellPower(xy[1], xy[2], serialNumber), Base.product(1:300, 1:300))
end

function findBlockPower(x::Int, y::Int, squareSize::Int, powers::Array{Int, 2})
    n = 0
    for i in x:x+squareSize-1
        for j in y:y+squareSize-1
            n += powers[i, j]
        end
    end
    return n
end

function findMaxPowerBlock(squareSize::Int, powers::Array{Int, 2})
    blockPowers = Dict()
    for x in 1:size(powers, 1)-squareSize+1
        for y in 1:size(powers, 2)-squareSize+1
            blockPowers[(x, y)] = findBlockPower(x, y, squareSize, powers)
        end
    end
    maxPower = -Inf
    maxPos = nothing
    for (pos, power) in blockPowers
        if power > maxPower
            maxPower = power
            maxPos = pos
        end
    end
    if maxPos == nothing
        return (nothing, nothing, squareSize), maxPower
    else
        return (maxPos[1], maxPos[2], squareSize), maxPower
    end
end

function solvePart1()
    findMaxPowerBlock(3, findAllCellPowers(serialNumber))[1][1:2]
end

"""I got this idea from someone on the AdventOfCode subreddit.
I've also seen it before, so I should have thought of it.

I wanted to see how fast it was. Much much much much much faster, as it turns out."""
function solvePart2()
    # Summed-area Table indexed starting with 0
    t = zeros(Int64, 301, 301)
    for x in 1:300
        for y in 1:300
            power = findCellPower(x, y, serialNumber)
            t[x+1, y+1] = power + t[x, y+1] + t[x+1, y] - t[x, y]
        end
    end
    best = nothing
    max = -Inf
    for s in 1:300
        for x in 1:300-s
            for y in 1:300-s
                total = t[x+s, y+s] - t[x, y+s] - t[x+s, y] + t[x, y]
                if total > max
                    best = (s, x, y)
                    max = total
                end
            end
        end
    end
    return best
end

"""This was my first attempt. It's very slow (~70 seconds)"""
function solvePart2Slow()
    powers = findAllCellPowers(6548)
    maxCoord = nothing
    maxPower = -Inf
    for size in 1:300
        coord, power = findMaxPowerBlock(size, powers)
        if power > maxPower
            maxCoord = coord
            maxPower = power
        end
    end
    x, y, s = maxCoord
    return s, x, y
end
end


text:
module Fairness

# ===================================================================
## IMPORTS

import Base
using Random
using DataFrames, CSV # For loading datasets
using MLJBase
using CategoricalArrays
using MLJModels, MLJModelInterface
using StatsBase # For reweighing algorithm
using JuMP, GLPK # For Equalized Odds Postprocessing algorithms
using Ipopt # For LinProgWrapper algorithm
using ForwardDiff # For MetaFair algorithm
using ZipFile # For downloading datasets
import Distributions
using HypothesisTests # For Benchmarking
using LinearAlgebra # For synthetic datasets
using Suppressor, Logging # To enable/disable unnecessary logs while evaluation
using ProgressMeter # To track progress while running benchmarking
# ===================================================================
## METHOD EXPORTS

export fair_tensor, fact
export _ftIdx
export istype

# Export the Boolean Metrics
export DemographicParity

export accuracy

#Export the metric instances from MLJBase to permit calculation of metrics without using MLJBase
export TruePositive, TrueNegative, FalsePositive, FalseNegative,
       TruePositiveRate, TrueNegativeRate, FalsePositiveRate,
       FalseNegativeRate, FalseDiscoveryRate, Precision, NPV,
       # standard synonyms
       TPR, TNR, FPR, FNR, FDR, PPV,
       # instances and their synonyms
       truepositive, truenegative, falsepositive, falsenegative,
       true_positive, true_negative, false_positive, false_negative,
       truepositive_rate, truenegative_rate, falsepositive_rate,
       true_positive_rate, true_negative_rate, false_positive_rate,
       falsenegative_rate, negativepredictive_value,
       false_negative_rate, negative_predictive_value,
       positivepredictive_value, positive_predictive_value,
       tpr, tnr, fpr, fnr,
       falsediscovery_rate, false_discovery_rate, fdr, npv, ppv,
       recall, sensitivity, hit_rate, miss_rate,
       specificity, selectivity, f1score, fallout

export PredictedPositiveRate, PP, predicted_positive_rate,
        ppr, FalseOmissionRate, FOR, false_omission_rate, foar,TruePositiveRateDifference, TPRD, true_positive_rate_difference,tprd,FalsePositiveRateDifference, FPRD, false_positive_rate_difference, fprd, FalseNegativeRateDifference, FNRD, false_negative_rate_difference, fnrd, FalseDiscoveryRateDifference, FDRD, false_discovery_rate_difference, fdr, FalsePositiveRateRatio, FPRR, false_positive_rate_ratio, fprr, FalseNegativeRateRatio, FNRR, false_negative_rate_ratio, fnrr, FalseOmissionRateRatio, FORR, false_omission_rate_ratio, forr, FalseDiscoveryRateRatio, FDRR, false_discovery_rate_ratio, fdrr, AverageOddsDifference, AOD, average_odds_difference, aod

# Export Fairness Metric Wrappers
export MetricWrapper, MetricWrappers

# Export the fairness metrics
export disparity, parity, Disparity, Disparities

export ReweighingWrapper, ReweighingSamplingWrapper
export EqOddsWrapper, LinProgWrapper, CalEqOddsWrapper
export PenaltyWrapper

# Export macros for datasets from datasets/
export @load_toydata, @load_toyfairtensor
export @load_compas, @load_adult, @load_german,
       @load_bank_marketing, @load_communities_crime, @load_student_performance

export genZafarData, genSubgroupData, genZafarData2, genBiasedSampleData, law_school
# -------------------------------------------------------------------
# re-export From CategoricalArrays
export categorical, levels, levels!

# re-export from MLJBase
export pretty

export Task, FairnessProblem
export fairevaluate, get_df, get_pareto_df
# export fairness problems of available datasets
export communities_crime, student, hmda, zafar, zafar2, subgroup, biased,
        adult, german, portuguese, framingham, loan, medical, fairness_problems
# ===================================================================
## CONSTANTS

const CategoricalElement = Union{CategoricalValue,CategoricalValue{String, R} where R}
const Vec = AbstractVector
const Measure =  MLJBase.Measure
const MMI =  MLJModelInterface

# the directory containing this file: (.../src/)
const MODULE_DIR = dirname(@__FILE__)

# This will be used to find whether the outcomes we will get will be probabilistic or deterministic
# This will can be defined overrided for any algorithm. This will specifically help in preprocessing algorithms
istype(x, y) = isa(x, y)
# ===================================================================
## Includes

include("utilities.jl")
include("fair_tensor.jl")
include("measures/measures.jl")
include("datasets/datasets.jl")
include("datasets/synthetic.jl")
include("algorithms/algorithms.jl")
include("benchmark/benchmark.jl")

end # module


text:
# Estimation of finite mixture model using EM

immutable FiniteMixtureEM
	maxiter::Int
	tol::Float64
	display::Symbol
	alpha::Float64
end

function fmm_em(;
	maxiter::Integer=100,
	tol::Real=1.0e-6,
	display::Symbol=:none,
	alpha::Float64=1.0)

	FiniteMixtureEM(maxiter, tol, display, alpha)
end

immutable FiniteMixtureEMResults{C}
	mixture::Mixture{C}
	Q::Matrix{Float64}
	L::Matrix{Float64}
	niters::Int
	converged::Bool
	objective::Float64
end

function fit_fmm!{C<:Distribution}(estimator::Estimator{C},
	data, Q::Matrix{Float64}, alg::FiniteMixtureEM)

	# basic numbers

	n = nsamples(estimator, data)
	@check_argdims size(Q, 1) == n

	K = size(Q, 2)
	maxiter = alg.maxiter
	tol = alg.tol
	verbose = verbosity_level(alg.display)

	# prepare storage

	components = Array(C, K)
	L = Array(Float64, n, K)
	E = Array(Float64, n, K)
	pi = Array(Float64, K)
	lpi = Array(Float64, K)
	qent = Array(Float64, n)
	objv = NaN

	# main loop

	t = 0
	converged = false

	if verbose >= VERBOSE_PROC
		@printf("%6s    %12s     %12s\n", "iter", "objective", "obj.change")
		println("------------------------------------------")
	end

	while !converged && t < maxiter
		t = t + 1

		# M-step
		sum!(pi, Q, 1)
		multiply!(pi, inv(n))
		for k in 1 : K
			lpi[k] = log(pi[k])
		end

		for k in 1 : K
			comp::C = estimate(estimator, data, refcolumn(Q, k))
			components[k] = comp
			logpdf!(refcolumn(L, k), comp, data)
		end

		# E-step
		vbroadcast!(Add(), E, L, lpi, 2)
		softmax!(Q, E, 2)

		# evaluate objective
		t_lpi = sum(Q * lpi)
		t_llik = dot(L, Q)
		t_qent = sum(entropy!(qent, Q, 2))
		t_lpri = 0.
		for k in 1 : K
			t_lpri += prior_score(estimator, components[k])
		end

		objv_pre = objv
		objv = t_lpi + t_llik + t_qent + t_lpri

		# decide convergence
		converged = t > 1 && abs(objv - objv_pre) < tol

		if verbose >= VERBOSE_ITER
			@printf("%6d    %12.4e     %12.4e\n", t, objv, objv - objv_pre)
		end
	end

	if verbose >= VERBOSE_PROC
		if converged
			println("EM procedure converged with $t iterations.")
		else
			println("EM procedure terminated after $t iterations without convergence.")
		end
	end

	mixture = Mixture(components, pi)
	return FiniteMixtureEMResults{C}(mixture, Q, L, t, converged, objv)
end

function fit_fmm!{C<:Distribution}(dty::Type{C}, data, Q::Matrix{Float64}, alg::FiniteMixtureEM)
	fit_fmm!(MLEstimator(C), data, Q, alg)
end

function fit_fmm{C<:Distribution}(dty::Type{C}, data, K::Int, alg::FiniteMixtureEM)
	n = nsamples(C, data)
	fit_fmm!(MLEstimator(C), data, qmatrix(n, K), alg)
end