text:
module GLRenderer

import PoseComposition: Pose
import Rotations
import Images as I
import Images: Color, RGBA
import Parameters: @with_kw
import GLFW
import PyCall
import FileIO
import MeshIO
import GeometryBasics as GB
import ColorSchemes
using ModernGL

@with_kw struct CameraIntrinsics
    width::Int = 640
    height::Int = 480
    fx::Float64 = width
    fy::Float64 = width
    cx::Float64 = width/2
    cy::Float64 = height/2
    near::Float64 = 0.001
    far::Float64 = 100.0
end

@with_kw mutable struct Mesh
    vertices::Any = nothing
    indices::Any = nothing
    normals::Any = nothing
    tex_coords::Any = nothing
    tex_path::Any = nothing
end

copy_mesh(m::Mesh) = Mesh(
                        vertices=m.vertices,
                        indices=m.indices,
                        normals=m.normals,
                        tex_coords=m.tex_coords,
                        tex_path=m.tex_path,
                    )

include("shaders.jl")
include("utils.jl")
include("mesh.jl")
include("point_cloud.jl")

abstract type RenderMode end
struct DepthMode <: RenderMode end
struct RGBBasicMode <: RenderMode end
struct RGBMode <: RenderMode end
struct TextureMode <: RenderMode end
struct TextureMixedMode <: RenderMode end

mutable struct Renderer{T <: RenderMode}
    window::Any
    camera_intrinsics::CameraIntrinsics
    shader_program::Any
    mesh_pointers::Any
    mesh_sizes::Any
    textures::Any
    perspective_matrix::Matrix
end

function setup_renderer(camera_intrinsics::CameraIntrinsics, mode::RenderMode; name="GLRenderer", gl_version=nothing)::Renderer
    if isnothing(gl_version)
        window = GLFW.CreateWindow(1,1, "dummy")
        GLFW.MakeContextCurrent(window)
        version_string = unsafe_string(glGetString(GL_VERSION))
        gl_version = (i -> parse(Int,i)).(split(split(version_string," ")[1],"."))
        GLFW.DestroyWindow(window)
    end

    perspective_matrix = get_perspective_matrix(
        camera_intrinsics.width, camera_intrinsics.height,
        camera_intrinsics.fx,
        camera_intrinsics.fy,
        camera_intrinsics.cx,
        camera_intrinsics.cy,
        camera_intrinsics.near,
        camera_intrinsics.far,
    )

    window_hint = [
        # (GLFW.SAMPLES,      0),
        # (GLFW.DEPTH_BITS,   24),
        # (GLFW.ALPHA_BITS,   8),
        # (GLFW.RED_BITS,     8),
        # (GLFW.GREEN_BITS,   8),
        # (GLFW.BLUE_BITS,    8),
        # (GLFW.STENCIL_BITS, 0),
        # (GLFW.AUX_BUFFERS,  0),
        (GLFW.CONTEXT_VERSION_MAJOR, gl_version[1]),
        (GLFW.CONTEXT_VERSION_MINOR, gl_version[2]),
        (GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE),
        (GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE),
        (GLFW.VISIBLE, GL_FALSE)
    ]
    for (key, value) in window_hint
        GLFW.WindowHint(key, value)
    end

    window = GLFW.CreateWindow(camera_intrinsics.width, camera_intrinsics.height, "GLRenderer")
    GLFW.MakeContextCurrent(window)
    glEnable(GL_DEPTH_TEST)
    glClear(GL_DEPTH_BUFFER_BIT)

    gl_version_for_shaders = "$(gl_version[1])$(gl_version[2])0"
    @show gl_version_for_shaders
    if typeof(mode) == DepthMode
        vertex_shader = createShader(vertex_source_depth(gl_version_for_shaders), GL_VERTEX_SHADER)
        fragment_shader = createShader(fragment_source_depth(gl_version_for_shaders), GL_FRAGMENT_SHADER)
    elseif typeof(mode) == RGBBasicMode
        vertex_shader = createShader(vertexShader_rgb_basic(gl_version_for_shaders), GL_VERTEX_SHADER)
        fragment_shader = createShader(fragmentShader_rgb_basic(gl_version_for_shaders), GL_FRAGMENT_SHADER)
    elseif typeof(mode) == RGBMode
        vertex_shader = createShader(vertexShader_rgb(gl_version_for_shaders), GL_VERTEX_SHADER)
        fragment_shader = createShader(fragmentShader_rgb(gl_version_for_shaders), GL_FRAGMENT_SHADER)
    elseif typeof(mode) == TextureMode
        vertex_shader = createShader(vertexShader_texture(gl_version_for_shaders), GL_VERTEX_SHADER)
        fragment_shader = createShader(fragmentShader_texture(gl_version_for_shaders), GL_FRAGMENT_SHADER)
    elseif typeof(mode) == TextureMixedMode
        vertex_shader = createShader(vertexShader_texture_mixed(gl_version_for_shaders), GL_VERTEX_SHADER)
        fragment_shader = createShader(fragmentShader_texture_mixed(gl_version_for_shaders), GL_FRAGMENT_SHADER)
    else
        error("unsupported mode $(mode)")
    end

    shader_program = glCreateProgram()
    glAttachShader(shader_program, vertex_shader)
    glAttachShader(shader_program, fragment_shader)
    glLinkProgram(shader_program)

    fbo = Ref(GLuint(0))
    glGenFramebuffers(1, fbo)

    depth_tex = Ref(GLuint(0))
    glGenTextures(1, depth_tex)
    color_tex = Ref(GLuint(0))
    glGenTextures(1, color_tex)

    glBindFramebuffer(GL_FRAMEBUFFER, fbo[])

    glBindTexture(GL_TEXTURE_2D, depth_tex[])
    glTexImage2D(
      GL_TEXTURE_2D, 0, GL_DEPTH24_STENCIL8, camera_intrinsics.width, camera_intrinsics.height, 0,
      GL_DEPTH_STENCIL, GL_UNSIGNED_INT_24_8, C_NULL
    );
    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_DEPTH_STENCIL_ATTACHMENT, GL_TEXTURE_2D, depth_tex[], 0);

    glBindTexture(GL_TEXTURE_2D, color_tex[])
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, camera_intrinsics.width, camera_intrinsics.height, 0,
                 GL_RGBA, GL_UNSIGNED_BYTE, C_NULL)
    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, color_tex[], 0)

    glViewport(0, 0, camera_intrinsics.width, camera_intrinsics.height)
    glDrawBuffers(2, [GL_DEPTH_STENCIL_ATTACHMENT, GL_COLOR_ATTACHMENT0])
    println(glCheckFramebufferStatus(GL_FRAMEBUFFER) == GL_FRAMEBUFFER_COMPLETE)

    Renderer{typeof(mode)}(window, camera_intrinsics, shader_program, [], [], [], perspective_matrix)
end

function activate_renderer(renderer::Renderer)
    GLFW.MakeContextCurrent(renderer.window)
end

function set_intrinsics!(renderer::Renderer, camera_intrinsics::CameraIntrinsics)
    renderer.camera_intrinsics = camera_intrinsics
    renderer.perspective_matrix = get_perspective_matrix(
        camera_intrinsics.width, camera_intrinsics.height,
        camera_intrinsics.fx,
        camera_intrinsics.fy,
        camera_intrinsics.cx,
        camera_intrinsics.cy,
        camera_intrinsics.near,
        camera_intrinsics.far,
    )

    fbo = Ref(GLuint(0))
    glGenFramebuffers(1, fbo)

    depth_tex = Ref(GLuint(0))
    glGenTextures(1, depth_tex)
    color_tex = Ref(GLuint(0))
    glGenTextures(1, color_tex)

    glBindFramebuffer(GL_FRAMEBUFFER, fbo[])

    glBindTexture(GL_TEXTURE_2D, depth_tex[])
    glTexImage2D(
      GL_TEXTURE_2D, 0, GL_DEPTH24_STENCIL8, camera_intrinsics.width, camera_intrinsics.height, 0,
      GL_DEPTH_STENCIL, GL_UNSIGNED_INT_24_8, C_NULL
    );
    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_DEPTH_STENCIL_ATTACHMENT, GL_TEXTURE_2D, depth_tex[], 0);

    glBindTexture(GL_TEXTURE_2D, color_tex[])
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, camera_intrinsics.width, camera_intrinsics.height, 0,
                 GL_RGBA, GL_UNSIGNED_BYTE, C_NULL)
    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, color_tex[], 0)

    glViewport(0, 0, camera_intrinsics.width, camera_intrinsics.height)
    glDrawBuffers(2, [GL_DEPTH_STENCIL_ATTACHMENT, GL_COLOR_ATTACHMENT0])
    println(glCheckFramebufferStatus(GL_FRAMEBUFFER) == GL_FRAMEBUFFER_COMPLETE)
end

function load_object!(renderer::Renderer{T}, mesh) where T <: RenderMode
    vao = Ref(GLuint(0))
    glGenVertexArrays(1, vao)
    glBindVertexArray(vao[])

    if T == DepthMode || T == RGBBasicMode
        @assert size(mesh.vertices)[1] == 3
        vertex_data = mesh.vertices
    elseif T == RGBMode
        @assert size(mesh.vertices)[1] == 3
        @assert size(mesh.normals)[1] == 3
        @assert size(mesh.vertices)[2] == size(mesh.normals)[2]
        vertex_data = vcat(mesh.vertices, mesh.normals)
    else
        @assert size(mesh.vertices)[1] == 3
        @assert size(mesh.normals)[1] == 3
        if T == TextureMixedMode && isnothing(mesh.tex_coords)
            tex_coords = zeros(2, size(mesh.vertices)[2])
        else
            @assert size(mesh.tex_coords)[1] == 2
            tex_coords = mesh.tex_coords
        end
        @assert size(mesh.vertices)[2] == size(mesh.normals)[2] == size(tex_coords)[2]
        vertex_data = vcat(mesh.vertices, mesh.normals, tex_coords)
    end

    vertex_data = Matrix{Float32}(vertex_data)
    indices = Matrix{UInt32}(mesh.indices)

    # copy vertex data into an OpenGL buffer
    vbo = Ref(GLuint(0))
    glGenBuffers(1, vbo)
    glBindBuffer(GL_ARRAY_BUFFER, vbo[])
    glBufferData(GL_ARRAY_BUFFER, sizeof(vertex_data), Ref(vertex_data, 1), GL_STATIC_DRAW)

    # element buffer object for indices
    ebo = Ref(GLuint(0))
    glGenBuffers(1, ebo)
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ebo[])
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(indices), Ref(indices, 1), GL_STATIC_DRAW)

    if T == DepthMode || T == RGBBasicMode
        glEnableVertexAttribArray(glGetAttribLocation(renderer.shader_program, "position"))

        # set vertex attribute pointers
        glVertexAttribPointer(glGetAttribLocation(renderer.shader_program, "position"),
            3, GL_FLOAT, GL_FALSE, 3 * sizeof(Float32), C_NULL)
        glEnableVertexAttribArray(glGetAttribLocation(renderer.shader_program, "position"))
    elseif T == RGBMode
        # set vertex attribute pointers
        glEnableVertexAttribArray(glGetAttribLocation(renderer.shader_program, "position"))
        glEnableVertexAttribArray(glGetAttribLocation(renderer.shader_program, "normal"))

        glVertexAttribPointer(
            glGetAttribLocation(renderer.shader_program, "position"),
            3, GL_FLOAT, GL_FALSE, 3 * sizeof(Float32) * 2, C_NULL)
        glVertexAttribPointer(
            glGetAttribLocation(renderer.shader_program, "normal"),
            3, GL_FLOAT, GL_FALSE, 3 * sizeof(Float32) * 2, C_NULL + (3 * sizeof(Float32)))
    else
        # set vertex attribute pointers
        glEnableVertexAttribArray(glGetAttribLocation(renderer.shader_program, "position"))
        glEnableVertexAttribArray(glGetAttribLocation(renderer.shader_program, "normal"))
        glEnableVertexAttribArray(glGetAttribLocation(renderer.shader_program, "vertTexCoord"))

        glVertexAttribPointer(
            glGetAttribLocation(renderer.shader_program, "position"),
            3, GL_FLOAT, GL_FALSE, 32, C_NULL)
        glVertexAttribPointer(
            glGetAttribLocation(renderer.shader_program, "normal"),
            3, GL_FLOAT, GL_FALSE, 32, C_NULL + 12)
        glVertexAttribPointer(
            glGetAttribLocation(renderer.shader_program, "vertTexCoord"),
            2, GL_FLOAT, GL_FALSE, 32, C_NULL + 24)

        PyCall.py"""
        from PIL import Image
        import numpy as np
        def load_texture_bytes(path):
            img = Image.open(path).transpose(Image.FLIP_TOP_BOTTOM)
            img_data = np.fromstring(img.tobytes(), np.uint8)
            return img_data, img.width, img.height
        """
        if isnothing(mesh.tex_path)
            push!(renderer.textures, nothing)
        else
            img_data, width, height = PyCall.py"load_texture_bytes"(mesh.tex_path);

            texture = Ref(GLuint(0))
            glGenTextures(1, texture)
            glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
            glBindTexture(GL_TEXTURE_2D, texture[])
            glTexParameterf(
                GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
            glTexParameterf(
                GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR)
            glTexParameterf(
                GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
            glTexParameterf(
                GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
            glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0,
                            GL_RGB, GL_UNSIGNED_BYTE, img_data)
            glGenerateMipmap(GL_TEXTURE_2D)

            push!(renderer.textures, texture)
        end

    end

    # unbind it
    glBindBuffer(GL_ARRAY_BUFFER, 0)
    glBindVertexArray(0)

    push!(renderer.mesh_pointers, vao[])
    push!(renderer.mesh_sizes, size(mesh.indices)[2] * 3)
    return true
end

function gl_render(
        renderer::Union{Renderer{T}}, mesh_ids::Vector{Int},
        poses::Vector{Pose}, camera_pose::Pose; colors=nothing
) where T <: RenderMode

    if isnothing(colors)
        colors = fill(I.colorant"black", length(poses))
    end
    colors = map(x -> convert(I.RGBA,x), colors)

    @assert length(mesh_ids) == length(poses)

    glClearColor(1.0, 1.0, 1.0, 1)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glEnable(GL_DEPTH_TEST)

    camera_pose_mat = pose_to_matrix(inv(camera_pose * Pose(zeros(3),Rotations.RotX(pi))))

    for (id, p, c) in zip(mesh_ids, poses, colors)
        if id < 1
            continue
        end
        vao = renderer.mesh_pointers[id]

        glUseProgram(renderer.shader_program)

        glUniformMatrix4fv(glGetUniformLocation(
                renderer.shader_program,"P"), 1, GL_FALSE, Ref(renderer.perspective_matrix, 1))
        glUniformMatrix4fv(glGetUniformLocation(
                renderer.shader_program,"V"), 1, GL_FALSE, Ref(camera_pose_mat, 1))
        glUniformMatrix4fv(glGetUniformLocation(
                renderer.shader_program,"pose_mat"), 1, GL_FALSE, Ref(pose_to_matrix(p), 1))

        if T == RGBBasicMode || T == RGBMode
            glUniform4fv(glGetUniformLocation(
                    renderer.shader_program,"color"), 1, Float32[c.r, c.g, c.b, c.alpha])
        elseif T == TextureMode || T == TextureMixedMode
            glUniformMatrix4fv(glGetUniformLocation(
                    renderer.shader_program,"pose_rot"), 1, GL_FALSE, Ref(pose_to_matrix(Pose(zeros(3),p.orientation)), 1))

            if isnothing(renderer.textures[id])
                glUniform1f(glGetUniformLocation(
                    renderer.shader_program,"textured"), 0.0)
                glUniform4fv(glGetUniformLocation(
                        renderer.shader_program,"color"), 1, Float32[c.r, c.g, c.b, c.alpha])
            else
                glUniform1f(glGetUniformLocation(
                    renderer.shader_program,"textured"), 1.0)
                glActiveTexture(GL_TEXTURE0)
                glBindTexture(GL_TEXTURE_2D, renderer.textures[id][])
                glUniform1i(glGetUniformLocation(
                        renderer.shader_program,"tex"), 0)
            end
        end

        glBindVertexArray(vao)
        glDrawElements(GL_TRIANGLES, renderer.mesh_sizes[id], GL_UNSIGNED_INT, C_NULL)
        glBindVertexArray(0)
    end

    cam = renderer.camera_intrinsics
    data = Matrix{Float32}(undef, cam.width, cam.height)
    glReadPixels(0, 0, cam.width, cam.height, GL_DEPTH_COMPONENT, GL_FLOAT, Ref(data, 1))
    near,far = cam.near, cam.far
    depth_image = far .* near ./ (far .- (far .- near) .* data)
    depth_image = permutedims(depth_image[:,end:-1:1])

    if T == DepthMode
        return depth_image
    end

    glReadBuffer(GL_COLOR_ATTACHMENT0)
    data1 = zeros(Float32, 4,cam.width, cam.height)
    glReadPixels(0, 0, cam.width, cam.height, GL_BGRA, GL_FLOAT, Ref(data1, 1))
    rgb = permutedims(data1,(3,2,1))[end:-1:1,:,:]
    rgb = cat(rgb[:,:,3],rgb[:,:,2],rgb[:,:,1],rgb[:,:,4],dims=3)
    return rgb, depth_image
end

# Viewing images
function view_depth_image(depth_image; scheme=nothing)
    max_val = maximum(depth_image)
    min_val = minimum(depth_image)
    # vals = depth_image[(depth_image .> (min_val + 1e-3)) .&  (depth_image .< (max_val - 1e-3))]
    # vals = depth_image[(depth_image .> (min_val + 1e-3)) .&  (depth_image .< (max_val - 1e-3))]
    # max_val = maximum(vals)
    # min_val = minimum(vals)
    d = clamp.(depth_image, min_val, max_val)
    img = get(ColorSchemes.hot, d, (min_val, max_val))
    I.convert.(I.RGBA, img)
end

function view_rgb_image(rgb_image; in_255=false)
    if in_255
        rgb_image = Float64.(rgb_image)./255.0
    end
    if size(rgb_image)[3] == 3
        img = I.colorview(I.RGB, permutedims(rgb_image,(3,1,2)))
        img = I.convert.(I.RGBA, img)
    else
        img = I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))
    end
    img
end

end


text:
struct NetworkOutput
    value::Matrix{Float32}
    reward::Matrix{Float32}
    policy_logits::Matrix{Float32}
    hidden_state::AbstractArray{Float32, 4}
end

function Base.:(==)(a::NetworkOutput, b::NetworkOutput)
    a.value == b.value &&
    a.reward == b.reward &&
    a.policy_logits == b.policy_logits &&
    a.hidden_state == b.hidden_state
end

function Base.show(io::IO, x::NetworkOutput)
    println(stdout, "value: $(x.value) reward: $(x.reward), policy: [$(join(x.policy_logits,","))]")
end

function NetworkOutput(tuple)
    value, reward, policy_logits, hidden_state = tuple
    NetworkOutput(value, reward, policy_logits, hidden_state)
end

tuple(x::NetworkOutput) = x.value, x.reward, x.policy_logits, x.hidden_state
value(x::NetworkOutput) = x.value
reward(x::NetworkOutput) = x.reward
policy_logits(x::NetworkOutput) = x.policy_logits
# @todo change this to get_policy
policy(x::NetworkOutput) = softmax(x.policy_logits)
hidden_state(x::NetworkOutput) = x.hidden_state

function next_action(x::NetworkOutput; greedy=false)
    if greedy
        return argmax(reshape(policy(x), :))
    end
    rand(Categorical(reshape(policy(x), :)))
end

struct SingleNetworkOutput
    output::NetworkOutput
    SingleNetworkOutput(tuple) = new(tuple |> NetworkOutput)
end

pred_values(x::NetworkOutput) = x.value
pred_policies(x::NetworkOutput) = softmax(x.policy_logits)
value(x::SingleNetworkOutput) = view(x.output.value, :, 1)[1]
reward(x::SingleNetworkOutput) = view(x.output.reward, :, 1)[1]
hidden_state(x::SingleNetworkOutput) = view(x.output.hidden_state, :, :, :, 1)
policy_logits(x::SingleNetworkOutput) = view(x.output.policy_logits, :, 1)
policy(x::SingleNetworkOutput) = view(policy(x.output), :, 1)
getpolicy(x::SingleNetworkOutput) = policy(x)
greedy_action(x::SingleNetworkOutput) = policy_action(x, greedy=true)

function policy_action(x::SingleNetworkOutput; greedy=false)
    greedy ? argmax(policy(x)) : rand(Categorical(policy(x)))
end

ext: jl


text:
type LSC
    n_clusters::Int64
    n_landmarks::Int64
    method::Symbol
    non_zero_landmarks::Int64
    bandwidth::Float64
    cluster_result::Vector
end

function LSC(;
             n_clusters::Int64 = 2,
             n_landmarks::Int64 = 150,
             method::Symbol = :Kmeans,
             non_zero_landmarks::Int64 = 4,
             bandwidth::Float64 = 0.4,
             cluster_result::Vector = zeros(10))

    return LSC(n_clusters, n_landmarks, method, non_zero_landmarks,
               bandwidth, cluster_result)

end

function gaussianKernel(distance, bandwidth)
    exp(-distance / (2*bandwidth^2));
end

function get_landmarks(X, p;method=:Kmeans)
    if(method == :Random)
        numberOfPoints = size(X,2);
        landmarks = X[:,randperm(numberOfPoints)[1:p]];
        return landmarks;
    end

    if(method == :Kmeans)
        kmeansResult = kmeans(X,p)
        return kmeansResult.centers;
    end
    throw(ArgumentError("method can only be :Kmeans or :Random"));
end

function compose_sparse_Z_hat_matrix(X, landmarks, bandwidth, r)

    distances = pairwise(Distances.Euclidean(), landmarks, X);
    similarities = map(x -> gaussianKernel(x, bandwidth), distances);
    ZHat = zeros(size(similarities));

    for i in 1:(size(similarities,2))
        topLandMarksIndices = selectperm(similarities[:,i], 1:r, rev=true);
        topLandMarksCoefficients = similarities[topLandMarksIndices, i];
        topLandMarksCoefficients = topLandMarksCoefficients / sum(topLandMarksCoefficients);
        ZHat[topLandMarksIndices,i] = topLandMarksCoefficients;
    end
    return diagm(sum(ZHat,2)[:])^(-1/2) * ZHat;

end

function train!(model::LSC, X::Matrix)
    if size(X,1) > size(X,2)
        X = transpose(X)
    end
    landmarks = get_landmarks(X, model.n_landmarks, method = model.method)
    Z_hat = compose_sparse_Z_hat_matrix(X, landmarks, model.bandwidth,
                                        model.non_zero_landmarks)
    svd_result = svd(transpose(Z_hat))
    temp = transpose(svd_result[1][:,1:model.n_clusters])
    model.cluster_result = kmeans(temp, model.n_clusters).assignments
end

function plot_in_2d(model::LSC, X::Matrix)
    if size(X,1) > size(X,2)
        X = transpose(X)
    end
    Gadfly.plot(x = X[1,:], y = X[2, :], color = model.cluster_result)
end

function test_LSC()
    path_ = "../data"
    data_name = ["smiley", "spirals", "shapes", "cassini1"]
    datasets = copy(data_name)
    for i in 1:4
        datasets[i] = datasets[i] * ".csv"
    end
    for i in 1:4
        datasets[i] = joinpath(path_, datasets[i])
    end
    clusters = [4, 2, 4, 3]
    n_landmarks = [50, 150, 50, 50]
    bandwidth = [0.4, 0.04, 0.04, 0.4]
    count = 0
    df = DataFrame()
    for i in 1:4
        data = readdlm(datasets[i])
        data = transpose(convert(Array{Float64,2}, data[:,1:2]))
        model = LSC(n_clusters = clusters[i], n_landmarks = n_landmarks[i], bandwidth = bandwidth[i])
        train!(model,data)
        dataframe = DataFrame(x=data[1,:],y=data[2,:],group=model.cluster_result,datasets=data_name[i])
        if count == 0
            df = dataframe
            count += 1
        else
            df = vcat(df, dataframe)
        end
        println("Progress: $(i/4*100)%....")
        #println("$(df)")
    end
    println("$(df)")
    println("computing finied, drawing the plot......")
    set_default_plot_size(25cm, 14cm)
    Gadfly.plot(df, x="x", y = "y", color = "group",xgroup = "datasets", Geom.subplot_grid(Geom.point), Guide.title("Large Scale Spectral Clustering"))
end


text:
using GLMakie
using FFTW
using SignalBase
using Statistics

export Scope

@enum Mode TIME FREQ TIMEFREQ

mutable struct Scope
  fap::Makie.FigureAxisPlot
  mode::Mode
  fs::Float32
  n::Int
  nfft::Int
  ch::Observable{Int}
  dc::Observable{Bool}
  y::Observable{Vector{Float32}}
  z::Observable{Matrix{Float32}}
  zmax::Observable{Float32}
  zrange::Observable{Float32}
  running::Observable{Bool}
  dirty::Bool
  ybuf::Vector{Float32}
  zbuf::Matrix{Float32}
  closed::Bool
  waitfordirty::Threads.Condition
  src::Any
end

function Scope(; fs=48000f0, bufsize=4096, nfft=1024, history=256)
  n = bufsize
  m = nfft ÷ 2 + 1
  t = range(0f0; length=n, step=1000f0/Float32(fs))
  y = Observable(zeros(Float32, n))
  ybuf = similar(y[])
  z = Observable(zeros(Float32, history, m))
  zbuf = similar(z[])
  fap = lines(t, y; axis=(xlabel="Time (ms)", xautolimitmargin=(0f0, 0f0)))
  display(fap)
  scope = Scope(fap, TIME, Float32(fs), n, nfft, Observable(1), Observable(true), y, z, Observable(20f0), Observable(50f0),
    Observable(true), false, ybuf, zbuf, false, Threads.Condition(), nothing)
  annotate(scope)
  on(events(fap.figure).keyboardbutton) do event
    Consume(keypress(scope, event.action, event.key))
  end
  task = @async monitor(scope)
  scope
end

function Scope(ai; nfft=1024, history=256)
  scope = Scope(; fs=inputframerate(ai), bufsize=inputblocksize(ai), nfft, history)
  scope.src = ai
  connect(ai)
  startinputstream(ai, (t, x) -> begin
    if scope.ch[] ≤ size(x,2)
      if scope.dc[]
        push!(scope, @view x[:,scope.ch[]])
      else
        x1 = x[:,scope.ch[]]
        x1 .-= mean(x1)
        push!(scope, x1)
      end
    end
  end)
  scope
end

Base.show(io::IO, scope::Scope) = print(io, "Scope(...)")

function Base.close(scope::Scope)
  scope.closed = true
  lock(() -> notify(scope.waitfordirty), scope.waitfordirty)
  nothing
end

Base.run(scope::Scope, v=true) = (scope.running[] = v)

function mode!(scope::Scope, mode::Mode)
  scope.mode == mode && return scope
  scope.mode = mode
  if mode == TIME
    t = range(0f0; length=scope.n, step=1000f0/scope.fs)
    scope.y = Observable(zeros(Float32, scope.n))
    scope.ybuf = similar(scope.y[])
    scope.fap = lines(t, scope.y; axis=(xlabel="Time (ms)", xautolimitmargin=(0f0, 0f0)))
  elseif mode == FREQ
    f = rfftfreq(scope.nfft, scope.fs) ./ 1000f0
    scope.y = Observable(zeros(Float32, length(f)))
    scope.ybuf = similar(scope.y[])
    scope.fap = lines(f, scope.y; axis=(xlabel="Frequency (kHz)", xautolimitmargin=(0f0, 0f0)))
    ylims!(scope.zmax[] - scope.zrange[], scope.zmax[])
  elseif mode == TIMEFREQ
    scope.z[] .= -Inf32
    scope.zbuf .= -Inf32
    n = size(scope.z[], 1)
    t = range(0f0; length=n, step=scope.n/scope.fs)
    f = rfftfreq(scope.nfft, scope.fs) ./ 1000f0
    scope.fap = heatmap(t, f, scope.z;
      colorrange=@lift(($(scope.zmax) - $(scope.zrange), $(scope.zmax))),
      axis=(xlabel="Time (s)", ylabel="Frequency (kHz)"))
    Colorbar(scope.fap.figure[1,2], scope.fap.plot)
  else
    throw(ArgumentError("Bad mode"))
  end
  annotate(scope)
  on(events(scope.fap.figure).keyboardbutton) do event
    Consume(keypress(scope, event.action, event.key))
  end
  display(scope.fap)
  scope
end

function psd(x, nfft)
  #w = 0.54 .+ 0.46 .* cos.(2π .* range(-0.5, 0.5; length=nfft))    # hamming
  w = cos.(π .* range(-0.5, 0.5; length=nfft)) .^ 2                 # hanning
  X = mean(abs.(rfft(x1 .* w)) for x1 ∈ Iterators.partition(x, nfft))
  20 .* log10.(X)
end

function Base.push!(scope::Scope, data)
  if scope.mode == TIME
    scope.ybuf .= data
  elseif scope.mode == FREQ
    scope.ybuf .= psd(data, scope.nfft)
  elseif scope.mode == TIMEFREQ
    scope.zbuf[1:end-1,:] .= scope.zbuf[2:end,:]
    scope.zbuf[end,:] .= psd(data, scope.nfft)
  end
  lock(scope.waitfordirty) do
    scope.dirty = true
    notify(scope.waitfordirty)
  end
  scope
end

function monitor(scope::Scope)
  try
    while true
      lock(scope.waitfordirty) do
        scope.dirty || wait(scope.waitfordirty)
        scope.dirty = false
      end
      scope.closed && break
      if scope.running[]
        if scope.mode == TIME || scope.mode == FREQ
          scope.y[] .= scope.ybuf
          scope.y[] = scope.y[]
        elseif scope.mode == TIMEFREQ
          scope.z[] .= scope.zbuf
          scope.z[] = scope.z[]
        end
      end
    end
  catch ex
    @error ex
  end
  try
    scope.src === nothing || close(scope.src)
  catch ex
    @warn ex
  end
  display(Figure())
end

function Makie.MakieLayout.autolimits!(scope::Scope)
  autolimits!(scope.fap.axis)
  if scope.mode == FREQ
    ymax = 5 * ceil(Int, maximum(scope.y[])/5)
    ylims!(scope.fap.axis, ymax - scope.zrange[], ymax)
  elseif scope.mode == TIMEFREQ
    scope.zmax[] = 5 * ceil(Int, maximum(scope.z[])/5)
  end
end

function keypress(scope::Scope, action, key)
  action == Keyboard.press || return false
  key == Keyboard._0 && reset_limits!(scope.fap.axis)
  key == Keyboard._1 && (scope.ch[] = 1)
  key == Keyboard._2 && (scope.ch[] = 2)
  key == Keyboard._3 && (scope.ch[] = 3)
  key == Keyboard._4 && (scope.ch[] = 4)
  key == Keyboard.a && autolimits!(scope)
  key == Keyboard.minus && (scope.zmax[] -= 5)
  key == Keyboard.equal && (scope.zmax[] += 5)
  key == Keyboard.t && mode!(scope, TIME)
  key == Keyboard.f && mode!(scope, FREQ)
  key == Keyboard.s && mode!(scope, TIMEFREQ)
  key == Keyboard.q && close(scope)
  key == Keyboard.space && (scope.running[] = !scope.running[])
  key == Keyboard.d && (scope.dc[] = !scope.dc[])
  false
end

annotation(ch, dc, running) = """
  Channel $ch
  $(dc ? "DC coupling" : "AC coupling")
  $(running ? "RUNNING" : "STOPPED")
  """

function annotate(scope::Scope)
  s = lift(annotation, scope.ch, scope.dc, scope.running)
  fig = scope.fap.figure
  fig[1,1] = Label(fig, s;
    textsize = 14.0f0,
    tellheight = false,
    tellwidth = false,
    padding = (10.0f0, 10.0f0, 10.0f0, 10.0f0),
    halign = :right,
    valign = :top)
end


text:
# This file was generated, do not modify it. # hide
Xnew = (Lag1 = [1.2, 1.5], Lag2 = [1.1, -0.8])
ŷ = predict(clf, Xnew)
ŷ |> pprint


text:
#in matlab:
#x=rand(512,1);
#[s,f,t,p]=spectrogram(x,ones(1,256),128,256,10);
#save
#
#in julia:
#using MAT
#
#matdata=matread("matlab.mat")
#
#for i in ("x", "f", "t", "p")
#  fid=open("spectrogram_$i.txt","w")
#  print(fid,matdata["$i"])
#  close(fid)
#end

using DSP, Test
using Statistics: mean
using FFTW: fftfreq

@testset "matlab ref" begin
    x0 = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t'))
    f0 = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_f.txt"),'\t'))
    t0 = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_t.txt"),'\t'))
    p0 = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_p.txt"),'\t')
    spec = spectrogram(x0, 256, 128; fs=10)
    p, f, t = power(spec), freq(spec), time(spec)

    # with real input matlab outputs a 1-sided PSD
    @test p0 ≈ p
    @test f0 ≈ f
    @test t0 ≈ t
end

@testset "0:7" begin
    #Matlab: p = pwelch(0:7, [1, 1, 1, 1, 1, 1, 1, 1], 0, 8, 1, 'twosided')
    data = 0:7
    data0 = Float64[98.0,
                    13.656854249492380,
                    4.0,
                    2.343145750507620,
                    2.0,
                    2.343145750507620,
                    4.0,
                    13.656854249492380]
    @test power(periodogram(data, onesided=false)) ≈ data0
    @test power(welch_pgram(data, length(data), 0, onesided=false)) ≈ data0
    @test power(spectrogram(data, length(data), 0, onesided=false)) ≈ data0
    @test power(periodogram(complex.([data;], [data;]), onesided=false)) ≈ data0*2
    @test power(welch_pgram(complex.([data;], [data;]), length(data), 0, onesided=false)) ≈ data0*2
    @test power(spectrogram(complex.([data;], [data;]), length(data), 0, onesided=false)) ≈ data0*2

    # # ~~~~~~~~ Tests with no window ~~~~~~~~~~~~~~~~~~~
    # Matlab: p = pwelch(0:7, [1, 1], 0, 2, 1, 'twosided')
    expected = Float64[34.5, 0.5]
    @test power(welch_pgram(data, 2, 0; onesided=false)) ≈ expected
    @test mean(power(spectrogram(data, 2, 0; onesided=false)), dims=2) ≈ expected

    # Matlab: p = pwelch(0:7, [1, 1, 1], 0, 3, 1, 'twosided')
    expected = Float64[25.5, 1.0, 1.0]
    @test power(welch_pgram(data, 3, 0; onesided=false)) ≈ expected
    @test mean(power(spectrogram(data, 3, 0; onesided=false)), dims=2) ≈ expected

    # Matlab: p = pwelch(0:7, [1, 1, 1], 1, 3, 1, 'twosided')
    expected = Float64[35.0, 1.0, 1.0]
    @test power(welch_pgram(data, 3, 1; onesided=false)) ≈ expected
    @test mean(power(spectrogram(data, 3, 1; onesided=false)), dims=2) ≈ expected

    # Matlab: p = pwelch(0:7, [1, 1, 1, 1], 1, 4, 1, 'twosided')
    expected = Float64[45, 2, 1, 2]
    @test power(welch_pgram(data, 4, 1; onesided=false)) ≈ expected
    @test mean(power(spectrogram(data, 4, 1; onesided=false)), dims=2) ≈ expected

    # ~~~~~~~~~~~ This one tests periodogram ~~~~~~~~~~~~
    # ~ If functionality of the other arguments has been
    # ~ tested above, we only test here that the correct
    # ~ value of the spectral density is obtained when
    # ~ using a window. More tests to be added if needed
    #Matlab: p = pwelch(0:7, window_func(8), 0, 8, 1, 'twosided')
    cases = (
        (hamming,  Float64[65.461623986801527,
                           20.556791795515764,
                            0.369313143650544,
                            0.022167446610882,
                            0.025502985564107,
                            0.022167446610882,
                            0.369313143650544,
                           20.556791795515764]),
        (bartlett, Float64[62.999999999999993,
                           21.981076052592442,
                            0.285714285714286,
                            0.161781090264695,
                            0.142857142857143,
                            0.161781090264695,
                            0.285714285714286,
                           21.981076052592442])
    )

    @testset for (window1, expected) in cases
        @test power(periodogram(data; window=window1, onesided=false)) ≈ expected
        @test power(welch_pgram(data, length(data), 0; window=window1, onesided=false)) ≈ expected
        @test power(spectrogram(data, length(data), 0; window=window1, onesided=false)) ≈ expected
        @test power(periodogram(data; window=window1(length(data)), onesided=false)) ≈ expected
        @test power(welch_pgram(data, length(data), 0; window=window1(length(data)), onesided=false)) ≈ expected
        @test power(spectrogram(data, length(data), 0; window=window1(length(data)), onesided=false)) ≈ expected
    end

    # Padded periodogram
    # MATLAB: a = periodogram(0:7, [], 32);
    expected = [
                      98
        174.463067389405
        121.968086934209
        65.4971744936088
        27.3137084989848
        12.1737815028909
        10.3755170959439
        10.4034038628775
                       8
        5.25810953219633
        4.47015397150535
        4.89522578856669
        4.68629150101524
        3.69370284475603
         3.1862419983415
        3.61553458569862
                       2
    ]
    @test power(periodogram(data; nfft=32)) ≈ expected
    @test power(welch_pgram(data, length(data), 0; nfft=32)) ≈ expected
    @test power(spectrogram(data, length(data), 0; nfft=32)) ≈ expected

    # Padded periodogram with window
    # MATLAB: a = periodogram(0:7, hamming(8), 32, 1)
    expected = [
          65.4616239868015
          122.101693164395
          98.8444689598445
           69.020252632913
          41.1135835910315
          20.5496474310966
          8.43291449161938
          2.78001620362588
         0.738626287301088
         0.174995741770789
        0.0501563022944516
        0.0327357460012861
        0.0443348932217643
        0.0553999745503552
        0.0561319901616643
        0.0526025934871384
        0.0255029855641069
    ]
    @test power(periodogram(data; window=hamming, nfft=32)) ≈ expected
    @test power(welch_pgram(data, length(data), 0; window=hamming, nfft=32)) ≈ expected
    @test power(spectrogram(data, length(data), 0; window=hamming, nfft=32)) ≈ expected

    # Test fftshift
    p = periodogram(data)
    @test power(p) == power(fftshift(p))
    @test freq(p) ≈ freq(fftshift(p))

    p = periodogram(data; onesided=false)
    @test fftshift(power(p)) == power(fftshift(p))
    @test fftshift(freq(p)) == freq(fftshift(p))
end

@testset "fftshift" begin
    data = 1:100

    p = spectrogram(data)
    @test power(p) == power(fftshift(p))
    @test freq(p) ≈ freq(fftshift(p))

    p = spectrogram(data; onesided=false)
    @test fftshift(power(p), 1) == power(fftshift(p))
    @test fftshift(freq(p)) == freq(fftshift(p))
end

@testset "2D" begin
    data2d = readdlm(joinpath(dirname(@__FILE__), "data", "per2dx.txt"),'\t')
    expectedsum = vec(readdlm(joinpath(dirname(@__FILE__), "data", "per2dsum.txt"),'\t'))
    expectedmean = vec(readdlm(joinpath(dirname(@__FILE__), "data", "per2dmean.txt"),'\t'))
    # 2-d periodgram (radialsum)
    # computed in octave with raPsd2d ((C) E. Ruzanski) replacing nanmean with nansum
    # P = raPsd2d(x,1)'*n^2
    @test power(periodogram(data2d,fs=1, radialsum=true)) ≈ expectedsum

    # 2-d periodgram (radialavg)
    # computed in octave with raPsd2d ((C) E. Ruzanski)
    # P = raPsd2d(x,1)'*n^2
    @test power(periodogram(data2d, fs=1, radialavg=true)) ≈ expectedmean

    # 2-d periodgram 2-d PSD
    @test power(periodogram(data2d, fs=1)) ≈ abs2.(fft(data2d))*1/prod(size(data2d))
    # 2-d periodgram 2-d PSD with padding
    pads = (size(data2d,1)+4,size(data2d,1)+7)
    data2dpad = zeros(Float64,pads...)
    data2dpad[1:size(data2d,1),1:size(data2d,2)] = data2d
    @test power(periodogram(data2d, fs=1, nfft=pads)) ≈ abs2.(fft(data2dpad))*1/prod(size(data2d))
    # 2-d periodgram radial freq
    @test freq(periodogram(data2d, fs=3.3, radialsum=true)) ≈ freq(periodogram(vec(data2d[1,:]), fs=3.3))
    # 2-d periodgram 2-d freq
    f1,f2 = freq(periodogram(data2d, fs=3.3))
    f1d = freq(periodogram(vec(data2d[1,:]), fs=3.3, onesided=false))
    @assert size(data2d,1)==size(data2d,2)
    for j=1:size(data2d,2)
        for i=1:size(data2d,1)
            @test  [f1[i],f2[j]] ≈ [f1d[i],f1d[j]]
        end
    end
    # Test fftshift
    p = periodogram(data2d)
    @test fftshift(power(p)) == power(fftshift(p))
    f = freq(p)
    @test (fftshift(f[1]),fftshift(f[2])) == freq(fftshift(p))
end

@testset "radial" begin
    # 2-d periodgram radial test for a non-square signal sparse in fft space
    n1 = 52
    n2 = 46  # assuming n1>n2
    nf = (22,7) # the non-zero location
    F = (fftfreq(n1,1),fftfreq(n2,1))
    a = [F[1][nf[1]],F[2][nf[2]]]
    FB = Bool[[F[1][i], F[2][j]]==a || [F[1][i], F[2][j]]==-a for i = 1:n1, j = 1:n2]
    ind = findall(FB)
    x = zeros(n1,n2)*0im;
    x[ind] = [1+2im,1-2im]
    y = real(ifft(x))

    fwn = round(Int, sqrt((a[1])^2+(a[2])^2)*n2)
    pe = zeros(n2>>1 + 1)
    pe[fwn+1] = 2*abs2(x[nf...])/n1/n2
    P = periodogram(y,nfft=(n1,n2),radialsum=true)
    @test power(P) ≈ pe
    @test freq(P)[fwn+1] ≈ fwn/n2
end

# Testing STFT function and comparing results with MATLAB
@testset "stft" begin
    fs = 16000
    nfft = 512
    nwin = 400
    nhop = 160
    s = vec(readdlm(joinpath(dirname(@__FILE__), "data", "stft_x.txt"),'\t'))

    Sjl = stft(s, nwin, nwin-nhop; nfft=nfft, fs=fs, window=hanning)
    Sml_re = readdlm(joinpath(dirname(@__FILE__), "data", "stft_S_real.txt"),'\t')
    Sml_im = readdlm(joinpath(dirname(@__FILE__), "data", "stft_S_imag.txt"),'\t')
    Sml = complex.(Sml_re, Sml_im)
    @test Sjl ≈ Sml
end

@testset "fft2oneortwosided!" begin
    n = 10
    floatt = Float64
    for onesided in (true, false),
            nfft in (n, n+2, n+3),
           atype in (floatt, Complex{floatt})
        nout = nout = onesided ? (nfft >> 1)+1 : nfft
        x = zeros(atype, nfft)
        if atype <: Real
            x[1:n] = rand(atype, n)
            xrcfft = rfft(x)
        else
            x[1:n] = rand(floatt, n)+im*rand(floatt, n)
            xrcfft = fft(x)
        end
        xfft = fft(x)
        out = zeros(fftouttype(atype),nout,3)
        if !(onesided == true && atype <: Complex)
            outft = DSP.Periodograms.fft2oneortwosided!(out, xrcfft, nfft, onesided, nout)
        end
        if onesided == true && atype <: Real
            @test out[:,2] ≈ xrcfft
            @test out[:,[1,3]] ≈ xrcfft.*[0 0]
        elseif onesided == false && atype <: Real
            @test out[:,2] ≈ xfft
            @test out[:,[1,3]] ≈ xfft.*[0 0]
        elseif onesided == false && atype <: Complex
            @test out[:,2] ≈ xfft
            @test out[:,[1,3]] ≈ xfft.*[0 0]
        else
            #onesided complex
        end
    end
end

@testset "mt_pgram" begin
    # MATLAB: x = pmtm(stft_x, 4, 5000, 16000, 'unity')
    s = vec(readdlm(joinpath(dirname(@__FILE__), "data", "stft_x.txt"),'\t'))
    mtdata = vec(readdlm(joinpath(dirname(@__FILE__), "data", "mt_pgram.txt")))
    @test power(mt_pgram(s; fs=16000)) ≈ mtdata
    @test power(mt_pgram(s; fs=16000, window=dpss(length(s), 4))) ≈ mtdata

    # error tests
    @test_throws ArgumentError periodogram([1 2 3])
    @test_throws ArgumentError periodogram(rand(2,3), nfft=(3,2))
    @test_throws ArgumentError periodogram([1 2;3 4],radialsum=true, radialavg=true)

    # #124
    q = arraysplit(ones(Float64, 1000),100,10);
    @test map(mean, q) == ones(Float64, 11)

    # test that iterating ArraySplit always yields the same Vector
    q = arraysplit(-10:10,4,2)
    for x in q
        @assert isa(x, Vector)
        @test pointer(x) == pointer(q.buf)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ ff81125a-086d-11eb-2226-eddf1955c141
begin
	# Activate environment. You do not need to do this if you are not using
	# separate environment
	push!(LOAD_PATH, "$(homedir())/pCloudDrive/julia_project")
	using Pkg
	Pkg.activate("/home/yusri/pCloudDrive/julia_project/MovingFrame/")
end

# ╔═╡ 5df5c832-ee72-11ea-20bf-7b42b7e25571
begin
	using GeometryBasics
	using FileIO
	using LinearAlgebra;
	using GaussianProcesses;
	using PlutoUI
	using MovingFrame;
	import Plots;
	const mva = MovingFrame
	const plt = Plots
end

# ╔═╡ 67beff68-ee73-11ea-3b77-fbf331f8ce9c
# Create MovingCells Type
struct MovingCells
    coordinates::Matrix{Float64}
    time::Vector{Float64}
    smooth_n_points::Int64
    smooth_coordinates::Matrix{Float64}
    smooth_time::Vector{Float64}
    T::Matrix{Float64}
    N::Matrix{Float64}
    B::Matrix{Float64}
    curvature::Vector{Float64}
    torsion::Vector{Float64}
    speed::Vector{Float64}
    original_shapes::Vector{GeometryBasics.Mesh}
    reoriented_shapes::Vector{GeometryBasics.Mesh}
    Spharm_Cx::Matrix{Float64}
    Spharm_Cy::Matrix{Float64}
    Spharm_Cz::Matrix{Float64}
end

# ╔═╡ 7139c988-ee73-11ea-1127-25bfb372000d
begin
	const base_dir = pwd();
	# You need to extract the cell_shape_final.zip first
	const input_data_dir = joinpath(base_dir, "cell_shape_final/")
	const input_files = readdir(input_data_dir);
	const n_points = length(input_files)
end

# ╔═╡ 9ebba1d6-ee73-11ea-1de2-d78152ef6b5e
begin
	# Calculating the cell trajectory
	coordinate = Matrix{Float64}(undef, n_points, 3)
	time = collect(1.0:n_points)
	shapes = Vector{GeometryBasics.Mesh}(undef, n_points)
	for i = 1:length(input_files)
		shape = load(joinpath(input_data_dir, input_files[i]))
		shapes[i] = shape
		verts = decompose(Point3f0, shape)
		# The trajectory is the mass centre of the cell
		coordinate[i, :] .= reduce(+, verts) / length(verts)
	end
	# Normalizing volume and translation
	shapes = (volume_normalizing ∘ mesh_centering).(shapes)
	new_coordinate = coordinate .- coordinate[1, :]'
end

# ╔═╡ 7d1cf63e-ee83-11ea-33e9-651efeecaf00
begin
	exp_lengthscale = 3.4 # lengthscale = e^3.4
	exp_variance= 0.0 #variance = e^0.0

	# Gaussian Processs using Zero mean function and
	# Squared exponential kernel (note that hyperparameters are on the log scale)
	smoother =  GaussianProcessSmoother(MeanZero(), SE(exp_lengthscale, exp_variance))
	smooth_coordinates = smoothing(new_coordinate, smoother, step=1.)

	smooth_time = range(1.0, stop=time[end], length=size(smooth_coordinates, 1)) |> collect
	T, N, B, curvature, torsion, speed = parallel_transport(smooth_coordinates)
end

# ╔═╡ 5b0d2618-0878-11eb-0fe0-69b9357b91e9
# Reoriented the cell shape to the Moving Frame basis
reoriented_shape = [
    mesh_reorientation(shapes[i], T[i, :], N[i, :], B[i, :])
    for i in eachindex(shapes)
]

# ╔═╡ 29c8a4b8-087a-11eb-3780-4fa5a2c12d0e
# Spherical parameterization
spheres = MovingFrame.spherization.(reoriented_shape, n_step=100)

# ╔═╡ ed0b2ee6-090d-11eb-3516-9b7788a1859b
begin
	# Calculating spherical harmonics coefficients
	spharm_coefs = [
	    MovingFrame.spharm_descriptor(
	        reoriented_shape[i],
	        GeometryBasics.coordinates(spheres[i]),
	    ) for i in eachindex(reoriented_shape)
	]
	Cx = MovingFrame.array_of_array_to_matrix([coef[:Cx] for coef in spharm_coefs])
	Cy = MovingFrame.array_of_array_to_matrix([coef[:Cy] for coef in spharm_coefs])
	Cz = MovingFrame.array_of_array_to_matrix([coef[:Cz] for coef in spharm_coefs])
end

# ╔═╡ 18462938-090e-11eb-2154-45881e7ef9f1
begin
	# Calculating shape eccentricity
	xy_ecc = Cx[:, 4] ./ Cy[:, 2]
	xz_ecc = Cx[:, 4] ./ Cz[:, 3]
	yz_ecc = Cy[:, 2] ./ Cz[:, 3]
end

# ╔═╡ cd0cd780-ee89-11ea-0e25-45fdbb85346d
begin
plt.plot(smooth_coordinates[:,1],smooth_coordinates[:,2],smooth_coordinates[:,3],label ="SMOOTH_TRAJECTORY")
plt.scatter!(new_coordinate[:,1],new_coordinate[:,2],new_coordinate[:,3],label ="OBSERVATION",markersize=1.0)

end

# ╔═╡ 46ff40b6-090e-11eb-1216-b1a6a3af5695
plt.plot(smooth_time,[xy_ecc xz_ecc],label=["XY Eccentricity" "XZ Eccentricity"])

# ╔═╡ 04cdf2de-3dcd-11eb-2123-9fef1cb5553e
begin
	p1 = plt.plot(smooth_time,speed,ylabel = "Speed")
	p2 = plt.plot(smooth_time,curvature,ylabel = "Curvature",)
	p3 = plt.plot(smooth_time,torsion,ylabel="Torsion",)
	plt.plot(p1,p2,p3,legend=false)
end

# ╔═╡ Cell order:
# ╠═ff81125a-086d-11eb-2226-eddf1955c141
# ╠═5df5c832-ee72-11ea-20bf-7b42b7e25571
# ╠═67beff68-ee73-11ea-3b77-fbf331f8ce9c
# ╠═7139c988-ee73-11ea-1127-25bfb372000d
# ╠═9ebba1d6-ee73-11ea-1de2-d78152ef6b5e
# ╠═7d1cf63e-ee83-11ea-33e9-651efeecaf00
# ╠═5b0d2618-0878-11eb-0fe0-69b9357b91e9
# ╠═29c8a4b8-087a-11eb-3780-4fa5a2c12d0e
# ╠═ed0b2ee6-090d-11eb-3516-9b7788a1859b
# ╠═18462938-090e-11eb-2154-45881e7ef9f1
# ╠═cd0cd780-ee89-11ea-0e25-45fdbb85346d
# ╠═46ff40b6-090e-11eb-1216-b1a6a3af5695
# ╠═04cdf2de-3dcd-11eb-2123-9fef1cb5553e


text:
#=============================#
#== put needed modules here ==#
#=============================#
import JuliaChem

#================================#
#== JuliaChem execution script ==#
#================================#
include("minimal-rhf-repl.jl")

JuliaChem.initialize()
minimal_rhf(ARGS[1])
JuliaChem.finalize()


text:
immutable TriangularDist <: ContinuousUnivariateDistribution
    a::Float64
    b::Float64
    c::Float64

    function TriangularDist(a::Real, b::Real, c::Real)
        a < b || error("TriangularDist: a < b must be true")
        a <= c <= b || error("a <= c <= b must be true")
        @compat new(Float64(a), Float64(b), Float64(c))
    end

    function TriangularDist(a::Real, b::Real)
        a < b || error("TriangularDist: a < b must be true")
        @compat a_ = Float64(a)
        @compat b_ = Float64(b)
        c_ = middle(a_, b_)
        new(a_, b_, c_)
    end
end

@distr_support TriangularDist d.a d.b

#### Parameters

params(d::TriangularDist) = (d.a, d.b, d.c)

#### Statistics

mode(d::TriangularDist) = d.c

mean(d::TriangularDist) = (d.a + d.b + d.c) / 3.0

function median(d::TriangularDist)
    (a, b, c) = params(d)
    m = middle(a, b)
    c >= m ? a + sqrt(0.5 * (b - a) * (c - a)) :
             b - sqrt(0.5 * (b - a) * (b - c))
end

_pretvar(a::Float64, b::Float64, c::Float64) = a*a + b*b + c*c - a*b - a*c - b*c

function var(d::TriangularDist)
    (a, b, c) = params(d)
    _pretvar(a, b, c) / 18.0
end

function skewness(d::TriangularDist)
    (a, b, c) = params(d)
    sqrt2 * (a + b - 2.0c) * (2.0a - b - c) * (a - 2.0b + c) / (5.0 * _pretvar(a, b, c)^1.5)
end

kurtosis(d::TriangularDist) = -0.6

entropy(d::TriangularDist) = 0.5 + log((d.b - d.a) / 2.0)

#### Evaluation

function pdf(d::TriangularDist, x::Float64)
    (a, b, c) = params(d)
    x <= a ? 0.0 :
    x <  c ? 2.0 * (x - a) / ((b - a) * (c - a)) :
    x == c ? 2.0 / (b - a) :
    x <= b ? 2.0 * (b - x) / ((b - a) * (b - c)) : 0.0
end

function cdf(d::TriangularDist, x::Float64)
    (a, b, c) = params(d)
    x <= a ? 0.0 :
    x <  c ? (x - a)^2 / ((b - a) * (c - a)) :
    x == c ? (c - a) / (b - a) :
    x <= b ? 1.0 - (b - x)^2 / ((b - a) * (b - c)) : 1.0
end

function quantile(d::TriangularDist, p::Float64)
    (a, b, c) = params(d)
    c_m_a = c - a
    b_m_a = b - a
    rl = c_m_a / b_m_a
    p <= rl ? a + sqrt(b_m_a * c_m_a * p) :
              b - sqrt(b_m_a * (b - c) * (1.0 - p))
end

function mgf(d::TriangularDist, t::Real)
    if t == zero(t)
        return one(t)
    else
        (a, b, c) = params(d)
        u = (b - c) * exp(a * t) - (b - a) * exp(c * t) + (c - a) * exp(b * t)
        v = (b - a) * (c - a) * (b - c) * t^2
        return 2.0 * u / v
    end
end

function cf(d::TriangularDist, t::Real)
    # Is this correct?
    if t == zero(t)
        return one(t)
    else
        (a, b, c) = params(d)
        u = (b - c) * cis(a * t) - (b - a) * cis(c * t) + (c - a) * cis(b * t)
        v = (b - a) * (c - a) * (b - c) * t^2
        return -2.0 * u / v
    end
end

#### Sampling

function rand(d::TriangularDist)
    (a, b, c) = params(d)
    b_m_a = b - a
    u = rand()
    b_m_a * u < (c - a) ? d.a + sqrt(u * b_m_a * (c - a)) :
                          d.b - sqrt((1.0 - u) * b_m_a * (b - c))
end


text:
module ModelTES
export
    transitionwidth,
    BiasedTES,
    TESParams,
    getlinearparams,
    NoiseModel,
    TESRecord,
    times,
    rk8,
    IrwinHiltonTES,
    stochastic,
    pulse

using Roots, ForwardDiff, DifferentialEquations
include("rk8.jl")
include("tes_models.jl")

const J_per_eV = 1.602177e-19 #unitless
const kb = 1.38064852e-23 #k boltzmann (J/K)

abstract type AbstractRIT end

# following Irwin-Hilton figure 3
mutable struct TESParams{RITType<:AbstractRIT}
    n       ::Float64   # thermal conductance exponent (unitless)
    Tc      ::Float64  # TES critical temperature (K)
    Tbath   ::Float64   # bath temperature (K)

    k       ::Float64   # Pbath pre-factor (W/K^n)
    C       ::Float64   # heat capacity of TES (J/K)

    L       ::Float64  # inductance of SQUID (H)
    Rl      ::Float64   # Thevenin-equivalent resistance Rload = (Rshunt+Rparasitic)(ohms)
                        # note Ibias = V/Rshunt
    Rp      ::Float64   # Rparastic; Rshunt = Rl-Rparasitic (ohms)
    Rn      ::Float64   # normal resistance of TES (ohms)

    RIT     ::RITType   # RIT surface
end

mutable struct ShankRIT <: AbstractRIT
    Tw::Float64 # transition width (K)
    A::Float64  # current dependence for R(T,I) (A/K^(3/2))
end
transitionwidth(RIT::ShankRIT)=RIT.Tw
transitionwidth(p::TESParams)=transitionwidth(p.RIT)

"Constructor that fixes `Tw` and `A` for a ShankRIT to have the given `alpha` and `beta`
parameters when biased at resistance `R0`."
function ShankRIT(alpha, beta, n, Tc, Tbath, k, R0, Rn)
    T0 = Tc /(1 + 3*beta/(2*alpha) - 2*(Rn-R0)/(Rn*alpha)*atanh(2*R0/Rn-1))
    I0 = sqrt(k*(T0^n-Tbath^n)/R0)
    Tw = T0*(Rn-R0)/(Rn*log(2)*alpha)
    A = I0*(2*alpha/(3*T0*beta))^(3/2)
    ShankRIT(Tw,A)
end

mutable struct BiasedTES{T}
    p::TESParams{T}
    I0::Float64 # intial current for diff eq, aka current through TES (A)
    T0::Float64 # initial temperature for diff equations, aka temperature of TES (K)
    V ::Float64 # Thévinen equivalent voltage V = I0*(p.Rl+p.R0), R0=quiescent resistance
                # also equal to Ibias*Rshunt. Careful! This V is a constant and is NOT
                # the voltage drop across the TES.
end

"For a given T0, the difference R-targetR. Use to solve numerically for T0."
function getR0error(T0, targetR, p::TESParams)
   I0 = sqrt(p.k*(T0^p.n-p.Tbath^p.n)/targetR)
   R(I0,T0,p)-targetR
end

"Find the initial conditions (I0, T0, V) that cause `p` to have resistance `targetR`."
function initialconditions(p::TESParams, targetR)
   T00 = fzero((t)->getR0error(t, targetR,p),p.Tc-10*transitionwidth(p), p.Tc+10*transitionwidth(p))
   I00 = sqrt(p.k*(T00^p.n-p.Tbath^p.n)/targetR)
   R00 = R(I00,T00,p)
   V00 = I00*(p.Rl+R00)
   # now evolve these conditions through integration to really lock them in.
   # shouldn't hard code step size here
   out = pulse(10,1e-1, BiasedTES(p, I00, T00, V00), 0)
   T0 = out.T[end]
   I0 = out.I[end]
   R0 = R(I0,T0,p)
   V = I0*(p.Rl+R0)
   I0,T0,V
end

"Created a biased tes with quiescent state resistance R0"
function BiasedTES(p::TESParams{T}, R0::Float64) where T
   @assert 0 < R0 < p.Rn
   I0, T0, V = initialconditions(p,R0)
   BiasedTES(p,I0,T0,V)
end

"iv_point(p::TESParams, V, I0, T0)
takes thevinin voltage `V`, and initial current `I0`, and intial temperature `T0`
evolves a pulse for 1 second, and takes the final values
returns I,T,V,R"
function iv_point(p::TESParams, V, I0, T0)
    # solve with an adapative algorithm that is fast for large time steps,
    # ask for very long time steps
    # we probably shouldn't hardcode the time, but 1 second is long for all TESs I know of
    out = pulse(2,1.0, BiasedTES(p, I0, T0, V), 0, method=DifferentialEquations.Rosenbrock23())
    T = out.T[end]
    I = out.I[end]
    R = out.R[end]
    V = I*(p.Rl+R)
    I,T,R,V
end

"iv_curve(p::TESParams, Vs)
takes a sorted array of V thevinin values `Vs`, calculates ivs points by
evolving a pulse for 1 second, and taking the last value
returns Is, Ts, Vs_out, Rs"
function iv_curve(p::TESParams, Vs)
    Is=Vector{Float64}(undef, length(Vs))
    Ts=Vector{Float64}(undef, length(Vs))
    Rs=Vector{Float64}(undef, length(Vs))
    Vs_out=Vector{Float64}(undef, length(Vs))
    @assert issorted(Vs)
    for i in length(Vs):-1:1
            if Vs[i]==0
                I,T,R,V=0.0, p.Tbath, 0.0, 0.0
            elseif i==length(Vs)
                # provide guesses that guaranteed to be resistive
                I,T,R,V = iv_point(p, Vs[i], Vs[i]/p.Rn, p.Tc)
            elseif i<length(Vs)
                #provide last solution as guesses
                # I got a speedup by providing nearby starting points, it was less than a factor of 2
                I,T,R,V = iv_point(p, Vs[i], Is[i+1], Ts[i+1])
            end
            Is[i]=I
            Ts[i]=T
            Rs[i]=R
            Vs_out[i]=V
    end
    Is,Ts,Rs,Vs_out
end

function dT_and_dI_iv_point(p,I,T,R,V)
    bt = BiasedTES(p, I, T, V)
    du = zeros(2)
    u = [T,I]
    bt(0.0, u, du)
    du
end

"Calculate `R0` the quiescent resistance of `tes`."
getR0(tes::BiasedTES) = R(tes.I0, tes.T0, tes.p)

"Calculate thermal conductivity from TES heat capacity to Tbath in small signal limit."
function getG0(tes::BiasedTES)
   p=tes.p
   p.k*p.n*tes.T0^(p.n-1)
end

"Calculate paramaters in Irwin-Hilton table 1."
function getlinearparams(bt::BiasedTES)
   f(x) = R(x[1], x[2], bt.p)
   p = bt.p
   R0 = getR0(bt)
   G0 = getG0(bt)
   drdi,drdt = ForwardDiff.gradient(f, [bt.I0, bt.T0])
   alpha = drdt*bt.T0/R0
   beta = drdi*bt.I0/R0
   PJ = bt.I0^2*R0
   loopgain = PJ*alpha/G0/bt.T0
   tauthermal = p.C/G0
   taucc = tauthermal/(1-loopgain) # constant current time constant
   r = p.Rl/R0
   taueff = (1+beta+r)/(1+beta+r+(1-r)*loopgain) # zero inductance effective thermal time constant
   tauelectrical = p.L/(p.Rl+R0*(1+beta))
   invtau = 1/(2*tauelectrical)+1/(2*taucc)
   a = (1/tauelectrical-1/taucc)^2
   b = -4*(R0/p.L)*loopgain*(2+beta)/tauthermal
   lcritical = -p.L*b/a
   invtaupm = 0.5*sqrt(a+b+0*im) # make it complex, so I can get a complex answer
   tauplus = 1/(invtau+invtaupm)
   tauminus = 1/(invtau-invtaupm)
   c = loopgain*(3+beta-r)+(1+beta+r)
   d = 2*sqrt(loopgain*(2+beta)*(loopgain*(1-r)+(1+beta+r)))
   f = R0*tauthermal/(loopgain-1)^2
   lcritplus = (c+d)*f
   lcritminus = (c-d)*f
   bt.I0, bt.T0, bt.V, p.Rl, p.Tbath, p.Tbath, p.L, R0, G0, p.C, alpha, beta, loopgain,
        tauthermal, taucc, taueff, tauelectrical, tauplus, tauminus, lcritplus, lcritminus, lcritical
end

"Paramters from Irwin-Hilton table one for modeling a linear TES. Defined in Table 1 of Irwin-Hilton chapter."
mutable struct IrwinHiltonTES
   I0::Float64
   T0::Float64
   V::Float64
   Rl::Float64
   Tl::Float64 # temperature of the load resistor, usually modeled as =Tbath, but really should come from modeling ep coupling in the load resistor
   Tbath::Float64
   L::Float64
   R0::Float64
   G0::Float64
   C0::Float64
   alpha::Float64
   beta::Float64 #βI
   loopgain::Float64 #ℒI
   tauthermal::Float64
   taucc::Float64 #τI
   taueff::Float64
   tauelectrical::Float64 #τel
   tauplus::Complex{Float64}
   tauminus::Complex{Float64}
   lcritplus::Float64
   lcritminus::Float64
   Lcritical::Float64
end

"Constructor for IrwinHiltonTES model"
function IrwinHiltonTES(tes::BiasedTES)
    linearparams = getlinearparams(tes)
    IrwinHiltonTES(linearparams...)
end

"`Z(tes::IrwinHiltonTES, f)`

Returns the impedance of the `tes` at frequency `f`.
Implements equation 42 of Irwin-Hilton chapter."
function Z(tes::IrwinHiltonTES, f)
  ω=2π*f
  tes.R0*(1+tes.beta) .+ tes.R0*tes.loopgain*(2+tes.beta)./((1-tes.loopgain)*(1 .+im*ω*tes.taucc))
end

"`Zcircuit(tes::IrwinHiltonTES, f)`

Returns impedance of complete circuit of `tes` at frequency `f`."
function Zcircuit(tes::IrwinHiltonTES, f)
  ω=2π*f
  tes.R0 .+ im*ω*tes.L .+ Z(tes,ω)
end

include("tes_noise.jl")
include("stochastic_integration.jl")

isoverdamped(tes::IrwinHiltonTES) = (isreal(tes.tauplus) && isreal(tes.tauminus) &&
                                    tes.tauplus<tes.tauminus)
isunderdamped(tes::IrwinHiltonTES) = (!isoverdamped(tes) && tes.tauplus!=tes.tauminus)

mutable struct TESRecord
    T::Vector{Float64}  # temperature (K)
    I::Vector{Float64}  # current (A)
    R::Vector{Float64}  # TES resistance (Ohm)
    dt::Float64   # seconds between samples (seconds)
end
times(r::TESRecord) = range(0, step=r.dt, length=length(r.I))
Base.length(r::TESRecord) = length(r.I)

# I want to re-write this code so that
# A and Tw are the TES params, alpha and beta are derived
# and R0 and be changed to change the bias conditions
# then I0 and T0 and derived from those 3 parametrs

"find bias parameters T0, I0, V for to make the TES have resistane R0 in steady state"
function findbiasparams(tes::TESParams{ShankRIT},R0)

    T0 = Tc /(1 + 3*beta/(2*alpha) - 2*(Rn-R0)/(Rn*alpha)*atanh(2*R0/Rn-1))

    # V = I0*(p.Rl+p.R0)
end

"TES resistance model (Shank et al. 2014)"
R(I, T, RIT::ShankRIT, Tc, Rn) = Rn/2*(1+tanh.((T-Tc+(max.(I,0.0)/RIT.A).^(2/3))/(2*log(2)*RIT.Tw)))
R(I,T, p::TESParams) = R(I,T, p.RIT, p.Tc, p.Rn)

"thermal TES equation
C*dT/dt = I^2*R - k*(T^n-Tbath^n)"
function dT(I, T, k, n, Tbath, C, R)
    # if you run into domain errors, try uncommenting the following line
    # but it probably represents a mistake in your timesteps that should be fixed.
    #T=max(T,0.0) # avoid domain errors when raising T to a power
    Q=-k*(T^n-Tbath^n)+I^2*R
    Q/C

end

"electrical TES equation
L*di/dt = (IBias-I)*Rs+I*Rs-I*Rtes"
function dI(I, T, V, Rl, L, R)
    (V-I*(Rl+R))/L
end

"Calling a BiasedTES gives the dI and dT terms for integration in an in place manner."
function (bt::BiasedTES)(du,u,p_,t) # use DifferentialEquation 4.0+ API
    T,I = u[1],u[2]
    p = bt.p
    r = R(I,T,p)
    # dT(I, T, p.k, p.n, p.Tbath, p.C, r)
    # dI(I,T, bt.V, p.Rl, p.L, r)
    du[1] = dT(I, T, p.k, p.n, p.Tbath, p.C, r)
    du[2] = dI(I,T, bt.V, p.Rl, p.L, r)
end
function (bt::BiasedTES)(t, u, du) # legacy API for rk8
    bt(du,u,nothing,t)
end

function rk8(nsample::Int, dt::Float64, bt::BiasedTES, E::Vector, npresamples::Int=0)
    out = Vector{TESRecord}(length(E))
    for i in 1:length(E)
        out[i] = rk8(nsample, dt, bt, E[i], npresamples)
    end
    out
end

function rk8(nsample::Int, dt::Float64, bt::BiasedTES, E::Number, npresamples::Int=0)
    # Pair of differential equations y' = f(t,y), where y=[T,I]
    p = bt.p
    # Integrate pair of ODEs for all energies EE
    T = Array{Float64}(undef, nsample)
    I = Array{Float64}(undef, nsample)
    T[1:npresamples] .= bt.T0
    I[1:npresamples] .= bt.I0 # set T0, I0 for presamples
    y = [bt.T0+E*J_per_eV/p.C, bt.I0]; ys = similar(y); work = Array{Float64}(undef, 14)
    T[npresamples+1]=y[1]
    I[npresamples+1]=y[2]
    # npresamples+1 is the point at which initial conditions hold (T differs from T0)
    # npresamples+2 is the first point at which I differs from I0
    for i = npresamples+2:nsample
        rk8!(bt, 0.0, dt, y, ys, work)
        y[:] = ys
        T[i] = y[1]
        I[i] = y[2]
    end
    Rout = [R(I[i],T[i],bt.p) for i=1:length(T)]

    TESRecord(T, I, Rout, dt)
end

# example of using the DifferentialEquations API to solve the relevant equations
function adaptive_solve(bt::BiasedTES, dt::Float64, tspan::Tuple{Float64,Float64}, E::Number, method, abstol, reltol, saveat)
    u0 = [bt.T0+E*ModelTES.J_per_eV/bt.p.C, bt.I0]
    prob = ODEProblem(bt, u0, tspan)
    sol = solve(prob,method,dt=dt,abstol=abstol,reltol=reltol, saveat=saveat, save_everystep=false, dense=false)
end
"pulse(nsample::Int, dt::Float64, bt::BiasedTES, E::Number, npresamples::Int=0; dtsolver=1e-9, method=DifferentialEquations.Tsit5(), abstol=1e-9, reltol=1e-9)"
function pulse(nsample::Int, dt::Float64, bt::BiasedTES, E::Number, npresamples::Int=0; dtsolver=1e-9, method=DifferentialEquations.Tsit5(), abstol=1e-9, reltol=1e-9)
    u0 = [bt.T0+E*ModelTES.J_per_eV/bt.p.C, bt.I0]
    saveat = range(0, step=dt, length=nsample-npresamples)
    prob = ODEProblem(bt, u0, (0.0, last(saveat)))
    sol = solve(prob,method,dt=dtsolver,abstol=abstol,reltol=reltol, saveat=saveat, save_everystep=false, dense=false)
    # npresamples+1 is the point at which initial conditions hold (T differs from T0) (sol[1])
    # npresamples+2 is the first point at which I differs from I0
    T = Vector{Float64}(undef, nsample)
    I = Vector{Float64}(undef, nsample)
    T[npresamples+1:end] = sol[1,:]
    I[npresamples+1:end] = sol[2,:]
    T[1:npresamples] .= bt.T0
    I[1:npresamples] .= bt.I0
    Rout = [R(I[i],T[i],bt.p) for i=1:length(T)]
    TESRecord(T,I, Rout,dt)
end

function pulses(nsample::Int, dt::Float64, bt::BiasedTES, Es::Vector, arrivaltimes::Vector; dtsolver=1e-9, method=DifferentialEquations.Tsit5(), abstol=1e-9, reltol=1e-9)
  u0 = [bt.T0, bt.I0]
  saveat = range(0, step=dt, length=nsample)
  prob = ODEProblem(bt, u0, (0.0, last(saveat)))
  Esdict = Dict([(at,E) for (at,E) in zip(arrivaltimes,Es)])
  # this defines a callback that is evaluated when t equals a value in arrival times
  # when evaluated it discontinuously changes the temperature (u[1])
  # the last (true,true) argument has to do with which points are saved
  # it doesn't appear to add extra points beyond saveat, so maybe save_everystep overrides it?
  function cbfun(integrator)
    integrator.u[1]+=Esdict[integrator.t]*ModelTES.J_per_eV/bt.p.C
    # modify the integrator timestep back to dtsolver, to take small steps on the rising edge of the pulse
    integrator.dtpropose=dtsolver # in future use modify_proposed_dt!, see http://docs.juliadiffeq.org/latest/basics/integrator.html#Stepping-Controls-1
  end
  cb = DiscreteCallback((u,t,integrator)->(t in arrivaltimes), cbfun; save_positions=(false,false))
  # tstops is used to make sure the integrator checks each time in arrivaltimes
  sol = solve(prob,method,dt=dtsolver,abstol=abstol,reltol=reltol, saveat=saveat, save_everystep=false, dense=false,callback=cb, tstops=arrivaltimes)

  T = sol[1,:]
  I = sol[2,:]
  Rout = [R(I[i],T[i],bt.p) for i=1:length(T)]
  TESRecord(T,I, Rout,dt)
end

end # module


text:
"""
    GLGM06{N, AM, S, D, NG, P, RM} <: AbstractContinuousPost

Implementation of Girard - Le Guernic - Maler algorithm for reachability of
linear systems using zonotopes.

## Fields

- `δ`                -- step-size of the discretization
- `approx_model`     -- (optional, default: `Forward`) approximation model;
                        see `Notes` below for possible options
- `max_order`        -- (optional, default: `5`) maximum zonotope order
- `static`           -- (optional, default: `false`) if `true`, convert the problem data
                        to statically sized arrays
- `dim`              -- (optional default: `missing`) ambient dimension
- `ngens`            -- (optional, default: `missing`) number of generators
- `preallocate`      -- (optional, default: `true`) if `true`, use the implementation
                        which preallocates the zonotopes prior to applying the update rule
- `reduction_method`    -- (optional, default: `GIR05()`) zonotope order reduction method used
- `disjointness_method` -- (optional, default: `NoEnclosure()`) method to check
                           disjointness between the reach-set and the invariant

## Notes

The type fields are:

- `N`  -- number type of the step-size
- `AM` -- approximation model
- `S`  -- value type associated to the `static` option
- `D`  -- value type associated to the dimension of the system
- `NG` -- value type associated to the number of generators
- `P`  -- value type associated to the `preallocate` option
- `RM` -- type associated to the reduction method

The sole parameter which doesn't have a default value is the step-size,
associated to the type parameter `N`. Parameters `D` and `NG` are optionally
specified (default to `Missing`). These parameters are needed for implementations
that require the size of the zonotopes to be known (fixed) at compile time, namely
the `static=true` version of this algorithm. Otherwise, the number of generators
is not necessarily fixed.

The default approximation model is

```julia
approx_model=Forward(sih=:concrete, exp=:base, setops=:lazy)
```
Here, `Forward` refers to the forward-time adaptation of the approximation model
from Lemma 3 in [[FRE11]](@ref). Some of the options to compute this approximation can be specified,
see the documentation of `Forward` for details.

## References

The main ideas behind this algorithm can be found in [[GIR05]](@ref) and [[GLGM06]](@ref).
These methods are discussed at length in the dissertation [[LG09]](@ref).

Regarding the zonotope order reduction methods, we refer to [[COMB03]](@ref),
[[GIR05]](@ref) and the review article [[YS18]](@ref).

Regarding the approximation model, we use an adaptation of a result in [[FRE11]](@ref).
"""
struct GLGM06{N, AM, S, D, NG, P, RM, DM} <: AbstractContinuousPost
    δ::N
    approx_model::AM
    max_order::Int
    static::S
    dim::D
    ngens::NG
    preallocate::P
    reduction_method::RM
    disjointness_method::DM
end

# TODO review setops "zonotope" / "lazy" options, used or ignored

# convenience constructor using symbols
function GLGM06(; δ::N,
               approx_model::AM=Forward(sih=:concrete, exp=:base, setops=:lazy), #CorrectionHull(order=10, exp=:base),
               max_order::Int=5,
               static::Bool=false,
               dim::Union{Int, Missing}=missing,
               ngens::Union{Int, Missing}=missing,
               preallocate::Bool=true,
               reduction_method::RM=GIR05(),
               disjointness_method::DM=NoEnclosure()) where {N, AM, RM<:AbstractReductionMethod, DM<:AbstractDisjointnessMethod}

    # algorithm with "preallocation" is only defined for the non-static case
    preallocate = !static
    n = ismissing(dim) ? missing : Val(dim)
    p = ismissing(ngens) ? missing : Val(ngens)
    return GLGM06(δ, approx_model, max_order, Val(static), n, p,
                  Val(preallocate), reduction_method, disjointness_method)
end

step_size(alg::GLGM06) = alg.δ
numtype(::GLGM06{N}) where {N} = N

@inline function setrep(alg::GLGM06{N, AM, Val{false}}) where {N, AM}
    ST = Zonotope{N, Vector{N}, Matrix{N}}
end

@inline function rsetrep(alg::GLGM06{N, AM, Val{false}}) where {N, AM}
    ST = setrep(alg)
    RT = ReachSet{N, ST}
end

@inline function setrep(alg::GLGM06{N, AM, Val{true}, Val{n}, Val{p}}) where {N, AM, n, p}
    VT = SVector{n, N}
    MT = SMatrix{n, p, N, n*p}
    ZT = Zonotope{N, VT, MT}
end

@inline function rsetrep(alg::GLGM06{N, AM, Val{true}, Val{n}, Val{p}}) where {N, AM, n, p}
    ZT = setrep(alg)
    RT = ReachSet{N, ZT}
end

@inline function rsetrep(alg::GLGM06{N, AM, Val{true}, Missing, Missing}) where {N, AM}
    throw(ArgumentError("the reach-set representation for a statically sized zonotope " *
                        "that the system's dimension argument has been defined; add `dim=...`" *
                        "to the algorithm constructor"))
end

include("post.jl")
include("reach_homog.jl")
include("reach_inhomog.jl")
#include("check.jl")


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 685479e8-1ad5-48d8-b9fe-f2cf8a672700
using AstroImages, PlutoUI

# ╔═╡ 59e1675f-9426-4bc4-88cc-e686ed90b6b5
md"""
Download a FITS image and open it.
Apply `restrict` to downscale 2x for faster rendering.
"""

# ╔═╡ d1e5947b-2c1a-46fc-ab8f-feeba03453e7
img = AstroImages.restrict(
	AstroImage(download("http://www.astro.uvic.ca/~wthompson/astroimages/fits/656nmos.fits"))
);

# ╔═╡ c9ebe984-4630-47c1-a941-795293f5b3c1
md"""
Display options
"""

# ╔═╡ a3e81f3f-203b-47b7-ac60-b4267eddfad4
md"""

| parameter | value |
|-----------|-------|
|`cmap` | $( @bind cmap  Select([:magma, :turbo, :ice, :viridis, :seaborn_icefire_gradient, "red"]) ) |
|`clims`| $( @bind clims Select([Percent(99.5), Percent(95), Percent(80), Zscale(), (0, 400)]) ) |
| `stretch` | $( @bind stretch  Select([identity, asinhstretch, logstretch, sqrtstretch, powstretch, powerdiststretch, squarestretch])) |
| `contrast` | $(@bind contrast Slider(0:0.1:2.0, default=1.0)) |
| `bias` | $(@bind bias Slider(0:0.1:1.0, default=0.5)) |
"""

# ╔═╡ 2315ffec-dc49-413a-b0d6-1bcce2addd76
imview(img; cmap, clims, stretch, contrast, bias)

# ╔═╡ d2bd2f13-ed23-42c5-9317-5b48ec3a8bb7
md"""
## `implot`
Uncomment the following cells to use `Plots` instead.
"""

# ╔═╡ fe6b5b76-8b77-4bfc-a2e8-bcc0b78ad764
# using Plots

# ╔═╡ f557784e-828c-415e-abb0-964b3a9fe8ef
# implot(img; cmap, clims, stretch, contrast, bias)

# ╔═╡ Cell order:
# ╠═685479e8-1ad5-48d8-b9fe-f2cf8a672700
# ╟─59e1675f-9426-4bc4-88cc-e686ed90b6b5
# ╠═d1e5947b-2c1a-46fc-ab8f-feeba03453e7
# ╟─c9ebe984-4630-47c1-a941-795293f5b3c1
# ╟─a3e81f3f-203b-47b7-ac60-b4267eddfad4
# ╠═2315ffec-dc49-413a-b0d6-1bcce2addd76
# ╟─d2bd2f13-ed23-42c5-9317-5b48ec3a8bb7
# ╠═fe6b5b76-8b77-4bfc-a2e8-bcc0b78ad764
# ╠═f557784e-828c-415e-abb0-964b3a9fe8ef


text:
# # Melt rates
#
# This example demonstrates the how to use melt rate parametrizations in WAVI.jl.
# We first demonstrate the parametrizations included in WAVI.jl, by produce a map of the melt rate for each of these in the MISMIP+ steady state geometry. (MISMIP+ is the latest ice sheet model intercomparison exercise, for more info see doi:10.5194/tc-14-2283-2020)
# Explicitly, these parametrizations are:
#   * Quadratic melt rates
#   * A plume melt emulator
#   * PICO melt rate model
#   * Binary file melt rate
# You can find more info on each of these models in the Physics --> Melt Rates tab.
#
# Secondly, we demonstrate how to add a simple melt rate model to WAVI.jl.

# ## Install dependencies
#
# First let's make sure we have all required packages installed.
# As well as WAVI and Plots for plotting, we're going to use the Downloads package to pull some data from a Github repository.

#using Pkg
#pkg.add("https://github.com/RJArthern/WAVI.jl"), Pkg.add(Plots), Pkg.add("Downloads")
using WAVI, Plots, Downloads

# ## Setting up the grid
# First we'll make the grid. The MISMIP+ grid has length 640km and 80km in the x and y directions, respectively. We're going to choose 2km resolution in both directions, i.e. 320 grid points in the x-direction and 40 grid points in the y-direction:
nx = 320;
ny = 40;
grid = Grid(nx=nx,ny=ny,dx=2000.,dy=2000.);

# For the bed, we'll cheat a little and use the "hard coded" form of the MISMIP+ bed (see the "overdeepened bed" example to see how to define this bed properly)
bed = WAVI.mismip_plus_bed; # function definition

# ## Getting the ice thickness
# To allow us to focus on the melt rate parametrizations, we're not going to run a simulation to define the steady state ice thickness, but rather pull it from Github, where it is stored as a binary file.
# (Note that this ice thickness is the result of the "overdeepened bed" example, but with 2km resolution.)
fname = Downloads.download("https://github.com/alextbradley/WAVI_example_data/raw/main/MISMIP_PLUS/WAVI_ice0_2km_thick.bin");
h = Array{Float64,2}(undef, nx, ny);
read!(fname, h);
h = ntoh.(h);
# (The final line simply converts the endianness from big-endian --  the format in which this file is stored -- to little-endian.)

# Now we make an InitialConditions object to store this ice thickness:
global initial_conditions = InitialConditions(initial_thickness = h);

# ## Melt Rate Models
# We'll loop over each of the melt rate models mentioned above. For each, we'll produce a map of the melt rate in the MISMIP+ geometry that we just pulled from online.

# First let's make the melt rate models in turn, starting with quadratic, which is imposed on the model via a `QuadraticMeltRate` object.
melt_quad = QuadraticMeltRate(γT = 0.745*1e-3)
# The parameter $\gamma_T$ is a normalization cofficient. The value here is chosen so that the mean melt rate on the shelf is approx 10m/a. (See the "WAVI Setup" -->  "Melt Rate Models" section for information on all keyword parameters, here and for the below models.)

# Next up: PICO. For this parametrization, we first have to make a mask defining where the ice front is, and then pass it as a keyword argument when we instantiate the `PICO` melt rate object.
ice_front_mask = zeros(nx,ny);
ice_front_mask[end,:] .= 1;
melt_PICO = PICO(ice_front_mask = ice_front_mask,
                    T0 = 1.2,
                    S0 = 34.6,
                    γT = 0.87e-5,
                    nbox = 5);
# Again, $\gamma_T$ is a normalization coefficient.

# Next up: plume model emulator:
melt_PME = PlumeEmulator(α=1.49);
# In this case, the normalization coefficent is named $\alpha$, and has a slightly different meaning to $\gamma_T$ above, but we use it in the same role: to the set the mean melt rate on the shelf.

# Finally, a binary file melt rate, in which the melt rate is read in from a binary file. First, we'll create such a file, which will set the melt rate to be uniform on the shelf.
isfloat = (h .< -918.0/1028.0 .* bed.(grid.xxh, grid.yyh)) #indices of floating elements
m = zeros(nx,ny);
m[isfloat] .= 10.0; #set everywhere floating to 10m/a
folder = joinpath(@__DIR__, "melt_rate_parametrizations");
isdir(folder) && rm(folder, force = true, recursive = true);
mkdir(folder) ;
out = open(joinpath(folder,"melt.bin"), "w");
write(out, m);
close(out);

# Now we instantiate out melt rate opject, pointing it to this file we just created:
binfile_melt = BinfileMeltRate(input_filename = joinpath(folder,"melt.bin"));

# It's useful to put these into a dictionary, so we can iterate over them:
global melt_rates = Dict("Quadratic" => melt_quad, "PME" => melt_PME, "PICO" => melt_PICO, "Binary file" => binfile_melt);

# ## Visualizing
# We'll loop over the melt rate models instantiated above. Each time, we make a WAVI.jl `Model` with the appropriate melt rate parametrization specified via the `melt_rate` keyword.
# We use the `update_state!` method to bring the melt (as well as all other quantities, such as the grounded fraction) in line with the specified thickness, and then plot the melt rate.
# (We alsoset the number of Picard iterations to 1 -- each time we do an `update_state!` we perform a velocity solve; in this example, we don't care about the velocity solve being completely accurate, so we'll do the minimum number of iterations in this process!)
for (key, melt) in melt_rates
    #instantiate the model with the appropriate melt rate
    local model = Model(grid = grid,
            bed_elevation = bed,
            initial_conditions = initial_conditions,
            solver_params = SolverParams(maxiter_picard=1),
            melt_rate = melt);

    #bring quantities in line with the thickness
    update_state!(model);

    #extract the melt rate, remove any grounded entries and saturate the melt rate to 50 m/a
    mcopy= deepcopy(model.fields.gh.basal_melt)
    mcopy[model.fields.gh.grounded_fraction .== 1.] .= NaN
    msat = deepcopy(mcopy)
    msat[msat .> 50] .= 50

    #plot the melt rate
    plt = Plots.heatmap(model.grid.xxh[:,1]/1e3, model.grid.yyh[1,:]/1e3, msat',
                        xlabel = "x (km)",
                        ylabel = "y (km)",
                        colorbar_title = "melt rate (m/yr)",
                        title = key,
                        framestyle = "box")
    xlims!((420, 640))
    plot!(size = (500,300))
    #display(plt) #uncomment to show in (e.g.) VSCode
end

# As expected, the binary file melt rate has the same (10m/a)melt rate over the whole shelf.
# The PICO parametrization, which divides the shelf up into discrete blocks, has a corresponding banded structure, with highest melt rates at the grounding line.
# The quadratic melt rate parametrizations similarly has the highest melt rate near the grounding line, but drops off with distance from the grounding line much quicker than the PICO parametrization.
# These plots can be compared to corresponding results for the Elmer ice sheet model  (Favier et al. 2019 doi:10.5194/gmd-12-2255-2019)

# ## Defining A New Melt Rate Model
# In this section, we show how to specify a new melt rate model. There are four steps:
# * Create a file to store code
# * Define the appropriate structure, which stores information required to prescribe the melt rate
# * Define a constructor of this structure
# * Write a function to update the melt rate appropriately, and export this
# This is quite abstract, so let's do an example. We'll create a melt rate model which sets the melt rate as it is specified in the MISMIP+ experiment:
# the melt rate on floating cells is $0.2 \tanh((z_d - z_b)/75) \max(-100 - z_d,0)$, where $z_d$ is the ice shelf draft and $z_d - z_b$ is the cavity thickness.
# We'll follow the steps above: first, we create a file to store the code. For this example, we've already create the file, you can see it at it at "src/MeltRate/mismip_melt_rate.jl".

# Next we define a structure, which stores parameters related to the melt rate model. Note that the melt rate model does not "own" the melt rate, the `model` does (and stores it in model.fields.gh.basal_melt, see below)
"""
struct MISMIPMeltRateOne{T <: Real} <: AbstractMeltRate
    α  :: T
    ρi :: T
    ρw :: T
end
"""
# In this case, a normalization coefficient $\alpha$, and the denisities of ice and ocean, $\rho_i$ and $\rho_w$, respectively.

# Now we define our "constructor", a function that defines how to create one of these structures:
"""
MISMIPMeltRateOne(; α = 1.0, ρi = 918.0, ρw = 1028.0) = MISMIPMeltRateOne(α,ρi, ρw)
"""
# In this case, the constructor simply sets the default values for the parameters $\alpha$, $\rho_i$, and $\rho_w$.
# For more complicated melt rate models, constructors might be more elaborate!

# The final step is to define a function `update_melt_rate!(melt_rate::TYPE, fields, grid)` which tells WAVI how to update the melt rate in this example.
# here, TYPE is the name of the structure we just made.
# Note that the arguments of this function must be as mentioned here, so that the multiple dispatch capability can be leveraged! This procedure defines another method named `update_melt_rate!`, which sets the melt rate according to this function when the input melt rate is of type "TYPE". One of these function is defined for each of the melt rate models mentioned above.
"""
function update_melt_rate!(melt_rate::MISMIPMeltRateOne, fields, grid)
    draft = -(melt_rate.ρi / melt_rate.ρw) .* fields.gh.h
    cavity_thickness = draft .- fields.gh.b
    cavity_thickness = max.(cavity_thickness, 0)
    m =  melt_rate.α .* 0.2*tanh.(cavity_thickness./75).*max.((-100 .- draft), 0)
    fields.gh.basal_melt[:] .= m[:]
end
"""

# Finally, we tell point WAVI to this code by adding `include("mismip_melt_rate.jl")` to the file "src/MeltRate/MeltRate.jl", and add this structure to the export section in "src/WAVI.jl" file.

# Now we can create a model which takes this melt rate andplot the result:
model = Model(grid = grid,
            bed_elevation = bed,
            initial_conditions = initial_conditions,
            solver_params = SolverParams(maxiter_picard=1),
            melt_rate = MISMIPMeltRateOne());

update_state!(model);
m = deepcopy(model.fields.gh.basal_melt)
m[model.fields.gh.grounded_fraction .== 1.] .= NaN
msat = deepcopy(m)
msat[msat .> 50] .= 50

plt = Plots.heatmap(model.grid.xxh[:,1]/1e3, model.grid.yyh[1,:]/1e3, msat',
            xlabel = "x (km)",
            ylabel = "y (km)",
            colorbar_title = "melt rate (m/yr)",
            title = "MISMIP melt rate",
            framestyle = "box")
xlims!((420, 640))
plot!(size = (500,300))
#display(plt)

# Hopefully this example demonstrates clearly the procedure for adding melt rate models to WAVI.jl. If there are any questions, don't hesistate to get in touch (see the "Contact Us" tab)

# Finally, let's clean up the files we just made
rm(joinpath(@__DIR__, "melt_rate_parametrizations"), force = true, recursive = true);


text:
"""
`AbstractGeoArray` wraps an array (or location of an array) and metadata
about its contents. It may be memory ([`GeoArray`](@ref)) or disk-backed
([`NCDarray`](@ref), [`GDALarray`](@ref), [`GRDarray`](@ref)).

`AbstractGeoArray`s inherit from [`AbstractDimArray`]($DDarraydocs)
from DimensionalData.jl. They can be indexed as regular Julia arrays or with
DimensionalData.jl [`Dimension`]($DDdimdocs)s. They will plot as a heatmap in
Plots.jl with correct coordinates and labels, even after slicing with
`getindex` or `view`. `getindex` on a `AbstractGeoArray` will always return
a memory-backed `GeoArray`.
"""
abstract type AbstractGeoArray{T,N,D,A} <: AbstractDimensionalArray{T,N,D,A} end

# Marker singleton for lazy loaded arrays, only used for broadcasting.
# Can be removed when DiskArrays.jl is used everywhere
struct LazyArray{T,N} <: AbstractArray{T,N} end

# Interface methods ###########################################################

"""
    missingval(x)

Returns the value representing missing data in the dataset
"""
function missingval end
missingval(x) = missing
missingval(A::AbstractGeoArray) = A.missingval

"""
    crs(x)

Get the projected coordinate reference system of a `Lat` or `Lon` `Dimension`,
or of the `Lat`/`Lon` dims of an `AbstractGeoArray`.

For [`Mapped`](@ref) mode this may be `nothing` as there may be not projected
coordinate reference system at all.
"""
function crs end
crs(A::AbstractGeoArray) =
    if hasdim(A, Lat)
        crs(dims(A, Lat))
    elseif hasdim(A, Lon)
        crs(dims(A, Lon))
    else
        error("No Lat or Lon dimension, crs not available")
    end
crs(dim::Dimension) = crs(mode(dim))

"""
    mappedcrs(x)

Get the mapped coordinate reference system for the `Lat`/`Lon` dims of an array.

In [`Projected`](@ref) mode this is used to convert [`Selector`]($DDselectordocs)
values form the mappedcrs defined projection to the underlying projection, and to
show plot axes in the mapped projection.

In `Mapped` mode this is the coordinate reference system of the index values.
"""
function mappedcrs end
mappedcrs(A::AbstractGeoArray) =
    if hasdim(A, Lat)
        mappedcrs(dims(A, Lat))
    elseif hasdim(A, Lon)
        mappedcrs(dims(A, Lon))
    else
        error("No Lat or Lon dimension, mappedcrs not available")
    end
mappedcrs(dim::Dimension) = mappedcrs(mode(dim))

# DimensionalData methods

units(A::AbstractGeoArray) = getmeta(A, :units, nothing)

for f in (:mappedbounds, :projectedbounds, :mappedindex, :projectedindex)
    @eval ($f)(A::AbstractGeoArray, dims_) = ($f)(dims(A, dims_))
    @eval ($f)(A::AbstractGeoArray) = ($f)(dims(A))
end

# Rebuild all types of AbstractGeoArray as GeoArray
rebuild(A::AbstractGeoArray, data, dims::Tuple, refdims, name, metadata, missingval=missingval(A)) =
    GeoArray(data, dims, refdims, name, metadata, missingval)
rebuild(A::AbstractGeoArray; data=data(A), dims=dims(A), refdims=refdims(A), name=name(A), metadata=metadata(A), missingval=missingval(A)) =
    GeoArray(data, dims, refdims, name, metadata, missingval)

Base.parent(A::AbstractGeoArray) = data(A)

"""
Abstract supertype for all memory-backed GeoArrays where the data is an array.
"""
abstract type MemGeoArray{T,N,D,A} <: AbstractGeoArray{T,N,D,A} end

"""
Abstract supertype for all disk-backed GeoArrays.
For these the data is lazyily loaded from disk.

To load a `DiskGeoArray` and operate on the data multiple times, use
[`Open`](@ref) and a `do` block.
"""
abstract type DiskGeoArray{T,N,D,A} <: AbstractGeoArray{T,N,D,A} end

filename(A::DiskGeoArray) = A.filename

"""
    data(f, A::DiskGeoArray)

Run method `f` on the data source object for `A`, as passed by the
`withdata` method for the array. The only requirement of the
object is that it has an `Array` method that returns the data as an array.
"""
data(A::DiskGeoArray) = withsourcedata(Array, A)

# Base methods

Base.size(A::DiskGeoArray) = A.size

Base.getindex(A::DiskGeoArray, i1::StandardIndices, i2::StandardIndices, I::StandardIndices...) =
    rebuildgetindex(A, i1, i2, I...)
Base.getindex(A::DiskGeoArray, i1::Integer, i2::Integer, I::Vararg{<:Integer}) =
    rawgetindex(A, i1, i2, I...)
# Linear indexing returns Array
Base.@propagate_inbounds Base.getindex(A::DiskGeoArray, i::Union{Colon,AbstractVector{<:Integer}}) =
    rawgetindex(A, i)
# Exempt 1D DimArrays
Base.@propagate_inbounds Base.getindex(A::DiskGeoArray{<:Any,1}, i::Union{Colon,AbstractVector{<:Integer}}) =
    rebuildgetindex(A, i)

rawgetindex(A, I...) =
    withsourcedata(A) do data
        readwindowed(data, I...)
    end

rebuildgetindex(A, I...) =
    withsourcedata(A) do data
        dims_, refdims_ = slicedims(dims(A), refdims(A), I)
        data = readwindowed(data, I...)
        rebuild(A, data, dims_, refdims_)
    end

Base.write(A::T) where T <: DiskGeoArray = write(filename(A), A)
Base.write(filename::AbstractString, A::T) where T <: DiskGeoArray =
    write(filename, T, A)

# Concrete implementation ######################################################

"""
    GeoArray(A::AbstractArray{T,N}, dims::Tuple;
             refdims=(), name=Symbol(""), metadata=NoMetadata(), missingval=missing)
    GeoArray(A::AbstractArray{T,N};
             dims, refdims=(), name=Symbol(""), metadata=NoMetadata(), missingval=missing)
    GeoArray(A::AbstractGeoArray; [data=data(A), dims=dims(A), refdims=refdims(A),
             name=name(A), metadata=metadata(A), missingval=missingval(A)]) =

A generic, memory-backed spatial array type. All [`AbstractGeoArray`](@ref) are
converted to `GeoArray` when indexed or otherwise transformed.

# Keyword Arguments

- `name`: `Symbol` name for the array.
- `dims`: `Tuple` of `Dimension`s for the array.
- `refdims`: `Tuple of` position `Dimension`s the array was sliced from,
  defaulting to `()`.
- `missingval`: Value reprsenting missing values, defaulting to `missing`.
  can be passed it.
- `metadata`: [`Metadata`](@ref) object for the array, or `NoMetadata()`.

## Example

```julia
A = GRDarray(gdalarray; name="surfacetemp")
# Select Australia using lat/lon coords, whatever the crs is underneath.
A[Lat(Between(-10, -43), Lon(Between(113, 153)))
```
"""
struct GeoArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Na,Me,Mi} <: MemGeoArray{T,N,D,A}
    data::A
    dims::D
    refdims::R
    name::Na
    metadata::Me
    missingval::Mi
end
@inline GeoArray(A::AbstractArray, dims::Tuple;
                 refdims=(), name=Symbol(""), metadata=NoMetadata(), missingval=missing) =
    GeoArray(A, formatdims(A, dims), refdims, name, metadata, missingval)
@inline GeoArray(A::AbstractArray; dims, refdims=(), name=Symbol(""), metadata=NoMetadata(), missingval=missing) =
    GeoArray(A, formatdims(A, dims), refdims, name, metadata, missingval)
@inline GeoArray(A::AbstractGeoArray; data=data(A), dims=dims(A), refdims=refdims(A),
                 name=name(A), metadata=metadata(A), missingval=missingval(A)) =
    GeoArray(data, dims, refdims, name, metadata, missingval)

Base.@propagate_inbounds Base.setindex!(A::GeoArray, x, I::Vararg{DimensionalData.StandardIndices}) =
    setindex!(data(A), x, I...)

Base.convert(::Type{GeoArray}, array::AbstractGeoArray) = GeoArray(array)


text:
const σ₀=I
const σ₁=[0 1;1 0]
const σ₂=[0 -im;im 0]
const σ₃=[1 0;0 -1]
const sigmax = [0 1;1 0]
const sigmay = [0 -im;im 0]
const sigmaz = [1 0;0 -1]
const hadamard = [1/sqrt(2) 1/sqrt(2);1/sqrt(2) -1/sqrt(2)]

export σ₀,σ₁,σ₂,σ₃,sigmax,sigmay,sigmaz,hadamard

max_stars_repo_name: Roger-luo/QuCmp

max_issues_repo_name: Roger-luo/QuCmp

max_forks_repo_name: Roger-luo/QuCmp


text:
arr = [1, 2, 23, 24, 6]
push!(arr, 7)
println(arr)

arr2 = [11, 12, 13]
append!(arr, arr2)
println(arr, ", ", arr2)

sort(arr) # Returns copy
println(arr)
sort!(arr) # Mutates array
println(arr)

println("sum(arr): $(sum(arr))")

function addone(x)
    return x + 1
end

a = map(addone, arr)
println(a)
b = zeros(3)
map!(addone, b, arr)
println(b)

c = [1, 2, 3, 4, 5]
filter(iseven, c)
println(c)
d = [1, 2, 3]
filter!(iseven, d)
println(d)

println(reduce(+, c))

println("c = $(c)")
println("c .+ 5 = $(c .+ 5)")
strings = ["this", "is", "an", "array"]
println("strings = $(strings)")
println("uppercase.(strings) = $(uppercase.(strings))")

println("splice!(strings, 3): $(splice!(strings, 3)),\tstrings: $strings")
println("pop!(strings): $(pop!(strings)),\tstrings: $strings")
println("arr: $(arr)")
println("popfirst!(arr): $(popfirst!(arr)),\tstrings: $arr")
println("pushfirst!(arr, 100, 101): $(pushfirst!(arr, 100, 101))")
println("deleteat!(arr, 3): $(deleteat!(arr, 3))")
println("insert!(arr, 5, 33): $(insert!(arr, 5, 33))")

println(collect(1:2:13))
println(collect("a b c d e f g"))
println(split("a b c d e f g"))
println(split("a_b_c_d_e_f_g", '_'))
t = ["I", "am", "groot"]
println(join(t, ' '))
println(join(t, '-'))

e = [1, 2, 3]
f = [4, 5, 6]
g = vcat(e, f)
println("e: $e,\tf: $f,\tg: $g")


text:
immutable NoiseModelGaussian{T<:AbstractFloat,V<:AbstractVector} <: AbstractNoiseModel
    variance::V
    distributions::Vector{Distributions.Normal{T}}
    NoiseModelGaussian(v::AbstractVector{T},d::Vector{Distributions.Normal{T}}) = new(v,d)
end

@inline _noise(::Type{Val{:gaussian}},v::AbstractVector) = NoiseModelGaussian{eltype(v),typeof(v)}(v,map(Distributions.Normal,zeros(v),sqrt(v)))

function loglikelihood(n::NoiseModelGaussian,measurements::AbstractArray,modeldata::AbstractArray)
    r = 0
    for j=1:length(n.distributions)
        @simd for i=1:size(measurements,1)
            @inbounds r+= logpdf(n.distributions[j],modeldata[i,j]-measurements[i,j])
        end
    end
    r
end

function applynoise!{T<:AbstractFloat}(n::NoiseModelGaussian,d::AbstractArray{T})
    for j=1:length(n.distributions)
        @simd for i=1:size(d,1)
            @inbounds d[i,j] += convert(T,rand(n.distributions[j]))
        end
    end
    d
end

noisemodelname(::NoiseModelGaussian) = "Gaussian Noise Model"


text:
#### Azure

export Azure

"""
    Azure

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Azure
    "Organization name"
    org::String
    "Package name (with extension)"
    pkg_name::String
    "Julia version"
    julia_version::String
    "Checks"
    checks::String="$(org).$(pkg_name)"
end

function toml(ci::Azure)
    return "trigger:
  batch: true
  branches:
    include:
    - staging
    - trying

jobs:
- job: Linux

  timeoutInMinutes: 0

  pool:
    vmImage: 'ubuntu-16.04'

  strategy:
    matrix:
      Julia $(ci.julia_version):
        JULIA_VERSION: '$(ci.julia_version)'

  steps:
  - bash: |
      set -o xtrace
      wget -nv https://julialang-s3.julialang.org/bin/linux/x64/\$(JULIA_VERSION)/julia-\$(JULIA_VERSION)-latest-linux-x86_64.tar.gz
      mkdir julia-\$(JULIA_VERSION)
      tar zxf julia-\$(JULIA_VERSION)-latest-linux-x86_64.tar.gz -C julia-\$(JULIA_VERSION) --strip-components 1
    displayName: 'Download and extract Julia'
  - bash: |
      set -o xtrace
      ./julia-\$(JULIA_VERSION)/bin/julia -e 'using InteractiveUtils; versioninfo()'
      ./julia-\$(JULIA_VERSION)/bin/julia --project=@. -e 'using Pkg; Pkg.instantiate()'
      ./julia-\$(JULIA_VERSION)/bin/julia --project=@. -e 'using Pkg; Pkg.test()'
    displayName: 'Run the tests'
    continueOnError: true

- job: macOS

  timeoutInMinutes: 0

  pool:
    vmImage: 'macOS-10.14'

  strategy:
    matrix:
      Julia $(ci.julia_version):
        JULIA_VERSION: '$(ci.julia_version)'

  steps:
  - bash: |
      set -o xtrace
      wget -nv https://julialang-s3.julialang.org/bin/mac/x64/\$(JULIA_VERSION)/julia-\$(JULIA_VERSION)-latest-mac64.dmg
      mkdir juliamnt
      hdiutil mount -readonly -mountpoint juliamnt julia-\$(JULIA_VERSION)-latest-mac64.dmg
      cp -a juliamnt/*.app/Contents/Resources/julia julia-\$(JULIA_VERSION)
    displayName: 'Download and extract Julia'
  - bash: |
      set -o xtrace
      ./julia-\$(JULIA_VERSION)/bin/julia -e 'using InteractiveUtils; versioninfo()'
      ./julia-\$(JULIA_VERSION)/bin/julia --project=@. -e 'using Pkg; Pkg.instantiate()'
      ./julia-\$(JULIA_VERSION)/bin/julia --project=@. -e 'using Pkg; Pkg.test()'
    env:
      MPICH_INTERFACE_HOSTNAME: localhost
    displayName: 'Run the tests'

- job: Windows

  timeoutInMinutes: 0

  pool:
    vmImage: 'VS2017-Win2016'

  strategy:
    matrix:
      Julia $(ci.julia_version):
        JULIA_VERSION: '$(ci.julia_version)'

  continueOnError: true

  steps:
  - powershell: |
      Set-PSDebug -Trace 1
      wget https://julialang-s3.julialang.org/bin/winnt/x64/\$(JULIA_VERSION)/julia-\$(JULIA_VERSION)-latest-win64.exe -OutFile julia-\$(JULIA_VERSION)-latest-win64.exe
      Start-Process -FilePath .\\julia-\$(JULIA_VERSION)-latest-win64.exe -ArgumentList \"/S /D=C:\\julia-$(JULIA_VERSION)\" -NoNewWindow -Wait
    displayName: 'Download and extract Julia'
  - powershell: |
      Set-PSDebug -Trace 1
      C:\\julia-$(JULIA_VERSION)\\bin\\julia.exe -e 'using InteractiveUtils; versioninfo()'
      C:\\julia-$(JULIA_VERSION)\\bin\\julia.exe --project=@. -e 'using Pkg; Pkg.instantiate()'
      C:\\julia-$(JULIA_VERSION)\\bin\\julia.exe --project=@. -e 'using Pkg; Pkg.test()'
    displayName: 'Run the tests'

"
end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates